Soil Formation and Mass Redistribution during the Holocene Using Meteoric 10 Be, Soil Chemistry and Mineralogy

: Soil development and erosion are important and opposing processes in the evolution of high-mountainous landscapes, though their dynamics are not fully understood. We compared soil development between a calcareous and a siliceous chronosequence in the central Swiss Alps at high altitudes, which both cover soil formation over the Holocene. We calculated element mass balances, long-term erosion rates based on meteoric 10 Be and we determined the rates of soil formation. We also analyzed the shifts in the mineralogical composition, weathering indices, the particle size distribution, carbon stocks and oxalate extractable Fe, Al, and Mn. The siliceous soils had high chemical weathering rates at the early stage of soil formation that strongly decreased after a few millennia. The development of calcareous soil was characterized by high carbonate losses and a shift to finer soil texture. Soil erosion hampered the upbuilding of soil horizons in the early stages of soil development, which led to a delay in soil and vegetation development. This study shows how soil formation drivers change over time. In the early stages of soil development, the parent material predominantly drives soil formation while at later stages the vegetation becomes more dominant as it influences surface stability, hydrological pathways, and that water drainage and retention. The pH was determined by mixing 5 g of fine earth with 12.5 mL of 0.01 M CaCl 2 , equating to a ratio of 1:2.5. We conducted the measurements on a Mehtrohm 692 pH/IonMeter with an absolute error of ±0.003. The carbonate content was calculated from total inorganic carbon (C inorg ) measurements using a MultiEA (Analytik Jena, Jena, Germany), which oper-ates with phosphoric acid.


Introduction
Soils and their formation, i.e., the physical and chemical alteration of the parent material, have been studied for more than a century. Dokuchaev [1] laid the foundation of soil science with his observations in the Russian steppe. He also laid the ground for the factorial model of soil development, which Jenny developed further in 1941 in his publication, Factors of Soil Formation [2]. Jenny described soil formation as the product of five independent factors: climate, organisms, relief, parent material, and time. Over the years, other models have been developed. Simonson [3] presented a model which looks at soils from a process-oriented point of view (i.e., additions, removals, translocations). Johnson and Watson-Stegner [4] introduced the concept of pathways. Their model allows soil formation to move in two directions and acknowledges soil formation as a highly dynamic process: soils experience progressive pedogenesis phases (characterized by horizonation, soil deepening, and developmental upbuilding) or regressive pedogenetic phases (halploidization, soil thinning, and retardant upbuilding). As the environmental conditions are never entirely stable, soil development may switch between progressive and regressive phases. The state of a soil at one point in time is merely the integration of all past soil development processes [4][5][6].
The formation of soils is in essence a balance between soil production (the conversion of bedrock or parent material to soil) and denudation (physical and chemical soil loss) [7]. This can be expressed as: Fsoil is the soil formation rate, Psoil the soil production rate and Dsoil the soil denudation rate.
Soil production (Psoil) is the sum of transformation of parent material to soil (TPsoil), atmospheric (dust) input (A), net organic matter input (O), and organic matter decay (G) (e.g., [8] Denudation (Dsoil) is the sum of chemical weathering (Wsoil) and erosion (Esoil): It has been shown that chemical weathering rates and physical erosion are closely linked [9]. Furthermore, the weathered layer acts as a protective layer which hinders weathering of the underlying parent material [10].
In glacial landscapes, glaciers shrink as a result of climate warming, exposing new surfaces and form fresh ground leading to significant amounts of sediment being eroded and accumulated downstream [11,12]. Periglacial areas are therefore useful to study processes related to early soil formation, chemical weathering, and erosion processes. They can provide good age constraints, and the till which generally forms the parent material can be regarded as well mixed. Based on the concept of Jenny [2] according to which the factor time and its influence on soil formation can be isolated when the other factors remain constant have become widely used approaches to trace and compare the transformation of soils over various timescales [13,14]. Chronosequences have been applied in various locations of the Alps for many decades. They show that most soil properties change very rapidly at early stages of soil development [15][16][17]. With time, this rate of alteration usually declines, reaching a plateau [14,[18][19][20][21][22]. This is due to the high abundance of easily weatherable material in proglacial areas, because the glacial till has already undergone considerable physical weathering and provides, based on its grain size distribution, large surface areas [14]. Thus, weathering fluxes are usually highest in the beginning and taper off, usually after several millennia [23,24]. The rates depend on the composition and properties of the parent material. Granite or gneiss-based parent materials are expected to lose mainly silicon, iron, and aluminum, as they are the major components of micas and feldspars [25]. Micas have been observed to weather easily and be transformed into vermiculite and smectites, i.e., clay minerals [26]. Limestone-based soils will mainly loose calcium and varying amounts of magnesium through the dissolution of carbonates [23]. These shifts in the mineralogy will alter the physical and chemical soil properties, leading to finer texture and a higher proportion of fine pores which has a generally positive effect on the water retention capabilities of the soil [27,28].
When discussing soil development, it is important to consider the development of the biosphere as it influences the substrate considerably. Although Jenny defined it as independent factor, it is not. As both the soils and the vegetation influence each other simultaneously, their connection is not easy to quantify [29,30]. In the early stages of soil development, the processes of physical and chemical weathering improve the conditions for plant growth by providing the substrate for root growth and plant stability, water storage, and a steady supply of nutrients. At first, a combination of pioneer plants, which reproduce rapidly but are not very resilient, will add organic matter to the soils at a high rate via roots [31][32][33]. Roots are able to release large quantities of C to the soils as root detritus, exudates, and by transferring C to root symbionts [34,35]. Organic carbon and nitrogen can therefore be expected to accumulate rapidly in the first centuries after deglaciation [27,36,37]. The plants stabilize the soil surfaces and thus reduce soil erosion, and promote soil deepening [20,29,30]. As erosion is reduced [38] and the soil conditions provide more favorable and stable habitats, the composition of plant communities shift towards more resistant and competitive species. If the conditions are favorable, primary successions from pioneer species to trees can develop in less than two centuries in the Alps [39]. The subsequent accumulation of soil organic matter (SOM) [40] plays a crucial role in increasing soil water retention (by increasing the microporosity, [41]) and nutrient supply [42] which further improve the soil's qualities as a habitat.
The extent to which the vegetation influences the substrate depends not only on the degree of coverage but also on the functional plant traits [43]. In a species-rich plant community, the diversity of different traits is assumed to positively influence ecosystem functioning, e.g., by the complementarity of nutrient recycling (diversity hypothesis). In order to quantify these relationships, scientists often use the concept of functional diversity, which is calculated as the combination of the abundance of plant species and their functional traits [44]. In principle, functional diversity cannot be described by one single index. One prominent measure of functional diversity often related to soil processes is the functional richness (FRic) meaning the total amount of niche space filled by the traits of all species in a certain community [45].
Meteoric 10 Be is a new and novel method to obtain information about the balance between soil production and soil erosion [46]. Meteoric 10 Be is, among others, produced in the upper atmosphere through spallation of nitrogen and oxygen by cosmic radiation [47]. Precipitation washes it out of the atmosphere and deposits it in soils (hence it is referred to as 'meteoric' 10 Be, [48]). It is adsorbed strongly to soil particles which makes it a useful tracer for long-term soil erosion and sedimentation processes [46], soil residence times (e.g., [49]), and soil production (e.g., [50]).This method enables us to calculate the net soil loss during soil evolution. In contrast, fallout radionuclides (e.g., 239+240 Pu) cover the last 60 years and enable the determination of erosion rates at a decadal scale. As soil erodibility generally decreases with (progressive) soil formation, the erosion rates are expected to decrease with increasing surface age. Vegetation and abiotic factors such as the topography and climatic conditions shape this temporal development. However, with increasing age, soils become less permeable and, thus, generate more surface runoff [51] which can induce more surface erosion. The roughness of the soil surface additionally influences the spatial distribution of the surface runoff. A rougher surface can lead to more channeling and, thus, increasing erodibility [52].
In order to improve our understanding of how high-mountainous (and, by extension, arctic) environments will develop under future climate warming, we must understand and quantify the processes that drive landscape evolution. This study explores soil formation in two proglacial areas with the aim to qualitatively and quantitatively assess the progression of chemical weathering and long-term erosion rates. We investigated soils on moraines of two periglacial chronosequences (one on granite/gneiss and the other one on carbonate parent material) which covered the Holocene and the Late Pleistocene. We used meteoric 10 Be to measure long-term soil erosion rates and related them with our data on weathering. We linked the results to short-term erosion rates [53], a vegetational survey [54], and hydrological properties of the soils [55,56], including a flow accumulation analysis, all of which were investigated in the same area as part of the same project.

Klausenpass
We studied the proglacial area of the Griess Glacier (2000-2200 m a.s.l.), a small glacier near the Klausenpass, central Switzerland ( Figure 1). The glacial sediment is mostly composed of limestone from the Early Cretaceous and the Late Jurassic period with some addition of Flysch [57]. The mean annual temperature in the area is about +2 °C (annual averages for the 1961-1990 period, [58]) and the mean annual precipitation is ca. 2000 mm  period, [59]). The snowy season typically lasts from mid-October to mid-June.
The chronosequence comprised four moraines: 110 a, 160 a, 4.9 ka, 13.5 ka (Table 1, Figures 2 and 3). We estimated the age of the 110 a site by tracking the glacier's retreat using historical maps [60]. The 160 a moraine corresponds to the Little Ice Age (LIA). We can identify LIA moraines with ease in the field and on aerial images as they are large and well preserved. Moreover, glaciers in the Alps advanced to almost the same extent as previously reached 3000 years ago, so that they form a prominent boundary between welldeveloped vegetation on the outside and the sparsely vegetated glacial foreland on the inside [61]. Approximations of the surface ages of the second oldest and the oldest moraines were obtained by radiocarbon dating of the H2O2 resistant fraction of soil organic matter in the top horizons [62]. Radiocarbon dating of soil organic matter is useful when no other dating methods are applicable, but there are limitations which need to be kept in mind. In general, the ages of the H2O2 resistant fraction indicate the oldest organic matter fraction in soils. These ages should rather be considered as minimum ages, because organic matter accumulation on moraines starts only then when plants are present which requires some time, especially when erosion rates are high (which may be the case in mountain soils, [53]). Musso et al. [62] obtained (2 σ) age ranges of 4969-4852 a cal BP for the second oldest and 14,001-13,760 and 13,287-13,120 a cal BP for the oldest moraine ( Table 2 in [62]). They used the 'averages' 4900 a cal BP and 13,500 a cal BP as the respective surface ages. For consistency and simplicity, this paper will use the same ages and henceforth refer to these two moraine ages as '4.9 ka' and '13.5 ka'.

Sustenpass
The study site at Sustenpass is located in the proglacial area of the Stein Glacier (Figure 1). The parent material originates from the Aarmassif and consists of metagranitoids (variscan intrusion), gneiss, amphibolite, and mica-rich schist [57]. The mean annual temperature is ca 0-2 °C (1961-1990 period, [58]) and the mean annual precipitation ca 1700-2000 mm  period, [59]) in the study area. The snowy season typically lasts from mid-October to mid-June.
The Sustenpass chronosequence had four moraines: 30 a, 160 a, 3 ka, and 10 ka (Figures 2 and 3, Table 1). The age and location of the 30 a site was estimated by tracking the glacier's retreat after its advance in the mid-1980s using maps as well as aerial images. The soil barely shows any weathering and the vegetation is sparse and consists mostly of an initial grassland vegetation. The 160 a moraine (corresponding to the LIA) is well preserved in the east side of the valley, but it was not possible to sample there for logistical reasons. We therefore studied historical maps and located a suitable area on the opposite side of the valley which had the same surface age (Figures 1 and 2, [63][64][65]).
The proglacial area of Steingletscher was previously studied by King [66] and Heikkinen and Fogelberg [67]. King identified three ridges ( Figure 1) which originated from different glacial advances. Two LIA moraines, one from the 17th century and one from the 19th century, are situated directly behind a distinctly older moraine. Schimmelpfennig et al. [68] provided 10 Be surface exposure ages on boulders and determined that the oldest of these three ridges was deposited nearly 3000 years ago. The 10 ka site was dated using the same method.   The profile name codes include the location (K/S = Klausen/Susten), the age category (A = youngest, D = oldest) and the whether the functional richness of the vegetation at those soil pits' surroundings was high (=1) or low (=2). The functional richness is a measure encompassing plant species and their traits of a certain community. It is used to predict ecological functions and processes. The bottom depth of each profile represents the depth to which we were able to dig the profile (y-axis = depth in cm). Note that in the 30 years profiles at Sustenpass, we found both rounded and unrounded rocks. In general, the depicted amount and sizes of the rocks represent the amount and size of the largest rocks encountered in the soil pits.  [62], aspect, slope (degrees), and elevation (m a.s.l.) and some basic soil parameters. The lowest depth of each profile refers to the depth that was reached during the dig (see also Figure 3). The soils were classified using the World Reference Base (WRB) classification [69]. The pH was measured with an absolute instrument error of ±0.003. LOI (Loss on ignition) represents the soil organic matter content of the fine earth fraction. The relative errors between lab replicates are 0.1%.

Sampling Strategy
This study is part of a series of studies which all used the same sampling scheme in order to ensure maximum comparability. Both the calcareous (Klausenpass) and the siliceous (Sustenpass) chronosequence comprised four moraines at which we excavated two soil pits each. We based the location of the pits on a vegetational mapping and the functional diversity of the vegetation. The functional diversity measures those components of plant diversity which influence the ecological functions and processes. The index used was the Functional Richness Index (FRI). It is based on the value, range, distribution, and the abundance of different plant traits in a vegetation community [70]. The parameters taken into account were: leaf dry matter content, specific leaf area, vegetative height [68], the plant life form, and the level of woodiness [71]. Soil pits were dug at two sites on each moraine where the FRI was highest and lowest, intending to cover the whole spectrum and see any effects, should they occur.

Physical and Chemical Properties
The skeleton (wt.%) is the weight percentage of the fraction > 2 mm and was obtained by dry-sieving 1 kg of soil. Note that, for logistical reasons, we did not include rocks larger than roughly 5 cm. Therefore, the skeleton content presented here underestimates the actual situation in the field as there were many large cobbles and boulders in the soil profiles (Figures 3 and 4). The bulk density was determined by sampling soil with an Eijkelkamp soil corer (volume: 100 cm 3 ) and measuring its dry weight. Due to the size of the steel cylinders, rocks larger than a few centimeters were not included in the sampling, thus leading to a slight underestimation of the bulk density. The values presented (Table 1) are the means of 2-3 field replicates. The organic matter (OM) content was obtained via loss on ignition (LOI). 2 g of fine earth were dry-ashed at 550 °C for 6 h. Corrected for the skeleton content, the lost material represents the organic matter content of the soil [72]. The pH was determined by mixing 5 g of fine earth with 12.5 mL of 0.01 M CaCl2, equating to a ratio of 1:2.5. We conducted the measurements on a Mehtrohm 692 pH/IonMeter with an absolute error of ±0.003. The carbonate content was calculated from total inorganic carbon (Cinorg) measurements using a MultiEA (Analytik Jena, Jena, Germany), which operates with phosphoric acid. The total element contents were measured by X-ray fluorescence (XRF, Spectro Xepos, Spectro Analytical Instruments GmbH, Kleve, Germany). The instrument is calibrated using two glass beads (FLX-SP1, FLX-SP2, Fluxana, Bedburg-Hau, Germany). We used the certified reference soil sample SO-4 (Canada Centre for Mineral and Energy Technology, Canada) as a control standard during measurements to detect drifts. To measure the samples, 5 g of fine earth were finely milled and measured as powder samples. The instrument performs 6 repeat measurements per sample and outputs mean and standard error.
We calculated the oxide forms of the most abundant elements in the soils (Na, Mg, Al, Si, P, K, Ca, Ti, Mn, Fe). In order to differentiate the two Ca species (CaCO3 and CaO), we subtracted Ca in carbonates (CaCO3 contents obtained from Cinorg) from the total Ca content. The difference represents Ca in silicates (CaO). Ctot and N contents of the fine earth were obtained with EA-IRMS (Elemental Analyzer Isotope Ratio Mass Spectrometry; Flash HT Plus CNSOH, Thermo Scientific, Bremen, Germany) measuring 2 lab replicates.
We calculated three weathering indices for all soil profiles. All indices compare the molar ratios of easily weathered ('mobile') elements to more 'immobile' elements: This ratio was introduced by Dorn [73] on rock varnishes and it has proven useful in soils as well, as Ti is mostly immobile in soils (alternatively, Zr can be used as the immobile element, [74]). Here, we calculated this ratio using the total Ca. It therefore also reflects the dissolution of carbonates. An alternative version to this ratio is : That only reflects weathering of silicates. We also used the Chemical Index of Alteration [75]. It is normally used in for igneous rocks or siliceous soils. There, an increase in the CIA (relative accumulation of Al) is indicative of stronger silicate weathering.
The oxalate extractable Fe, Al, and Mn contents of the fine earth fraction were determined according to [76] with 2-3 replicates The Fe and Al measurements were conducted on a flame-AAS (atomic absorption spectrometer; ContrAA 700, Analytik Jena, Jena, Germany). The values presented are the averages of 2-3 replicates.

Mineralogical Analysis
The mineralogy of the samples was determined on randomly oriented powder specimens with X-ray diffraction (XRD) analysis. First, the samples were air-dried. The fine earth fraction (<2 mm) of the bulk soil was comminuted by hand using a mortar and sieved (400 µm). This fraction was subsequently Ca saturated. A representatively split aliquot of about 2 g was then milled in ethanol to a grain size below 20 µm with a McCrone micronizing mill and dried afterwards at 65 °C. For frontloading preparation, about 1 g of the powdered material was gently pressed in a sample holder for packing, sample-height adjustment, and forming a flat surface. Preferred orientation was avoided by using a blade for surface treatment [77]. For an improved identification of the clay mineralogy, appr. 1 g of the untreated <2 mm sample was sieved with a 20 µm sieve to enrich the clay minerals. From this <20 µm-fraction, we produced oriented specimens by carefully smearing a thin layer of the sample onto a glass slide. They were prepared for enhancement of the basal reflexes of layer silicates, thereby facilitating their identification. The textured samples were measured before and after treatment. The changes in the reflex positions in the XRD pattern by intercalation of different organic compounds (e.g., ethylene glycol) and after heating were used for the identification of expandable clay minerals.
X-ray diffraction measurements were made using a Bragg-Brentano X-ray diffractometer (D8 Advance, Bruker AXS, Karlsruhe, Germany) using CoKα (35 kV, 40 mA) radiation. The instrument worked with automatic beam optimization (equipped with an automatic theta compensating divergence slit and an automatic air scatter screen). The detector was an energy dispersive Lynx-Eye XE-T detector. The powder samples were step-scanned at room temperature from 2 to 80°2Theta (step width 0.02°2Theta, counting time 2 s per step).
The qualitative phase analysis was carried out with the software package DIFFRACplus (Bruker AXS, Karlsruhe, Germany). The phases were identified on the basis of the peak positions and relative intensities in the comparison to the PDF-2 database (International Centre for Diffraction Data). The quantitative amount of the mineral phases was determined with Rietveld analysis using the Rietveld program Profex/BGMN [78]. This full pattern-fitting method consists of the calculation of the X-ray diffraction pattern and its iterative adjustment to the measured diffractogram. In the refinements phase, specific parameters and the phase content were adapted to minimize the difference between the calculated and the measured X-ray diffractogram.

Meteoric 10 Be Sample Preparation and Analysis
Our method of extraction of 10 Be from soil is based on the methods described in [10,79]. A fine-milled soil sample is spiked with 9 Be standard solution and leached overnight twice using 16% HCl. The leachate undergoes a series of pH adjustments in order to separate metals, mainly iron, by precipitating as hydroxides. Thereafter, a fraction containing beryllium is isolated on ion-exchange resin Bio-Rad AG 50W-X8. A special grade of chemical reagents and water are used through all steps of sample preparation to ensure an ultra-low level of 10 B, which is an isobar of 10 Be.
Prior to the accelerator mass spectrometry (AMS) measurement, all the resulting BeO samples were mixed with Nb powder and pressed in a cathode insert. The measurements were performed with MILEA AMS system at ETH Zurich, Switzerland. The results were normalized to in-house standards S2007N and S2010N with nominal values of 10 Be/ 9 Be = 28.1 × 10 −12 and 10 Be/ 9 Be = 3.3 × 10 −12 and errors of 2.7% and 2.2%, respectively [80]. The inhouse standards are calibrated to a primary standard ICN 01-5-1 [81]. The final 10 Be concentrations were corrected to preparation blank, while the final errors include a preparation blank error and an error of the in-house AMS standards.

Calculation of Weathering and Erosion Rates
The progression of chemical weathering was determined by comparing the differences in the chemical composition of the weathered and unweathered soil material. Using an immobile element (Ti), we calculated element-specific mass balances. As the soil ages were known, it was also possible to calculate long-term weathering rates (mass balances) [80,81]. We calculated these for Na, Mg, Al, Si, P, K, Ca, Mn, and Fe. We will refer to the sum of these elements as the 'total' mass balance, henceforth.
Following [82], we first calculated the volume change between the unweathered parent material and the weathered material, or strain coefficient ε, for each weathered horizon: where ρp and ρw are the bulk densities (kg m −3 ) of the parent material (=lowermost horizon of each profile) and the weathered material, Ci,p (g kg −1 ) is the concentration of the immobile element (Ti was used) of the parent material, and Ci,w (g kg −1 ) is the concentration of the immobile element (Ti) of the weathered soil. The open-system mass transport function τj,w [82] is defined as: , , with Cj,w (g/kg) as the concentration of element j in the weathered material, and Cj,p for the unweathered parent material (=lowermost horizon) and the strain coefficient ɛw. Ci,w and Ci,p (g kg −1 ), the concentrations of the immobile element, were calculated with Ti. The change in mass for an element j (mj, flux(z) in kg m −2 ) for n soil layers is defined by Egli and Fitze [83]: The long-term erosion rates were calculated using meteoric 10 Be in the soils. When the surface ages are known, soil erosion rates can be calculated by comparing the measured concentrations of 10 Be in the soils with the expected concentrations based on their age, as demonstrated by [84][85][86].
With no erosion, the surface age on a soil is defined as: with λ = the decay constant of 10 Be (4.997 × 10 −7 y −1 ), q (atoms cm −2 year −1 ) = annual deposition rate of 10 Be (calculated according to [46,87,88]), and Nexp (atoms cm −2 ) = the expected 10 Be inventory in the soil profile (i.e., the accumulated 10 Be in the soil due to atmospheric deposition, assuming no erosion or other disturbances). The deposition rates of 10 Be will have to be estimated for a specific area, as they are not known precisely. However, Maejima et al. [84] have shown that they are mostly dependent on the amount of precipitation.
With soil erosion, Equation (9) is rewritten as: where C10Be (atoms g −1 ) = concentration of 10 Be in the top eroding horizons (ca. the top 20 cm), Esoil (cm year −1 ) = soil erosion rate, f = fine earth fraction, and ρ (g cm −3 ) = bulk density of the top horizons, λ = decay constant of 10 Be, N = inventory of 10 Be. As C10Be changes over time, it can be approximated by using the average between t = 0 and t, i.e., 0.5 × C10Be(today).
We also assume that the erosion occurs mostly in the top 20 cm of the soil profile. We calculate the 10 Be inventory N using [86]: for n horizons, and with the horizon thickness z, the bulk density ρ, the 10 Be concentration C and the fine earth fraction f of the weathered (w) material.
For comparison, we also used the method of Lal et al. [47]: KE is the first order rate constant for removal of soil, with NS = 10 Be inventory in the S layer (comprising O and A horizons), ND = 10 Be inventory in the D layer (below 20 cm), Q = flux of atmospheric 10 Be into topsoil (atoms cm −2 year −1 ), qa = flux of meteoric 10 Be (atoms cm −2 year −1 ), and λ = decay constant of 10 Be. The soil erosion rate Esoil (cm year −1 ) is then calculated using z0 = thickness of topsoil horizons (comprising O and A horizons) multiplied by KE.
In addition, the rates of soil formation were calculated. Most approaches are based on the assumption that the soil in question is in equilibrium, i.e., that the soil thickness remains stable because soil production and soil denudation are equal. As this is clearly not the case for such young and still intensely developing soils (see Figure 3), we used the changes in the profile thickness as an approximation for soil formation rates. The soil thickness function was first developed by [89] and further developed by Johnson et al. [90]. According to this model, the soil thickness (T) is a function of soil deepening factors (SD), upbuilding factors (U), and removals (R): We calculated the soil formation rates (F = SD + U) based on the changes in thickness (z) [91]:

Flow Accumulation Analysis
Digital surface models (DSM) were created using centimeter-resolution aerial drone images. The drone flights and the post-processing of the images were conducted by WWL Umweltplanung und Geoinformatik GbR [92]. Flow direction and flow accumulation analyses were conducted in ArcMap (ESRI). The flow direction analysis uses the DSM to compute the most likely flow path that water would take according to the steepness when moving from one grid cell to the next. The D-∞ ('D-infinity') method was used [93]. The flow accumulation shows the spatial distribution and intensity of the flow pathways, i.e., how much water is directed to each pathway. We used this output in combination with the microtopography measurements of Maier et al. [56] to better evaluate the results of our erosion analysis, as more bundled and therefore more intense flow paths have a higher potential for surface erosion than dispersed ones.

Uncertainty Calculation and Propagation
We defined the standard uncertainties for the different used parameters by calculating the standard deviation between field replicates (bulk density) and lab replicates (oxalate extractable Fe, Al, Mn; Corg, N contents) and by using the instrument measurement errors (elemental contents, mineralogy).
When using these datasets in further calculations (weathering indices, Corg, Corg stocks, C/N, element mass balances, 10 Be erosion rates), the standard uncertainties were accounted for by propagating them following the recommendations of the Guide of Expression of Uncertainty [94].

Profile Descriptions
The soils of the 110 a moraine were classified as Hyperskeletic Leptosol. The moraine was covered in boulders and cobbles of varying sizes and was vegetated with initial snowbed communities of the type Arabidion caeruleae and Saxifraga aizoides (Figures 2-4, Table  1). The 160 a moraine was classified as Hyperskeletic Leptosol and had already accumulated a humus-rich A horizon. The vegetation was dominated by Thlaspietum rotundifolii and Dryadetum octopetalae. The 4.9 ka soils were Calcaric Skeletic Cambisols which reached a depth of approximately 60 cm. They were vegetated by dense grass, mainly Caricetum ferruginei and Rhododendretum hirsuti with roots to 40-50 cm. The 13.5 ka moraine had Calcaric Skeletic Cambisols as well and was also densely vegetated. The dominating species were Seslerio-Caricetum sempervirentis and Caricetum ferruginei.
The two soils at the 30 a moraine Hyperskeletic Leptosols, one of which had accumulated an A horizon (Figures 2-4). The surface was covered by pioneer species (mainly Epilobium fleischeri) and initial grassland vegetation rich in Trifolium badium/pallescense and Luzula alpino-pilosa ( Table 1). The 160 a soils were also Hyperskeletic Leptosols with 10-20 cm thick A horizons. They were covered by a grassland vegetation of the type Poion aplinae. There were also creeping dwarf shrubs with a high coverage of Salix retusa present. The 3 ka soils were Skeletic Cambisols vegetated by grasses (mainly Agrostio rupestris-Sempervivetum montani and Carici sempervirentis). The 10 ka soils were Entic Podzols with thick, humus-rich topsoils that reached a depth of 30-40 cm. The surface was largely covered by shrubs (Rhododendro ferruginei-Vaccinietum) and some grasses (Geo montani-Nardetum) in between and underneath the shrubs.

Soil Properties
The bulk density is similar in both locations. At Klausenpass, it varied from 0.7 to 1.8 g cm −3 in the younger soils and from 0.3 to 1.9 g cm −3 in the older soils. The soils at Sustenpass display clear age trends, ranging from 1.2 to1.7 g cm −3 in the 30 a, 160 a, and 3 ka soils and 0.7-1.4 in the 10 ka soils.
The soils had a high skeleton content throughout all sites, often over 50% (Table 1). At Klausenpass, it decreased with age, most strongly in the topsoil. There is also a strong increase in the skeleton fraction with depth. At Sustenpass, the skeleton content varied more strongly from moraine to moraine and the depth trend was less clear compared to Klausenpass.
The highest organic matter (OM) contents (measured as loss on ignition, Table 1) were found in the topsoil of the oldest soils, ranging from 19% to42% at Klausenpass and 22% to38% at Sustenpass. The OM in the young topsoils at Klausenpass (110 a, 160 a: 1.8-8.9%) was distinctly smaller compared to the young soils at Sustenpass (30 a, 160 a: 0.8-38.4%). The differences between the profiles of one moraine is also higher at Sustenpass than at Klausenpass.
The pH at Klausenpass was mostly in the neutral to weakly alkaline range (pH 6.7-7.8) except for the top horizons in the older soils where the pH was below 6. At Sustenpass, the pH was in the neutral to strongly acidic range (pH 3.4-6.7).

Klausenpass
The elemental composition is presented in Table 2. The calculated strain coefficients and open-system mass transport functions for each Klausenpass soil profile are shown in Tables 3 and 4. The element mass fluxes of each moraine are summarized in Figure 5a. High losses were recorded for the mobile elements Mg and Ca, as well as P, while lessmobile elements passively accumulated. The total losses are highest at the 13.5 ka soil (−152 ± 49 kg m −2 ), Figure 5b. When including the soil age to calculate rates of chemical weathering, the rates do not change very strongly over time (+28 ± 2 to −7 ± 1 m −2 kyear −1 in the youngest soils to +1 ± 0 to −11 ± 4 m −2 kyear −1 in the oldest soils).
The weathering indices (Figure 7a) show strong alterations in the uppermost 20 cm of the soils. The variability between the soils at each moraine becomes slightly larger with age.
Both the (Ca + K)/Ti and the (Na + K)/Ti ratios (Figure 7b) show similar patterns among the Sustenpass soils. The soils of each moraine show a homogenous profile and the soils are clearly grouped by age. The Chemical Index of Alteration (CIA) mostly increases with soil age. The variability is however high among the top horizons, particularly at the 160 a moraine site.
The ratio of Feo/Fetot only increased slowly; fluctuating at around 5-6 ± 0.5% in the first 3 ka of soil development. Between 3 ka and 10 ka, the ratio sharply increased to 29 ± 1.9-35 ± 1.2% in the topsoil (Figure 8b). Similarly, the Alo/Altot ratio changed little during the first 3 ka of soil development (0.4-2.3 ± 0.2% in the topsoil) and then increased to 9-10% in the 10 ka soils. The Mno/Mntot ratio increased too with age, although the variability is very high.

Flow Accumulation
The GIS analysis to model the flow accumulation ( Figure 11) produced a homogeneous pattern on the Klausenpass slopes. There, the surface runoff is dispersed with a low intensity. The Sustenpass slopes have a runoff pattern that indicates channeling and more intense flow accumulation (dark blue) of the surface water runoff. The same patterns also emerged in the younger moraines.    [86].

Soil Development in Calcareous Parent Material
Weathering in calcareous soils is primarily characterized by carbonate leaching [95]. The values of carbonate leaching in the study area generally fit to the range of observations from other calcareous soils in the Alps [24,96]. The weathering front seems to advance in the soil column 'stepwise'. This leads to a divergence of top-and subsoil properties, where the topsoil is strongly altered, while the subsoil has changed very little or might have accumulated elements that were leached above (such as Ca). We observed this in the pH (drop of pH to <6 in the topsoil), mineralogy (dissolution of calcite), weathering indices and τ values (high element losses in the topsoil) as well as the oxalate extractable fractions (accumulation of Feo, Alo and Mno). Soil texture data (Figures 12 and 13, obtained from [54] at the same moraines) showed that the topsoils at this location were altered very strongly with time, shifting from sandy loam to a silty clay loam texture. Despite 110 and 160 years of rapid soil carbonate leaching, the profiles of the youngest two moraines showed almost no carbonate losses. We calculated Ca mass balances between +2 and −2 kg m −2 which would translate to +6 to −6 kg m −2 CaCO3. In comparison, [24] calculated losses of 15-30 kg m −2 in 193 a soils. The 4.9 and 13.5 ka soils had mass balances between +7 and −154 kg m −2 . Indeed, the rates of weathering (based in mass balances) remained roughly the same over the course of 13,500 years. Because the mass balances were calculated by comparing the soil column with the lowermost horizon, soil erosion may account for the small carbonate losses in the young soils. As soil erosion removes the weathered residue of the soil and exposes fresh material, it effectively leads to a rejuvenation of the soil.

Soil Development in Siliceous Parent Material
Soil development in siliceous parent materials is characterized by soil acidification and weathering of primary minerals. With weathering, elements are either rearranged into secondary minerals (i.e., clays) and remain in the soil profile [97], are used as nutrients by the vegetation, or lost through leaching into the groundwater. Overall, the progression of chemical weathering and the accumulation of soil organic carbon at Sustenpass is in agreement with other observations of mountainous soils of similar age and parent material [37,[98][99][100]. Siliceous soils from other regions of the world show similar weathering and carbon accumulation rates [100][101][102][103]. The chemical weathering rates were high in the <160 a soils (losses of up to 69 kg m −2 kyear −1 ) and decreased thereafter rapidly. These findings match the hydrological data on the same moraines by Maier et al. [51] and Hartmann et al [104]. The saturated hydraulic conductivity (Ksat) decreased significantly with soil age [51], leading to a less efficient draining of the soil column in the older moraines. Blue dye experiments by Hartmann et al. [104] also showed that the flow paths become more heterogeneous with age and that the macropore flow (via root channels) becomes more important in the older moraines. This concurs with the shift in soil texture from rather coarse, sandy soils to more silty compositions. The overall texture at Sustenpass is coarser compared to Klausenpass.
The elemental composition (Table 2) revealed some irregularities in the parent material composition. We found the Fe contents increased strongly with soil age and with depth. The concentrations of Fe in the oldest soil are roughly twice as high as in the youngest soils. This can be linked to the high chlorite concentrations in the 3 ka and 10 ka soils ( Figure 9). As pedogenic chlorite forms very slowly, this must be primary chlorite, which suggests a fluctuation in the mineral composition of the glacial sediment. The distribution of the chlorite in the profiles of the 10 ka and 3 ka moraines is not uniform either. There is an accumulation of chlorite in lower horizons and an apparent depletion in the topsoil.

Soil Production vs. Erosion
The most challenging part in calculating erosion rates using meteoric 10 Be is the estimation of the meteoric 10 Be deposition rates. We applied a variety of approaches given in the literature [48,87,105]. Independent of the approach used, we found the expected inverse relationship between soil age and erosion at both locations. In addition, lower rates of long-term soil loss at the calcareous site (Klausenpass) were found when compared to the siliceous site (Sustenpass). This is in contrast to findings of the short-term erosion rates (averaged over decades). Short-term erosion rates using fallout plutonium showed similar erosion rates at both locations [53]. Overall, the long-term erosion rates measured at these two proglacial areas are in a similar range to findings of [106] for soils of the eastern Swiss Alps (up to 0.44 t ha −1 year −1 of soil loss). Interestingly, ref. [107] found denudation rates that are twice as high at formerly glaciated slopes, which are now very active geomorphologically (760-2100 mm kyear −1 ). They found comparable rates to the Sustenpass soils at slopes which were not influenced by glaciers (60-560 mm kyear −1 ). In Southern Alpine soils (New Zealand), a range of 1.2-14 t ha −1 year −1 soil erosion was found, which is on the higher end compared to the 0.1-10 ha −1 year −1 for global soils [108].
The differences in the long-term erosion rates between the calcareous and the siliceous chronosequence are surprising. The sites at Klausenpass have a steeper slope and the material is siltier. These conditions usually would give rise to higher erosion rates [109,110]. Ref. [54] showed that these soils all have similar root tissue densities in the topsoil, from which we can derive that the soil surfaces experience similar stabilization by the plant roots. The SOM, which also greatly influences soil erodibility by acting as a binding agent for aggregation and protective layer for the mineral soil underneath [111,112], is higher at the 4.9 ka soil (187 g kg −1 vs. 37 g kg −1 at Sustenpass in the top horizon), and lower in the 13.5 ka soil (424 g kg −1 vs. 439 g kg −1 ). The most important reason why the Klausenpass soils have lower erosion rates than expected may be related to the small-scale topography. Surface roughness determines the pathways of the surface water runoff and can, thus, have a large influence on the erodibility of a slope [52,[113][114][115]. A flow direction and accumulation analysis revealed that the Sustenpass and Klausenpass slopes differ in their flow patterns ( Figure 11). The Sustenpass slopes produced patterns that rather lead to a channeling of the surface water runoff, indicating a higher terrain roughness compared to Klausenpass, where the surface runoff is more spatially homogeneous and less intense. This is especially pronounced at the 10 ka moraine (Sustenpass), where the Rhododendron bushes promoted the formation of small-scale channels on the soil surface. The Klausenpass slopes, covered with grasses and having only few boulders interspersed on the surface (Figure 2), were able to retain a smoother surface. This GIS analysis has some limitations, as it is based on a digital surface model (DSM) which includes the volume and shapes of the vegetation, as it was produced from drone images, and not by using LiDAR (light detection and ranging, 'laser scanning'). However, we can support this hypothesis with the results of a detailed investigation of the terrain microtopography on the same slopes by Maier and van Meerveld [56]. The authors calculated a tortuosity index based on high-resolution transect measurements of the surface roughness of the moraine slopes. The tortuosity index at Sustenpass (10 ka: 0.64-1.14; 3 ka: 0.33-0.65) is indeed higher than at Klausenpass (13.5 ka: 0.25-0.30; 4.9 ka: 0.26-0.40), supporting the argument that the microtopography may be the factor causing lower erosion rates at Klausenpass [114,116].
We estimated the soil production rates (Equation (1)) using soil formation, weathering and erosion rates ( Table 7). The soil production rates were considerably higher at the siliceous soils compared to the calcareous soils. Furthermore, the Klausenpass soils show clear signs of regressive soil formation. The denudation at the 110 a moraine greatly outstrips the soil production, preventing it from building up an organic top horizon ( Figures  3 and 4). Such a delay in the soil development was also mirrored in the Corg stocks, which increased at a much smaller pace compared to the Sustenpass soils ( Figure 6). It is further accompanied by a slower development of the vegetation: the Klausenpass moraines took more time to be covered in vegetation and to inhabit more complex species (i.e., grasses and woody species) compared to Sustenpass [54].
The soil production rates of the siliceous soils showed a strong decline with increasing age, which was to be expected. The 30 a moraine is estimated to have soil production rates of at least 34 t ha −1 year −1 while the 10 ka moraine has 3.5 t ha −1 year −1 . This is three times as high as other Alpine soils (ca 0.4-10 t ha −1 year −1 [107,117]). It is also high compared to global soil production and formation rates measured in soils and river sediments which are in the range of 1-30 t ha −1 year −1 [8,108,[118][119][120]. Table 7. Estimation of soil production rates. Profile thickening (Figure 3) was used for soil formation (Fsoil), soil erosion rates from 239+240 Pu and 10 Be were used as denudation rates (Dsoil), and the total mass balances served as chemical weathering rates (Wsoil). The meteoric 10 Be erosion rates (erosion = negative values) were used for the 3 ka, 4.9 ka, 10 ka, and 13.5 ka soils. The denudation rates of the young soils (30 a, 110 a, 160 a) were derived from the short-term ( 239+240 Pu) erosion rates.  1 Short-term erosion rate (erosion = negative values) derived from 239+240 Plutonium. It is the mean yearly soil loss over the past 54 years (data: [53]). 2 Equation (1): Psoil = Fsoil + Dsoil = Fsoil + Wsoil + Esoil. * These are minimum values as the depth of the soil pits did not reach the bottom of the BC horizon (see Figure 3). ** Approximated using the erosion rate of S-B.

Interaction with Vegetation and Hydrology
In many areas of the Alps, the snow cover lasts for roughly six months or more and in summer, the day and night temperatures on the ground can vary strongly, as the shallow and debris-rich soils in the proglacial areas offer little moisture and little vegetation to buffer these fluctuations [121]. Therefore, in the earlier stages of soil development, the habitability for plants would be primarily given by the physical properties of the soil substrate. A loamy soil texture would have a higher water retention [122] but may also lead to a higher erodibility [123], which would slow down if not stop a deepening of the soil profile, thus inhibiting the progression of soil development. If pioneer plant species manage to inhabit the substrate, they will stabilize its surface with their roots and enhance soil aggregation by accumulating organic matter in the soils (via litter and root exudates).
These improved conditions for the soils would enable new species to inhabit these soils, thus setting in motion a succession of plant species and communities as well as a progressive soil development [124].
Overall, the texture at Klausenpass is consistently finer compared to the Sustenpass location, which remains low in clay-sized particles (Figures 12 and 13). Due to the dissolution of carbonates, the clay and oxide residues accumulate passively giving rise to finegrained soils. Concurrent to the decrease in particle sizes, the water retention also increases with age at both sites [54]. What may have reinforced this effect of decreasing particle sizes in the topsoil, are loess inputs, which are common in the Alps [125].
Some data seem to indicate a noticeable influence of the vegetation on soil development. For example, the SOC stocks of the siliceous chronosequence increased consistently faster under the more functionally rich vegetation, while we could only see differences between some of the calcareous soil profiles in SOC accumulation ( Figure 6). The total mass balances (Figure 5b) of the older soils, showed higher losses in the low FRI soils at both locations.
The weathering indices ( Figure 7) showed stronger weathering in the low FRI profile at 3 ka. At the 10 ka moraine, this relationship is only detectable with the CIA. The oxalateextractable fractions do not exhibit a clear and consistent pattern ( Figure 8). Only the 10 ka soils show a difference between the two FRI types. Figure 13. Particle size distribution on top of the USGS soil texture classification for each depth increment and for all depths together (same dataset as Figure 12 [55], plots were drawn in R [126]). The blue dots represent the Klausenpass and the red dots the Sustenpass soils. Light hues represent younger and the dark hues the older ages. Per age/moraine and depth, there are three replicates. (Data reproduced with permission of [55]).It therefore seems that the functional diversity of the vegetation has an effect on soil formation. The variability of functional diversity gives rise to a different litter input, root activity, and the microbial activity which in turn influence soil weathering and stabilization processes. How the FRI precisely relates to soil formation is not clear yet. The FRI, however, might be a useful parameter to better describe the natural variability of soils. Both chronosequences demonstrate that the FRI only started to have an effect on soil evolution when the vegetation was already well established, i.e., sometime between 160 and 3000 years of soil development. Before that, the physical site properties and the climate (mainly the precipitation, [127]) are the stronger drivers. The small-scale topography of the moraines (indentations offer protection from wind and erosion and have a higher retention of moisture [128]), inhomogeneity in the mineralogical and chemical composition of the parent material, and soil texture drive the soil formation and are therefore responsible for the high spatial variability of mountainous soils. These properties also influence the hydrological conditions of the soils [129] and have therefore a strong control on the early stages of vegetation settlement [90,130,131,39].
At later stages, however, the vegetation influences the site properties by increasing chemical weathering, reducing surface erosion and altering the hydrological properties of the substrate as well as the vertical water transport by causing heterogeneous infiltration patterns, leading to an increase in preferential flow paths [51,55,104]. The functional richness of the vegetation affects soil parameters that are directly related to the vegetation (e.g., SOC).

Conclusions
Calcareous and siliceous soil evolution in two proglacial areas in the Swiss Alps follow, not surprisingly, different traits of chemical weathering. Chemical weathering of the calcareous soils was characterized by high CaCO3 losses over time, whereas the siliceous soils lost predominantly Si, Fe, and Al. By using meteoric 10 Be, we found higher long-term erosion rates along the siliceous chronosequence than along the calcareous chronosequence. A flow direction and flow accumulation analysis revealed that the lower terrain roughness at Klausenpass likely led to a lower erodibility. The temporal evolution of the erosion rates was completed by a combined approach of short-term ( 239+240 Pu) and longterm ( 10 Be) rates. The soil erosion rates counterbalanced the soil production rates at the calcareous location, delaying the development of both soil and vegetation in comparison to the siliceous soils. Overall, we found that the relative importance of the parent material decreased with increasing soil age and vegetation coverage. The vegetation is an important driver of soil development as it accelerates chemical weathering, increases surface stabilization, and changes the hydrological properties. Differences in the functional diversity are, however, only detectable with parameters that strongly depend on the vegetation such as SOC.

Data Availability Statement:
The data that support the findings of this study will be made available in Pangaea (https://pangaea.de/, accessed on: 12 February 2022).