Topsoil organic matter build-up in glacier forelands around the world.

Since the last glacial maximum, soil formation related to ice-cover shrinkage has been one major sink of carbon accumulating as soil organic matter (SOM), a phenomenon accelerated by the ongoing global warming. In recently deglacierized forelands, processes of SOM accumulation, including those that control carbon and nitrogen sequestration rates and biogeochemical stability of newly sequestered carbon, remain poorly understood. Here, we investigate the build-up of SOM during the initial stages (up to 410 years) of topsoil development in 10 glacier forelands distributed on four continents. We test whether the net accumulation of SOM on glacier forelands (i) depends on the time since deglacierization and local climatic conditions (temperature and precipitation); (ii) is accompanied by a decrease in its stability and (iii) is mostly due to an increasing contribution of organic matter from plant origin. We measured total SOM concentration (carbon, nitrogen), its relative hydrogen/oxygen enrichment, stable isotopic (13 C, 15 N) and carbon functional groups (C-H, C=O, C=C) compositions, and its distribution in carbon pools of different thermal stability. We show that SOM content increases with time and is faster on forelands experiencing warmer climates. The build-up of SOM pools shows consistent trends across the studied soil chronosequences. During the first decades of soil development, the low amount of SOM is dominated by a thermally stable carbon pool with a small and highly thermolabile pool. The stability of SOM decreases with soil age at all sites, indicating that SOM storage is dominated by the accumulation of labile SOM during the first centuries of soil development, and suggesting plant carbon inputs to soil (SOM depleted in nitrogen, enriched in hydrogen and in aromatic carbon). Our findings highlight the potential vulnerability of SOM stocks from proglacial areas to decomposition and suggest that their durability largely depends on the relative contribution of carbon inputs from plants.

INTRODUCTION

Since the Last Glacial Maximum (ca. 20 kyr), more than 10% of the Earth's land surface area has been freed from its ice cover (Adams & Faure, 1998). In 2010, the total glacierized area (ca. 200 000 glaciers excluding the Greenland and Antarctic ice sheets) was estimated at 0.72 Mkm2 (Pfeffer et al., 2014; RGI n class="Chemical">Consortium, 2017), but ongoing climate change is accelerating glacier ice loss (Hock et al., 2019). As an example, in the European Alps, 25%–30% of ice cover disappeared over the past 60 years (Gardent et al., 2014; Smiraglia & Azzoni, 2015). Based on the Representative Concentration Pathway 8.5 scenario, mountain glaciers are expected to lose 37%–57% of their mass by 2100, and many will disappear even for lower emission scenarios (Hock et al., 2019). This accelerating ice shrinkage will lead to the expansion of novel terrestrial ecosystems conditioned by ecological succession (Clements, 1916), landform changes and soil development (Eichel et al., 2016).

Globally, soils are a major terrestrial reservoir of carbon (e.g., Gherardi & Sala, 2020; Jobbágy & Jackson, 2000). Soil development has led to the accumulation of significant amounts of n class="Chemical">carbon as organic matter in formerly glaciated areas (e.g., Albrecht, 1938; Schlesinger, 1995). Despite high uncertainties, up to one‐third (i.e., 490 GtC) of the current global soil organic carbon (SOC) has been sequestered in soils since the Last Glacial Maximum due to a large increase in the areas of forests and other SOC‐rich ecosystems (Adams & Faure, 1998; Adams et al., 1990). The rate of postglacial carbon sequestration as soil organic matter (SOM) strongly varies with time after barren substrate exposure. Net SOC accumulation rates are usually greatest during the early stages of soil development in proglacial areas, and also depend on climatic conditions as some authors have pointed out without drawing up a general law (Bockheim et al., 2000; Egli et al., 2012; Harden et al., 1992; Schlesinger, 1990).

By modeling postglacial SOC sequestration, Harden et al. (1992) emphasized the strong sensitivity of their simulations to SOC decomposition dynamics while regretting the lack of knowledge on the relative n class="Chemical">contributions of labile (fast cycling) and stable (slow cycling) SOC kinetic pools within soil chronosequences. Accumulation of SOM in glacier foreland soils with low net primary production and low rate of SOM degradation is also a key element of ecosystem development in the early stages of successions (Schulz et al., 2013; Wietrzyk et al., 2018; Yoshitake et al., 2018). During the past 15 years, researchers have examined the distribution of SOC in different kinetic pools as soils form in glacier forelands (e.g., Bardgett et al., 2007; Egli et al., 2012; Schweizer et al., 2018) using chemical (e.g., Egli et al., 2012), physical (e.g., Conen et al., 2007) or thermal (Bardgett et al., 2007) SOM fractionation methods, in some cases in combination with elemental isotopic signatures of SOM (15N, 13C, 14C; Bardgett et al., 2007; Smittenberg et al., 2012) or other SOM characterization techniques such as 13C nuclear magnetic resonance or Fourier transform mid‐infrared spectroscopy (Dümig et al., 2012; Egli et al., 2010). The main assumptions from these works are as follows: (i) a stable and ancient SOC fraction—whose origin is still debated—may be the primary source of carbon and energy for microbial communities during the first decades of soil development (Bardgett et al., 2007; Guelland et al., 2013), while nitrogen would mostly originate from atmospheric deposition (Smittenberg et al., 2012); (ii) organo‐mineral associations evolve in early stages of soil development, with SOC sequestration occurring at a faster rate than soil mineral weathering (Dümig et al., 2012; Schweizer et al., 2018). However, discrepancies remain among studies regarding the evolution of SOC stability during the development of proglacial soils and the build‐up of their SOM stocks. Indeed, Bardgett et al. (2007) showed that SOC stability decreased during soil development, whereas Egli et al. (2010) showed that the proportion of stable SOC increased and other researchers did not identify clear temporal changes in SOC stability (Conen et al., 2007; Egli et al., 2012). Such inconsistencies may relate to local environmental conditions. Moreover, the methods used for partitioning stable from labile SOC may also affect study outcomes, as the currently available SOM fractionation methods cannot reliably isolate SOC fractions that are unique and non‐composite in terms of carbon turnover times in soil (von Lützow et al., 2007; Smith et al., 2002). The same is true for thermal methods (e.g., Balesdent, 1996), which have, however, demonstrated that biogeochemically stable SOC was thermally stable (Barré et al., 2016; Gregorich et al., 2015; Plante et al., 2013; Sanderman & Grandy, 2020).

To date, most studies of the dynamics of SOM and the fate of SOC fractions in glacier forelands have focused on one or a few glacier forelands within the same area. In addition, most studies have been performed in the European Alps (e.g., Damma, Morteratsch and Ödenwinkelkees glaciers), the U.S. Rocky Mountains (e.g., Wind River Range) and the high Arctic region (Ny‐Ålesund) (e.g., Burga et al., 2010; Dümig et al., 2012; Eckmeier et al., 2013; Egli et al., 2012; Nakatsubo et al., 2005; Schurig et al., 2013; Schweizer et al., 2018), with very limited information from tropical areas and the Southern Hemisphere. As a result, we are constrained in our ability to draw general n class="Chemical">conclusions about the drivers of net SOM accumulation rates and the build‐up of SOC kinetic pools in proglacial environments.

In this paper, we use a global dataset on SOM dynamics during the initial stage (i.e., up to 410 years) of topsoil and ecosystem development in mountain proglacial areas. Specifically, we test three hypotheses: (i) the accumulation of SOM in glacier forelands at the early stage of topsoil development is affected by time and is accelerated by a warmer climate; (ii) irrespective of local n class="Chemical">conditions, there is a common pattern of decreasing SOM stability along soil chronosequences in glacier forelands (i.e., newly accumulated SOM in glacier forelands soils is mostly labile) and (iii) the accumulation of SOM in glacier forelands is mainly driven by an increasing contribution of organic matter from plant origin. To test these hypotheses, we examine the relationships between local climatic variables and the dynamics of SOM accumulation in soil chronosequences in 10 glacier forelands around the world, and evaluate the qualitative evolution of newly accumulated SOM over time by studying the thermal stability and the functional groups of organic carbon, the elemental stoichiometry of SOM and the stable isotopes of nitrogen and carbon in topsoils.

MATERIALS AND METHODS

Study sites, soil sampling and climatic data

This study is based on data from 10 glacier forelands spread over four continents (Figure 1): Glaciers Noir/Blanc (France), Forni (Italy), Gergeti (Georgia), Tiedemann (Canada), Charquini and Zongo (Bolivia), Antisana (Ecuador), Perito Moreno (Argentina), Lobuche (Nepal) and Apusinikajik (Greenland). Study sites enn class="Chemical">compass a wide range of geographic and climatic variables (i.e., latitude, elevation, mean annual temperature and precipitation). Soil sampling was performed between October 2014 and July 2017. In each glacier foreland, we sampled the topsoil from three to eight dated plots to obtain a chronosequence of soil development (see kmz file in Supplementary Information). Well‐dated soil chronosequences were established using dendrochronology, lichenometry, radiocarbon, optically stimulated luminescence techniques, photogrammetry and time‐series reconstruction based on satellite images and old maps (Table 1). The short timescale studied (<410 years after bedrock exposure) and the similar soil parent material within each glacier foreland (Table 1) ensure that the chronosequences are a suitable space‐for‐time substitution tool (Johnson & Miyanishi, 2008; Walker et al., 2010). If possible, plots were 100 × 2 m strips of glacier foreland, corresponding to the known position and shape of the glacier forefront at the specific dates. In a few cases, plots were narrower because of accessibility problems and/or if the undisturbed portion of the foreland was narrower.

FIGURE 1

Locations of the 10 glacier forelands of this study. Background map is modified from Randolph Glacier Inventory under an Attribution 4.0 International license (RGI Consortium, 2017)

TABLE 1

Glacier forelands sampled and bibliographic sources used to date moraines

Name	Mountain region	Country	Lat. (°E)	Long. (°N)	Elevation of glacier front (m a.s.l.)	Elevation of oldest sample site (m a.s.l.)	Lithology	Studied time period since deglacierization (years) and number of study sites	Data sources for soil chronology
Apusinikajik	Renland	Greenland	71.26	−25.82	75	55	Granite and gneiss	10–150; n = 3	Medford (2013)
Perito Moreno	South Andes	Argentina	−50.5	−73.04	180	340	Granite and granodiorite	100–410; n = 3	Aniya and Skvarca (2012)
Tiedemann	North Pacific Range	Canada	51.32	−124.923	950	815	Granodiorite and orthogneiss	36–116; n = 3	Larocque and Smith (2003)
Forni	Central European Alps	Italy	46.41	10.57	2600	2200	Granite	10–150; n = 7	Pelfini et al. (2014)
Glaciers Noir/Blanc	Western European Alps	France	44.92	6.41	2670	1890	Granite	14–166; n = 8	Cossart et al. (2006); Rabatel et al. (2008)
Gergeti	Greater Caucasus	Georgia	42.66	44.55	3220	2770	Andesite and dacite	15–150; n = 5	Tielidze et al. (2019)
Lobuche	Central Himalaya	Nepal	27.96	86.81	5100	5020	Black gneiss, metapelite and quartzite	20–300; n = 5	Richards et al. (2000)
Charquini	Central Andes	Bolivia	−16.31	−68.11	5070	4830	Granodiorite and granite	9–350; n = 7	Rabatel et al. (2005)
Zongo	Central Andes	Bolivia	−16.27	−68.13	4940	4830	Granite	9–351; n = 7	Rabatel (2005)
Antisana	Northern Andes	Ecuador	−0.47	−78.15	4870	4780	Andesite and volcanic ash	17–150; n = 5	Collet (2010)

For each plot, five topsoil samples (0–15 cm; 15 g each) were randomly collected at a distance of more than 20 m from one another on loose materials (glacial till) with no evident geomorphological activity. Sampling soils “by depth” rather than “by horizons” is indeed common practice procedure in soil monitoring networks and in other works focusing on soil organic matter build‐up in glacier forelands (e.g., Orgiazzi et al., 2018, Schweizer et al., 2018) even though with this approach, different thicknesses of the A layer between samples might impact the analyses of SOM. The five samples were mixed together to form a composite topsoil sample. For all topsoil samples, any pure organic surface layer, litter and vegetation was excluded from the sampling, but the eventual biological soil crust (i.e., the association between soil particles and soil microorganisms such as cyanobacteria, algae or fungi; e.g., Belnap et al., 2003) was included in the topsoil sample when present. All soils sampled in the 10 glacier forelands were classified as Leptosols (FAO, 2014) with mostly sandy texture, single grained structure for the youngest soils and granular structure for the most evolved soils. This approach was supported by field data showing that topsoil samples (0–15 cm) corresponded to C horizons (weathered soil parent material) on the most recent soils, and to the grouping of A (surface organo‐mineral soil horizon) and C horizons on the most evolved soils (see examples of soil profiles sampled in Figure S1). The C horizon was thus reached within the first 15 cm for all the studied soil profiles (aged up to 410 years). We did not describe any clearly defined B horizon in the studied soil profiles (Figure S1), similarly to Egli et al. (2012, 2001) and Dümig et al. (2011) in soils of glacier forelands of similar ages.

For each glacier foreland, we extracted two climate variables from the CHELSA dataset (Climatologies at High resolution for the Earth's Land Surface Areas; Karger et al., 2017), with a resolution of 30 arc‐seconds (approximately 1 km at the equator) averaged over the period 1979–2013: mean air temperature of warmest quarter (T, °C) and precipitation of warmest quarter temperature (P, mm month−1). The CHELSA dataset uses ERA‐Interim reanalyzes to downscale climate surface variables acn class="Chemical">counting for topography (Karger et al., 2017). Other global climate gridded datasets such as TerraClimate database (Abatzoglou et al., 2018) yield similar results.

Soil organic matter analysis

Total soil carbon and n class="Chemical">nitrogen concentrations were measured in each composite soil sample by elemental analysis (OEA Flash2000; ThermoFisher). As none of the foreland soils has a carbonate parent material, total soil carbon concentration corresponds to the total SOC concentration. 13C and 15N content were measured by isotope ratio mass spectrometry (ELEMENTAR Isoprime) and the results expressed in δ13C and δ15N abundance ratios (i.e., in parts per thousand (‰) relative to international standards. These indexes provide information about the origin of the SOM. In proglacial soils, high δ13C values have been linked to the presence of ancient carbon (Bardgett et al., 2007), although it also depends on the type of photosynthesis of autotrophic organisms that are present. Persistent SOC has been shown to have higher δ13C values in comparison to bulk SOC (Balesdent & Mariotti, 1996; Brüggemann et al., 2011; Menichetti et al., 2015). Negative δ15N values are generally associated with biological nitrogen fixation or atmospheric nitrogen deposition (Smittenberg et al., 2012).

The composition of n class="Chemical">organic carbon functional groups in topsoil samples was assessed by Fourier transform mid‐infrared spectroscopy (FTIR) in the attenuated total reflectance (ATR) mode. FTIR measurements were performed using a Nicolet iS10 spectrometer (Thermo Fischer Scientific) equipped with a diamond crystal ATR device. Finely ground samples were dried overnight (40°C) to standardize their water content prior to analysis. FTIR spectra of bulk soil were acquired over the 4000–650 cm−1 spectral range, with a spectral resolution of 4 cm−1 and 16 scans. Spectra were corrected for atmospheric interferences (H2O and CO2), and absorbance values were calculated as the inverse logarithm of the measured reflectance values. FTIR spectra were pre‐processed (offset‐correction using the 4000–3950 cm−1 spectral region as a reference, and setting the minimum absorbance value to zero). Three FTIR waveband regions were selected, corresponding to specific organic carbon functional groups, and relatively free from artifact absorptions by soil minerals in soils free of carbonates (adapted from Soucémarianadin et al., 2019): (i) the aliphatic (CH2, CH3 stretch) waveband region between 2930 and 2900 cm−1; (ii) the carbonyl and carboxyl (C=O stretch) waveband region between 1750 and 1670 cm−1; (iii) the aromatic carbon (C=C bond) waveband region between 1610 and 1590 cm−1. Relative ratios for these three regions were then calculated (FTIR C‐H; FTIR C=O; FTIR C=C) using Equation (1):

Areas of waveband‐regions were divided by their respective spectral region width (30 cm−1 for C‐H, 80 cm−1 for C=O, 20 cm−1 for C=C) prior to the calculation of relative ratios.

The organic matter bulk chemistry and thermal stability was characterized by Rock‐Eval® thermal analysis using a Rock‐Eval® 6 Turbo device (Vinci Technologies). This technique, which does not require any chemical pre‐treatment of the soil sample, involves the measurement of n class="Chemical">carbon as gaseous effluent during two phases (Behar et al., 2001; Cécillon et al., 2018; Disnar et al., 2003). The first phase is a pyrolysis stage (200–650°C) in a N2 atmosphere, during which CO and CO2 gases are quantified with an infrared detector and volatile hydrocarbon effluents (CH) are quantified using a flame ionization detector. The second phase is an oxidation stage (300–850°C) in a laboratory air atmosphere, during which CO and CO2 gases are quantified with an infrared detector.

Using the combination of the five Rock‐Eval® thermograms (Behar et al., 2001), we defined three SOC pools of different thermal stability (Figure S2): (i) a small and highly thermolabile pyrolyzable n class="Chemical">organic carbon pool (hereafter termed POC 1), corresponding to the carbon evacuated as CH, CO or CO2 at 200°C over 3 min during the pyrolysis phase; (ii) an intermediate pyrolyzable organic carbon pool (hereafter termed POC 2), corresponding to the carbon evacuated as CH, CO or CO2 during the temperature ramp‐up (30°C min−1) of the pyrolysis phase between 200 and 650°C and (iii) a pool of organic carbon resistant to pyrolysis (hereafter termed ROC), corresponding to the carbon evacuated during the oxidation phase as CO or CO2. As is the case with all SOM fractionation methods, thermal analysis does not isolate unique and non‐composite SOC kinetic pools (von Lützow et al., 2007), yet SOC thermal stability is positively correlated to SOC biogeochemical stability (Barré et al., 2016; Sanderman & Grandy, 2020). Therefore, the SOC fractions POC 1, POC 2 and ROC correspond to three SOC fractions containing increasing proportions of biogeochemically stable SOC. They are expressed as concentration (g C kg−1) or as percentage of total SOC. The thermal stability of the SOC fractions that are pyrolyzable or resistant to pyrolysis has been shown positively correlated to the proportion of persistent SOC (Cécillon et al., 2018). We thus additionally determined the temperature corresponding to the release of 50% and 90% of the carbon from the SOC fractions POC 2 and ROC (i.e., T50‐POC 2, T90‐POC 2, T50‐ROC, T90‐ROC; expressed in °C) to characterize more accurately the thermal stability of the carbon included in these two RE thermal SOC fractions. Higher values of these indexes indicate higher thermal stability of the considered thermal SOC fraction. For the two pyrolyzable SOC fractions (POC 1 and POC 2), we also quantified the proportion of SOC carbon released as volatile hydrocarbon effluents (CH) to estimate the SOM enrichment in hydrogen for POC 1 (CH‐POC 1, unitless) and POC 2 fractions (CH‐POC 2, unitless). These two indexes are negatively correlated to the oxygen content of SOM in the POC 1 and POC 2 fractions.

Low SOC concentrations of topsoil samples from glacier forelands present a challenge to the quality of the Rock‐Eval® signals. Soil samples with a SOC n class="Chemical">concentration below 1 g C kg−1 were discarded, as we considered those values below the detection limit of the Rock‐Eval®6 Turbo device. Similarly, we discarded soil samples with carbon yields below 75% and above 125% of the SOC concentration determined with elemental analysis (EA). With these two constraints, we retained Rock‐Eval® results for 37 topsoil samples, representing eight glacier forelands (all Rock‐Eval® results from soil samples of the Lobuche Glacier foreland in the central Himalaya and the Charquini glacier foreland in the Central Andes were discarded).

Statistical analysis

We used linear mixed models to assess the factors related to the evolution of SOM characteristics in glacier forelands (SOC, total nitrogen n class="Chemical">concentration, carbon and nitrogen stable isotope signatures [δ13C and δ15N], organic carbon functional groups [FTIR C‐H, FTIR C=O, FTIR C=C], Rock‐Eval® POC 1, POC 2 and ROC fractions [% of total SOC and g C kg−1], thermal stability of POC 2 and ROC fractions [T50‐POC 2, T90‐POC 2, T50‐ROC, T90‐ROC] and the CH proportion of the POC 1 and POC 2 fractions [CH‐POC 1 and CH‐POC 2]). Explanatory variables include time since barren substrate/debris exposure (soil age, year), mean air temperature (T, °C) and precipitation (P, mm) of warmest quarter, averaged over the period 1979–2013. Because the effect of climate on SOM might change through time, we tested the relationship between age and climatic variables. Non‐significant effects were removed from the final models. To take into account the heterogeneity between glacier forelands, the identity of the foreland was included as a random factor. In preliminary mixed models, we also tested a quadratic term (age2) to verify the linearity of the SOM/age relationship. Significant quadratic terms would indicate that relationships between SOM and age were nonlinear. However, the quadratic term was not significant in any the models, indicating a lack of significant deviation from the linear patterns. Therefore, in the final models, we only retained the linear relationship. For each dependent variable, we tested mixed models with both random intercept (RI) and random slope (RS), and selected the one with lowest Akaike's Information Criterion (AIC; Burnham et al., 2002; Schielzeth & Forstmeier, 2009). The RI model assumes that the relationship between age and SOM has the same slope, but different intercept across sites, whereas the RS model assumes that both intercept and slope vary across sites.

Additionally, we used mixed models to assess the evolution of SOM characteristics with SOC and total nitrogen n class="Chemical">concentrations. Statistical analyses were performed with R v.3.5 (R Core Team, 2015) and the MuMIn and lme4 libraries. Before running analyses, we log‐transformed soil age, temperature, precipitation, SOC and total nitrogen concentrations, the concentration of the different RE thermal SOC fractions, and T50 and T90 to improve normality and reduce skewness. The CH‐POC 1, CH‐POC 2 indexes, the proportions of the different RE thermal SOC fractions and the FTIR indexes were logit‐transformed for the same reasons. The variables C/N, δ15N, δ13C were normally distributed.

RESULTS

Organic carbon, total nitrogen concentrations and C/N ratio of topsoils

Soil organic carbon n class="Chemical">concentrations are very low (generally <2 g C kg−1) in all soils less than 25 years old, but thereafter increased linearly with age without reaching a plateau (Figure 2a; Table 2). Statistical analysis of the relationship between soil age and temperature showed that the rate of SOC accumulation through time is significantly different among glacier forelands. SOC concentration is highest in forelands with high temperature during the warmest quarter. Furthermore, a significant interaction between age and temperature indicates faster SOC accumulation in forelands with high temperatures during the warmest quarter of the year (Table 2). We found no significant relationship between SOC concentration and precipitation of the warmest quarter of the year (Table 2).

FIGURE 2

Plots of (a) soil organic carbon (SOC) concentration, (b) total N (Ntot) concentration, and (c) C/N ratio of topsoil samples versus time for the 10 soil chronosequences

TABLE 2

Results of general mixed models that assess relationships between the SOM characteristics and soil age, mean temperature of warmest quarter (T) and precipitation of warmest quarter (P). Two types of mixed models were tested: models with random intercept (RI) and random slope (RS). The table includes results only for the mixed models with the lowest AICs values (cells in gray correspond to the variables discarded). Symbols for p values: ***p < 0.001; **p < 0.01; *p < 0.05; † p < 0.1; NS >0.1. B values indicating the direction of the relationships are in brackets. Other models tested and detailed results are presented in Tables S1 and S3

	Model	Soil age	T	P	Age:T	Age:P	R2
SOC (g C kg−1)	RS	*** (0.85)	* (0.70)		* (0.26)		0.51
Ntot (g N kg−1)	RI	*** (0.66)	† (0.47)		** (0.27)	† (0.14)	0.41
C/N	RI	*** (2.39)	* (1.88)		† (0.56)	† (−0.46)	0.54
δ13C (‰)	RS		** (−0.82)				0.20
δ15N (‰)	RI	** (0.71)	* (−0.82)		† (0.49)		0.26
POC 1 (g C kg−1)	RI	** (0.32)	* (0.38)	* (0.37)			0.36
POC 2 (g C kg−1)	RI	*** (0.90)	* (0.67)	* (0.64)	* (0.43)	† (0.41)	0.57
ROC (g C kg−1)	RI	*** (0.95)	* (0.68)	* (0.60)	* (0.45)	NS	0.62
POC 1 (% of total SOC)	RI	*** (−0.48)	† (−0.31)	* (−0.35)			0.50
POC 2 (% of total SOC)	RI	† (−0.06)				* (0.09)	0.17
ROC (% of total SOC)	RI	*** (0.13)	† (0.05)			** (−0.10)	0.38
T50‐POC 2 (°C)	RI	*** (−0.03)					0.33
T50‐ROC (°C)	RI	*** (−0.02)					0.08
CH‐POC 1	RI	*** (−0.39)	NS		* (−0.37)	NS	0.30
CH‐POC 2	RI	* (0.18)		† (0.13)	* (0.29)	** (0.36)	0.38
FTIR C=C	RS	** (0.04)	* (0.02)		* (0.03)		0.36
FTIR C=O	RS	* (−0.02)	† (0.02)		* (−0.01)		0.29
FTIR C‐H	RI	*** (−0.02)			** (−0.02)		0.29

Plots of (a) soil organic carbon (SOC) n class="Chemical">concentration, (b) total N (Ntot) concentration, and (c) C/N ratio of topsoil samples versus time for the 10 soil chronosequences

Results of general mixed models that assess relationships between the SOM characteristics and soil age, mean temperature of warmest quarter (T) and precipitation of warmest quarter (P). Two types of mixed models were tested: models with random intercept (RI) and random slope (RS). The table includes results only for the mixed models with the lowest AICs values (cells in gray n class="Chemical">correspond to the variables discarded). Symbols for p values: ***p < 0.001; **p < 0.01; *p < 0.05; † p < 0.1; NS >0.1. B values indicating the direction of the relationships are in brackets. Other models tested and detailed results are presented in Tables S1 and S3

Nitrogen n class="Chemical">concentration is very low on soils younger than 25 years (between 0.04 and 0.46 g N kg−1; Figure 2b). Nitrogen exceeded 1 g N kg−1 in topsoils about 50 years post‐deglacierization for mountain soil chronosequences. The Zongo, Gergeti and Charquini forelands show soil nitrogen concentrations above 1 g N kg−1 only after about 150 years since deglacierization. Soil nitrogen concentration is strongly related to soil age, but with a less steep slope than for SOC concentration (Table 2). We did not detect a significant slowdown in soil nitrogen accumulation over the first 410 years of soil development. Moreover, the rate of soil nitrogen accumulation is significantly faster in forelands with warmer climates (Table 2).

The C/N ratio of topsoils is generally very low (<5) at the earliest stage of soil development, but it increases with soil age (Figure 2c). At Perito Moreno, Apusinikajik, Zongo and Charquini, the C/n class="Chemical">N ratio of topsoils reaches a plateau about 100 years after deglacierization, with values of about 10, whereas in the others forelands the C/N ratio continues to increase to 14–21 after 100 years (Tiedemann, Forni, Gergeti and Noir/Blanc). The model testing a nonlinear relationship between C/N and soil age is not significant (p = 0.69), suggesting a linear increase through time. Glacier forelands with warmer temperatures show significantly higher soil C/N ratios and a faster increase in soil C/N ratio with soil age than those with cooler temperatures (Table 2).

Stable carbon and nitrogen isotopes and carbon functional groups in soil organic matter

The δ13C signature of SOM differs among the glacier forelands, especially in topsoils with SOC n class="Chemical">concentration below 3 g C kg−1 where it is ranging from −30‰ to −21‰ (Figure 3a). Soil δ13C is significantly lower in warmer than cooler forelands (Table 2). We found no significant relationship between δ13C and soil age (Table 2), but there is a significant negative relationship between the δ13C and SOC concentration (Figure 3a): soil δ13C is higher on low‐carbon soils and decreases significantly as SOC increases (B = −0.37, F = −2.50, df =52.20, p = 0.016).

FIGURE 3

Plots of (a) soil δ13C versus soil organic carbon (SOC), and (b) soil δ15N evolution versus Ntot

Plots of (a) soil δ13C versus soil n class="Chemical">organic carbon (SOC), and (b) soil δ15N evolution versus Ntot

The δ15N signature of SOM also differs n class="Chemical">considerably among the glacier forelands (Figure 3b), ranging from −4‰ to +4‰ on average. Soils from the Gergeti, Noir/Blanc, Tiedemann, Forni and Apusinikajik forelands have the lowest δ15N values (around −4‰); those from the Charquini, Zongo, Perito Moreno, Lobuche and Antisana forelands are greater than −2‰. Soil δ15N values increase significantly with soil age, and the signature is lowest in forelands with warmer temperatures (Table 2). There is also a weak positive relationship between δ15N values and the total concentration of nitrogen (B = 0.42, F = 1.76, df =52.33, p = 0.08).

The FTIR signature of SOM changes with soil age and SOC concentration. The aromaticity (FTIR C=C) of SOM increases significantly with soil age and SOC n class="Chemical">concentration, whereas the SOM content in aliphatic carbon (FTIR C‐H), carbonyl/carboxyl functional groups (FTIR C=O) decrease (Table 2; Table S2). SOM content in aromatic carbon (FTIR C=C) is positively related to the temperature of the warmest quarter. The analysis of the relationship between soil age and temperature showed that the rate of increase in the aromatic carbon content of SOM through time is significantly different among glacier forelands, with the fastest increase in SOM aromaticity in the warmest forelands (Table 2).

Soil organic carbon thermal stability

The summary statistics of the Rock‐Eval® carbon yield (% of SOC n class="Chemical">concentration determined by elemental analysis) for the 37 retained topsoil samples are as follows: mean =96%; median =94%; minimum =77%; maximum =124% (Table S3). The average concentrations of thermal SOC fractions POC 1, POC 2 and ROC are, respectively, 0.1, 1.2 and 1.8 g C kg−1 between 9 and 50 years post‐deglacierization and, respectively, 0.1, 5.6 and 8.5 g C kg−1 between 150 and 220 years post‐deglacierization. The absolute amounts of all these pools increase significantly with soil age and are highest in forelands with highest temperatures and precipitation (Table 2). The proportions of POC 1, POC 2 and ROC thermal SOC fractions are 3.2%, 40.8% and 56.0% of total SOC, respectively, between 9 and 50 years post‐deglacierization while they represent 1.3%, 39.2% and 59.5% of total SOC, respectively, from 150 to 220 years post‐deglacierization. The proportion of the POC 1 fraction significantly decreases with soil age, whereas the contribution of the ROC fraction significantly increases with soil age (Table 2; Figure 4a,b). The proportion of POC 1 decreases with SOC accumulation, whereas the other fractions did not change significantly (Table S2).

FIGURE 4

Relationships between soil age and (a) POC 1 (% of total SOC), (b) ROC (% of total SOC), (c) T50‐POC 2, (d) T50‐ROC, (e) CH‐POC 1 index, and (f) CH‐POC 2 index for eight proglacial soil chronosequences. CH, hydrocarbon effluents; POC, pyrolyzable organic carbon; SOC, soil organic carbon

Relationships between soil age and (a) POC 1 (% of total SOC), (b) ROC (% of total SOC), (c) T50‐POC 2, (d) T50‐ROC, (e) CH‐POC 1 index, and (f) CH‐POC 2 index for eight proglacial soil chronosequences. CH, hydrocarbon effluents; POC, pyrolyzable n class="Chemical">organic carbon; SOC, soil organic carbon

We detected a significant decrease in Rock‐Eval® T50 and T90 indexes with soil age and with SOC concentration, which indicates a general decrease in the thermal stability of the POC 2 and ROC fractions (Table 2; Table S2; Figure 4c,d). Furthermore, we found a significant decrease in the proportion of CH with soil age and SOC in the POC 1 fraction. In n class="Chemical">contrast, we found that the proportion of CH in the POC 2 fraction increases with soil age and SOC (Table 2; Table S2; Figure 4e,f).

DISCUSSION

In all the glacier forelands we studied, the amount of organic matter in topsoils increased with time. This SOM build‐up seems significantly modulated by climate: a warmer climate accelerates SOM accumulation. We also found a general pattern of decreasing SOM stability along the soil chronosequences, observed over all the glacier forelands studied. Then, the observed changes in SOM elemental stoichiometry, aromaticity and stable isotope signature with SOM accumulation, suggest an increasing n class="Chemical">contribution of organic matter from plant origin during the first centuries of topsoil development.

Temporal and climatic trends of organic matter accumulation in topsoils of glacier forelands

The measured concentrations of SOC and total n class="Chemical">nitrogen in the 10 glacier forelands are in good agreement with previous observations. For example, similar SOC concentrations were reported for topsoils in the Morteratsch glacier foreland in Switzerland (i.e., about 3 g C kg−1 around 30 years; Eckmeier et al., 2013), and similar very low nitrogen concentrations have been reported for the Zongo and Charquini forelands in Andes (i.e., <1 g N kg−1 for soils less than 80‐year old; Schmidt et al., 2008). For all the glacier forelands we studied, time since exposure of the substrate was the major driver of the variation in the concentrations of SOC and total nitrogen over the first 410 years of topsoil development (Table 2). Some studies analyzing the evolution of organic matter quantity in glacier foreland soils showed a nonlinear relationship between carbon and nitrogen concentrations or stocks and soil age (e.g., Darmody et al., 2005; Harden et al., 1992), with a slowdown of organic matter accumulation after several hundred years (300–700 years), and a plateau after several hundred or thousand years (Bockheim et al., 2000; Dümig et al., 2011; Egli et al., 2001, 2012; Harden et al., 1992; He & Tang, 2008; Mavris et al., 2010). Our results show linear, non‐saturating, SOC accumulation during the early stages of the soil development. The short timescale of our study probably precluded the detection of later slowdowns of SOM accumulation, but it suggests a more regular pattern of organic matter accumulation during the early stages (first 2–4 centuries) of establishment of the new ecosystem.

Contrary to previous reports (Göransson et al., 2011; Smittenberg et al., 2012), the C/N ratio of SOM increases with soil age at all sites. This result suggests that SOC storage does not require n class="Chemical">nitrogen storage equivalent to the initial C/N of SOM, as shown by Erktan (2013) in the early stage of topsoil development in different ecosystems, with young topsoils starting at a similar C/N ratio of ca. 5–6. Therefore, the low amount of nitrogen in topsoil is not a limiting factor for accumulation of SOM, nor probably for SOC storage in the early stages of glacier foreland topsoil development.

We found a positive relationship between temperature of the warmest quarter and SOM accumulation dynamics. The relationship between SOM accumulation and climate in glacier forelands has previously been observed and attributed to the effect of temperature and moisture (Bockheim et al., 2000; Egli et al., 2012). However, to date, the lack of in situ measurements has hampered a more in‐depth analysis of the local climate drivers of SOM accumulation in glacier forelands. Even if global gridded datasets cannot precisely capture the local climates, our findings provide evidence that temperature has a major role on SOM accumulation in glacier forelands. This inference is consistent with results of studies of high‐elevation en class="Chemical">cosystems, showing that primary productivity is primarily linked to temperature (Körner, 2003). Further analyses should investigate the interplay between temperature, precipitation and soil moisture on SOM accumulation, preferably using locally observed values rather than global gridded data sets.

The effect of temperature during the warmest quarter on SOM accumulation can be explained by an increase in SOM input. Notably, this positive effect is not offset by the increased heterotrophic activity and SOM mineralization expected under a warmer climate (e.g., Guelland et al., 2013). Nevertheless, a better quantification of soil carbon fluxes (aboveground and belowground inputs, and outputs by dissolved n class="Chemical">organic carbon leaching and heterotrophic respiration) along each soil chronosequence would be required to refine our understanding of SOM accumulation processes. To overcome the difficulties of assembling field data at the global scale, further studies could investigate the potential of using high‐resolution remote sensing data to develop proxies of aboveground and belowground primary productivity in glacier forelands (Fischer et al., 2019; Gherardi & Sala, 2020). Finally, it is important to note that we focused on SOM concentration in topsoils and that patterns could be different for SOM stocks.

Changes in SOC stability during topsoil development and SOM accumulation

Despite the diversity of climates in our studied glacier forelands, we detected consistent and general patterns in the build‐up of SOC fractions with different stability during topsoil development and the associated accumulation of SOM (Figure 5). First, the thermal stability of SOC decreased with soil age and with SOC n class="Chemical">concentration, as shown for the POC 2 and the ROC fractions (representing ca. 97–99% of total SOC; Figures 4c,d and 5). This suggests that the overall biogeochemical stability of SOC decreases as the soil age increases and SOC is sequestered in the soil. The accumulation of SOC during the first centuries of topsoil development in glacier forelands is thus driven preferentially by the accumulation of labile SOC. The decrease in SOC biogeochemical stability with increasing SOC concentration is consistent with the associated δ13C decrease (Figure 3a), as previously reported for glacier forelands in the Alps and Alaska (Bardgett et al., 2007; Guelland et al., 2013; Malone et al., 2018). It has been shown that SOC biogeochemical stability increases with δ13C signature of SOM (Balesdent & Mariotti, 1996; Brüggemann et al., 2011; Menichetti et al., 2015). However, we observed no correlation between δ13C and soil age (Table 2), as previously reported in an Alpine foreland (Smittenberg et al., 2012; but see Guelland et al., 2013, for contradictory results from the same foreland), highlighting that the stable carbon isotope signature of SOM is complex and does not solely reflect its biogeochemical stability.

FIGURE 5

Overall trends in the biogeochemical signatures of soil organic matter during its accumulation in recently deglacierized (up to four centuries) topsoils

The decrease in SOC stability with age and SOC concentration is also in line with previous observations of a decrease in SOC thermal stability during the first 145 years of soil development in a glacier foreland in the European Alps (Bardgett et al., 2007). In n class="Chemical">contrast, Egli et al. (2010), who used a chemical fractionation technique in their study, did not confirm this relationship. These inconsistencies are probably due to the use of different analytical methods. It should be noted, however, that chemical oxidation (e.g., hydrogen peroxide treatment) is now not recommended to isolate the stable SOC kinetic pool (Lutfalla et al., 2014; Poeplau et al., 2019). All SOM fractions isolated by physical, chemical or thermal methods are a mixture of labile and stable carbon (Balesdent, 1996; von Lützow et al., 2007). In our study, the additional information on thermal stability (T50 and T90 indexes) of the two main isolated thermal SOC fractions (POC 2 and ROC thermal fractions) helped reducing biases inherent to SOM thermal fractionation methods. Indeed, the slight increase in the proportion of the thermally stable ROC fraction (ca. +3.5%; observed with soil age, not with SOC concentration; Table 2; Table S2) does not mean that SOC accumulation is associated with an increase in SOC biogeochemical stability, as evidenced by the consistent decrease in Rock‐Eval® T50 and T90 indexes observed with soil age and SOC concentration (Table 2; Table S2). This slight increase in the proportion of the ROC fraction might be explained by the progressive accumulation of SOM of plant origin, which is rich in cellulose and lignin, two compounds with a higher ROC fraction compared to other SOM compounds such as lipids and proteins (Carrie et al., 2012). The increase in CH‐POC 2 with soil age and SOC concentration also highlights the accumulation of labile organic matter enriched in hydrogen moieties over time (Barré et al., 2016; Gregorich et al., 2015; Poeplau et al., 2019; Saenger et al., 2015; Soucémarianadin et al., 2018). Finally, the consistent decrease in SOC stability with SOC accumulation in glacier forelands supports the increasing body of evidence that SOC storage in most terrestrial ecosystem is driven by the accumulation of labile SOC (see e.g., Barré et al., 2017, and references therein). This labile SOC fraction is considered potentially more vulnerable to microbial decomposition and thus to loss (Cotrufo et al., 2019; Viscarra Rossel et al., 2019).

Second, we detected a decreasing proportion of the thermally highly labile POC 1 fraction in the early stages of soil development (from ~3.2% of total SOC, 9–50 years post‐deglacierization to ~1.3% in the oldest soils; Figures 2a and 5). In comparison, in "mature" alpine grassland topsoils (0 to 5–10 cm) in the Grand Galibier massif (SW French Alps), the average proportion of the POC 1 fraction was just 0.006% (unpublished results). The high proportion of labile SOC in the initial stages of new soils compared to "mature" ones has been previously reported and interpreted in different ways. It might be due to a particularly high soil surface community consisting of algae, fungi and bacteria, called biological soil crust (Breen & Lévesque, 2008; Yoshitake et al., 2010), or to an important necromass consisting of dead cell‐envelope fragments from autotrophic or heterotrophic microbes (Bardgett et al., 2007; Schurig et al., 2013) that could contribute to SOM during the initial stage of soil development. It might also be due to the higher proportion of an easily mineralizable SOC fraction, such as dissolved organic carbon produced within the glacier (Guelland et al., 2013; Yoshitake et al., 2018), or to an easily mineralizable fraction of fossil organic carbon (hydrogen‐rich and thermally labile; Copard et al., 2006; Graz et al., 2011). The decrease in the CH‐POC 1 index with soil age and with SOC concentration corresponds to a decrease in thermally labile hydrogen‐rich, and a relative increase in oxygen‐rich, organic compounds with age in the POC 1 fraction. It has been shown that the microbial biomass of some taxa (e.g., cyanobacteria) is enriched in hydrogen moieties (Carrie et al., 2012); thus, the decrease in the CH‐POC 1 index with soil age might reflect a temporal decrease in the specific soil microbial biomass, or the mineralization of the above‐mentioned fraction of fossil organic carbon (Copard et al., 2006; Graz et al., 2011). However, due to the limited chemical information provided by Rock‐Eval®, we cannot assess the exact chemical nature or origin of the POC 1 thermally labile SOC pool.

Origin of soil organic matter

The C/N ratio and the relative n class="Chemical">composition of SOM in organic carbon functional groups provide information on the origin of SOM, regarding the contribution of plant material and microorganisms. The C/N ratio significantly increased with soil age (Figures 2c and 5), just as the SOM content in aromatic carbon as previously shown by Egli et al. (2010) in the Swiss Alps (Table 2; Figure 5). A SOM C/N ratio of about 6, observed in recently deglacierized topsoils, is characteristic of SOM of microbial origin (Paul & Clark, 1996) but does not provide information on its age (see e.g., Graz et al., 2011, for observed low C/N ratios of fossil SOM). After 100 years of soil development, the C/N ratio of SOM reached a value around 10 for certain soil chronosequences, matching the typical C/N of surface soils (e.g., 11.6; Kirkby et al., 2011). At older sites, the higher C/N ratio, the higher aromaticity (with FTIR C=C ratio approaching values typical from the particulate organic matter fraction of topsoils; Soucémarianadin et al., 2019), and the CH‐POC 2 index of SOM suggest a transition toward a SOM largely derived from plants (Figure 5). This transition occurred faster in glacier forelands with warm climates, where vegetation colonization and growth tend to be more rapid than in forelands located in cold climates (Fickert et al., 2017).

Stable isotopes of carbon and n class="Chemical">nitrogen can also provide useful, although not unambiguous, information on the origin and cycling of the SOM (Brüggemann et al., 2011; Craine et al., 2015; Malone et al., 2018; Whiticar, 1996). High values of δ13C in glacier foreland topsoils are generally interpreted, when C4 plants are absent, as the presence of ancient carbon that could come from relics of soils of previous interglacial cycles, or cryoconites soot and dust that have been deposited on the glacier (Baccolo et al., 2017; Bajerski & Wagner, 2013; Guelland et al., 2013; Sattin et al., 2009; Stubbins et al., 2012). High δ13C values of SOC and topsoil CO2 effluxes have been correlated to low 14C signals corresponding to ancient carbon in soils from two glacier forelands in the Alps (Bardgett et al., 2007; Guelland et al., 2013). Our results showed a high δ13C signature, especially in young topsoils (Figures 3a and 5). Along with this 13C enrichment, we note that the thermal stability of the POC 2 and ROC fractions is higher in younger topsoils. Although ancient organic carbon's Rock‐Eval® signature is quite variable (containing both thermally labile and thermally stable carbon), some studies have demonstrated the relevance of this analysis for detecting ancient organic carbon (e.g., Copard et al., 2006; Vindušková et al., 2015). Therefore, our results suggest that a small (from ca. 0.5 to 2 g C kg−1; Figure 2a), ancient and thermally stable SOC fraction is present in most of the young proglacial topsoils that developed in the early stages of glacier retreat. The use of this stable SOM fraction as a substrate for heterotrophic soil microorganisms is likely (Bardgett et al., 2007), but not strictly demonstrated here. Finally, the significant decrease in δ13C with increased SOC concentration (Figure 3a) is an additional evidence suggesting that newly accumulated SOM in topsoils is mostly of plant origin (C4 plants being absent from glacier forelands ecosystems).

The relative depletion of 15N in SOM indicates that a significant proportion of n class="Chemical">nitrogen can come from atmospheric deposition (Handley et al., 1999; Lehmann et al., 2004) or microbial fixation (Boddey et al., 2000). Both mechanisms can coexist, but more negative values of the δ15N signature were systematically found in forelands from the northern hemisphere. Similarly, Smittenberg et al. (2012) observed comparable rather low δ15N values (i.e., −4‰) in topsoils from a glacier foreland in the Alps. These observations are consistent with inorganic nitrogen deposition documented over large areas of Europe and North America (Galloway et al., 2004; Holland et al., 1997), and support the hypothesis of a significant nitrogen contribution from atmospheric depositions. This nitrogen deposition is probably progressively buffered as soils age through the arrival of plant‐derived nitrogen, including isotope fractionation during SOM decomposition, which explains the δ15N increase with time (Figures 3b and 5; Amundson et al., 2003; Craine et al., 2015; Malone et al., 2018).

In conclusion, our results indicated highly n class="Chemical">consistent patterns of SOM build‐up in topsoils of glacier forelands at the global scale. The rate of SOM accumulation in topsoils is accelerated by higher temperature during the warmest quarter, suggesting that climate and time are key drivers of SOM build‐up during the initial stages of topsoil development after glacier retreat. The increase in the C/N ratio in SOM with soil age at all sites indicates that SOC accumulation occurred in spite of a slower, yet significant nitrogen storage in topsoil. Furthermore, a highly stable and possibly ancient SOC fraction can act as starting point for the initial SOM build‐up, providing a key source of energy for early soil food webs. Finally, the general decrease in SOC biogeochemical stability and the general increase in SOM aromaticity indicate that the SOM newly accumulated in glacier forelands topsoils is mostly labile and of plant origin. This highlights the potential vulnerability of SOC stocks from proglacial areas to decomposition, and suggests that their maintenance in a warmer world largely depends on increased soil carbon inputs from plants.

Supplementary Material

Click here for additional data file.