The limited effect of deforestation on stabilized subsoil organic carbon in a subtropical catchment

Predicting the quantity of soil organic carbon (SOC) requires understanding how different factors control the amount of SOC. Land use has a major influence on the function of the soil as a carbon sink, as shown by substantial organic carbon (OC) losses from the soil upon deforestation. However, predicting the degree to which land use change affects the OC content in soils and the depth down to which this occurs requires context-specific information related to, for example, climate, geochemistry, and land use history. In this study, 266 samples from forests and agricultural fields were collected from 94 soil profiles down to 300 cm depth in a subtropical region (Arvorezinha, southern Brazil) to study the impact of land use on the amount of stabilized OC along the soil profile. We found that the stabilized OC content was not affected by land use below a depth of 90 cm. Along the soil profile, the amount of stabilized OC was predominantly controlled by land use and depth in addition to the silt and clay content and aluminium ion concentrations. Below 100 cm, none of the soil profiles reached a concentration of stabilized OC above 50 % of the stabilized OC saturation point (i.e. the maximum OC concentration that can physically be stabilized in these soils). Based on these results, we argue that it is unlikely that deeper soil layers can serve as an OC sink over a timescale relevant to global climate change due to the limited OC input in these deeper layers. Furthermore, we found that the soil weathering degree was not a relevant control on the amount of stabilized OC in our profiles because of the high weathering degree of the studied soils. It is therefore suggested that, while the soil weathering degree might be an effective controlling factor of OC stabilization over a large spatial scale, it is not an informative measure for this process at regional and local scales (with similar climate, bedrock, and weathering history) in highly weathered soils.

1 Introduction

Soils contain an estimated amount of 1.5 Gt of organic carbon (OC) in the upper metre (excluding permafrost) (Balesdent et al., 2018; Jobbágy and Jackson, 2000; Scharlemann et al., 2014). It has been estimated that ca. 44 % of global OC present in soils is stored in the tropics (Veldkamp et al., 2020), which are undergoing the greatest rate of land use (LU) change globally, mainly as a consequence of the rapid conversion of natural land uses to agricultural land (Hansen et al., 2013). This rapid land use change has considerable impacts on the amount of soil organic carbon (SOC) stored in tropical ecosystems. Indeed, the conversion of tropical forests to agricultural fields leads to an average decrease in the OC concentration in the topsoil (< 30 cm depth) of up to 50 % after 25 years (Veldkamp et al., 2020).

Nevertheless, estimates of changes in OC vary a lot between studies (de Blécourt et al., 2019; Bonini et al., 2018; Comte et al., 2012; Detwiler, 1986; Don et al., 2011; Powers et al., 2011) because this depends on different factors, such as the time since land use conversion and agricultural management practices (e.g. residue inputs reduce OC loss after conversion) (Detwiler, 1986; Veldkamp et al., 2020). While it is commonly accepted that land use has an impact on the amount of SOC, it is unclear down to which depth SOC storage is affected by land use change for specific climatic regions. Therefore, estimates of the effect of land use change on the OC stock below a depth of 30 cm are uncertain. Although some studies only observed a change in the topsoil (e.g. Guillaume et al., 2015), Strey et al. (2016) argue that soils should be studied at least down to 100 cm depth to fully assess the impact of land use change on the OC stock. For example, Veldkamp et al. (2020) found that forest soils below 0.5 m depth contained, on average, 35 % more OC compared to crop plantations that had been converted more than 25 years previously. This illustrates that the subsoil needs to be accounted for when assessing the impact of land use change on the OC content.

The SOC pool consists of a wide variety of organic molecules, ranging from unprocessed plant material to organic molecules with a wide range of oxidation states (Kögel-Knabner, 2002; Lehmann and Kleber, 2015). The residence time of OC in soils is not uniform, with respired CO₂ generally being substantially younger than bulk SOC (Trumbore et al., 1995), which can have an age on the order of thousands of years (Mathieu et al., 2015). As unprocessed organic residues decompose, they break down into smaller OC-containing particles that can form chemical bonds with soil minerals, a process that at least partially protects OC from microbial degradation and thus contributes to long-term carbon storage (Lehmann, J. & Kleber, 2015; Cotrufo et al., 2013; Kögel-Knaber et al., 2002; Schrumpf et al., 2013; von Lützow et al., 2008). It is estimated that this mineral-associated OC (MAOC) constitutes more than 65 % of total SOC globally (Georgiou et al., 2022).

To quantify the proportion of SOC pools with different properties and a different turnover rate, soil organic matter (SOM) is generally separated into particulate organic matter (POM), which consists of relatively unprocessed plant material, and mineral-associated organic matter (MAOM) (Cotrufo et al., 2019; Lavallee et al., 2020; Lehmann and Kleber, 2015; Kögel-Knabner et al., 2008) using a variety of soil fractionation methods (Poeplau et al., 2018). POM is characterized by a low density, a high C:N ratio, and a relatively young age, while MAOM has a high density, a low C:N ratio, and a relatively old age (Kögel-Knabner et al., 2008; Lavallee et al., 2020; von Lützow et al., 2008). Information on how SOC is distributed among these pools is useful, for example, to assess how agricultural management practices can store atmospheric CO₂ in SOC pools with a long residence time (Chenu et al., 2019; Kallenbach et al., 2015; Kögel-Knabner et al., 2022; Tiefenbacher et al., 2021; Zhang et al., 2021) or to inform SOC models that simulate the proportion of SOC in pools with different turnover rates (e.g. Abramoff et al., 2018; Ahrens et al., 2015; Tang and Riley, 2015).

The dynamics of MAOC have mainly been studied in the topsoil (here defined as the upper 30 cm of the soil) in temperate regions (Cotrufo et al., 2019; Georgiou et al., 2022; Lugato et al., 2021; Kleber et al., 2015; Kramer and Chadwick, 2018; Rocci et al., 2021; Sokol et al., 2022). At the scale of the European continent, for example, it has been shown that between ca. 25 % and 100 % of topsoil SOC consists of MAOC, depending on land use and the type of mycorrhizal association of the dominant vegetation (Cotrufo et al., 2019; Lugato et al., 2021). However, the subsoil (≥ 0.30 m depth) contains a large portion of the total SOC stock (ca. 55 % of total SOC stock globally; Lal, 2018). This OC has average residence times ranging from decades to millennia (Mathieu et al., 2015), making it an important long-term carbon sink in terrestrial ecosystems (Rumpel and Kögel-Knabner, 2011). Subsoil OC consists of more processed and microbially derived molecules with a lower $\mathrm{C}/\mathrm{N}$ ratio compared to the topsoil (Schrumpf et al., 2013). These OC molecules are small and abundant in polar groups (Kleber et al., 2021). Thereby, they are highly reactive towards the mineral matrix, causing the share of MAOC to total organic carbon (TOC) to generally increase with increasing soil depth (Cotrufo et al., 2013; Schrumpf et al., 2013).

Improving global projections of SOC stabilization requires predictive models that allow us to assess the amount of MAOC that currently is and potentially can be stored in soils under specific climate and land use regimes. Building such models requires that the factors controlling MAOC be understood. For example, soils with high silt (soil particles in 2–53 µm) and clay (soil particles < 2 µm) contents generally store more SOC because of their higher specific surface area as compared to coarse-textured soils (Amato and Lass, 1992; Feller and Beare, 1997; Hassink, 1994; Kleber et al., 2015). However, there is an upper limit in the capacity of soil minerals to stabilize organic molecules (Hassink, 1997; Six et al., 2002c), generally referred to as the maximum sorption capacity (Guggenberger and Kaiser, 2003; Kothawala et al., 2009). The latter has been shown to be controlled by edaphic factors such as mineralogy and grain size (Abramoff et al., 2021; Feng et al., 2013; Rasmussen et al., 2018). Across the soil profile, the ratio of MAOC to maximum sorption capacity (i.e. the MAOC saturation) generally decreases with soil depth, with subsoils having the largest MAOC deficit and, thus, at least theoretically, the largest potential to stabilize additional OC (Chen et al., 2018; Georgiou et al., 2022). Therefore, it has been proposed that increasing the amount of subsoil MAOC can lead to a net transfer of atmospheric CO₂ to the SOC pool for timescales of decades to millennia (Rumpel et al., 2020), e.g. through growing deep-rooting crops (Kell, 2012), deep tillage (Alcántara et al., 2016), or deep soil flipping (Schiedung et al., 2019). While the potential for additional OC storage in the subsoil has been relatively well studied for different soil types and ecosystems in temperate ecosystems, data for tropical ecosystems are scarce (Georgiou et al., 2022).

In addition to soil texture, mineralogical and morphological characteristics of the soil have been identified as controls on SOC content (Rasmussen et al., 2018), while polyvalent cations such as Al³⁺ enable negatively charged OC molecules to bind to the mineral surface (i.e. cation bridging) (von Lützow et al., 2006). Furthermore, mineral soil weathering, which is the degree of degradation of bedrock and minerals due to physical, chemical, and biological factors, can have either a positive or a negative control on SOC stabilization. On the one hand, it promotes the formation of OC–mineral associations by increasing the specific surface area (SSA) of soil particles or by releasing Al and Fe oxides and cations (Depetris et al., 2014; Kleber et al., 2015). On the other hand, a high weathering degree might alter the structure of secondary minerals, reducing their reactivity and thus the potential for OC stabilization (Doetterl et al., 2018). However, all these parameters might influence each other, and their interaction depends on the climate and geochemistry (Doetterl et al., 2015a). For example, Six et al. (2002b) observed a lower association between silt and clay particles and OC in tropical soils in comparison to temperate soils. Therefore, it is necessary to study specific environments, such as highly weathered tropical soils, in order to disentangle the complex relationship between SOC storage and abiotic factors. Without a proper understanding of the controls of SOC stabilization in tropical soils, the extrapolation of the findings from a limited number of field observations may lead to erroneous results (Powers et al., 2011).

Here, we studied the effect of forest conversion to agriculture on stabilized OC along the soil profile in a subtropical environment. We analysed how deforestation in a catchment in southern Brazil has affected OC associated with the silt and clay fraction down to 3 m depth. The contents of OC associated with silt and clay within the bulk soil are referred to as the silt and clay organic carbon (OC_S+C) and represent OC stabilized on soil minerals (i.e. MAOC). The main factors controlling the amount of OC_S+C along the soil profile were assessed using a general additive mixed-effect model (GAMM). The silt and clay contents were used as limiting parameters to estimate the stabilized OC saturation point (i.e. the maximum OC concentration that can physically be stabilized in these soils). We hypothesized that OC_S+C concentrations differ between the two land use types in the topsoil, but this difference is no longer present below 100 cm depth. We furthermore hypothesized that soil texture and weathering degree are important controlling factors of the OC_S+C concentration and that OC_S+C in the subsoils of both environments is well below the maximum sorption capacity of the soil.

Figure 1 Map showing the sampling locations. Location of the study area in South America (left) and the municipalities of Arvorezinha, Itapuca, and Iliópolis (right), with an indication of the sampling locations (46 under forest and 48 under arable land). The elevation map was obtained from the R package elevatr (Hollister et al., 2022).

2 Material and methods

2.1 Study design and soil sampling

The study area is located in southern Brazil at the southern edge of the Paraná Basin in the Serra Gerral region (state of Rio Grande do Sul) (Fig. 1). The climate is classified as warm and humid with no dry season (Köppen classification Cfb) (Alvares et al., 2013; Peel et al., 2007).), with precipitation being evenly distributed throughout the year with an annual mean of 1900 mm, while the mean monthly temperature varies from 12 °C in July to 22 °C in January (Robinet et al., 2018a, b). The elevation is between 370 and 790 m above sea level, and the slopes vary between 2.5 and 30° (Brosens et al., 2021). The dominant soils in the study area are Acrisols, Leptosols, and Cambisols, overlying a basalt and dacite–rhyolite bedrock (Turner, Simon et al., 1994; Caner et al., 2014; Minella et al., 2009). The natural vegetation of the study area consists of a mixed ombrophilous forest (Araucaria forest) that has been subject to intensive land use change from forest to agriculture since the beginning of the 20th century (Morellato and Haddad, 2000; Lopes, 2006). The cultivated crops on the agricultural fields where samples were collected were tobacco (Nicotina tabacum), erva-maté (Ilex paraguariensis), maize (Zea mays), and soybean (Glycine max).

In order to assess how deforestation affected the distribution of OC_S+C down the soil profile, soil samples collected by Brosens et al. (2020, 2021) were used. The soil samples were originally collected to investigate the influence of slope gradient on soil thickness and chemical weathering. The study area was selected for its relatively uniform lithology, ensuring that variations in soil depth and weathering would primarily reflect differences in topography. A total of 226 soil samples were collected from mid-slope positions in 46 forested locations and 48 agricultural fields (94 soil profiles in total) using an Edelman auger. Sample collection is briefly described here; for more details, we refer the reader to the original studies. Sampling locations were selected using a stratified random approach to cover a wide range of slopes in a circular region of 250 km² (Fig. 1) using an approach that aimed to sample the widest variation in hillslope gradients in the study area. Only agricultural soils that were converted from forest at least 30 years prior were selected (Table S3 in the Supplement). Sampling was restricted to mid-slope positions. At each sampling location, one soil sample was taken from each soil horizon down to the bedrock or saprolite, with a maximum sampling depth of 2.95 m. The number of samples collected per site depended on the number of soil horizons, which ranged between one and six. All samples were collected from the centre of the extracted soil to avoid contamination from more shallow soil layers.

2.2 Measurements and laboratory analysis

2.2.1 Sample selection

All laboratory analyses were performed on a subset of samples from 40 selected sites (20 sites for each land use), with the exception of the mid-infrared spectra (see below), which were measured on all 229 samples from the 94 sites. The sites from which we obtained the samples analysed in the laboratory (the 40 selected sites) were selected according to three criteria: location (evenly distributed throughout the study area), soil depth (representing all possible maximum soil depths evenly), and the weathering degree (all levels of weathering degree are represented evenly), which were previously measured in a separate study (see Sect. 2.2.3).

2.2.2 Measuring mid-infrared spectra

After drying, sieving, and grinding, the soil samples were measured with mid-infrared reflectance spectroscopy (MIR, Frontier™, with an auto-sampler from PerkinElmer, Waltham, MA, USA) by Brosens et al. (2021). The measured wavenumbers ranged between 400 and 4000 cm⁻¹, with a resolution of 1 cm⁻¹. Signal scattering was removed, and data quality and calibration were improved using various pre-processing treatments. More details about these analyses can be found in Brosens et al. (2021).

2.2.3 Weathering indices

Brosens et al. (2021) used mid-infrared spectroscopy to predict multiple weathering indices of the collected soil samples using the samples of Vanacker et al. (2019) that were collected in the same catchments to calibrate their spectral model. Two soil weathering indices measured by Brosens et al. (2021) are used in this study: the total reserve in bases (TRB, cmol_c kg⁻¹) and the chemical index of alteration (CIA, dimensionless). The chemical index of alteration reflects the degree of weathering as the relative proportion of Al₂O₃ (i.e. a conservative oxide) to the sum of Al₂O₃, CaO, Na₂O, and K₂O (i.e. major oxides), which can be interpreted as the extent of the conversion of feldspar to clay minerals by evaluating the mobility of the cations Ca²⁺, K⁺, and Na²⁺ (i.e. Eq. 1) (Burke et al., 2007; Nesbitt and Young, 1982). The higher the CIA, the more the labile cations were set free and the higher the weathering degree.

$$\begin{array}{}\text{(1)}& \mathrm{CIA}=\mathrm{100}\cdot \left[{\displaystyle \frac{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}}{{\mathrm{Al}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{3}}+\mathrm{CaO}+{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}+{\mathrm{K}}_{\mathrm{2}}\mathrm{O}}}\right]\end{array}$$

The total reserve in bases (TRB, cmol_c kg⁻¹) is the sum of the alkaline and alkaline-earth element (i.e. Ca²⁺, Mg⁺, K⁺, and Na²⁺) cations present in the soil (basic cations) (Eq. 2). The lower the TRB, the higher the weathering degree (Delvaux et al., 1989).

$$\begin{array}{}\text{(2)}& \mathrm{TRB}={\mathrm{Ca}}^{\mathrm{2}}+{\mathrm{Na}}^{+}+{\mathrm{K}}^{+}+{\mathrm{Mg}}^{\mathrm{2}+}\left[{\mathrm{cmol}}_{\mathrm{c}}{\mathrm{kg}}^{-\mathrm{1}}\right]\end{array}$$

2.2.4 Soil organic carbon fractionation

To quantify the portion of SOC associated with soil minerals, OC fractionation was performed on the soil samples collected by Brosens et al. (2020). The OC associated with fine soil particles (< 53 µm, referred to here as silt and clay organic carbon, OC_S+C) was separated from the bulk soil (Del Galdo et al., 2003; Six et al., 2002a). The resulting OC_S+C represents the OC protected by minerals in the soils (i.e. MAOC) (Cotrufo et al., 2019; Lavallee et al., 2020).

To fractionate the soil samples, 10 g of dried and 2 mm sieved bulk soil were weighed in a 50 mL Falcon tube. Next, 30 mL of 0.5 % sodium hexametaphosphate (NaHMP) and 20 glass beads were added, and the vial was shaken on a reciprocal at 150 rpm for 18 h. The dispersed sample was poured over a 53 µm sieve in a basin (30 cm diameter and 8 cm depth). The soil remaining on the sieve was thoroughly rinsed with deionized water until the water passing through the sieve was clear. The soil remaining on the sieve (i.e. the soil fraction > 53 µm, referred to as the POM fraction) was poured into a pre-weighed aluminium tray. The silt and clay soil fraction (i.e. the fraction < 53 µm, referred to as the S&C fraction) was transferred into a pre-weighed aluminium tray. Both soil fractions were dried at 90 °C for 48 h. Afterwards, the aluminium trays were weighed to obtain the mass of both soil fractions.

2.2.5 Total and stabilized soil organic carbon

Prior to analysing the soil samples for OC% and total N, large organic particles were manually removed from the samples using tweezers. Afterwards, 200 mg of dried and finely ground soil was analysed for OC% and total N using a CHN628 Series Elemental Determinator (LECO Corporation). The OC concentration of the non-fractionated bulk soil was also determined. As the pH was lower than 7 for all soil samples, it was assumed that the soil did not contain any carbonates. The carbon concentration of the S+C fraction was obtained from the fraction < 53 µm. The OC concentration of the S+C fraction in the bulk soil is referred to here as silt and clay organic carbon (OC_S+C). It is defined as follows: the OC_S+C is the mass of OC in the S+C fraction divided by the mass of bulk soil (Fig. S1 in the Supplement). When the ${\mathrm{OC}}_{\mathrm{S}+\mathrm{C}}/\mathrm{TOC}$ ratio is > 1, it was assumed that 100 % of the OC was OC_S+C. As bulk density was not measured on the samples, the OC stocks could not be calculated. Variations in the OC content of multiple samples are expressed as the standard deviation.

2.2.6 pH

To obtain ${\mathrm{pH}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$, 10 g of dried soil was weighed in a Falcon tube, and 25 mL of deionized water was added. The vials were shaken on a reciprocal shaker at 150 rpm for 2 h. The slurry was allowed to settle for 24 h, after which the pH was measured using a pH metre (Thermo Scientific™ Eutech™ 150 series waterproof handheld meters).

2.2.7 Cation exchange capacity and aluminium cations

Before measuring the cation exchange capacity (CEC) and aluminium cations (Al³⁺), visible organic particles were manually removed from the soil sample using tweezers. Next, 2 g of dried and ground soil was weighed in a centrifugation tube, 25 mL 0.1 M BaCl₂ was added, and the samples were shaken on a reciprocal shaker at 150 rpm for 2 h. Next, the sample was centrifuged for 10 min at 2500 rpm, and the supernatant was passed through a Whatman 41 filter. After filtering, 1 mL of the solution was diluted in 4 mL water, which was analysed with an ICP-OES (G8010A Agilent 5100 SVDV ICP-OES, Parent Asset SYS-10-5100). The concentration of Al, Ca, K, Mg, Mn, and Na cations in the solution was converted to cmol_c kg⁻¹ soil before calculating the effective cation exchange capacity by summing the above-mentioned cations (Hendershot and Duquette, 1986).

2.2.8 Soil texture

Prior to the analysis of soil texture, visible organic particles were manually removed using tweezers. Aggregates were gently crushed manually, and the soil was sieved to 2 mm. Next, between 200 and 300 mg of soil was put in 5 mL of 5 % (NaPO₃)₆ and shaken for 4 h, after which the solution was sonicated for 1 min. The soil grain size was determined using an LS13320 (Beckman Coulter). Particle size classes were as follows: clay (< 2 µm), silt (2–53 µm), and sand (53–2000 µm). To assess if the presence of organic matter in the samples affected the measured grain size, 20 samples were measured twice: (1) by treating them with 40 mL 35 % H₂O₂ before analysis and (2) without H₂O₂. As the results of both treatments showed no significant differences, the remaining samples were not treated with H₂O₂ prior to measurement, and dispersion with (NaPO₃)₆ was the only pre-treatment.

Figure 2 Silt and clay organic carbon (OC_S+C) concentration along the depth profile for agricultural and forest soils (dots), with trends indicated by a loess smoother (lines). Vertical bars on the right represent the soil depth layers over which the OC_S+C concentrations of both land uses were compared. Blue vertical bars indicate that the difference between the land uses is significant over the considered depth layer, while yellow bars indicate no significant difference.

2.2.9 Mid-infrared estimations of total and stabilized organic carbon concentrations

MIR spectra were used to determine the OC concentrations of the samples that were not analysed with the dry-combustion method. To that end, a spectral model of the concentration of total organic carbon (i.e. TOC) and stabilized organic carbon (i.e. OC_S+C) was constructed using the MIR spectra of selected samples (n=116). This model was used to predict the TOC and OC_S+C concentration of the remaining samples (n=113). The spectral processing and modelling for both TOC and OC_S+C were done with the R package simplerspec (Baumann, 2019). The spectra were averaged for three replicate measurements for each soil sample. Then, a partial least-square regression (PLSR) calibration model was developed. For the final model, the absorbance spectra were pre-processed using a Savitzky–Golay first derivative with a window of 21 points (Savitzky and Golay, 1964). The models were tuned and evaluated using 10-fold cross-validation. The resulting models for TOC and OC_S+C had an R² of, respectively, 0.94 and 0.88 and an RMSE of, respectively, 0.22 (% OC) and 0.26 (% OC_S+C), with a slight tendency for upper values to be overestimated (Fig. S2 in the Supplement). Predicted values of TOC and OC_S+C smaller than 0 were set to 0. The OC_S+C-to-TOC ratio was calculated for all measured and modelled data. Some of the modelled data had an unrealistically low OC_S+C-to-TOC ratio (Fig. S3 in the Supplement). Therefore, we removed all samples with an OC_S+C-to-TOC ratio below 0.5 (n=7). The resulting dataset was only used in the statistical analysis to determine the impact of land use on OC_S+C along the soil profile (Sect. 2.6.1).

2.3 Statistical analyses

All statistical analyses were performed using R (R Core Team, 2022).

2.3.1 Effect of land use on the concentration of stabilized organic carbon along the soil profile

The depth profiles of OC_S+C were compared for soils under forest and agricultural land use. The dataset was tested for outliers using the boxplot method with the R function identify_outliers, and one sample was removed from the analysis (Kassambara, 2022). In order to get a normal distribution, the OC_S+C data were transformed using Box–Cox transformation with λ=0.4 using the R function bcPower (Fig. S4 in the Supplement) (Fox and Weisberg, 2018). Homogeneity of variance was assessed, and no evident relationship between residuals and fitted values was observed (Fig. S5 in the Supplement). The difference in OC_S+C between both land uses along the depth profile was assessed by performing multiple analyses of variance (ANOVA) over 20 different depth intervals. As the number of data points decreased with depth, the depth layers were chosen in such a way that they all contained 40 ± 5 samples (i.e. 20 ± 5 forest samples and 20 ± 5 agriculture samples). Consequently, the depth intervals to which the ANOVAs were applied had different thicknesses (between 10 cm in the topsoil and 160 cm for the deepest layer) and overlapped with each other (Fig. 2; Table S1 in the Supplement). For each depth layer, an ANOVA was performed. If the confidence interval contained 0, the difference between forest and agriculture was considered not to be significant (Table S1). A power analysis was performed using the R function power.anova.test to determine whether the number of replicates was optimal. Known between-group and within-group variances were used, along with a significance level of α=0.05. This analysis was performed for each depth layer to determine the power of each ANOVA (Table S1).

2.3.2 Soil characteristics controlling the concentration of stabilized organic carbon

A statistical model was created to determine the most important variables controlling the concentration of stabilized OC along the sampled depth profiles under forest and agriculture. For this purpose, a generalized additive mixed-effect model (GAMM) was constructed. The model included OC_S+C as the response variable and selected characteristics of the samples as predictor variables (i.e. land use, texture, weathering indices, local hillslope gradient, Al³⁺ content, and pH), as well as a smoother and a power variance function for depth (more details in Supplementary information 1; Fig. S5; Table S2 in the Supplement). To construct the GAMM, only the 116 samples on which soil properties were analysed in the lab used as MIR predictions were not available for all the soil characteristics.

The predictor variables in the GAMM were tested for correlation and multicollinearity. As the two weathering indicators, the total reserve in bases (TRB) and the chemical index of alteration (CIA), were correlated (Pearson correlation = 0.72), only the CIA was kept to fit the model (Fig. S7 in the Supplement). Because the variance inflation factor (VIF) of Al and CEC was above 5, CEC was removed from the analysis (Fig. S8 in the Supplement). An the initial GAMM, including all selected predictor variables and realistic interactions between these variables, was fitted using the gamm function from the package mgcv (Wood, 2011). Using the summary statistics of the LMM part of the GAMM, the variable with the lowest p value was removed, and the model was fitted again. This stepwise removal was repeated until the p values of all the remaining variables were significant (i.e. p value < 0.05). The model was fitted again with standardized variables to allow the comparison among the regression coefficients.

Before validation, the model was tested for the following assumptions: the linear relationship between explanatory and response variables, the homogeneity of the variance, and the normal distribution of the residuals (Fig. S10 in the Supplement). All these assumptions were confirmed. The final GAMM was validated by comparing fitted values (i.e. the fitted value of the linear mixed-effect (lme) part of the model) with the measured values of OC_S+C (Fig. S9 in the Supplement). The model performed well overall (i.e. R² = 0.78 %, RMSE = 0.29 % OC_S+C) but tended to underestimate the higher values (> 2.5 %) of OC_S+C. The final GAMM was also compared to fitted values of a more complex model (i.e. a model with the same structure but including all initial variables and interactions) (Fig. S9). The values of both models were similar (i.e. R² = 0.84 %, RMSE = 0.26 % OC_S+C), which confirms that the final model performed nearly as well as the more complex model.

2.3.3 Degree of OC saturation of soil minerals as a function of the silt and clay content

To determine the influence of the portion of < 53 µm soil particles (i.e. the S+C fraction) on the maximum OC_S+C, the boundary line method was used (Fig. 4) (Feng et al., 2013). For different levels of S+C fraction, the samples with OC_S+C below a certain threshold (i.e. the boundary line) can be considered to have an OC_S+C saturation deficit while having a physical potential to stabilize more OC. For performing the boundary line analysis, the samples were grouped according to the proportion of S+C fraction in the bulk soil. The first category included the three samples with a proportion of S+C fraction lower than 50 %. The samples with more than 50 % S+C were separated into five categories with intervals of 10 %, resulting in six S+C fraction categories: 0 %–50 %, 50 %–60 %, 60 %–70 %, 70 %–80 %, 80 %–90 %, and 90 %–100 %. The top 20 % of samples with the highest OC_S+C of each group were selected (except for the first group, where all data were selected) and used to fit a linear regression, of which the intercept was forced to 0. The regression line was considered to be the upper boundary of OC_S+C concentration in the soil (Feng et al. 2013). This boundary line represents the maximum OC_S+C concentration at each level of S+C fraction. The saturation level of OC_S+C (%) is defined as the ratio of measured OC_S+C concentration to the maximum OC_S+C concentration at the corresponding S+C fraction level.

The saturation levels of OC_S+C under forest and agricultural land use were compared for different depth layers (Fig. 5b, Table S6 in the Supplement). The dataset was tested for outliers with a boxplot method using the R function identify_outliers, and two samples were removed from the analysis (Kassambara, 2022). In order to get a normal distribution, the stabilized OC saturation data were transformed using Box–Cox transformation with λ=0.38 using the R function bcPower (Fig. S17 in the Supplement) (Fox and Weisberg, 2018). Homogeneity of variance was assessed, and no evident relationship between residuals and fitted values was observed (Fig. S18 in the Supplement).

Figure 3 (a) Total organic carbon (TOC) concentration along the depth profile for agricultural and forest soils (dots), with trends indicated by a loess smoother (lines). (b) Ratio of the concentration of silt and clay organic carbon (OC_S+C) over the concentration of TOC along the depth profile. Measured data have been obtained by direct analysis, while modelled data have been predicted using mid-infrared spectroscopy.

3 Results

3.1 Soil characteristics and depth profiles

The collected data showed no specific pattern with depth for most soil characteristics, except for Al³⁺, which increases with depth (Fig. S12 in the Supplement). The grain size varied between the different profiles and tended to be slightly finer under agriculture. The variability of soil ${\mathrm{pH}}_{{\mathrm{H}}_{\mathrm{2}}\mathrm{O}}$ was largest in the topsoil, and in the overall profiles, it was between 4.4 and 6.3.

3.2 Stabilized soil organic carbon for different land uses

Our results showed that the difference in the concentration of stabilized OC between forest and agricultural soils was only significantly different in the top 90 cm of the soil (Fig. 2), indicating that deforestation did not significantly affect deeper soil OC_S+C. The OC_S+C concentration in the forest topsoil (i.e. 0–30 cm depth) was 2.05 ± 0.64 %, and it was 1.27 ± 0.56 % in the arable topsoil. The average TOC concentration of the topsoil under forest was 2.50 ± 0.80 %, and it was 1.50 ± 0.56 % under arable land (Figs. 3a and S11 in the Supplement). The difference between the two land uses was statistically significant down to 90 cm depth. Below this depth, the TOC% under arable land and forest was 0.39 ± 0.25 % and was, as for OC_S+C, not statistically significantly different between the two land uses (Tables S4 an S5 in the Supplement).

The impact of land use on OC_S+C concentrations below 90 cm could only be tested for soils that were sampled deeper than 90 cm. Therefore, it was necessary to confirm that the above-mentioned result was due to the depth effect as it could be an artifact of site subsetting (i.e. if the sites with deeper soils do not have a significant land use effect in their topsoil). To test this, an ANOVA was performed on the land use effect on OC_S+C concentration using a subset (36 sites) containing only the data points from above 90 cm and excluding the profiles that did not have any measurements below 90 cm. The resulting p value was significant (i.e. < 0.05), reflecting the fact that there is no reason to suspect such a bias. However, the power of the tests for the layers that were not significantly different were all below 0.8, indicating that the number of replicates was not optimal for assessing the difference in land use effect below 90 cm (Table S1).

3.3 Controlling factors in terms of the amount of stabilized soil organic carbon

The parameters selected for the final GAMM model were depth, land use, clay, silt, Al³⁺, the interaction of land use and pH, and the interaction of silt and Al³⁺ (Table 1). These parameters are therefore the important controlling factors in terms of OC_S+C concentration in these soils. The other parameters (i.e. CIA, pH, and slope) and most interactions were not significant to the description of the OC_S+C concentration. The lme part of the model had an R² of 0.78 and an RMSE of 0.29 (% OC_S+C) and tended to underestimate values with higher OC_S+C concentration (Fig. S9). The contribution of the smoother to the fitted value is very high as it ranges between ca. −1 and +1 (% OC_S+C) and decreases with depth; it is positive above 66 cm and becomes negative below 66 cm (Fig. S13 in the Supplement). This means that the depth smoother increases the slope of the regression curve for depths below 66 cm and decreases the slope of the regression curve above this depth threshold. Therefore, the parameters that influenced the OC_S+C concentration the most were depth and land use.

Table 1 Statistics of independent variables selected by the GAMM to predict OC_S+C concentrations.

^∗ The effect of the independent variable is significant.

3.4 Saturation of stabilized OC as a function of clay and silt content

The saturation of stabilized OC on soil minerals was higher in the topsoil, and the degree of stabilized OC saturation decreased with soil depth under both land uses (Fig. 5a). The saturation was significantly different between the land uses above 50 cm depth but not below this depth (Fig. 5b, Table S6). The silt and clay content were selected by the GAMM as good predictors for OC_S+C concentration, and the S+C fraction (i.e. soil fraction < 53 µm) was well adapted for the boundary method as there seems to be a linear correlation between the S+C fraction and the OC_S+C for the selected points (Fig. 4). Therefore, the S+C fraction was used as an independent variable to create a linear model that predicts the maximum OC_S+C concentration. The statistical model showed a boundary line (with a slope of 2.8 %, describing the maximum amount of OC_S+C at any S+C fraction level) that enabled us to estimate the saturation of stabilized OC for each soil sample (Fig. 5a).

Figure 4 The silt and clay organic carbon (OC_S+C) concentration as a function of the fraction of the soil particles < 53 µm (S+C fraction). The boundary line (blue) indicates the estimated maximum amount of OC_S+C that can be stabilized for any given portion of soil particles in the S+C fraction within the 95 % confidence interval. The purple dots were used to construct the linear regression for the boundary line.

Figure 5 (a) Soil profiles of the saturation level of stabilized OC (i.e. OC_S+C), expressed as the percentage of the maximum of stabilized OC_S+C and indicated by the size of the dots for soils under forest and agriculture. (b) Boxplots of the saturation level of OC_S+C per depth layer and land use. The asterisks above the plots indicate statistically significant differences between both land uses.

4 Discussion

4.1 The effect of land use on stabilized organic carbon along the soil profile

Our results show that deforestation over at least 30 years on the studied subtropical soils led to a significant decrease in stabilized OC in the top 90 cm, while no significant effect on stabilized OC deeper in the soil was detected. Most studies on the effect of deforestation on SOC stocks focus on the top 10 to 30 cm of the soil and consider mainly total organic carbon (TOC). On average, observed relative losses of TOC in tropical topsoils after deforestation range between 10 % and 58 % (de Blécourt et al., 2019; Detwiler, 1986; Don et al., 2011; Gurmessa et al., 2016; Kassa et al., 2017; Powers et al., 2011; Rittl et al., 2017; Veldkamp et al., 2020). According to our results, deforestation led to a decrease of 45 % in the OC_S+C concentration in the top 20 cm, a decrease of 42 % between 20 and 30 cm depth, and a decrease of 27 % to 38 % in depth layers between 30 and 90 cm (Table S4). The decrease in TOC concentration was very similar to the decrease in the OC_S+C concentration (Table S5).

These results indicate that stabilized OC in the studied subsoils is much less dynamic than in the topsoil. This implies that increasing the amount of OC stored below 90 cm depth in the studied agricultural soils would be a very slow process. This is evident from the lack of significant differences in OC_S+C between both land uses below this depth, showing that the OC_S+C concentration at these depths has a slow turnover rate. Hence, reforesting arable land is not likely to lead to an increase in stabilized SOC below this depth over timescales relative to ongoing climate change, limiting the potential atmospheric CO₂ withdrawal by reforestation into the subsoil in this ecosystem. Indeed, we observed that, decades after deforestation, the concentration of OC stabilized in the deeper agricultural soil layers remains largely unchanged, and the lowest soil layer with a significant difference in terms of stabilized OC between forest and crop plantation is located at 45–90 cm depth. Our observation of a significant effect of deforestation on SOC down to 90 cm depth is a deeper threshold compared to similar studies in the tropics, in which the difference in SOC content between these land uses was found not to be significant below 30 cm depth (de Blécourt et al., 2019; Dechert et al., 2004; Gurmessa et al., 2016; Kirsten et al., 2019; Nagy et al., 2018). There are several possible reasons for these different findings. One factor that could explain this difference is the study design – for example, in Dechert et al. (2004), the agricultural fields were much younger. Kirsten et al. (2019) collected far fewer samples (i.e. they sampled soil profiles from four forest sites and three cropland sites) for their study, which might not be sufficient to detect a significant difference. Also, these studies all took place in different locations and soil types, where climate and/or vegetation will certainly influence the effect of land use on SOC. Nevertheless, our findings are consistent with the findings of other studies with a similar design (Kassa et al., 2017; Tesfaye et al., 2016) and emphasize that SOC changes after land use conversion need to be monitored at least down to 100 cm depth (Strey et al., 2016). Overall, the results are consistent with the observations of Balesdent et al. (2018); according to their study, OC inputs in deeper soil layers are hardly affected by land use change as OC below 100 cm is very old (> 1000 years).

Current knowledge of SOC dynamics allows us to hypothesize which factors shaped the effect of land use change on OC_S+C concentrations along the depth profile in the studied soils. Concerning the significant difference in OC_S+C concentration in the top 90 cm of the two land uses, a possible hypothesis could be that the difference in OC_S+C was due to soil erosion of the arable soils. During soil erosion, the topsoil is removed, and the new upper layers at the eroded site consist of former subsoil that is typically less saturated in OC_S+C compared to the topsoil (Doetterl et al., 2015b). This phenomenon would have an impact on our results as our samples were collected at mid-slope positions. However, soil erosion rates in our study area are low (0.9–1.4 $\mathrm{t}{\mathrm{km}}^{-\mathrm{2}}{\mathrm{yr}}^{-\mathrm{1}}$) and are therefore unlikely to substantially affect our observations on the considered timescale of a century (Brosens et al., 2020). A more probable hypothesis is that the different soil OC_S+C concentration after deforestation is the result of a change in the quantity and quality of OC inputs. For instance, trees have a much higher root density than crops. Therefore, higher rates of root exudation and root decomposition in forests lead to higher OC inputs to the soil than in arable land (Jackson et al., 1996; Sokol et al., 2019). Furthermore, belowground carbon inputs are more likely to form stabilized SOC than aboveground inputs (Jackson et al., 2017; Sokol et al., 2019; Sokol and Bradford, 2019). This is especially noticeable in tropical soils, where POM accumulation from litter deposition is limited by the quick decomposition of the organic matter (Sokol et al., 2022). Therefore, the difference in OC_S+C concentration over the 0–90 cm depth layer is likely due to the higher belowground OC inputs in forest soils. Below 90 cm depth, OC_S+C is likely to originate primarily from root-derived dissolved OC that slowly leaches downwards along the soil profile while experiencing cycles of sorption and desorption (Uselman et al., 2007; Kaiser and Kalbitz, 2012). Our results indicate that there is no significant difference in soil properties below a depth of 90 cm. This lack of significance suggests that the variance within this depth range is too large to detect meaningful differences. A potential limitation of our study is that our sampling design was not tailored to account for potential controlling factors of soil organic carbon content. As a result, spatial variability likely contributed to a larger variance than if our sampling design had been specifically designed to detect changes in land use (Blécourt et al., 2017). In addition, the deeper soil layers were thicker than the upper soil layers, and the greater the depth intervals between samples of the same layer, the greater the variance. These limitations may have hindered our ability to detect differences in soil layers below 90 cm. We performed a power analysis for each depth layer (see Table S1), which showed that the power fell below 0.8 for depth layers below 40–80 cm. This indicates that the variability within each land use was too high compared to the variability between land uses, suggesting that more samples are required to avoid the risk of a type-II error. A more appropriate sampling design, such as paired plots or sampling at stable landscape locations, would likely have increased the power of the test. However, soil carbon concentration decreases with depth, resulting in small and similar values for both land uses. Consequently, the impact of the sampling design on the power of the test might be minimal at deeper depths. This is because the variance within land use increases relative to the variance between land uses as depth decreases. To enhance statistical power at lower depths, it is crucial to not only optimize the sampling design but also increase the number of samples. Nevertheless, the required sample size increases rapidly as carbon content decreases. For example, to achieve a power of 0.8 in the deepest layer (115–280 cm) of our study, at least 7562 samples would be required, which is an unrealistic number of samples to collect and analyse.

4.2 Factors controlling the concentration of stabilized organic carbon along the soil profile

The general additive mixed model (GAMM) showed that information on land use and soil characteristics that are generally used as indicators of OC_S+C in tropical regions, such as texture and Al³⁺, can also be used to predict the concentration of OC_S+C at soil depths below 1 m. However, our results showed no indication of an important control of the weathering degree of soil minerals on the OC_S+C concentration. This contrasts with Doetterl et al. (2018), who identified long-term weathering as a dominant controlling factor for OC storage in soils. Nevertheless, their results were obtained by studying geochemically different soils across a spatial gradient, while the soil samples collected for our work were highly weathered along the soil profile (the interquartile of CIA ranges from 84 to 96) and were geochemically similar. This suggests that, while the weathering degree of soils is an important control on potential SOC storage across spatial scales (Doetterl et al., 2018), it does not provide insights into the SOC storage potential along the depth profile of highly weathered (tropical) soils, highlighting the scale specificity of SOC stabilization indicators (Manning et al., 2015).

The obtained controlling factors of the OC_S+C concentration are consistent with a large number of studies that emphasize the role of fine soil texture and cations in stabilizing OC in soils (Hassink, 1997; Six et al., 2002c; Kunhi Mouvenchery et al., 2012; von Fromm et al., 2021; Zinn et al., 2007). The relative importance of these factors depends on clay mineralogy (Quesada et al., 2020). Mineralogy and, more broadly, soil type are frequently pointed out as important controlling factors of SOC (Powers et al., 2011; Singh et al., 2018; Wiesmeier et al., 2019). In the literature, climatic variables, in contrast to soil properties, are frequently highlighted as the most important controlling factors of the amount of SOC (Don et al., 2011; Haddix et al., 2020; Hombegowda et al., 2016; Marín-Spiotta and Sharma, 2013; von Fromm et al., 2021). Climatic variables were not considered in this study because samples were all collected in the same region, where no major climatic differences between the sampling locations are expected (Fig. 1). Besides, studies focusing on subsoils have observed that climatic variables impact SOC only in the topsoil and not in the subsoil (van Straaten et al., 2015; Luo et al., 2019). Based on our results, we suggest that variables related to climate and soil weathering are important at large spatial scales, while at regional and local scales (with similar climate and bedrock characteristics), soil texture becomes more important in controlling the amount of stabilized SOC.

4.3 Maximum amount of stabilized organic carbon as a function of clay and silt content

To estimate the potential of the studied soils in stabilizing OC, the S+C fraction (i.e. the soil fraction < 53 µm) of the samples was used as an independent variable to estimate the upper limit of OC_S+C concentration using the boundary line method (Feng et al., 2013). The slope of the boundary line is 2.8 %OC_S+C per unit silt and clay fraction (Fig. 4). This is a gentle slope when compared to the result of the global analysis including multiple land uses from Georgiou et al. (2022), who used a similar approach. These authors found a slope between 4.3 and 5.1 for low-activity minerals (i.e. kaolinite, gibbsite), which are dominant at our study site (Ito and Wagai, 2017). This suggests that the increase in OC_S+C concentration with S+C fraction is lower for the studied soils compared to the average soil with a similar mineralogy. However, land use also plays a role in the magnitude of this slope. For example, Georgiou et al. (2022) considered multiple land use types in their study (e.g. grassland, savanna, shrubland), while Six et al. (2002c) showed that the increase in the concentration of OC associated with soil particles of size < 50 µm was faster for grassland soils than for agricultural or forest soils.

The S+C fraction was selected among other parameters because both silt and clay contents were important controlling factors for OC_S+C according to the GAMM. Doing the same analysis using Al³⁺ (which was also an important controlling factor of OC_S+C concentration) did not result in a similar relation (Fig. S14 in the Supplement). This is due to the increasing concentration of Al³⁺ with decreasing depth, which is common in acidic soils (Kalbitz and Kaiser, 2008). Therefore, although Al³⁺ is a positive controlling factor of OC_S+C concentration, these are negatively correlated because of the opposite pattern with soil depth (Fig. S12 in the Supplement). Using the boundary line method with the chemical index of alteration (CIA) as the independent variable confirmed that this parameter was also not a good indicator of the potential for soils to stabilize carbon for the studied soils. Concerning the relationship between CIA and OC_S+C concentration, it would be expected that the higher the CIA, the higher the potential of the soil to store OC. However, the relationship between the CIA with the highest OC_S+C concentration and the OC_S+C concentration is negative (Fig. S15 in the Supplement). Therefore, using the boundary line method with the CIA did not work for the soil samples in our study. This is in line with the result of the GAMM that the CIA does not control the amount of OC_S+C along soil profiles.

Our results show that only a small proportion of the studied soils had a topsoil OC_S+C concentration close to saturation, while most topsoil and all subsoil samples were substantially undersaturated in OC (Figs. 5a, b and S16 in the Supplement). This indicates that deeper soil layers do not reach their OC storage potential, even under natural vegetation. Therefore, we argue that, while there might be a potential to increase the concentration of stabilized OC in the top 90 cm of the studied agricultural soils by converting them back to forest, it does not seem possible to increase the OC_S+C concentration in the deeper soil layers over a timescale of decades. This supports multiple studies suggesting that adapted agricultural management techniques and strategic crop breeding have the potential to increase the SOC concentration only in the topsoil (Angers et al., 2022; Chenu et al., 2019; Kell, 2012; Hombegowda et al., 2016; Lynch and Wojciechowski, 2015; Poffenbarger et al., 2020; Sayer et al., 2019). Our data, together with other studies, show that the potential to use the subsoil as an atmospheric carbon sink might be limited to the top 100 cm of the soil (Kirschbaum et al., 2021; Lorenz and Lal, 2005; Mathieu et al., 2015; Rumpel and Kögel-Knabner, 2011). In fact, below 50 and 100 cm depth, none of the forest soils reached a degree of OC_S+C saturation above 75 % and 50 %, respectively (Figs. 5b and S16). This indicates that these soil layers are limited by C inputs with respect to OC_S+C saturation and not by the potential of minerals to associate with more OC. Therefore, we argue that using the maximum OC storage potential under natural vegetation as the maximum attainable OC storage below 50 cm depth would not be an appropriate measure to assess potential increases in OC concentration after the conversion of arable land to natural vegetation.

5 Conclusion

This study shows that the conversion of forest to agriculture in a subtropical region affected the concentration of stabilized OC (i.e. OC_S+C) down to 90 cm depth, while no significant differences in the OC_S+C between 90 and 300 cm were detected. We found a difference in the OC_S+C concentration of, respectively, 44.2 %, 41.2 %, and 27.7 % at the depth layers of 0–20, 25–50, and 45–90 cm over a time period of at least 30 years. The most important factors controlling the concentration of stabilized OC along the studied soil profiles were land use, texture (i.e. silt and clay content), and aluminium cations (Al³⁺). In these highly weathered soils, the silt and clay contents were the best predictors of the maximum potential for SOC storage along the soil profile, while no effect of soil weathering degree was detected. This shows that, while soil mineralogy and weathering status may control maximum mineral-associated OC stocks at large spatial scales, this is not the case for the pedon scale in highly weathered subtropical soils. Last, it was shown that, while soil profiles below 90 cm were highly undersaturated in stabilized OC under both forest and agriculture, it is unlikely that subsoils below this depth have the potential to store additional stabilized OC by reforestation over decadal timescales as the OC_S+C content was not affected below 90 cm by deforestation.