Initial carbonate weathering is linked with vegetation development along a 127-year glacial retreat chronosequence in the subtropical high mountainous Hailuogou region (SW China)


 
 The retreat of glaciers is exposing new terrains to primary plant succession around the globe. To improve the understanding of vegetation development along a glacier retreat chronosequence, we (i) evaluated a possible link between base metal (Ca, Mg, K, Na) supply and vegetation establishment, (ii) determined the rates of the establishment of soil and plant base metal stocks, and (iii) estimated the size of the main base metal fluxes.
 
 
 We determined base metal stocks in the soil organic layer, the mineral topsoil (0–10 cm), and in leaves/needles, trunk, bark, branches and roots of the dominating shrub and tree species and estimated fluxes of atmospheric deposition, plant uptake and leaching losses along the 127-yr Hailuogou chronosequence.
 
 
 Total ecosystem Ca and Mg stocks decreased along the chronosequence, while those of K and Na were unrelated with ecosystem age. Fortyfour and 30% of the initial stocks of Ca and Mg, respectively, were leached during the first 47 years, at rates of 130 ± 10.6 g m−2 year−1 Ca and 35 ± 3.1 g m−2 year−1 Mg. The organic layer accumulated at a mean rate of 288 g m−2 year−1 providing a bioavailable base metal stock, which was especially important for K cycling.
 
 
 We suggest that the initial high Ca bioavailability because of a moderately alkaline soil pH and carbonate depletion in 47 years, together with the dissolution of easily-weatherable silicates providing enough Mg and K to the pioneer vegetation, contributed to the establishment of the mature forest in ca. 80 years.


) and tropical regions (Cullen et al. 2013;Seehaus et al. 2019). Following glacial retreat, glacial debris is continuously deposited in the retreat areas serving as new substrate for soil and ecosystem development. Particularly during the time between deposition of the glacial debris and establishment of a vegetation cover, there is an enhanced risk of natural hazards typical for mountain areas, such as landslides, debris flows, erosion, and flooding. The main causes of this are the abundance of unconsolidated material and the lack of a vegetation cover which reduces surface runoff and wind and water erosion (Crozier 2010;Haeberli et al. 2017;Keiler et al. 2010;Portes et al. 2018;Richardson and Reynolds 2000). The similar climatic and geological conditions of a glacial retreat area allow us to apply space-for-time substitution approaches (Pickett 1989). Studies on primary succession on glacial moraines have been performed since the end of the nineteenth century (e.g., Coaz 1887; Crocker and Major 1955;Dickson and Crocker 1953;Friedel 1938) and continue to date (see compilation of studies in Supplementary Table 1). Glacial retreat has accelerated during the second half of the twentieth century (Zemp et al. 2008;Zhou et al. 2013) and is expected to continue for global mountain systems during the twentyfirst century (Hock et al. 2019;Nogués-Bravo et al. 2007). Thus, the knowledge of the processes that drive vegetation succession is crucial to cope with the risks originating from the globally increasing glacier retreat areas.
Because of the harsh, high-elevation site conditions, glacier forelands are challenging environments for vegetation development. Although pioneer plants usually start colonizing young surfaces a few years after glacier melt, it needs various decades, if not centuries to develop a full vegetation cover (Burga et al. 2010;Conen et al. 2007;Vilmundardóttir et al. 2015). Many young primary successions (<200 yr) developing on different parent materials in different climate zones are dominated by pioneer mosses, dwarf shrubs and shrubs but lack fully developed forests (see compilation of studies in Supplementary Table 1). The primary succession developing in the Hailuogou glacial retreat area has, however, surprisingly achieved the average biomass estimated for mature temperate coniferous forests (307 tones ha −1 , Cole and Rapp 1981) within a period of ~80 years ( Supplementary  Fig. 1, Luo et al. 2004).
Besides nitrogen (N) supply by microbial fixation from the atmosphere, the release of nutrients from parent materials by chemical weathering also plays a key role in vegetation succession (Chadwick et al. 1999). During chemical weathering processes, nutrients are released from minerals to soil solution as ions, which can readily move among the various geochemical and biological reservoirs and be lost via deep infiltration and stream water (Likens 2013). Plants can influence chemical weathering rates by a factor of two to five compared to the rates when plants are absent (Berner et al. 2004;Drever and Stillings 1997). Plants influence the biogeochemical cycle of nutrients by affecting their bioavailability in soils, taking up nutrients from the soil solution via their roots and their associations with mycorrhiza, storing them in tissues, and finally returning them to the soil via canopy leaching, litterfall, and organic matter decomposition (Vogt et al. 1986). In this context, the formation of a soil organic layer on top of the mineral soil plays an important role, because it stores nutrients and slowly releases them for re-use by plants (Lilienfein et al. 2001;Wilcke et al. 2002). Different tree species can further influence nutrient cycling in soil, because of their different nutrient requirements (Larcher 2003). Deciduous tree species usually have higher Ca, Mg and K concentration in leaves and greater foliar biomass than coniferous tree species (Ericsson 1994), which causes larger metal fluxes with litterfall in deciduous than coniferous forests (Carnol and Bazgir 2013). Although N and phosphorous (P) are considered the most commonly limiting nutrients affecting biomass production in terrestrial ecosystems (Elser et al. 2007;LeBauer and Treseder 2008;Vitousek and Howarth 1991;Lei et al. 2021;Yang et al. 2021), previous studies have shown that Ca, Mg and K can have a great influence on tree growth and vegetation development (Baribault et al. 2012;Burstrom 1968;Epron et al. 2012;Federer et al. 1989;McLaughlin and Wimmer 1999;Paoli and Curran 2007;Wright et al. 2011). Sodium (Na) can also be taken up and used by vegetation (Amtmann and Sanders 1999). Sodium regulates the cellular osmotic balance of microorganisms and thereby influences the decomposition of soil organic matter (Jia et al. 2015;Kaspari et al. 2009Kaspari et al. , 2014 and the resupply of bioavailable mineral nutrients. However, when high concentrations accumulate in the cytoplasm, Na interferes with K and Ca nutrition and disturbs efficient stomatal regulation (Marschner 2012;Tavakkoli et al. 2010). Moreover, high concentrations of any of these elements (i.e., Ca, Mg, K and/or Na) can lower or even inhibit the uptake of the others, which is known as cationic antagonistic effect (Diem and Godbold 1993;Ertiftik and Zengin 2017;Fageria 1983;Rhodes et al. 2018). Thus, at high Ca availability, e.g., in soils on limestone, Mg and K uptake can be limited (Ertiftik and Zengin 2017).
Our overall goal was to improve the understanding of the drivers of vegetation development on glacial debris. We hypothesized that (i) the vegetation establishment in the Hailuogou area, which surprisingly reaches the stage of a mature forest in only ~80 years, is linked with changes in total base metal stocks of the ecosystems attributable to a high initial bioavailability of Ca, Mg, K and Na, because the carbonic acid/carbonate buffer system maintains a near-neutral soil pH with a high base saturation in the early phase of the vegetation establishment. (ii) This high base metal availability allowed for the faster buildup of large base metal stocks in plants and in the soil organic layers than at other <200 yr-old glacier retreat chronosequences (Supplementary Table 1), driven by particularly high base metal fluxes between soils and plants. (iii) After the carbonates and easily weatherable silicates have been dissolved, the base metal supply decreases considerably. However, this is compensated by the presence of an organic layer which provides plant-available nutrients to the vegetation and by a shift in the vegetation composition from deciduous to coniferous trees which reduces the base metal demand of the vegetation. The Hailuogou chronosequence has not been strongly disturbed by human activities, making it an ideal area to study biogeochemical element cycling in the early stage of soil and vegetation development at decadal scale.

Study area
The Hailuogou Glacier flows down the eastern slope of the Gongga Mountain, in the transition zone of the Sichuan Basin and the Tibetan Plateau, southwest China (Fig. 1). Because of a temperature increase, the Hailuogou Glacier has markedly retreated since the late nineteenth century (Li et al. 2010), developing an approximately 2 km-long, 50-200 m wide chronosequence, with an elevational difference of ~130 m. The climate is mainly controlled by the Southeast Monsoon. The mean annual temperature and precipitation at the Gongga Mountain Alpine Ecosystem Observation Station, Chinese Academy of Sciences, which is located in the same Hailuogou catchment (Fig. 1) are 4.2 °C and 1947 mm, respectively. For the eastern slope of the Gongga Mountain, where our study site is located, Zhou et al. (2016a) have reported a mean annual temperature of 5.6 °C and 3.9 °C at 2772 m.a.s.l. and 3060 m.a.s.l., i.e., an elevational temperature gradient of ca. 0.006 °C m −1 or 0.77 °C for the ca. 130 m elevational difference. In the same study, the mean annual precipitation was 1923 mm at 2772 m.a.s.l. and 1933 mm at 3060 m.a.s.l. resulting in a difference of ca. 4.5 mm (or 0.23% of the precipitation at 2772 m a.s.l.) between the youngest and oldest study sites. Rain falls mainly during the vegetation period (May to September; Wu et al. 2013). The parent material of soil formation is moraine and is mainly composed of silicates, including plagioclase (285 mg g −1 ), quartz (245 mg g −1 ), biotite (120 mg g −1 ), hornblende (120 mg g −1 ) and K-feldspar (100 mg g −1 ) and non-silicatic carbonates (<103 mg g −1 ) and apatite (<21 mg g −1 ) Zhou et al. 2016b). The short time of pedogenesis (<130 years) resulted in the formation of initial soils without B horizon classified from the youngest to the oldest soils as Leptic Calcaric to Folic Dystric Regosols according to the World Reference Base for Soil Resources (IUSS Working Group WRB 2014). The intense carbonate weathering and the rapid establishment of the vegetation decreased the pH of the mineral soil from 8.0 to 5.8 in only 47 years (Table 1, Supplementary Table 2).
Exposure ages of the different sites and trees ages at the Hailuogou Chronosequence have been determined based on historic records including aerial fotos and tree rings, respectively (Zhong et al. 1999). This study includes seven dated sites, de-glaciated between 0 (Site 1; 2982 m a.s.l.) and 127 years before our sampling (Site 7; 2855 m), which were previously described by Zhou et al. (2013). A complete primary vegetation successional sequence has developed along the retreat area, from bare land with pioneer mosses via a bush stage to forests successively dominated by the shrub species Hippophae rhamnoides L., which had a tree-like habit (i.e., it showed  Luo et al. (2004).
Previous research on biologically relevant elements in the Hailuogou region focused on the change in C, N and P mass ratios and stocks (e.g., Bing et al. 2016b;He and Tang 2008;Wu et al. 2015;Yang et al. 2021;Zhang et al. 2021;Zhou et al. 2013) or contamination by trace metals (e.g. Bing et al. 2016aBing et al. , 2019Luo et al. 2015;Wang et al. 2020). The latter studies described a rapid accumulation of organic matter and N in the topsoil, while the soil formation rate along the Hailuogou chronosequence decreased with increasing ecosystem age (He and Tang 2008). Zhou et al. (2016b) reported that weathering processes in the Hailuogou are initially dominated by the weathering of carbonates followed by a more intense biogeochemical weathering (>80 yr) because of lower pH values. In a recent study, Yang et al. (2021) suggested that N is the main growth-limiting element in the Hailuogou area based on N:P stoichiometric ratios in leaves. Wang et al. (2021) recently reported a positive effect of pioneer N 2 -fixing plants on the establishment of other non-N 2 -fixing species after glacier retreat. However, the relationship between the bioavailability of base metals (Ca, Mg, K, Na) and vegetation development in this area have not yet been evaluated.

Field sampling
Soil samples were collected in August 2017 from seven ecosystem succession stages (Sites 1-7), which were exposed for 0, 5, 37, 47, 59, 87 and 127 years since glacier retreat (Fig. 1). Each ecosystem succession stage was sampled in triplicate. The distance between the sampled soil profiles was at least 20 m, except at Sites 1 and 2, where the distance was reduced to 10 m because the studied valley is narrower in the proximity of the glacier. Soil profiles were hand-dug and five soil horizons sampled: Oi (fresh litter), Oe (shredded litter), Oa (dark layer of decomposed humus), A (surface mineral soil with humus enrichment), and C (weathered soil parent material). No organic layers were present at Sites 1 and 2 (0 and 5 years old, respectively).
Freshly cut leaves, 3-year-old needles, bark, branches, trunk, and roots of the dominant tree and shrub species were collected in triplicate for elemental analysis, in the surrounding of our replicate soil profiles, between August and October 2017 at each study site. To achieve a more representative sample of each replicate, a minimum of three individuals were randomly sampled from the trees of the same area. Four species were sampled along the chronosequence: H. rhamnoides (Site 3, 37 years), P. purdomii (Site 4, 47 years), A. fabri (Sites 5-6, 59 and 87 years), and P. brachytyla (Site 7, 127 years). Tree branches were randomly sampled from the tree canopy using pole shears. Bark samples were collected using an outdoor knife. The depth of cut was adjusted depending on the bark thickness in order to take a representative sample without reaching the cambium. Trunk samples equivalent to the radius of the sampled tree were collected using a tree corer with a diameter of 5.15 cm. Roots were manually washed out of a soil sample with defined volume.
Soil samples (mineral and organic horizons) were air-dried to constant weight in a drying room located in the research station. Mineral samples were sieved to collect the two fractions fine earth (<2 mm) and stones (soil >2 mm). Plant residues were removed from the soil. Leaf, needle, branch, bark, trunk and root samples were oven-dried to constant weight at 40 °C during 72 h and stored in sealed bags. Trunk and root samples were manually homogenized using an agate pestle and mortar. Bark and branch samples were homogenized using a blade grinder (Polymix PX-MFC 90 D) equipped with a 3.0 mm mesh size sieve. Aliquots of all samples were ground using a ball mill equipped with a zirconium oxide jar.

Chemical analyses
Soil pH was determined from the air-dried samples (<2 mm) by a glass electrode (WTW SenTix 81) in a 1:5 (v/v) soil:water suspension. Total element concentrations in the organic horizon, fine earth and stones were determined after digestion with concentrated HNO 3 /HF/H 2 O 2 (4:1.5:1, v:v:v) and in the plant compartments after digestion with 8 mL of concentrated HNO 3 and 2 mL of concentrated H 2 O 2 in a microwave oven (MARS6Xpress, CEM, Kamp-Lintfort, Germany). Elemental concentrations in the digests were analyzed using an inductively-coupled plasma optical-emission spectrometer (5100 VDV ICP-OES, Agilent, Waldbronn, Germany). Accuracy was assessed by the analysis of certified reference materials: SRM1547 (peach leaves) and SRM1515 (apple leaves) for organic samples and BCR2 (Columbia River Basalt 2) for soil samples. Average recoveries ± single standard deviation were 100 ± 10% for all certified element concentrations. Precision determined by duplicate measurements was <5%. Carbonate concentrations (CO 3 2− ) in the fine earth was determined by measuring the volume of emitted carbon dioxide after reaction of the sample with 10% HCl in a Scheibler calcimeter. Inorganic C concentrations in the stones were measured with an elemental analyser (EuroEA Elemental Analyzer, HEKAtech, Wegberg, Germany), after muffling at 550 °C.
Exchangeable cations (Ca, Mg, K, Na, and Al) in the mineral soil (2 g) were extracted with 1 M NH 4 NO 3 in a 1:25 soil to solution ratio following the method described by Zeien and Brümmer (1989) and measured by ICP-OES. All our element concentrations and stocks refer to the air-dry mass of the samples.

Stocks
Stocks (E plant compartment ) of the respective element E (i. e., Ca, Mg, K or Na) in the different tree compartments (leaves/needles, branches, trunk, bark and roots) of the dominant shrub or tree species at each site were calculated by multiplying the corresponding compartment biomass (B c ) with the element concentrations measured in the different compartments (C E ) (Eq.1).
The B c was calculated by multiplying the total biomass at each plot (Supplementary Table 3), estimated for our study site from Luo et al. (2004) after applying a logistic fitting (Supplementary Fig. 1), with the proportional contribution of the tree compartments to total biomass taken from Zhou (2013) and for partitioning between trunk and bark from Wilcke and Lilienfein (2004) (Supplementary Table 4). Tree biomass published by Luo et al. (2004) was calculated with the allometric equations reported by Zhong et al. (1997), which are based on tree height and diameter at breast height. The biomass of the shrub species was obtained by Luo et al. (2004) through destructive sampling after harvesting and drying. With this approach, we assumed that the mean base metal concentrations of the dominant species was representative for the whole vegetation on the considered plot. To support this assumption, we calculated according to the mass contributions of the various tree compartments (leaves/needles, branches, bark, trunk, roots) weighted mean base metal concentrations of the dominant species ( Supplementary Fig. 2). The results illustrate that the nutrient concentrations in the shrub and trees were low at Sites 3, 6, and 7 and higher at Sites 4 and 5. The coniferous tree species A. fabri showed significantly higher base metal concentrations at Site 5 than at Site 6. Moreover, at Site 3 the tree-like shrub H. rhamnoides contributed more to the total biomass than can be expected from its cover ( Table 1), because of its multiple branches and much taller stature than the second most abundant herb A. adsurgens. At Sites 6 and 7, the cover of the two coniferous tree species A. fabri and P. brachytyla Plant Soil together was >80% of the total cover, which we took as approximate measure of the contribution of the individual species to the total biomass on the respective plot. The two coniferous tree species did not show significantly different mass-weighted mean base metal concentrations ( Supplementary Fig. 2). With respect to our extrapolation, the most critical study sites were Sites 4 and 5 with a mixed vegetation. To assess the validity of our extrapolation, we compared the base metal stocks of an assumed 100% cover of A. fabri with an assumed 100% cover of P. purdomii and did not find significant differences in the Ca, Mg and K stocks in the biomass ( Supplementary Fig. 3). The estimate of the Ca, Mg and K stock in P. purdomii on Site 5 is based on the assumption that the lack of a difference in the concentrations of these elements in the leaves between Sites 4 and 5 implies that the concentrations in all other plant compartments are also similar.
Element stock in each soil organic horizon (E Oa/Oe/Oi ) was calculated considering the elemental concentration (C E Oa/Oe/Oi ), the bulk density (BD i ) and the thickness (T i ) of the respective horizon (Oa, Oe, Oi) (Eq. 2). Total organic layer stock at each plot was calculated by summing up all individual horizons (Eq. 3).
Soil mineral stocks (E A/C ) were calculated considering the contribution of the stones (χ i stones , vol.%) in the soils along the chronosequence, T i , BD i , and the elemental concentration in the stones (>2 mm, C Ei stones ) and in the fine earth (fe <2 mm, C Ei fs ), for the A and C horizons, respectively (Eq. 4-6). where: Bulk densities (BD i ) and stone volume were taken from Zhou et al. (2016b) and Wang et al. (2019), who E A∕C stones g m −2 = C Ei stones g g −1 x 2.65 x 10 6 g m −3 x T i (m) studied the same sites as in our study. The soil particle density was assumed as 2.65 g cm −3 (Eq. 6; Zhou et al. 2016b). To calculate stocks of the extractable cations in the mineral soil, we only considered the fine earth, assuming that no extractable cations were present in the stones.
To account for the decrease in bulk density of the A horizon related with the loosening of soil during pedogenesis by weathering and ecosystem succession, the thickness of the C horizon (T C ) underlying the A horizon was adjusted in our stock calculation using Eq. 7.
Where T A is the thickness of the A horizon, 0.1 m is the initial mineral soil depth and BD A and BD C are the field bulk densities of the local A and C horizons, respectively. Total mineral soil stock (E mineral horizons ) was calculated by summing up the individual A and C horizon stocks. Additionally, we calculated the stocks of the mineral soil assuming that the density of the C horizon was initially the same at all the sites (BD C = BD C at Site 1, along the chronosequence) ( Supplementary Fig. 4). There were no significant differences between the two approaches so that we only present the results of Eq. 7.
Total stock (g m −2 ) at each plot was calculated by summing up the individual stocks of the different plant compartments, (E plant compartment ; leaves/needles, branch, trunk, bark and root) and organic and mineral horizons (Eq. 8).

Fluxes
Annual plant uptake of base metals was calculated by multiplying the annual net primary productivity (NPP) (Supplementary Table 3) with the massweighted mean element concentration of the plant (C E plant ) (Eq. 9). NPP was taken from Luo et al. (2004), who derived it from measurements of the tree diameter growth at breast height of all individual trees or shrubs on a study plot, which was translated into biomass growth with the help of the allometric equations of Zhong et al. (1997).
E total stock = E plant compartment + E organic horizons + E mineral horizons Total deposition (TD E , g m −2 year −1 ) for any element E was calculated with Eq. 10.
where BDep E represents the bulk deposition measured with a Hellmann-type rain collector and DD E is the estimation of the fine particulate dry deposition. Data of BDep E were obtained from the Gongga Mountain Alpine Ecosystem Observation Station (3000 m above sea level), Chinese Academy of Sciences ), which is located in the Hailuogou river catchment, where our study sites are. Dry deposition (DD E ) was roughly estimated by multiplying the relationship DD/BDep E published by Wilcke et al. (2017), calculated using the canopy budget method from Ulrich (1983) in a tropical forest with similar precipitation than that in the Hailuogou area. Because the scavenging of fine particulates from the atmosphere is related with the size of the plant surface of the forest canopy (Freer-Smith et al. 2005;Song et al. 2015), we scaled the DD/BDep E ratio by multiplying with the ratio of biomass at each individual site to that at Site 7 (127 years). This assumes that only the forest at Site 7 had a similar DD/BDep E ratio than the native montane forest in Ecuador used as reference and that the scavenging surface of the canopy is proportional to the biomass.
Leaching of the element E (L E ) from the soil to the stream (= element export with the stream) was estimated from the concentration in the stream water (C E SW , data from Zhou et al. 2016b), the mean annual precipitation (Pr, data from Wu et al. 2013) and the fraction of rainfall lost via the stream (α) (Eq. 11). The runoff coefficient α was individually adjusted to each of our sites by multiplying the transpiration of a mature forest taken from Sun et al. (2020) with the ratio of biomass at each individual site to that at Site 7. The canopy interception loss was taken from Sun et al. (2013), who provided interception loss values for different succession stages of our chronosequence allowing for interpolating data for two sites for which no data were available (Supplementary Table 3). The fraction of rainfall lost via the stream (9) E plant uptake g m −2 year −1 = NPP g m −2 year −1 x C E plant g g −1 (10) TD E = BDep E + DD E near the glacier forefront was estimated at 86% (Brock et al. 2010;Small et al. 2018).
Our estimate of the base metal leaching losses is based on the assumptions that all water passes through the deep subsoil, which consists of the unchanged glacial debris at all study sites and that during this passage chemical equilibrium of the soil solution with the substrate is reached via almost instantaneous cation exchange and mineral precipitation processes because any water flow through the debris and soils takes long enough to reach these fast equilibria. This assumption is reasonable because the length of the water path from the study soils to the stream is at least in the range of >1 m so that even at a high water conductivity of 10 −4 m s −1 , it would take the water >160 h to pass from the soil to the stream and thus long enough for equilibration. As a consequence, the base metal output of our study sites only depended on the fraction of water percolating through the subsoil to the stream described by α and the mean base metal concentration in stream water.
Mean depletion rates (g m −2 year −1 ) of the Ca and Mg stocks in the total ecosystem along the first stage of the chronosequence (from 0 to 47 years) were approximated by fitting a linear relationship of the total stocks versus time. Mg and K release from the depletion of biotite and hornblende was estimated from mineral composition data taken from Zhou et al. (2016b) for the same study sites.
In order to compare the budget of the independently estimated ecosystem fluxes (input: atmospheric deposition; outputs: loss with stream flow, accumulation in the organic layer and biomass; all in g m −2 year −1 ) with the cumulative stock change in the mineral topsoil (uppermost 10 cm, in g m −2 ), we multiplied the respective flux by the site age.

Statistical evaluation
Normal distribution of the data was checked with the Shapiro-Wilk test. Pearson correlation, linear regression, and non-linear function fitting were used to evaluate relationships between variables. Calculation of total stocks using (i) the local BD C (Eq. 7) or (ii) assuming that the density of the C horizon was initially the same at all sites (i.e., equal to the glacial (11) L E g m −2 year −1 = C E SW g L −1 x Pr L m −2 year −1 x α debris, 0 yr old site) were tested for significant differences with independent two-sample t-tests at each site. One-way analysis of variance (ANOVA) and a Tukey's HSD post-hoc test were applied to detect significant differences between total stocks among the various study sites. Normal distribution of residuals was visually inspected. Homoscedasticity was confirmed for our data after applying Levene's test. Statistical analyses were conducted with the statistical software R (R Core Team 2019). Tukey's HSD posthoc tests following ANOVA were performed with the 'agricolae' package (de Mendiburu 2021). Significance was set at p < 0.05.

Base metal depletion rates in the entire ecosystem
Total ecosystem stocks of Ca and Mg decreased along the chronosequence (Fig. 2a, b). Most of the depletion of Ca and Mg occurred during the first 47 years of soil exposure, paralleling the pH decrease from 8.0 to 5.8 in the mineral topsoil and the loss of carbonates (Table 1), likely mainly calcite (CaCO 3 ) ). Ca and Mg stocks decreased at rates of 130 ± 10.6 g m −2 and 35 ± 3.1 g m −2 over the first 47 years of soil development, respectively. This entailed that 44% of the initial stock of Ca was leached during the first 47 years after glacier retreat, accompanied by 30% of the initial Mg stock. After most of the carbonates were dissolved and the pH decreased to the strong acid/aluminum oxide buffer range, the annual depletion of Ca and Mg slowed down (Fig. 2a, b; Table 2).
Contrary to Ca and Mg, there was no clear relationship between K and Na total ecosystem stocks and development time (Fig. 2c, d). Instead, K and Na total stocks remained similar at on average 5162 ± 710 g m −2 K and 3530 ± 429 g m −2 Na along the chronosequence.

Base metal stocks in the soils
Most of the total element (i.e., Ca, Mg, K and Na) stock was stored in the mineral topsoil (0-10 cm;  Table 5). The initial weathering of carbonates and subsequent loss of Ca (Tables 1 and  2) accounted for 40-60% of the total loss of Ca in the mineral topsoil along the chronosequence, compared to the immediate glacier forefield (0 yr-old site). Mg stocks in the mineral topsoil soil decreased from 4990 ± 196 g m −2 to 2620 ± 202 g m −2 over the first 47 years, and further decreased only slightly from 47 to 127 years (Supplementary Table 5). K and Na mineral soil stocks remained similar at on average 5075 ± 700 g m −2 K and 3459 ± 420 g m −2 Na along the whole chronosequence.
The exchangeable Mg stock in the mineral soil consistently increased along the chronosequence from 0.7 ± 0.0 g m −2 to 2.3 ± 0.5 g m −2 , while no trend was found in the exchangeable K and Na concentrations with time (Supplementary Table 5). It was not possible to quantify exchangeable Ca in the carbonatecontaining soils, because of the partial dissolution of calcite in our extract (Tessier et al. 1979;Klimova et al. 2011).
The accumulation of organic matter on top of the mineral soil, which increased at a mean annual rate of 288 g m −2 over the studied time span, promoted the formation of an organic layer of increasing thickness (Supplementary Table 2) that stores a large amount of nutrients (Supplementary Table 5). The strongest element accumulation in the combined organic horizons occurred from 37 to 47 years after glacier retreat (Fig. 3), accompanied by the establishment of a young P. purdommii-dominated forest. After the establishment of the coniferous forest, 87 years after deglaciation, the annual accumulation sharply decreased for all the studied elements ( Table 2). The mean concentration of the base metals in the organic horizons increased with depth (i.e., Oi < Oe < Oa) (Fig. 3).

Base metal fluxes
Annual atmospheric deposition to the whole study area varied from 1.3-3.2 g m −2 Ca, 0.3-0.7 g m −2 Mg, 3.6-8.0 g m −2 K, and 0.3-0.8 g m −2 Na ( Table 2). The estimated dry deposition increased along the chronosequence, since the scavenging of fine particulates from the atmosphere by the vegetation increases with increasing surface area of the canopy (Freer-Smith et al. 2005;Song et al. 2015), which we approximated by biomass (Table 2). Annual losses of Ca and Mg via the stream were between one and two orders of magnitude higher than the estimated bulk deposition, while those of K and Na were in a similar range in bulk deposition and stream export (Table 2; Fig. 4).
Total aboveground and belowground plant biomass increased with site age along the  Supplementary Fig. 5) at rates of 380 ± 33 g m −2 year −1 . Annual plant uptake ranged from 22 to 64 g m −2 Ca, 2.1-4.8 g m −2 Mg and 29-44 g m −2 K, but did not show a consistent temporal trend along the chronosequence (Table 2). Because Ca, Mg, and K concentrations were considerably higher in P. purdomii foliage than in all other sampled tree species, sites where P. purdomii was growing usually showed the highest plant uptake of base metals. Annual Na uptake ranged from not detected to 0.1 g m −2 . The accumulation of base metals in the organic layers along the chronosequence generally decreased in the order Ca > K > Mg > Na (Table 2).
To determine whether the individual fluxes matched the total base metal loss, we compared the release of base metals in the mineral topsoil (0-10 cm), which represented the part of the mineral soil where most of the roots were located, with the budgets of the independently determined individual fluxes along the chronosequence (Fig. 5). Moreover, we included the part of the change in the Ca, Mg and K stocks, which was attributable to the loss of easily weatherable minerals (calcite, biotite and hornblende). The CO 3 2− concentrations decreased in the topsoil at a rate of 113 ± 15 g m −2 yr −1 CO 3 2− during the first 47 years after glacier retreat and was then gone (Table 1). This resulted in a Ca depletion of 79 g m −2 yr −1 . The Mg depletion in the mineral topsoil (0-10 cm) attributed to the loss of hornblende + biotite was 21 g m −2 yr −1 in the first 47 years of the chronosequence, during which the pH dropped markedly (Table 1) and 8 g m −2 yr −1 in the later 80 years. The K depletion attributed to the loss of biotite was 19 g m −2 yr −1 (0-47 yr) and 6 g m −2 yr −1 (47-127 yr).
Other minerals included in the glacial debris and containing Mg or K, such as chlorite or K-feldspar, have not been included in this estimation because their concentration hardly changed after 127 years. The chlorite concentration decreased from 40 to 30 mg g −1 in the topsoil while that of K-feldspar  Fig. 4 only explained a small part of the changes in the stocks of Ca and Mg, while the loss of the considered easily weatherable minerals explained a larger part.

Change in the total and soil base metal stocks and fluxes
The most rapidly depleted element from the mineral topsoil along the chronosequence was Ca, followed by Mg, while K and Na stocks remained similar (Fig. 2). We attribute the higher loss of Ca relative to the other elements to the depletion of carbonates (Supplementary Tables 2 and 5), dominated by calcite (CaCO 3 , Zhou et al. 2016b). The CO 3 2− decrease rate in the topsoil of 113 ± 15 g m −2 yr −1 CO 3 2− during the first 47 years after glacier retreat was higher than in other soil chronosequences evolving from parent materials containing a broad range of CaCO 3 concentrations (20-340 mg g −1 ), where the loss of CO 3 2− ranged from 14 to 65 g m −2 year −1 (Van Breemen and Protz 1988). However, our CO 3 2− dissolution rate was in the range of that found in older soils (>10.000 years) in Switzerland (91-136 g m −2 year −1 CO 3 2− ; Egli and Fitze 2001) and slightly lower than in ~193 yr old Calcaric Regosols soils in the Alps, which, in contrast to our study site, originally contained more carbonates (70-180 mg g −1 of CO 3 2− ) in the topsoil (Egli and Fitze 2001). The higher temperatures and higher precipitation together with the influence of the vegetation succession may have contributed to the more rapid CO 3 2− depletion at the Hailuogou chronosequence compared to other proglacial areas (Supplementary Table 1). Because the dolomite concentrations remained constant at ca. 20 mg g −1 along the chronosequence Zhou Fig. 5 Release of Ca (a), Mg (b), K (c) and Na (d) from the mineral topsoil (−10 cm) compared with the budget of the independently determined individual fluxes including the inputs atmospheric deposition and the outputs accumulation in the different tree compartments and organic layers and export with the stream (sum of the fluxes shown in Fig. 4). We also included the release of Ca, Mg and K attributed to calcite (Ca), biotite (Mg and K) and hornblende (Mg) depletion (CaCO3 + Bt + Hbl) relative to the 0 yr-old site. Error bars represent single standard deviations of three replicate plots at each site, which was in the case of the flux budget calculated by Gaussian error propagation from the standard errors of the individual fluxes et al. 2016), the Mg release cannot be attributed to the carbonate weathering but must be related with other easily weatherable minerals (such as e.g., biotite and hornblende).
The order of the element release at the beginning of soil development was consistent with the order of mineral dissolution (i.e., carbonate minerals > ferromagnesian minerals > feldspars; Lichter 1998;White et al. 1996;Zhou et al. 2016b). However, our findings are not consistent with the sequence of the ion release by weathering inferred from open-system element budgets at other chronosequences, where Na and K are the most easily leached ions (Bain et al. 1993;Harden 1988). We attribute this discrepancy to the low concentration of K and Na in easily weatherable minerals which release high quantities of Ca and Mg in the early stage of the chronosequence but little K and Na. Mg-rich minerals, such as hornblende and biotite (which also contain K), are more readily weathered than those that contain high K and Na concentrations but little Mg, e.g., orthoclase and albite (Clow and Drever 1996;Wilson 2004). At the Hailuogou chronosequence, the dissolution of Mg-rich silicate minerals contributed to the faster decrease of Mg (Fig. 2b) than K stocks (Fig. 2c), which is in line with the findings of Zhou et al. (2016). Moreover, the release of Mg from hornblende and biotite had a similar size as the total Mg loss from the studied ecosystems (Fig. 5). The strong mineral soil acidification indicated by the drop in soil pH along the chronosequence (Table 1 and Supplementary Table 2), accompanied by the establishment of a coniferous forest, may have further intensified the weathering of silicate minerals in the older stages of the chronosequence, promoting Mg and K release.
The one order of magnitude higher K than Mg atmospheric deposition (Table 2) but similar Mg and K export with the stream (Fig. 4b, c) kept K stocks nearly constant along the chronosequence (Fig. 2c), while those of Mg decreased (Fig. 2b). The accumulation of K in biomass and in the organic layers (Fig. 4c) contributed to the slight decrease of K in the mineral soil along the chronosequence (Fig. 5c), but compensated the total ecosystem stocks which remained unrelated with ecosystem age (Fig. 2c). In line with our observations, Lichter (1998) reported that K showed a negligible depletion during the first 1000 years of soil development but the K stock linearly decreased in the longer term (> 4000 years) along a sand dune chronosequence in Wilderness State Park (Michigan, US), in which the dissolution of carbonates also promoted an initially rapid decrease in pH from 8.5 to 4.3.
Because of the low Ca and Mg deposition from the atmosphere, the accumulation of deposited Ca and Mg contributed little to the total stocks of Ca and Mg (Table 2). Annual losses of Ca and Mg via the stream were between one and two orders of magnitude higher than the estimated bulk deposition (Table 2), which contributed to the stock decrease of both elements along the chronosequence (Fig. 2a,b). Na inputs and outputs were low (Fig. 4d) and fluxes hardly changed along the chronosequence (Table 2) explaining that the Na stocks remained nearly unchanged (Fig. 2d). We attribute this to the prevalence of Na in little weatherable minerals such as e.g., albite.
The finding that the flux budget did not explain the total losses of Ca, Mg and K from the ecosystems of the Hailuogou chronosquence (Fig. 5) indicates that we missed output fluxes. We speculate that an important missing flux might be leaching into the subsoil (i.e. the C horizons), where the elements are stored in the exchangeable pool or as re-precipitated minerals.
We explain the low K and Na stock 37 years after glacial retreat and the related apparent strong depletion of K and Na between 5 and 37 years and later accumulation between 37 and 47 years to the spatial heterogeneity of the glacial debris (Fig. 2c, d), which might have resulted in a locally lower initial stock of K and Na in the parent material of the 37 yr-old site (Supplementary Table 5). Because Ca and Mg are associated with other minerals (and rocks) than K and Na, the heterogeneity of the glacial debris does not necessarily affect all elements in the same way.

Biological stocks and fluxes of base metals
The high release of Ca from the intense initial carbonate weathering enabled the build-up of a large biomass stock of Ca (Supplementary Table 5), which was supported by the release of sufficient Mg and K to soil solution via weathering of easily weatherable silicate minerals, such as biotite and hornblende. Ca and K stocks in the foliage of the tree species developing at Hailuogou ranged from 5.4 to 13.3 g m −2 and 7.2 to 18.8 g m −2 , respectively (Supplementary Table 5), and were mostly higher than observed in different similarly aged forests of the temperate zone, for which we found comparison values, e.g., the 70 year-old forests in the Walker Branch watershed (Tennessee, US) on acidic soils developed from paleozoic dolomitic bedrock (5.8-7.5 g m -2 Ca, 3.6-5.3 g m -2 K, Johnson and Henderson 1989). In contrast, Mg stocks ranged from 0.6 to 2.1 g m −2 and were slightly below those in the Walker Branch watershed or in the same range (1.2-2.1 g m −2 Mg, Johnson and Henderson 1989). Along the Hailuogou chronosequence, plant biomass reached up to 40,630 g m −2 within 127 years of ecosystem development (calculated from Luo et al. 2004; Supplementary Fig. 5). Total biomass estimated at the different 70 year-old forests in the Walker Branch watershed ranged between 16,380 and 19,770 g m −2 (Edwards et al. 1989), which is less than the total biomass which developed at the Hailuogou chronosequence within only 47 years (Supplementary Fig. 5). The higher biomass production and metal uptake at the Hailuogou chronosequence occurs in spite of less favorable climatic conditions for forest growth than at the Walker Branch watershed (mean annual precipitation of 1947 mm and mean annual temperature of 4.2 °C at Hailuogou  vs. 1410 mm and 14 °C at the Walker Branch watershed (Edwards et al. 1989). Cool temperature and a higher risk of waterlogging at the Hailuogou chronosequence should have reduced biomass production (Eamus 2003). The aboveground biomass in another young (55 year-old) forest developed on glacial till originating from the local metasedimentary and magmatic rocks (Barton et al. 1997) in the Hubbard Brook Valley with a mean annual temperature of 3.7-6.7 °C and a mean annual precipitation of 1326 to 1607 mm (Campbell 2021), which should be again more favorable for forest growth than at the Hailuogou chronosequence showed also considerably lower element stocks (38.3 g m −2 Ca, 3.6 g m −2 Mg and 15.5 g m −2 K; Likens 2013) and, consequently, lower annual vegetation uptake (6.2 g m −2 year −1 Ca, 0.9 g m −2 year −1 Mg and 6.4 g m −2 year −1 K; Likens 2013) than the >47 year-old Hailuogou sites. This comparison indicates that the vegetation along the Hailuogou glacial retreat chronosequence took up nutrients more rapidly than other young forests, which is strongly associated with the faster rate of biomass production. Moreover, nutrient release by weathering of the young carbonatic glacial debris was likely faster than on the older substrates at the forest sites used for comparison.
Coincidence of carbonate depletion and shift from deciduous to coniferous forest After the reactive carbonates were leached as indicated by the drop of the mineral topsoil pH from 8.0 to 5.8 over the first 47 years of the chronosequence (Table 1), net weathering rates decreased because less soluble minerals remained in the soil. This occurred in line with the vegetation change from deciduous to coniferous trees and their associated lower nutrient demand and slower nutrient cycling (Table 2; Carnol and Bazgir 2013;Ericsson 1994). As a consequence, nutrient supply by weathering and nutrient demand by the vegetation seemed to match well. Because the differences in mean annual climatic conditions between the youngest and the oldest study site separated by ca. 130 m elevation difference were small (0.77 °C, 4.5 mm), we assume that climate played a minor role in explaining the changing vegetation composition.
The most rapid elemental accumulation in the organic layers occurred between 37 and 47 years after glacier retreat, accompanied by the stabilization of a half-mature P. purdommii forest and sharply decreased after the stabilization of a coniferous forest, 87 years after deglaciation (Fig. 3), likely because of the lower element concentration in foliage of coniferous than deciduous trees (Richardson and Friedland 2016). The rapid initial accumulation of organic matter in the organic layer on top of the mineral soil created an important reservoir of nutrients for the ecosystem, because the decomposition of organic matter mineralizes nutrients which can be reused for plant uptake (Lilienfein et al. 2001;Wilcke et al. 2002). Calcium availability in the forest soil has been shown to be directly related with tree growth (Baribault et al. 2010;Gradowski and Thomas 2008;Long et al. 1997) and K availability contributes to a higher N uptake (Stevens et al. 1993) and water-use efficiency (Bradbury and Malcolm 1977), promoting faster vegetation growth. The high K uptake in the early stage of the vegetation succession (0-47 yr) and the low K export with the stream at the older sites (47-127 yr) indicated that this nutrient was efficiently internally recycled in our study ecosystems. This together with the release of other essential nutrients from easily weatherable minerals, that we did not study, might have contributed to allow the likely N-limited forest to make full use of the fixed N 2 facilitating the fast establishment of a deciduous forest in the early stage of the chronosequence. The faster vegetation development than at many other glacier retreat chronosequences (Supplementary Table 1) might have been further supported by the higher mean annual temperature at the Hailuogou sites.

Conclusions
The development from bare soil to mature forest in < 80 years along the 127-year old Hailuogou glacier retreat chronosequence was linked with a strong soil acidification and an intense loss of carbonates in the first 47 years of vegetation succession. This development was associated with a decrease in total stocks of Ca and Mg, mainly because of the dissolution of carbonates and easily weatherable Mg-containing minerals, while the stocks of K and Na changed little.
Part of the released Ca and Mg was incorporated into the growing aboveground biomass and thereby retained in the ecosystem. The smaller K release than that of Ca and Mg in the early stage of the chronosequence was compensated by higher K inputs with atmospheric deposition, which were in the same order of magnitude than K losses through the stream. We suggest that the strong supply of Ca and Mg and a closed K cycling contributed to the establishment of a deciduous forest allowing the likely N-limited vegetation to make full use of the fixed N 2 .
The slower element release by weathering after leaching of carbonates occurred synchronously with the vegetation change from deciduous to coniferous forest with its lower nutrient demand, while climatic differences along the chronosequence were minor. The vegetation development resulted in the formation of an organic layer on top of the mineral soils with increasing thickness along the chronosequence that stored a large quantity of nutrients. The storage of K in the organic layers was particularly important, since this element was tightly cycled between soil and vegetation (i.e., via plant uptake and return to the soil via throughfall and litterfall). As a consequence, nutrient supply by weathering and nutrient demand by vegetation combined to promote the vegetation development along the Hailuogou chronosequence.