Depletion of soil carbon and aggregation after strong warming of a subarctic Andosol under forest and grassland cover

The net loss of soil organic carbon (SOC) from terrestrial ecosystems is a likely consequence of global warming and may affect key soil functions. The strongest changes in temperature are expected to occur at high northern latitudes, with forest and tundra as prevailing land cover types. However, specific soil responses to warming in different ecosystems are currently understudied. In this study, we used a natural geothermal soil warming gradient (0–17.5 ^∘C warming intensity) in an Icelandic spruce forest on Andosol to assess changes in the SOC content between 0 and 10 cm (topsoil) and between 20 and 30 cm (subsoil) after 10 years of soil warming. Five different SOC fractions were isolated, and their redistribution and the amount of stable aggregates were assessed to link SOC to changes in the soil structure. The results were compared to an adjacent, previously investigated warmed grassland. Soil warming depleted the SOC content in the forest soil by −2.7 g kg^−1 ∘C⁻¹ (−3.6 % ^∘C⁻¹) in the topsoil and −1.6 g kg^−1 ∘C⁻¹ (−4.5 % ^∘C⁻¹) in the subsoil. The distribution of SOC in different fractions was significantly altered, with particulate organic matter and SOC in sand and stable aggregates being relatively depleted and SOC attached to silt and clay being relatively enriched in warmed soils. The major reason for this shift was aggregate breakdown: the topsoil aggregate mass proportion was reduced from 60.7±2.2 % in the unwarmed reference to 28.9±4.6 % in the most warmed soil. Across both depths, the loss of one unit of SOC caused a depletion of 4.5 units of aggregated soil, which strongly affected the bulk density (an R² value of 0.91 and p<0.001 when correlated with SOC, and an R² value of 0.51 and p<0.001 when correlated with soil mass in stable aggregates). The proportion of water-extractable carbon increased with decreasing aggregation, which might indicate an indirect protective effect of aggregates larger than 63 µm on SOC. Topsoil changes in the total SOC content and fraction distribution were more pronounced in the forest than in the adjacent warmed grassland soils, due to higher and more labile initial SOC. However, no ecosystem effect was observed on the warming response of the subsoil SOC content and fraction distribution. Thus, whole profile differences across ecosystems might be small. Changes in the soil structure upon warming should be studied more deeply and taken into consideration when interpreting or modelling biotic responses to warming.

1 Introduction

Global warming is inexorably progressing, with the largest expected changes to occur in the high northern latitudes (Diffenbaugh and Giorgi, 2012). The Intergovernmental Panel on Climate Change (IPCC) worst-case scenario (RCP8.5) predicts an air temperature increase of up to 11 ^∘C in areas north of 60^∘ latitude by the end of this century (IPCC, 2013). This will lead to strong responses from ecosystems, one of which will be increased microbial activity and, thus, oxidation of carbon (Melillo et al., 2002). Predicted alterations in soil organic carbon (SOC), as the largest terrestrial carbon (C) pool (Scharlemann et al., 2014), are inducing a positive climate–carbon cycle feedback loop. The highest SOC stocks are located in high northern ecosystems (Tarnocai et al., 2009). This spatial coherence of the strongest warming and the highest SOC stocks is expected to turn the vast land masses at high northern latitudes into a major C source. The simple extrapolation of short-term soil warming experiments has predicted a global SOC loss of up to 203±161 Pg C with a 1 ^∘C warming by 2050 (Crowther et al., 2016), which equals one-fourth of the current atmospheric C pool. More conservative estimates from the same authors still predict losses of 55±50 Pg C. This range of possible SOC changes, as well as the large standard errors associated with each of the estimates, points towards the high uncertainty of potential changes in carbon fluxes from terrestrial ecosystems to the atmosphere (van Gestel et al., 2018).

One of the major uncertainties in predicting SOC responses to warming is due to an incomplete mechanistic understanding of the temperature sensitivity of different functional SOC pools. For example, owing to different methodological approaches and partly also to misinterpretations (Conant et al., 2011), slow-cycling SOC is found to be more (Lefevre et al., 2014) or equally (Fang et al., 2005) sensitive to warming compared with fast-cycling SOC. As a consequence, SOC models frequently use the same temperature sensitivity for all SOC functional pools. However, it has been recently suggested that the implementation of carbon turnover and stabilisation in many models is outdated (Bradford et al., 2016) and that more wholistic experimental knowledge on warming-induced mechanisms related to carbon turnover in soils is necessary (Conant et al., 2011). Thus, isolated quantifications of CO₂ fluxes, bulk SOC or even SOC fractions might not yield enough insights to understand and predict SOC dynamics under global warming. Furthermore, individual soil warming experiments are mostly restricted to one ecosystem type and differ strongly with respect to methodology, i.e. the type and degree of warming. Therefore, comparisons across ecosystems are hampered (Crowther et al., 2016), but they might be critically important to (i) foster the understanding of underlying processes driving the SOC responses to warming and (ii) inform land surface models to increase their accuracy.

Apart from its significant role in the global carbon cycle, soil organic matter has numerous functions related to soil fertility and soil health: it is an important food source for soil biota (Barrios, 2007), contains and binds major plant nutrients and trace elements, has a large water storage capacity and is directly linked to soil structure, i.e. the three-dimensional arrangement of soil particles and pore space (Larsbo et al., 2016). The soil structure drives water and gaseous fluxes through the soil matrix, root growth and nutrient uptake, and the susceptibility of soils to compaction and erosion (Johnston et al., 2009; Chepil, 1951; Horn et al., 1994). Thus, in addition to the enrichment of atmospheric CO₂, soil carbon loss upon warming might also deteriorate the soil quality, with potential consequences for net primary production. To date, such effects, and the mechanisms involved, have been shown little attention, which might be related to the fact that most warming experiments have only been run for a relatively short period of time and with moderate warming treatments (Rustad, 2001; Conant et al., 2011).

In essence, long-term multi-ecosystem warming studies with strong soil warming gradients that might even exceed realistic temperature changes are ideal for advancing our understanding of carbon cycling and related changes in soil function under global change (Kreyling et al., 2014). An experiment such as this has been established in southern Iceland, where an earthquake shifted geothermal channels within the bedrock in 2008, resulting in strong gradients in soil warming (up to ∼80 ^∘C) in previously unwarmed grassland and forest soils. A growing community of scientists has been investigating the warming effects in permanent monitoring plots on virtually all ecosystem aspects since 2013 (https://www.forhot.is, last access: 3 February 2020). In a previous study, Poeplau et al. (2017) quantified the effect of soil warming on bulk SOC and five different SOC fractions with distinct turnover rates in the unmanaged grassland soil. The authors found a strong decline in soil mass and C in the stable aggregate fraction, indicating that either (i) warming-induced SOC depletion led to a destabilisation of aggregates or (ii) warming-induced aggregate breakdown led to a destabilisation of SOC.

In this study, we isolated the identical SOC fractions from an equally warmed adjacent forest soil to (i) advance our understanding of the temperature response of different SOC fractions representing kinetic pools, (ii) assess the role of the ecosystem type in the temperature response of SOC, and (iii) investigate potential links between SOC loss and soil structure changes.

2 Materials and methods

2.1 Study site and experimental design

In May 2008, a major earthquake in southern Iceland affected geothermal channels close to its epicenter (Halldórsson and Sigbjörnsson, 2009). Due to this event, a geothermal system in Reykir, close to the village of Hveragerði (64.008^∘ N, 21.178^∘ W) was moved to a previously unwarmed area, which is now constantly warmed with strong temperature gradients of up to ∼80 ^∘C (O'Gorman et al., 2014). This recently warmed area is covered by a Sitka spruce forest (Picea sitchensis (Bong.) Carr.) that was planted in 1966 and an adjacent unmanaged treeless grassland dominated by common bent (Agrostis capillaris, L.). These two ecosystems are located on a southwest sloping hillside (83–163 m a.s.l.). The mean annual temperature and precipitation between 2003 and 2015, as measured at the closest weather station, were 5.2 ^∘C and 1457 mm respectively (Sigurdsson et al., 2016). According to the world reference base, the soil is characterised as a silandic Andosol with a silt loam texture (clay : silt : sand ratio of 8 : 61 : 31 in the forest and 6 : 53 : 41 in the grassland) (Sigurdsson et al., 2016). The soil pH is slightly acidic (5.3 in the forest and 5.7 in the grassland), and the average SOC content at a soil depth of between 0 and 10 cm in the unwarmed soils is 75 g C kg⁻¹ in the forest (present study) and 54 g C kg⁻¹ in the grassland (Poeplau et al., 2017). Between autumn 2012 and spring 2014, a total of 30 permanent plots were installed in each ecosystem, comprising six different degrees of warming along five different transects. In 2014, the permanently monitored average soil temperature changes due to geothermal warming were 0, 1.0, 1.9, 2.7, 5.8 and 17.5 ^∘C in the forest and 0, 0.5, 2.1, 3.9, 10.5 and 17.3 ^∘C in the grassland (Sigurdsson et al., 2016).

2.2 Soil sampling, fractionation and analysis

In late April 2018, i.e. almost exactly 10 years after the warming was initiated, mineral soils of all permanent forest plots (six warming intensities and five replicates each) were sampled. Before sampling, the litter layer was carefully removed. Sampling was carried out with a thin auger (3 cm diameter) to a depth of 30 cm in direct proximity of the plot. For each plot, three individual soil cores were taken, split into 0–10, 10–20 and 20–30 cm depth increments, and pooled per depth. In case of soil compaction within the auger, the increment depth was adjusted linearly. For example, a compaction of 3 cm over the whole soil core resulted in the sampling of 0–9, 9–18 and 18–27 cm increments. For this study, only the 0–10 and 20–30 cm depth increments were used, which will hereafter be referred to as “topsoil” and “subsoil”. After sampling, soils were oven dried at 40 ^∘C and sieved to less than 2 mm.

Fractionation of SOC was performed as initially described by Zimmermann et al. (2007) and refined by Poeplau et al. (2013). A scheme can be found at https://www.somfractionation.org/combined-meth/part-dens-oxid-zimmermann/ (last access: 3 February 2020). The procedure involves chemical (oxidation) and physical (size and density separation) fractionation steps, based on the current understanding of prevailing SOC stabilisation mechanisms in soils. In a recent comprehensive method comparison, this method was among the most efficient for isolating SOC fractions with varying turnover rates (Poeplau et al., 2018). In brief, 20 g of sieved soil was suspended in 150 mL of deionised water and subjected to a light ultrasonic treatment of 21 J mL⁻¹ at 30 W to disperse the most unstable aggregates and associations. Subsequently, the soil was wet-sieved with a fixed amount of water over 63 µm to separate silt- and clay-sized particles from sand-sized particles. Several pretests with the most extreme warming treatments and the unwarmed reference revealed that 1400 mL of deionised water was sufficient for a complete separation of the coarse (>63 µm) and fine fraction (<63 µm) particles, as indicated by clear rinsing water. The fine fraction containing silt- and clay-sized particles (SC) was centrifuged for 15 min at 1000 g, and an aliquot of the supernatant was filtered over 0.45 µm to derive the dissolved organic carbon fraction (DOC). Fine and coarse fractions were oven-dried at 40 ^∘C and weighed. Sodium polytungstate (SPT) with a density of 1.8 g cm⁻³ was used to separate the coarse light fraction, i.e. particulate organic matter (POM), from the coarse heavy fraction, i.e. the sand and stable aggregates fraction (SA). To do so, about 40 mL of SPT was added to the coarse fraction in a centrifuge tube and stirred gently. The stirred samples were left standing for several hours at room temperature so that particles could float or sink and were subsequently centrifuged for 15 min at 1000 g for the complete separation of the light and heavy fractions. The supernatant was decanted into a sieve bag with a 50 µm mesh size. The density fractionation procedure was repeated once to ensure the complete separation of the light and heavy fractions. After the second SPT treatment, the remaining heavy fraction was transferred to a sieve bag with a 50 µm mesh size, and both the heavy and light fractions were washed thoroughly to remove all SPT, dried at 40 ^∘C and weighed. Based on this procedure, we use the term aggregates in the following for the 63–2000 µm aggregate size fraction, which comprises larger microaggregates as well as macroaggregates (Totsche et al., 2018). Finally, the SC fraction was subjected to sodium hypochlorite (NaOCl) oxidation, which is undertaken to mimic strong enzymatic decay and isolate the oxidation-resistant SOC fraction (rSOC). To do so, NaOCl with 6 % Cl was first adjusted to a pH of 8 using concentrated HCl. A 1 g aliquot of the SC fraction was then mixed with 40 mL NaOCl. After a 17 h reaction time, samples were centrifuged, decanted and washed once with deionised water. The whole procedure was repeated twice to ensure complete oxidation of NaOCl-oxidizable SOC (SC-rSOC). Thereafter, soil was dried at 40 ^∘C and weighed to determine the mass loss caused by oxidation. All solid fractions and the bulk soil were ball-milled and measured for the respective C and N content via dry combustion (LECO TruMac, St Joseph, MI, USA). The DOC fraction was measured using a liquid analyser (DIMATOC, Dimatec, Essen, Germany). The average mass recovery was 97±2 %, and the average C recovery was 99±21 %. In the following, two different measures of SOC in the isolated fractions will be used, depending on the context: (i) the SOC concentration, which indicates the amount of SOC in each fraction per fraction mass (g C kg fraction⁻¹), and (ii) the SOC content, which indicates the amount of SOC in each fraction per bulk soil mass (g C kg soil⁻¹).

To determine the total amount of soil in stable aggregates, i.e. to separate the SA fraction into sand and stable aggregates, another 4 g of each bulk soil sample was subsequently used. Instead of the soft ultrasonic treatment of 21 J mL⁻¹, we applied 500 J mL⁻¹ at a high amplitude (70 %) to completely disperse all aggregates (Schmidt et al., 1999). After subsequent wet sieving, the mass proportion of the coarse fraction (>63 µm) containing POM and pure sand grains was determined and subtracted from the earlier coarse fraction to determine the mass proportion of stable aggregates.

To evaluate the effect of bulk SOC and SOC fractions on soil structure, we determined the poured bulk density in the bulk soil as well as the coarse (SA + POM) and fine (SC) fractions of each sample. The poured bulk density is also known as the aerated bulk density and is a measure of the structural strength of loose material (Abdullah and Geldart, 1999). This was done by pouring material of known weight into a scaled cylindric flask to measure the volume of the sample. The poured bulk density of each individual sample (ρ_i g cm⁻³) was then calculated as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\rho }}_i={\displaystyle \frac{{\mathrm{mass}}_i}{{\mathrm{volume}}_i}},\end{array}$$

where mass_i is the total soil mass of the individual fraction (g) and volume_i is the volume of the individual fraction (cm⁻³). We assumed that a higher poured bulk density would indicate less structure and hypothesised that ρ_i would be negatively correlated with the SOC content in the SA fraction in particular.

Soil sampling of the adjacent grassland SOC (data from previous study) was carried out in December 2014, 6 years after the warming was initiated, and involved the same experimental design and analyses as that undertaken on the forest soil (Poeplau et al., 2017).

2.3 Statistics

The balanced design of the experiment, i.e. the six warming intensities, five transects (replicates) and two different sampling depths, allowed for the use of an analysis of variance (ANOVA) to test the differences between warming intensities in the bulk SOC and SOC fractions for significance. Moreover, a nonparametric analysis of similarity (ANOSIM), as implemented in the vegan R package (Oksanen et al., 2019), was used to test if warming significantly altered the SOC composition, i.e. its distribution in different fractions. Finally, an analysis of covariance was used to assess whether forest SOC (data from this study) and grassland SOC (data from a previous study) would differ in their response to soil warming. This was done using an ANOVA including ecosystem, warming intensity and their interaction. Linear or logarithmic regression models were used to describe the warming response of bulk SOC and SOC fractions. The Akaike information criterion (AIC) was used to select the most suitable model for each individual case. Despite the fact that some temperature responses were non-linear, we used linear regressions to derive absolute and relative changes in SOC concentration per degree Celsius as a proxy in order to compare the different fractions. Whenever necessary, data were log-transformed to an approximate normal distribution, which was visually assessed using histograms. Significance was assessed at a level of p<0.05. All statistical tests and plots were carried out in R (R Development Core Team, 2010). The ggplot2 package was used for plots (Wickham, 2016).

3 Results

3.1 Warming-induced changes in forest soil organic carbon

After 10 years of soil warming, the bulk SOC content in the forest soil had dropped severely in all of the investigated warming treatments. In the forest soil, warming-induced SOC losses increased linearly with the degree of warming (Fig. 1a, b, Table 1) in both depth increments. Absolute losses in the topsoil (−2.7 g kg^−1 ∘C⁻¹, Table 1) were more pronounced than absolute losses in the subsoil (−1.6 g kg^−1 ∘C⁻¹, Table 1). In the topsoil, SOC dropped from 75.1 g kg⁻¹ in the unwarmed soil to 26.5 g kg⁻¹ in the most warmed soil; in the subsoil it dropped from 36.2 to 4.0 g kg⁻¹. Thus, relative losses were even more pronounced in the subsoil (−4.5 % SOC ^∘C⁻¹) compared with the topsoil (−3.6 % SOC ^∘C⁻¹). Despite these strong linear trends, the SOC content in the bulk soil was only significantly different from the unwarmed reference at warming intensities of 5.8 and 17.5 ^∘C (topsoil) as well as 17.5 ^∘C (subsoil) (Table 1). The same was true for the SOC content in SA and POM; however, for SC and rSOC, only a warming intensity of 17.5 ^∘C was enough to significantly decrease the SOC content at both depths after 10 years. For DOC, significant changes with warming were only observed in the subsoil. In the topsoil, relative changes in the SOC content were in the following order: POM > SA > bulk soil > DOC > SC > rSOC. This order is in agreement with the concept of the fractionation method, i.e. a stronger decline in the most labile fractions and a slower decline in the more stable fractions. However, this was not the case for the subsoil, in which the order of relative SOC changes almost reversed to rSOC > SC > POM > bulk soil > SA > DOC (Table 1). Nevertheless, the strong changes in rSOC and SC were mainly driven by the 17.5 ^∘C warming intensity.

Figure 1 Areal plots of (a) the soil organic carbon (SOC) content in the topsoil and (b) the SOC content in the subsoil as well as (c) the SOC proportion in each fraction of the topsoil and (d) the SOC proportion in each fraction of the subsoil as a function of warming intensity. Fractions were dissolved organic carbon (DOC), particulate organic matter (POM), SOC in sand and aggregates (SA), non-oxidation-resistant silt- and clay-sized SOC (SC-rSOC), and oxidation-resistant silt- and clay-sized SOC (rSOC).

Table 1 The average soil organic carbon (SOC) content (n=5) of all fractions and the bulk soil including their standard errors. The letters following the numbers in the table indicate significant differences (p<0.05) across warming intensities (^∘C) within one soil depth. Absolute and relative changes in the SOC content as derived from linear regression models are also displayed for both of the investigated soil depths. Although this was not the best model in all cases, we used this value as a proxy to compare the warming response among fractions. Fractions were dissolved organic carbon (DOC), particulate organic matter (POM), SOC in sand and aggregates (SA), total silt- and clay-sized SOC (SC), and oxidation-resistant silt- and clay-sized SOC (rSOC).

The depletion of the bulk SOC content led to altered relative distributions of SOC in the isolated fractions (Fig. 1c, d). The ANOSIM revealed that warming intensities of 5.8 and 17.5 ^∘C were necessary to significantly change the topsoil SOC distribution (Table 2). In the subsoil, the fraction distribution was significantly different from the unwarmed reference at warming intensities of 2.7 and 5.8 ^∘C. In the topsoil, the unwarmed reference soil was strongly dominated by SOC in the POM and SA fractions (∼90 % combined), which were strongly depleted with warming (Fig. 1). This led to a relative increase in SOC stored in the fine fractions (SC-rSOC and rSOC). In the topsoil, even an absolute increase of SOC in these fractions was observed upon warming (Fig. 1a), which strongly indicated a redistribution of fraction masses. Indeed, the soil mass of the SA fraction decreased with warming, while the mass of the SC fraction increased (Fig. 2). This was true for both of the investigated soil depths, with the mass distribution of the subsoil at the 17.5 ^∘C warming intensity being an exception. As expected, the second ultrasonic step revealed that only the aggregates depleted within the SA fraction, while the proportion of sand-sized mineral particles remained stable across warming levels (Fig. 2). Therefore, the aggregate mass proportion in the topsoil decreased from 60.7±2.2 % in the unwarmed reference to 28.9±4.6 % in the 17.5 ^∘C warmed soil. In the subsoil, it decreased from 43.7±3.8 % in the unwarmed reference to 17.7±2.9 % in the 5.8 ^∘C warmed soil, while at a warming intensity of 17.5 ^∘C, the mass proportion of aggregates amounted to 32.9±4.9 %. The average sand content of 28 % determined after the second ultrasonic treatment (Fig. 2) concurred well with the 31 % sand content of the texture analysis.

Figure 2 Areal plots of the soil mass distribution in the particulate organic matter (POM), sand and stable aggregates (SA), and silt and clay (SC) fractions as a function of warming intensity.

Table 2 Summary of the analysis of similarity (ANOSIM) that tested the differences in the distribution of SOC in the investigated fractions for all warming intensities against the unwarmed reference. p values < 0.05 indicate significant differences, whereas “n.s.” indicates nonsignificant differences. An R value close to 1 suggests dissimilarity between groups.

Within the fine fraction, the relative mass proportion of rSOC was expected to increase with warming due to its proposed higher biogeochemical stability compared with the NaOCl-oxidised part of the SC fraction. However, this was not the case: across all warming intensities and both soil depths, we found a significant linear correlation between rSOC and total SOC in the SC fraction (y=0.319x, R²=0.92, p<0.001). Thus, the NaOCl treatment constantly oxidised two-thirds of the SC fraction across all warming intensities, indicating that no relative accumulation of rSOC occurred within the silt- and clay-sized soil fraction.

Interestingly, the proportion of SOC that was water soluble (DOC) tended to increase with warming in both of the investigated depth increments (Fig. 1c, d), which was not significant. However, for the topsoil, we detected a significantly negative relationship of the percentage of total SOC in SA and the percentage of total SOC in DOC (Fig. 3), which might point towards the SOC-stabilising function of aggregates.

Figure 3 Correlation between the proportion of soil organic carbon (SOC) in the sand and aggregates (SA) and dissolved organic carbon (DOC) fractions in the topsoil with the 95 % confidence interval.

3.2 Forest vs. grassland soil carbon responses to warming

The observed changes in the bulk and fraction SOC in the forest soil were generally comparable with those in the adjacent grassland soils (Figs. 4, 5). Especially in the subsoil, the interaction effect of ecosystem and warming on SOC was not significant for four out of five fractions and the bulk soil, indicating the same SOC response to warming in both ecosystem types (Fig. 5, Table 3). Moreover, the difference between ecosystems with respect to the subsoil SOC content was less pronounced than in the topsoil. This might be partly related to the fact that the forest was planted on an unmanaged grassland and that the forest subsoil SOC was still grassland-derived to a high extent. However, for the topsoil, we found significant interactive effects of ecosystem and warming for four out of five fractions and the bulk soil (Table 3). The forest soil, which had a considerably higher bulk SOC content in the unwarmed reference than in the grassland, showed a stronger response to warming. The predominant SOC fraction in the forest topsoil was the SA fraction, which responded most strongly to warming (Fig. 1). This was generally observed in both ecosystems. However, the stronger redistribution of soil mass across fractions in the forest soil compared with the grassland soil led to very distinct responses from SC-rSOC and rSOC, with stronger warming-induced increases of these fractions in the forest soil (Fig. 4). Furthermore, the POM fraction of the forest soil responded more negatively to warming than that in the grassland soil. With respect to the warming response of DOC alone, we did not detect any differences between the ecosystems in the topsoil. Interestingly, despite the differences in the initial SOC and warming duration, i.e. 10 years for the forest and 6 years for the grassland, the SOC in both ecosystems approached an almost equal SOC content in the most extreme warming intensities (Fig. 4).

Figure 4 Soil organic carbon (SOC) mass in the bulk soil and fractions of the forest and grassland topsoils (0–10 cm) as a function of warming intensity with linear and logarithmic fits and the 95 % confidence intervals. Fractions were dissolved organic carbon (DOC), particulate organic matter (POM), SOC in sand and aggregates (SA), non-oxidation-resistant silt- and clay-sized SOC (SC-rSOC), and oxidation-resistant silt- and clay-sized SOC (rSOC).

Figure 5 Scatter plots showing the soil organic carbon (SOC) content in the bulk soil and fractions of the forest and grassland subsoils (20–30 cm) as a function of warming intensity with linear and logarithmic fits and the 95 % confidence intervals. Fractions were dissolved organic carbon (DOC), particulate organic matter (POM), SOC in sand and aggregates (SA), non-oxidation-resistant silt- and clay-sized SOC (SC-rSOC), and oxidation-resistant silt- and clay-sized SOC (rSOC).

Table 3 Summary of the linear regression models (p values) assessing the effects of warming, ecosystem (grassland vs. forest) and their interaction on soil organic carbon (SOC) for the bulk soil and all isolated fractions.

3.3 Structural changes following soil carbon loss

As expected, we found a strong negative correlation between the SOC content and the poured bulk density (Fig. 6a, R²=0.92, p<0.001). A very similar relationship with an identical slope was observed for the coarse (>63 µm) soil fraction, comprising SA and POM (Fig. 6b). In contrast, we did not detect any correlation between the SOC content and the poured bulk density in the silt and clay fraction (data not shown). A direct link between the poured bulk density and aggregates is given in Fig. 6c. Finally, in agreement with the strong decline in SOC and soil mass in the SA fraction with warming intensity (Figs. 1, 2), we found a strong positive correlation between the SOC mass and soil mass in the coarse soil fraction comprising SA and POM (Fig. 6d). The slope of the regression was 4.5, indicating that one unit of SOC was causing the aggregation of 4.5 units of soil. The effects of SOC on soil structure were equally observed in the topsoil and subsoil. Furthermore, for all of the structure-related parameters shown in Fig. 6, observations of both of the investigated soil depths scattered approximately around the same regression line. This might indicate that SOC depletion, rather than soil warming, induced the breakdown of aggregates.

Figure 6 Poured bulk density as a function of the soil organic carbon (SOC) content in (a) the bulk soil and (b) the coarse (>63 µm) soil fraction (sand and stable aggregates, SA, and particulate organic matter, POM); (c) the poured bulk density as a function of soil mass in aggregates and (d) soil mass in the coarse soil fraction as a function of SOC mass in the coarse soil fraction with regression models fitted to all observations (p<0.001 for all models).

4 Discussion

4.1 Warming effects on forest soil organic carbon and its fractions

A total of 10 years of forest soil warming caused a strong decline in the SOC content. Along the temperature gradient, SOC changes followed a linear response, with a −3.6 % ^∘C⁻¹ and −4.5 % ^∘C⁻¹ change in topsoil and subsoil respectively. Thus, under the most extreme warming intensity treatment of 17.5 ^∘C, SOC was depleted by 65 % and 89 % in the topsoil and subsoil respectively. Considering that an air temperature increase of up to 11 ^∘C by the end of the century is within the possible range of IPCC climate change projections (IPCC, 2013), we assume that a soil warming intensity of up to 5.8 ^∘C can be considered realistic. For example, Zhang et al. (2005) showed that soil temperature increase (+0.6 ^∘C) generally followed the air temperature increase (+1 ^∘C) in Canada during the 20th century. At a warming intensity of 5.8 ^∘C, the investigated soil lost 29 % (topsoil) and 37 % (subsoil) of the SOC in 10 years. This is in line with other studies, which also reported significant losses of SOC upon warming (Crowther et al., 2016, and papers cited therein). In the investigated experiment, there is no doubt that potential warming-induced changes in net primary productivity (NPP; Sigurdsson et al., 2014) did not offset increased soil microbial activity. In fact, root biomass between 0 and 10 cm decreased in both ecosystems (data not shown), leading to weak positive correlations (R²=0.37 for forest and R²=0.29 for grasslands) between SOC and root biomass. Moreover, aboveground plant litter tended to decline in both ecosystems. This suggests that SOC losses were partly driven by decreasing C input with warming and not by increased microbial activity alone. However, a clear picture of absolute C inputs in the experimental plots is not available yet, as NPP and biomass turnover need to be simultaneously taken into account.

Similar or relatively even more pronounced losses of SOC from the subsoil compared with the topsoil are confirmed by the results of a recent whole profile forest soil warming study, which concluded that subsoils will be an important source of CO₂ under climate change (Hicks Pries et al., 2017). Higher relative losses of SOC in the subsoil could potentially be driven by warming-induced changes in C input patterns. Indeed, especially fine root production and turnover of trees in the boreal zone has previously been found to increase with moderate warming (Leppälammi-Kujansuu et al., 2014; Majdi and Öhrvik, 2004), and fine roots are primarily located in the uppermost centimetre of forest soils (Hansson et al., 2013; Leppälammi-Kujansuu et al., 2013). However, at the investigated site, the amount of fine roots and mycorrhizal production has been found to decrease at the more extreme warming levels (Parts et al., 2019; Rosenstock et al., 2019). In addition, in this geothermal warming experiment, heat was coming from below, leading to slightly more intense soil warming in the subsoil. This is likely to explain the stronger relative SOC depletion in the subsoil to a certain extent. However, except for the highest warming level, the vertical gradients within the top 30 cm of soil were not substantial (Sigurdsson et al., 2016).

A major strength of a warming gradient approach is the identification of potential tipping points, which may mark abrupt changes in ecosystem functionality (Kreyling et al., 2014). However, the present study did not reveal such tipping points for bulk SOC content, which changed in a surprisingly linear manner with increasing temperature in both of the investigated depth increments. Despite certain methodological drawbacks of the geothermal (or any other manipulated) soil warming experiment, such as very abrupt initial temperature changes and soil warming from below instead of whole ecosystem warming from above, it can be inferred that climate change is likely to strongly affect SOC stocks of subarctic forests. These forests cover an area of approximately 15×10⁶ km² or one-third of the global forest area (Bonan, 2008). The analysis of the soil warming gradient also revealed detection limits for warming effects on SOC that are intrinsically very heterogeneous in space and respond slowly to environmental change (Smith, 2004): even after 10 years of chronic soil warming, changes in topsoil SOC were only significant at a warming intensity of at least 5.8 ^∘C, when assessed using the ANOVA approach. An ANOVA, instead of a regression analysis, is required when only one warming treatment is investigated (e.g. Schnecker et al., 2016). If this treatment is relatively mild, e.g. below 4 ^∘C, changes might easily be undetectable against the background heterogeneity of SOC. This is an important insight considering the ongoing debate regarding whether SOC is lost upon warming or not (Crowther et al., 2016; van Gestel et al., 2018). The majority of currently available datasets are based on such experiments with relatively short, mild and singular warming treatments (van Gestel et al., 2018). However, the transferability of the results in this study to the SOC response to global warming is still rather limited and can only slightly reduce the following given uncertainties: (i) we studied the soil temperature, not the air temperature increase; (ii) the warming occurred abruptly and not gradually; and (iii) we studied an Andosol. Thus, extrapolations to larger areas or longer time periods should be undertaken carefully and were not intended in this study.

The fractionation method used in this study isolates SOC pools of different biogeochemical stabilities (Zimmermann et al., 2007). Turnover rates are estimated to range from several years in the POM fraction to centuries in the oxidation-resistant rSOC fraction that is associated with silt and clay particles (von Lützow et al., 2007). Such differences are mainly related to different degrees of physicochemical stabilisation in the soil, such as the interaction with the mineral phase or occlusion into aggregates (von Lützow et al., 2007). Due to differences in composition and the bioavailability of these SOC fractions, distinct responses to warming were expected in the following order: POM > DOC > SA > bulk > soil > SC-rSOC > rSOC. Indeed, the average relative decrease in the SOC content, which might be the best indicator to describe a fraction's sensitivity to warming, was observed to follow a similar order in the topsoil: POM > SA > bulk soil > DOC > SC-rSOC > rSOC. This concurs well with the sensitivity of these fractions to land use change as observed across different land use changes by Poeplau and Don (2013). However, the difference in the warming response between SC-rSOC (−2.14 % ^∘C⁻¹) and rSOC (−2.05 % ^∘C⁻¹) was negligible, which was also reflected in the stable proportion of rSOC in the total SC fraction throughout the warming gradient. This indicated that NaOCl oxidation did not yield a meaningful fraction with regard to biogeochemical resistance. This has been observed before and calls the notion that this oxidation-resistant pool can be linked to a centennially persistent or even inert SOC pool into question (Lutfalla et al., 2014; Poeplau et al., 2019, 2017; Zimmermann et al., 2007). At the same time, NaOCl-resistant SOC has often been described as substantially older and, thus, slower cycling as bulk SOC (Helfrich et al., 2007) and has also been found to correlate with the abundance of Al and Fe oxides in the soil (Mikutta et al., 2005). Therefore, the strong warming response of this fraction is somewhat in contrast with the slow responses observed in other treatments, such as C₃–C₄ vegetation changes (Poeplau et al., 2018). In the subsoil, the average relative depletion in rSOC was even strongest across all fractions and the bulk soil. However, this was related to the very low carbon content of the highest warming intensity (17.5 ^∘C) driving the slope of the regression. Only when the highest warming intensity was excluded, the sensitivity of fractions followed the observed order in the topsoil, with DOC being an exception: POM > SA > bulk soil > SC-rSOC > rSOC > DOC.

4.2 Is aggregate breakdown induced by soil organic carbon losses or vice versa?

The most significant warming effect on the distribution of SOC in the isolated fractions was the strong decrease in SA. In the unwarmed reference soil, it accounted for the highest proportion of soil mass and SOC content. However, with warming, aggregates collapsed, leading to strong mass increases in the fine SC fractions, which even increased the carbon mass upon warming. The second ultrasonic step, which was used to distinguish sand from aggregates in the SA fraction, provided evidence that the aggregate size fraction investigated (63–2000 µm) was strongly reduced. A tipping point for aggregate breakdown appears to be located between the warming intensities of 2.7 and 5.8 ^∘C. The same mechanism, although less pronounced, was observed for the adjacent grassland (Poeplau et al., 2017). Observing SOC depletion and aggregate breakdown at the same time raises the question of cause and effect: aggregates – at least microaggregates smaller than 250 µm – are acknowledged to protect organic matter from microbial decomposition (Six et al., 2002). At the same time, organic matter, especially mucilage, polysaccharides and fungal hyphae, acts as an aggregate binding agent (Tisdall and Oades, 1982). Answering the question of whether warming has intrinsically fostered aggregate breakdown via changes in biotic and abiotic environmental conditions might be of critical importance for conceptualising and modelling warming effects on SOC dynamics. However, the results of the present study suggest that the major cause of aggregate breakdown was not necessarily warming and could instead be well described by loss of SOC: we found a very strong positive correlation of SOC mass and total soil mass in the coarse soil fraction (comprising POM and SA) – 1 g kg⁻¹ of SOC was keeping 4.5 g kg⁻¹ of soil aggregated. Topsoil and subsoil samples scattered approximately around the same regression line. This indicates that the abundance of young and coarse SOC, rather than the degree of soil warming, is driving the amount of stable aggregates in the soil. This is well known and, thus, in accordance with the literature (Franzluebbers, 2002; Oades, 1984; Shepherd et al., 2002). Another reason to doubt that warming-induced aggregate breakdown caused the destabilisation of SOC is the fact that the SOC protection capacity of macroaggregates is debatable (Six et al., 2004). For example, Bischoff et al. (2017) found higher heterotrophic respiration in uncrushed soil compared with the same soil with crushed macroaggregates. To some extent, a positive feedback loop, i.e. SOC depletion causing aggregate breakdown which in turn causes the mineralisation of accessible C, might indeed be possible. The fact that the proportion of water soluble SOC in the topsoil increased with decreasing aggregation points in this direction. The desorption of carbon compounds from the mineral phase is likely to be fostered by increased surface area, which is the case when aggregates disintegrate. However, soil pH is also acknowledged to affect DOC formation (Kalbitz et al., 2000), which might be another possible explanation for the observed increase in the proportion of DOC: in both ecosystems, the soil pH increased by up to 0.5 units under the highest warming intensity treatment (Sigurdsson et al., 2016).

4.3 Linking losses in soil organic carbon to changes in soil structure

As a consequence of SOC loss, the total pore space decreased strongly as indicated by the poured bulk density. The poured bulk density was used as a proxy for the in situ bulk density in the undisturbed soil, which was unfortunately not determined in the present study. However, the relationship between SOC and the poured bulk density was in the range of established pedotransfer functions (PTFs) for field bulk density estimation using SOC content. In a literature review comparing different PTFs (De Vos et al., 2005), the slopes of regression models using SOC content (g C kg⁻¹) to predict soil bulk density (g cm⁻³) ranged from −0.003 to −0.011, while the slope in the present study was −0.005 for both the bulk soil and the SA fraction. This negative correlation is due to the much lower specific gravity of organic matter compared with mineral particles and also to the effect of organic matter on aggregation (De Vos et al., 2005). The variation in slopes, i.e. the effect of SOC on bulk density, is most likely related to the soil's capability to form aggregates. In very sandy soils with a single grain structure, even high organic matter contents do not lead to the considerable formation of aggregates, and the organic matter effect on the bulk density is mainly restricted to a gravity effect. Using a two-pool mixing model of mineral particles with a density of 2.5 g cm⁻³ and soil organic matter with a density of 1, i.e. ignoring the structural effect of organic matter, we found a slope of −0.0026. Accordingly, Callesen et al. (2003) reported a PTF for sandy forest soils with a slope of approximately −0.0028 in the range of 0–80 g SOC kg⁻¹ (non-linear function). Thus, the slope of −0.005 found in this study might indicate that approximately 50 % of the SOC effect on the poured bulk density can be assigned to a structural effect. Indeed, we also found a strong negative correlation between the soil mass stored in aggregates and the poured bulk density. To conclude, the slope of the regression between SOC and bulk density, at least in unmanaged soils, might be a good indicator of the aggregation affinity of a soil. However, the poured bulk density of disturbed and sieved soil can only express a potential and should be treated as such. Conversely, factors like position in the soil profile that strongly influence the packing density of the soil are cancelled out, enabling a direct comparison of topsoil and subsoil samples.

Strong systematic gradients in the SOC content in the same soil, such as those created by the soil warming in our study, are rare and extremely valuable for improving our understanding of organic matter functions. Larsbo et al. (2016) used a natural SOC gradient to evaluate its effect on pore networks, influencing solute and gaseous transport in the soil. Changes in soil structure induced by large SOC loss might also affect other key ecosystem properties, such as NPP (Oldfield et al., 2019), microbial biomass (Walker et al., 2018) or other soil biota. For example, in the adjacent warmed grassland plots, Holmstrup et al. (2018) detected a warming-induced shift in Collembola species abundance towards species with smaller body size. An increase in the bulk density with the associated decrease in pore space might have fostered this physiological response, although this was not explicitly mentioned by the authors. Moreover, a positive correlation between pore volume and microbial and nematode biomass was found by Hassink et al. (1993). In the present study, aggregation and poured bulk density were assessed on sieved soils, which provided valuable initial information on warming-induced changes in basic soil structural parameters. A follow-up study should investigate soil structure and other physical parameters in undisturbed soil samples for two major reasons: (i) the gradient in the SOC content is unique and can be used to improve the general understanding of the link between organic matter and soil functions; (ii) the warming responses of many ecosystem aspects are studied along the investigated warming gradients and knowledge on changes in soil physical properties might be central to interpret such responses. Moreover, these structural changes most likely led to a certain sampling bias and, thus, a slight overestimation of SOC losses: sampling fixed depth increments ignores the fact that depth increments change with changes in bulk density. Therefore, the depth increments sampled under the higher warming intensity treatments do not exactly match the depth increments sampled under the lower warming intensity treatments. However, this effect is expected to be more pronounced in the topsoil, where the SOC depth gradient is largest and a shift in reference soil depth would, therefore, have the strongest impact on the bulk SOC content. However, relative losses in SOC were even more pronounced in the subsoil, indicating that the sampling bias might have been small. However, it should be mentioned that mass-based instead of depth-based sampling (Don et al., 2019) or at least an a posteriori soil mass correction (Ellert and Bettany, 1995) would be indispensable to accurately estimate SOC stock changes.

4.4 Comparing forest and grassland soil carbon responses to warming

To date, warming experiments have mostly focused on one single type of ecosystem. However, the warming response could be ecosystem specific (Shaver et al., 2000), which can only be investigated using a paired ecosystem approach. In the present study, we investigated a small stretch of forest located directly adjacent to a similarly warmed grassland. Changes in the SOC content and the relative distribution of fraction masses in the grassland soils have been previously investigated (Poeplau et al., 2017). Both ecosystems showed a similarly strong response to warming. The fact that no difference in subsoil SOC dynamics in the bulk soil or any isolated fractions were observed might indicate that the same mechanisms of SOC depletion were involved in both ecosystems. For example, aggregate breakdown and an equal decrease in rSOC and SC-rSOC were also observed in the grassland. However, the initial SOC content and fraction distribution in the topsoil differed across ecosystems, leading to distinct responses to warming: the unwarmed forest had about 50 % more SOC in the topsoil compared with the grassland, and about 150 % more SOC was stored in the SA fraction. Furthermore, the POM fraction in the forest was almost twice that in the grassland, with proportionally less SOC stored in more stable fractions. The shift in the fraction mass distribution, i.e. aggregate breakdown, was more pronounced in the forest topsoil, leading to the increase in the fine-fraction SOC with warming, which was not observed in the grassland. Crowther et al. (2016) reported that SOC loss upon warming is a function of initial SOC, and the present study confirms this. In fact, to some extent, the explanation for this might be the higher proportion of labile SOC in soils with higher SOC stocks (Besnard et al., 1996). It has been previously reported that forest SOC is more labile than grassland SOC (Poeplau and Don, 2013). The forest was sampled after 10 years of warming, and the grassland was sampled after 6 years. However, (i) subsoils showed an almost identical response to warming and (ii) there were indications that the grassland had already reached a new SOC steady state after 6 years of warming (Walker et al., 2018). Therefore, it seems likely that amount and fraction distribution of SOC drove the ecosystem-specific warming response in the topsoil. This difference in the topsoil SOC and fraction distribution has been found before and is related to the different sources and qualities of fresh organic matter inputs (Poeplau and Don, 2013; Huang et al., 2011). In particular, needle litter is acknowledged to decompose slowly (Prescott et al., 2000). Differences in POM as well as total SOC stocks are observed to level off with increasing soil depth (Davis and Condron, 2002; Poeplau and Don, 2013). This might also be true for the response to warming, as indicated by the present study. Finally, the SOC content in both ecosystems approaches a similar baseline under the highest warming intensity. This might indicate that the specific amount of biogeochemical persistent SOC does not depend on land cover or vegetation type but is rather controlled by mineralogy.

5 Conclusion

Using a strong geothermal warming gradient, we found a clear link between SOC losses and soil structural changes. A total of 10 years of soil warming created a steep gradient in the SOC content that is rare and should be used to study the links between organic matter and soil structure and soil function more deeply. The results of the present study reveal that the effects of warming on biogeochemical cycles are most likely not restricted to direct effects on biotic processes and that changes in the microbial habitat and possibly abiotic soil properties should be considered. These factors are likely to exert a strong indirect influence on any biotic response. Differences in the warming response of the bulk SOC and SOC fractions between ecosystems have only been found in the topsoil, which might be related to the fact that the forest was planted on unmanaged grassland half a century ago. In the forest, the depletion of SOC was more pronounced in the subsoil, which calls for more whole soil profile warming studies.