Root density drives aggregate stability of soils of different moraine ages in the Swiss Alps

The stability of hillslopes is an essential ecosystem service, especially in alpine regions with soils prone to erosion. One key variable controlling hillslope stability is soil aggregate stability. We aimed at identifying dominant controls of vegetation parameters on aggregate stability and analysed their importance for soil aggregate stability during landscape development. We quantified the aggregate stability coefficient (ASC) and measured plant cover, diversity, root mass and root length, density (RMD, RLD) along two chronosequences with contrasting bedrocks (siliceous, calcareous) in the Swiss Alps. We found that ASC developed slower along the calcareous chronosequence. Furthermore, we observed a significant positive effect of vegetation cover and diversity on ASC that was mediated via root density. These relationships developed in a time-depended manner: At young terrain ages, vegetation parameters had a strong effect on aggregate stability compared to older stages. Moreover, RLD was the most powerful predictor of ASC on young terrain, whereas on older moraines RMD became more important. We highlight that root density plays a major role in governing ASC for soils differing in moraine ages. The changing importances of RLD and RMD for ASC development suggest different mechanistic linkages between vegetation and hillsope stability during landscape development.

areas, soil formation and development are slow and soil degradation is common [Food and Agriculture Organization of the United Nations (FAO) 2015]. Mountainous landscapes consist largely of steep hillslopes that are prone to erosion, especially above the timberline (Pohl et al. 2009). This high risk of soil erosion is further enhanced by climate change, which is increasing the frequency and intensity of rainfall events in the alpine region (Anache et al. 2018;Nearing et al. 2004) and accelerating glacial retreat, leaving fresh and unconsolidated glacial debris behind (Gobiet et al. 2014). Together with anthropogenic pressure induced by tourism and land-use changes, this leads to an increased risk of natural hazards provoked by soil erosion such as landslides and debris flows (Gobiet et al. 2014;Hudek et al. 2017b). Hence, understanding the development of hillslope stability and identifying its drivers in alpine habitats is key to avoiding further soil degradation and connected risks to human safety in these areas (Hudek et al. 2017b).
Vegetation structures are known to improve hillslope stability and are frequently used in ecological restoration programs (e.g. Bardgett et al. 2014;Liu et al. 2018;Stokes et al. 2009Stokes et al. , 2014Vannoppen et al. 2015). The aboveground compartments of vegetation enhance hillslope stability by reducing the kinetic energy of rainfall and surface runoff through canopy structure and leaf traits, surface runoff through an increased surface roughness, soil water content via transpiration, and sediment transport since vegetation functions as a sediment trap (Geißler et al. 2012;Goebes et al. 2015;Kervroëdan et al. 2019Kervroëdan et al. , 2021Poesen 2018;Vannoppen et al. 2015). Aboveground vegetation characteristics, such as cover and diversity, are known to play a key role for hillslope stability. That is because a densely vegetated plant community with different growth forms and functional traits will be more effective for maintaining hillslope stability than a vegetation cover with few species (Körner und Spehn 2002). The belowground root system of plants controls hillslope stability by increasing the infiltration capacity via root channels, modifying soil texture, organic content as well as chemical composition, and affecting soil shear strength (Amezketa 1999;Burylo et al. 2012;Stokes et al. 2014;Vannoppen et al. 2015). In the context of hillslope stability, belowground plant characteristics have been studied less than the aboveground compartments and are referred to as the "hidden half" (Bardgett et al. 2014;Eshel und Beeckman 2013). In most cases, high aboveground vegetation cover is positively correlated with root biomass (Martin et al. 2010;Pohl et al. 2009). However, since roots have a more direct impact on many soil stabilizing processes than vegetation cover, it is of relevance to quantify root characteristics. As they are relatively easy to assess, Root Mass Density (RMD = root dry mass per soil volume) and Root Length Density (RLD = root length per soil volume) are normally used in studies dealing with hillslope stability (Pohl et al. 2012;Erktan et al. 2015;Hudek et al. 2017b). However, in alpine habitats, there are only few studies dealing with roots, and even less is known about their role in maintaining ecosystem services above the timberline (Bast et al. 2016;Hudek et al. 2017a, b;Pohl et al. 2009Pohl et al. , 2011. One key variable reflecting hillslope stability is soil aggregate stability. Soil aggregates are groups of soil particles that cohere due to chemical and physical forces (Amezketa 1999). An improved soil aggregate stability will increase soil permeability, lead to an increased infiltration rate, and reduce surface runoff (Vannoppen et al. 2015). Stable aggregates withstand slaking when they are wetted rapidly as they can tolerate the pressure due to swelling (Tisdall and Oades 1982). Consequently, aggregate stability is measured best by testing their water stability, which simulates the natural forces during precipitation. The classical method to determine soil aggregate stability is to measure the aggregate size distribution with a wetsieving approach and then to relate the specific weight of large, stable aggregates to the same weight of small aggregates (Tisdall and Oades 1982). Alternatively, for alpine habitats, a new method [Aggregate Stability Coefficient-(ASC)] has been developed for soils with high stone contents and reducing work effort to a minimum (Bast et al. 2015).
Soil aggregate stability is connected to various biotic (soil organic matter, plant roots, soil fauna, and microorganisms), abiotic (clay minerals, sesquioxides, exchangeable cations), and environmental parameters (soil temperature, moisture and lithology) (Amezketa 1999;Márquez et al. 2004, Pérès et al. 2013. For example, plant litter increases the activity of soil microbes (Blume et al. 2010), which release enzymes that mineralise high molecular weight compounds and polysaccharides (Amezketa 1999). Both factors contribute to soil aggregation. Furthermore, to support microorganisms at the rhizosphere, increase nutrient availability, and reduce element toxicity, roots release exsudates, which positively affect soil aggregate stability (Bardgett et al. 2014). This is because exudates consist of polysaccharides and proteins, which bond mineral particles together. The flux of primary metabolites through root exsudation is mostly located at the root tip (Canarini et al. 2019) and thus, the amount of fine roots present is of special significance for soil aggregate stability. Another actor at the rhizosphere level secreting exsudates and thereby influencing soil aggregate stability is the hyphal network created by mycorrhizal fungi. As plant species differ in the amounts as well as in the composition of their exsudates and because the occurrence of mycorrhizal fungi is highly plant species-dependent, species composition of plant communities may have a strong impact on soil aggregate stability (Bardgett et al. 2014). One prominent abiotic factor controlling the formation of soil aggregates is lithology (Duan et al. 2021) as it is affecting soil mineralogy and texture as crucial factors for soil aggregation (Bronick and Lal 2005). In soils developed on calcareous bedrock, high clay and calcium concentrations drive the stabilisation of organic matter and the formation of aggregates. In more sandy or acidic soils, the absence of alkaline cations and the coarse soil texture impede the formation of coarse aggregates (Bronick and Lal 2005;Yang et al. 2020).
The most relevant previous studies dealing with the impact of vegetation dynamics on hillslope functioning above the timberline have mainly dealt with narrow environmental gradients (Bast et al. 2014(Bast et al. , 2016Martin et al. 2010;Pohl et al. 2009Pohl et al. , 2012. Therefore, little is known about the long-term development of hillslope stability during landscape evolution. We, especially, lack a holistic understanding of its influencing factors, which applies to the underlying interactions, mechanisms, and their relation to the resulting pattern, also in terms of scale linkages. The expanding terrain of glacier forelands is a field laboratory for studying the development of such features as a response to vegetation dynamics. In these habitats, former glacier positions can be identified by moraines representing distinct ages of substrate exposure (e.g. Egli et al. 2001;Musso et al. 2019Musso et al. , 2020. This space for time substitution, called chronosequence approach, enables scientists to test hypotheses related to terrain age and ecologists have frequently used it to study plant succession or soil development. In this study, we used two proglacial chronosequences with contrasting bedrocks (siliceous versus calcareous) to study the development of soil aggregate stability for soils differing in moraine ages and to reveal its interactions with above-and belowground vegetation parameters. Chiefly, we tested the following hypotheses ( Fig. 1): We expect that (1) vegetation cover and diversity are positively correlated to RMD and RLD because of density and complementarity effects, respectively (e.g. Ravenek et al. 2014). In addition, we assume that (2) soil aggregate stability correlates positively to RMD and RLD. Consequently, if hypothesis 1 and 2 are true, (3) vegetation cover and diversity should be positively correlated with soil aggregate stability. Because Pohl et al. (2012) found that the impact of vegetation on aggregate stability was strongest on slopes with high disturbance, we assume that (4) the relationships described in hypothesis 1-3 vary along the chronosequences, covering a gradient of young moraines with high disturbance regime to old moraines with lower disturbance levels. Specifically, we expect that vegetation parameters play a more important role in governing soil aggregate stability on young than on old moraines. Finally, bedrock is one of the key factors determining rates and speed of vegetation dynamics (Vítovcová et al. 2021), whereby succession is thought to be decelerated in calcareous glacier forelands (Lüdi 1958;Oechslin 1935). Because the vegetation is critical for determining aggregate stability, we assume (5) a slower development of the latter variable at the calcareous glacier foreland.

Study sites
The study was conducted in two proglacial areas of Central Switzerland ( Fig. 2; Table 1). The first site is located west of the Susten Pass, in the foreland of the Stein Glacier, Canton Berne (47°43′N, 8°25′E). The local bedrock at the Stein Glacier foreland is covered by glacial till on pre-Mesozoic silicate parent material, comprised of metamorphosed metagranitoids, gneisses, and amphibolites (Musso et al. 2019;Schimmelpfennig et al. 2014). The second site is situated 33 km northeast of Susten Pass, close to the Klausen Pass, in the glacier foreland of the Griess Glacier, Canton Uri (46°50′N, 8°49′E). Here, the geological substrates are characterized by calcareous material (Globigerinen limestone), schists (Pectinien schist), marl (Globigerinen marl), and some quartzites (Musso et al. 2019;Oechslin 1935). The soils of both chronosequences developed from glacial till; an overview of the soil conditions at the study sites is given in Musso et al. (2019).
The two sites are located above the timberline at an elevation between 1800 and 2300 m above sea level.  46°50′N, 8°49′E). The coloured areas represent the studied moraines forming two chronosequences. Images were taken from DigitalGlobe/Google EarthTMmapping service According to the closest weather station, at the Grimsel Hospiz (1980 m a.s.l.; 18 km away from Stein Glacier, and 49 km from Griess Glacier), the annual air temperature is 1.9 °C, the mean annual precipitation is 1856 mm year −1 and snowfall events are common even during summer (Meteo Swiss 2016). At both sites, we established one chronosequence consisting of four moraines each. At the Stein Glacier site, the terrain ages of the moraines were 30 a (a = years), 160 a (Little Ice Age), 3 ka (ka = 1000s of years), and 10 ka. At Griess Glacier, the ages were 110 a, 160 a (Little Ice Age), 4.9 ka, and 13.5 ka (for details of terrain age determination see Musso et al. 2019).
The plant species composition differed strongly between the study sites as well as across the moraines. At the Stein Glacier chronosequence, the vegetation was clustered into four distinct successional stages. At the Griess Glacier site, the two younger and the two older moraines had a similar species composition ( Fig. S1a). In terms of plant growth forms, vegetation was mainly dominated by forbs and grasses (Fig. S1b). Typical colonizers of the youngest two moraines on both study areas were creeping dwarf shrubs, e.g. Dryas octopetala, Salix retusa, and pioneer species, e.g. Saxifraga aizoides, Poa alpina. The vegetation of the old moraines was characterized by alpine grassland species, e.g. Nardus stricta, Trifolium pratensis ssp. nivale (Stein Glacier) and Carex ferruginea, Sesleria albicans (Griess Glacier; Table S1).

Sampling design
On all moraines, three plots (4 m × 6 m each) were installed, thus establishing a total number of 24 plots. Detailed information on the single plots is provided in Table S1. Plot selection was done based on a structural vegetation complexity measure. Briefly, structural vegetation complexity was defined as an index based on vegetation cover and functional diversity data (Maier et al. 2020). Functional diversity was calculated based on the following eight traits characterizing species along the main axes of plant performance (Garnier et al. 2016): specific leaf area, nitrogen content, leaf dry matter content, Raunkiaeŕs life form, seed mass, clonal growth organ, root type, and stem growth form. For the plot selection, we conducted a vegetation mapping differentiating between vegetation units classified according to the characteristic, recurrent combination of species. In each unit, we recorded all vascular plant species with proportion cover and calculated the vegetation complexity index. On every moraine, the three plots were placed within the surface units with the lowest, intermediate, and highest vegetation complexity (Maier et al. 2020).
In July 2018 (Stein Glacier) and July 2019 (Griess Glacier), a total of 144 soil cores (six per plot) were extracted. Since the study was performed in the scope of the larger HILLSCAPE project including rainfall simulations (see Lustenberger 2019), the plots needed to be undisturbed, and hence sampling took place outside of the plots but as close to them as possible ("off-plot locations", Fig. 3). To ensure that the data gathered was representative for the original HILLSCAPE plots, we visually selected off-plot locations that represented the two dominant vegetation strata of each plot, characterized by growth forms and species identity of its dominant plants. The number of soil cores taken per stratum was determined based on the percentage proportion of the plot area covered by the same stratum. For instance, assuming that a hypothetical plot would have two-thirds of the cover dominated by shrubs (such as Rhododendron ferrugineum and Vaccinium myrtillus), and one third dominated by grassland species (e.g. Nardus stricta and Calamagrostis villosa), four soil cores would be taken in the stratum dominated by shrubs and two soil cores in the stratum dominated by the grassland species in the corresponding off-plot locations. The cores were taken at six randomly selected sampling points that fitted the predefined vegetation strata.

Aboveground characteristics
Vegetation surveys (1 m × 1 m each) were conducted at two of the sampling points of each off-plot location (each one per vegetation stratum; in total 48 vegetation surveys). The vegetation surveys were conducted within 1 week in mid-July 2018 (Stein Glacier) and 2019 (Griess Glacier), recording the percentage of plant coverage by visual estimation. The species nomenclature followed the latest edition of Flora Helvetica (Lauber et al. 2018). We calculated the Shannon diversity index (H) as a basic diversity parameter (Pielou 1966). The vegetation cover was determined by summing up all single species coverages in each vegetation survey.

Soil cores
The soil samples were taken within three consecutive days without rainfall in July 2018 (Stein Glacier) and 2019 (Griess Glacier) (Fig. 3). The soil cores were extracted carefully, ensuring that at least the uppermost 10 cm of the soil below the humus layer (O horizon) was captured. To extract the soil cores a HUMAX probe and a sledgehammer were used (for a detailed description see Bast et al. 2015). When hammering the probe into the ground the soil core was  pushed directly into a plastic cylinder (25 cm × 5 cm), thereby minimizing disturbance of the sample. Due to large quantities of stones at several sites sampling was limited to 10 cm depth. During the coring process, the samples were carefully placed into a box with the topsoil part up and kept in the shade. Afterwards, the samples were stored in a refrigerator at approximately 4 °C. 1 week later they were reloaded into a padded box and transported by car to the laboratory. There they were kept at 4 °C until soil aggregate stability measurements were completed.

Soil aggregate stability
To quantify soil aggregate stability we calculated the ASC as proposed by Bast et al. (2015). The method has been previously shown to be successfully applied in alpine environments and was designed for coarsegrained soil with gravel (> 2 mm) content above 50% (Bast et al. 2015(Bast et al. , 2016. For measuring ASC, the soil samples (soil cores) were trimmed to a length of 10 cm, excluding the litter layer above the topsoil. Subsequently, the samples were placed on a sieve with a mesh width of 20 mm standing in a drainable bucket. The bucket was filled with tap water until the sample was covered. After 5 min, a ball valve at the bottom of the bucket was opened to release the water. The soil material trapped within the sieve with a mesh size of 20 mm as well as the material that passed through the sieve was flushed into two separate bowls. The two fractions were cleaned from roots and other plant material. Stones larger than 2 mm were removed from both fractions to determine soil volume. The different fractions were dried at 105 °C for at least 24 h and then weighed separately. The different fractions were used to calculate the ASC, which is defined as: represents the fraction of soil remaining on the 20 mm mesh sized sieve, m stones [g] equals stones larger than 20 mm, and m total [g] means the entire soil sample before water treatment (Bast et al. 2015). The ASC ranges from zero to one, zero suggesting a completely disaggregated sample and one implying a stable, fully aggregated sample.

Root mass and root length density
The roots were sorted during the measurement of soil aggregate stability and then stored in a freezer at − 15 °C until further processing. The roots were then cleaned thoroughly from any remaining soil and prepared for scanning. For this purpose, according to their diameter roots were separated into fine (< 1 mm), and coarse roots (> 1 mm). The roots were put into water-filled petri dishes and scanned with a flatbed scanner. The scanned petri dishes were put into a compartment drier at 60 °C. After the samples were completely dried the weight of the roots was determined.
To reduce the workload, fine roots were subsampled when more than five scans would have been necessary to scan all fine roots of one soil sample. For subsampling, over the mesh grid of a sieve (diameter 20 cm; mesh width 1 mm), the roots were uniformly distributed on a paper towel with a defined size. With a round hollow puncher, five randomly selected pieces of the towel were punched out covering 1% of the area of the sieve mesh grid. The roots attached to the five punched-out pieces were washed into petri dishes and processed analogously to the other roots. The roots still attached to the rest of the paper towel were put into a compartment drier and weighed after the drying process. In this way the root length of the entire sample could be estimated using the dry weight and root length of the scanned subsample as well as the dry weight of the remaining sample.
The root length was measured using WinRHIZO (Régent Instruments Inc 2013). Images were analyzed at 800 dpi. RMD was determined by dividing the dry weight of the sample by the soil volume of the corresponding sample. RLD was obtained by dividing the root length by the soil volume. Importantly, we used all root fractions (fine and coarse roots) to calculate RMD and RLD. Because coarse roots considerably affect root biomass, RMD was interpreted as a measure reflecting the accumulation of soil organic matter and related decomposition processes. RLD is instead predominantly driven by fine root and therefore RLD is more related to direct interactions at the root-soil interface, such as the secretion of exudates. To determine the soil volume of each sample, the volume of gravel (> 2 mm) was subtracted from the volume of the corresponding soil cylinder. The volume of stones was determined using the water displacement method.

Statistical analyses
For testing the study hypotheses, we omputed Linear Mixed-Effect Models (LMERs) (Bates et al. 2007). We used "plot", nested within "moraine" as random effects to account for the spatial nestedness of our data. To detect if the effects on the response varied depending on the terrain age, the interaction between terrain age and a relevant explanatory variable was inspected. Terrain age was implemented as a continuous predictor and was transformed with the decadic logarithm due to the large scaling difference between the younger and older moraines. The study site (Stein Glacier, Griess Glacier) was used as a covariable in the models. To improve model quality, all continuous predictor variables (terrain age, RMD, RLD, Shannon index, vegetation cover) were scaled. The models were built according to the three study hypotheses. As vegetation characteristics (Shannon index, cover) and root density (RMD, RLD) were correlated significantly (ρ > 0.6, p < 0.05) we set up a single model of every parameter. Further details on the single models are listed in Table 2. All models were checked for outliers and assessed in terms of homogeneous residual distribution and approximate normality (Fig. S2). We calculated post-hoc multiple comparisons of group-specific means to get p-values for the linear mixed effect models (Kuznetsova et al. 2017). Collinearity was tested with variance inflation factors (Zuur et al. 2010). Marginal R-squared values, which are interpreted as the variance explained by the fixed effects, were calculated after the method of Nakagawa and Schielzeth (2013). To visualize the models, effect plots were generated (Fox and Weisberg 2018). The plots show the predictions of model interaction terms by displaying the effects of the root and vegetation predictors at distinct terrain age steps. For this purpose, we chose four different terrain ages (70 a, 160 a, 4 ka, 12 ka) that represented mean terrain ages for each pair of moraines of the two study sites (e.g. mean of the two youngest moraines, Stein Glacier (30 a) and Griess Glacier (110 a), which is 70 a).
In addition to the LMERs, we used Principal Component Analysis (PCA) and Structural Equation Modelling (SEM) to investigate how root densities, Shannon index and vegetation cover directly and indirectly affected ASC. SEM is used to test direct and indirect relationships between variables in a multivariate approach. Our study design served as the hypothetical base for the initial SEMs. The PCA and SEM analyses were conducted based on the results of the LMERs, which showed a clear separation for young and for old moraines concerning the relationships between vegetation parameters and ASC. Therefore, we subsetted the whole dataset into two groups: the first group included the two young moraines of each glacier foreland (data of moraines ≤ 160a: Stein Glacier: 30a, 160a; Griess Glacier: 110a, 160a) and the second group included the two old moraines of each glacier foreland (data of moraines ≥ 3 ka: Stein Glacier: 3 ka, 10 ka; Griess Glacier: 4.9 ka, 13.5 ka). In order to account for the terrain-age dependent effects observed from the LMERs, we then ran the PCAs and SEMs separately on each group. The initial two SEM models were built based on the significant relationships found in the LMERs (f1: ASC ~ RLD + RMD, f2: RLD ~ Shannon Index, f3: RMD ~ Cover). The models were then simplified until an adequate fit was achieved. The adequacy of the models was determined via χ 2 -tests, AIC and RMSEA. Adequate model fits are indicated by a non-significant χ 2tests (p > 0.05), low AIC and low RMSEA (< 0.05) (Grace 2006). All statistical analyses were performed with R version 3.5.1 (R Core Team 2018).

Effects of terrain age on vegetation characteristics, root density, and ASC
Along both chronosequences, the vegetation characteristics varied across the moraines (Table 3). At the Stein Glacier foreland, the plant cover reached its peak on mid-successional moraines (137% at 160 a and 125% at 3 ka moraines; Table 3). At the Griess Glacier foreland, the plant cover was increasing significantly with moraine age (Tukey Test, p < 0.05). A similar pattern was observed for the Shannon diversity index (Table 3). The highest levels of diversity (H≈ 2.5 at 160 a and 3 ka moraines) were observed at the Stein Glacier foreland for the midsuccessional moraines. Conversely, at the Griess Glacier site, the Shannon index was significantly the highest at old moraines (H ≈ 2.3 at 4.9 ka and 13.5 ka moraines) compared to young moraines (H ≈ 1.7 at 110 a and 160 a moraine). At the Stein Glacier foreland, RMD significantly varied across the moraines ranging from 2.5 kg m −3 at the 30 a moraine to 17.7 kg m −3 at the 10 ka moraine (Table 3). At the Griess Glacier foreland, old moraines (4.9 ka, 13.5 ka) had significantly higher levels of RMD (≈ 11.5 kg m −3 ) than young moraines (3.1 kg m −3 at the 110 a, 5.3 kg m −3 at the 160 a moraine; Table 3). The amount of fine roots (< 1 mm) contributing to total root biomass ranged between 81 and 24%. Hence, RMD was considerably influenced by coarse roots (> 1 mm) ( Table 3).
At Stein Glacier foreland, RLD was highest at mid-successional stages (951 km m −3 at the 160 a, 845 km m −3 at the 3 ka moraines; Table 3). At the Griess Glacier site, RLD values of young moraines (110 a, 160 a ≈ 300 km m −3 ) were significantly lower compared to old moraines (4.9 ka, 13.5 ka ≈ 900 km m −3 ). Root length was almost exclusively driven by the fine root fraction (Table 3). ASC values increased with terrain age at both study areas, ranging from 0.33 to 0.97 at Stein Glacier and from 0.23 to 0.92 at Griess Glacier, respectively (Table 3). At the Griess Glacier foreland, the two oldest moraines had significantly higher ASC levels compared to the two youngest moraines. At the 160 a and 3 ka moraines of the Stein Glacier foreland, ASC was significantly higher compared to the 30 a moraine (Tukey Test, p < 0.05), and significantly lower compared to the 10 ka moraine (Table 3). Table 3 Development of plant cover, Shannon index, Root Mass Density (RMD), Root Length Density (RLD), and Aggregate Stability Coefficient (ASC) along the chronosequences of the two study sites (Stein Glacier, Griess Glacier).
Mean values of raw data and standard errors (n = 6) are given for all moraines. A post-hoc Tukey statistic was used to conduct multiple pairwise comparisons between the moraines of each of the two sites. Significant differences are highlighted with different letters. The fine root fractions of RMD and RLD are given in parentheses  age. For models 2 and 4 illustrating the relationship between vegetation cover, terrain age, and RMD, the interaction term was not significant.
RLD and RMD were both significantly positively affecting ASC (models 5 and 6, Table 4, Fig. 5). Moreover, the interaction term including terrain age was a  Table 2 are shown. The predictor variables were displayed for four distinct terrain ages [70 a (= years), 160 a, 4 000a, 12 000a] that represented mean ter-rain ages for each pair of moraines of the two study sites [e.g. mean of the two youngest moraines, Stein Glacier (30 a) and Griess Glacier (110 a), which is 70 a]. The black lines represent the model predictions and the black shades are the 95% confidence intervals. The black points are the partial residuals of the models. significant predictor in both models showing decreasing effect sizes with terrain age (model 5: p = 0.002, R 2 = 0.70; model 6: p = 0.001, R 2 = 0.67).  Table 2 are shown. The predictor variables were displayed for four distinct terrain ages [70 a (= years), 160 a, 4 000 a, 12 000a] that represented mean terrain ages for each pair of moraines of the two study sites (e.g. mean of the two youngest moraines, Stein Glacier Pass (30 a) and Griess Glacier (110 a), which is 70 a). The black lines are the model predictions and the blue shades are the 95% confidence intervals. The black points are the partial residuals of the models.
Testing the independent effects of plant cover and diversity on ASC, no significant relationship could be detected (Table 4; Fig. 5). However, the interaction effect of cover and terrain age was indicated resulting in decreasing effect sizes with terrain age (model 8: p = 0.17, R 2 = 0.67).
Because LMERs suggested different relationships between ASC and vegetation parameters (Fig. 5), we grouped the dataset into young (≤ 160 a) and old (≥ 3 ka) moraines and performed separate PCAs on each subset. Both PCAs showed two principal components with an eigenvalue larger than 1. The first two axes accounted for 83.27% (Fig. 6a) and 59.58% (Fig. 6b) to total inertia. The multivariate analyses revealed that ASC was more closely related to RLD on young terrain (Fig. 6a), whereas later during landscape evolution RMD seemed to be more important for ASC development (Fig. 6b).
As with the PCAs, SEMs were performed on young and old moraines separately. Both SEMs adequately fit the data on aggregate stability for the moraines with terrain ages ≤ 160a (χ 2 = 1.17, P = 0.28; AIC = 265.03; RMSEA ≤ 0.05, Fig. 7a) and for the moraines ≥ 3 ka (χ 2 = 0.25, P = 0.62; AIC = 314.49; RMSEA < 0.001, Fig. 7b) (standardized path coefficients are given in Fig. 7). Similarly, to the results of the PCAs we found that RLD increased ASC on young terrain, whereas RMD was the more important variable for ASC on older moraines (Fig. 7).

Discussion
Effects of parent material and terrain age on vegetation characteristics and soil aggregate stability As ASC generally increased positively with terrain age, this study supports the hypothesis that soil aggregate stability increases with progressing vegetation and soil development (Amezketa 1999). Soil development (i.e. the weathering of primary minerals, accumulation of organic matter, and the formation of secondary minerals) generally progresses at high rates in young soils and slows down after several centuries or millennia (Egli et al. 2001;Sauer 2015). Due to an increasing vegetation cover and increased above-and belowground biomass production with terrain age, the organic matter content of the topsoil increases as well (Burga et al. 2010;D'Amico et al. 2014). This was also confirmed by Musso et al. (2019) who found strong changes in However, such processes in soil development are expected to differ with geology. Carbonate-rich parent material weathers mainly by the dissolution of solid carbonates, which buffers soil pH in the neutral range, thereby slowing down chemical weathering processes (Egli et al. 2008). As a consequence, the accumulation of aggregate-forming clay minerals (Amezketa 1999) might be reduced, which could be an explanation for the lower ASC levels found at the youngest two moraines of the Griess Glacier site compared to the Stein Glacier site.
Furthermore, higher levels of vegetation cover prevent soil weathering products from being washed away by rainfall, thereby increasing the relative amount of fine soil particles (Burri et al. 2009). The organic compounds and the fine soil particles form micro-aggregates, mainly shaped and stabilized by plant roots and microorganisms (Amezketa 1999;Rillig and Mummey 2006;Tisdall and Oades 1982). The soil aggregates modify the soil structure and enhance the water and nutrient retention capacity of the soil, thereby again benefitting root growth and plant performance (Miller et al. 2002;Pohl et al. 2012).

The effect of vegetation characteristics on ASC is mediated via root density
Shannon diversity index and plant cover showed a positive correlation with root density, and thus we could confirm our first hypothesis assuming that plant community characteristics affect root density. A high species richness might increase the probability that plant species with deeper or denser root systems or higher biomass production occur ("selection effect" sensu Loreau and Hector 2001, see also Tilman et al. 2014, Garnier et al. 2016, which would then largely control the effects of vegetation on ASC and soil stability. In addition, higher plant diversity offers the possibility that the available soil volume might be more densely occupied by roots due to niche partitioning and subsequent resource use complementarity ("complementarity effect"). We found evidence that the effect of plant diversity on root density was getting stronger with terrain age. Similar results were found by Ravenek et al. (2014) who detected an increasing correlation between plant species richness and belowground productivity with time. A reason might be that during succession not only the total species number but also the number of highly productive and competitive species increases (Walker and del Moral 2003). This potentially enhances belowground Fig. 7 Structural equation models of vegetation effects on aggregate stability (ASC) in alpine glacier forelands for moraines ≤ 160 a (a) and for moraines ≥ 3 ka (b). Numbers on arrows are standardized path coefficients (equivalent to correlation coefficients). Asteriscs mark significant relationships (***: p-value < 0.001, *: p-value < 0.05). The initial models were build based on the study hypotheses and the results of the LMERs (f1: ASC ~ RLD + RMD, f2: RLD ~ Shannon Index, f3: RMD ~ Cover). The results were the simplified until adequate fit was achieved. a χ 2 = 1.17, P = 0.28; AIC = 265.03; RMSEA = < 0.05; b χ 2 = 0.25, P = 0.62; AIC = 314.49; RMSEA < 0.001. Acronyms: ASC = Aggregate Stability Coefficient, RLD = Root Length Density, RMD = Root Mass Density complementarity effects leading to an increase of root productivity with terrain age.
Testing the second hypothesis showed that RLD and RMD had a significant positive effect on the ASC, whereby the correlation weakened with increasing age. Such positive effects support the current view that roots are one of the main drivers of soil aggregate stability (Bardgett et al. 2014;Pohl et al. 2009;Reubens et al. 2007). Roots are a crucial part of the soil-forming processes mentioned in the previous section such as weathering and accumulation of organic matter. They contribute to increasing stability by mechanisms such as particle binding through root exudates and also physically stabilize the soil in the form of cohesive strength. Thereby, fine roots are more important than coarse roots for soil processes that directly depend on the physiological interactions between plants and soil, such as the formation of soil aggregates and nutrient cycling (Bardgett et al. 2014;Loades et al. 2010;McCormack et al. 2015). Our results highlight that the measured root characteristics are useful for predicting soil aggregate stability.
The third hypothesis addressed the direct impact of aboveground vegetation and diversity on ASC. Although our data did not provide statistically significant evidence for a direct effect of plant cover and diversity on ASC, we could clearly show that these parameters influenced root density, and those, in turn, drove ASC. We thus conclude that the effect of cover and diversity on ASC is mediated via root density. However, on young moraines, plant cover positively influenced ASC. This finding is in accordance with Pohl et al. (2012) who showed that species richness, vegetation cover, and root density were positively related to ASC at disturbed alpine slopes. Therefore, plant cover might not only have an indirect impact on ASC through roots but it could be also beneficial in itself (Cerdà 1998;Pohl et al. 2012). It reduces hydrological sediment transport and soil erosion on hillslopes induced by raindrop splash and surface runoff (Geißler et al. 2012;Goebes et al. 2015;Poesen 2018;Vannoppen et al. 2015). Furthermore, plant cover protects the surface from wind erosion and collects wind-blown silt, which enhances soil development as well as aggregate formation (Gyssels et al. 2005;Walker and del Moral 2003). ASC predictions change with terrain age: link to the concept of geomorphic succession The importance of the different predictor variables for ASC changed with terrain age (hypothesis 4). On older moraines, the impact of root density on ASC was weakened. Instead, species diversity was positively related to root biomass at those habitats. Such diverse alpine grassland biotopes are known to produce very stable hillslopes as described by Körner and Spehn (2002), picturing root systems of speciesrich slopes as the screws and nails of various forms and sizes in alpine soils. On young surfaces, patches with high plant cover and high root density had the highest ASC levels. The results of the LMERs indicate that rather cover than diversity affected ASC at early stages of soil formation suggesting that species identity was a stronger driver than species diversity. Moreover, we found evidence that RLD, which is linked to direct interactions at the root-soil interface, was the more important variable early during succession. Instead, on older moraines, RMD was a stronger driver of ASC suggesting that the accumulation of SOM is an important process for hillslope stability in the long term. On young surfaces, matforming species such as the creeping dwarf shrubs Salix retusa and Dryas octopetala build microhabitats with enhanced organic matter accumulation and microbial activity due to their dense rootstocks (Bjorbaekmo et al. 2010). Thus, in early-successional glacier foreland habitats, soil formation seems to develop initially from separate patches of ecosystemengineering species (Eichel et al. 2013;Eichel 2019). This applies to the concept of biogeomorphic succession, which is a sequence of three biogeomorphic phases forming a gradient of decreasing geomorphic and increasing ecological importance (Corenblit et al. 2007). In the pioneer phase, high geomorphic activity is assumed and only pioneer species adapted to the dominant geomorphic processes occur. In the biogeomorphic phase, geomorphic activity is lower and community composition is dominated by dwarf shrubs acting as geomorphic-engineering species. In the ecological phase, autecological processes such as competition prevail and geomorphic-engineering species are missing (Eichel et al. 2013). In light of this concept, the major changes in ASC development of this study could be attributed to the biogeomorphic phase where the strongest interactions of geomorphic and ecological processes are expected. The results of the present study suggest that this phase is crucial for the development of hillslope functioning and that important biogeomorphological processes, such as biostabilisation or bioaccumulation (Corenblit et al. 2011), change during this phase.

Conclusion
To our knowledge, this is the first study to investigate the development of ASC along alpine chronosequences by quantifying its interrelations with various vegetation characteristics. The results highlight that root density plays a major role in governing ASC on moraines of different substrate exposure. We show evidence that the effect of vegetation cover and Shannon index on ASC is mediated via root density. These relationships are most pronounced on young surfaces where engineering species create patches with high biogeomorphic activity.
As our study focused on the biotic factors, we encourage further research to include the impact of soil and hydrological parameters on ASC development. Especially for modelling ASC development, future research should take into account the grain size distribution as patches with high silt or clay content might be crucial for the initial aggregate formation on young surfaces. Furthermore, for a better understanding of the development of hillslope stability, it is important to measure further soil aggregate related parameters such as aggregate size and aggregate amount. Finally, in order to more mechanistically unravel the direct plantsoil interactions affecting aggregate stability at the hillslope scale the in-situ quantification of mycorrhiza and root exudates should be considered.