Edaphic controls of soil organic carbon in tropical agricultural landscapes

Study area and soil collection

Twenty NRCS map units were selected across Hawaii Commercial & Sugar Company (HC&S) in central Maui that represented seven soil orders, 10 NRCS soil series, and approximately 77% of the total plantation area (Fig. 1). Soil heterogeneity across the landscape allowed for the comparison of a continuum of soil and soil properties that have experienced the same C₄ grass inputs and agricultural treatment under sugarcane production for over 100 years. Conventional sugarcane production involved 2-year growth followed by harvest burn, collection of remaining stalks by mechanical ripper, deep tillage to 40 cm, no crop rotations, and little to no residue return. The sampled soils, collected from September-August 2015, thus represent a baseline of SOC after input-intensive tropical agriculture and long-term soil disturbance. Elemental analyses from this work show consistent agricultural disturbances led to degraded SOC content ranging from 0.23 to 2.91% SOC of soil mass with an average of only 1.16% SOC across all locations and depths.

The homogenized land use history allowed focused investigation of soil property effects on SOC storage across heterogenous soils (Table 1). Though soil inputs (e.g. water, nutrients, root inputs, residue removal) and disturbance regimes (e.g. burn, rip, till, compaction, no crop rotation) were consistent across the 20 field locations, average annual surface temperatures varied from 22.9 to 25.1 °C with a mean of 24.4 °C, average annual relative humidity varied from 70.4 to 79.2% with a mean of 73.4%, and average annual rainfall varied from 306 to 1493 mm with a mean of 575 mm. Gradients of rainfall, relative humidity, and elevation across the site generally increase in an east/north-east direction towards the prevailing winds and up the western slope of Haleakalā. In contrast, surface temperatures increase in the opposite direction towards Kihei and the southern tip of the West Maui Mountains.

^aSoil descriptions²⁶.

^bInterpolated estimates from Ref.²⁵.

Soil sampling and analysis

Pit locations were identified with a handheld GPS and were sampled using NRCS Rapid Carbon Assessment methods²⁷. A total of 75 horizons were identified from the 20 selected map units to a depth of 1 m^28,29. The central depth of each horizon was sampled using volumetric bulk density cores up to 50 cm. After 50 cm, a hand auger was used to check for any further horizon changes. The bulk density of horizons past 50 cm were estimated using collected soil mass and known auger size. Collected soils were air dried, processed through a 2 mm sieve, and analyzed for total C and nitrogen percent, SOC percent, soil texture, iron (Fe) and aluminum (Al) minerals, pH, cation and anion exchange capacity, extractable cations, wet and dry size classes, aggregate stability, and soil water potential at -15 kPa. Total C and nitrogen were measured by elemental analysis (Costech, ECS 4010, Valencia, CA), with SOC content determined by elemental analysis after hydrochloric acid digestion to remove carbonates. Soil texture was measured using sedimentary separation, while a 10:1 soil slurry in water was used to test soil pH. Soil pressure plates were used to measure soil water potential at -15 kPa.

Fe and Al oxides were quantified in mineral phases using selective dissolutions of collected soils, including: (1) a 20:1 sodium citrate to sodium dithionite extraction, shaken 16 h, to quantify total free minerals³⁰, (2) 0.25 M hydroxylamine hydrochloride and hydrochloric acid extraction, shaken 16 h, to quantify amorphous minerals³¹, and (3) 0.1 M sodium pyrophosphate (pH 10), shaken 16 h and centrifuged at 20,000g, to quantify organo-bound metals³⁰. Extracted Fe, Al, and Si from al extractions were measured by inductively coupled plasma analysis (PerkinElmer, Optima ICP-OES, Norwalk, CT). Exploratory ratios of Fe/Al, Fe/Si, and Al/Si for the citrate/dithionite (c), hydroxylamine (h), and pyrophosphate (p) extractions were calculated. Crystalline Fe, operationally-defined as the difference between the citrate dithionite and hydroxylamine extraction, and Al + ½ Fe³² were calculated for each extraction.

Plant-available phosphorus was extracted by the Olsen method using 0.5 M sodium bicarbonate adjusted to pH 8.5 and measured by continuous flow colorimetry (Hach, Lachat Quickchem 8500, Loveland, CO). Exchangeable cations (i.e. calcium, magnesium, potassium, and sodium), effective cation exchange capacity, and anion exchange capacity were measured by compulsive exchange using barium chloride and magnesium sulfate³³. Cations were quantified by continuous flow colorimetry and flame-spectroscopy (Hach, Lachat Quickchem 8500, Loveland, CO). Field soils were air dried and initially passed through a 2 mm sieve before size classes of macroaggregate (2 mm – 250 µm) and microaggregate (< 250 µm) were separated using a 250 µm sieve. Further wet sieving at 250 µm was conducted on 4 g of dry sieved macroaggregates using a wet sieving apparatus (Eijkelkamp, Wet Sieving Apparatus, Morrisville, NC). In short, dry sieved macroaggregates were wet sieved at 60 oscillations a minute for 45 min in distilled water. Soil particles passing through the 250 µm sieve were classified as non-water-stable aggregates, while soil particles retained were dispersed by 16 h of shaking in 2% sodium hexametaphosphate. After soil dispersion, roots and rocks were separated from remaining soil on a 250 µm sieve using distilled water and light agitation with a rubber policeman. Roots were further separated from rocks using a combination of density separation in water and visual removal with tweezers. Though the Roots metric technically includes particulate organic matter (POM), most recovered material was visually identifiable as root biomass.

Soil incubation

The equivalent of 50 dry grams of soil were incubated for 90 days at 25 °C in 500 mL septa capped cell culture bottles after moisture adjustment to -15 kPa moisture content using soil pressure plates^28,29. Evolved gases were sampled using 10 mL syringes, stored in evacuated 3 mL Exetainers (Labco Limited, Lampeter, UK), and analyzed by GC (Shimadzu, GC-2010 Green House Gas analyzer, Kyoto, Japan) using a 4-point calibration curve of mixed CO₂, N₂O, and CH₄ standard gas. After each sampling event, incubation bottles were flushed with air, resealed, sampled for initial concentration, and returned to controlled environment chambers (Caron, 7000-33). All incubations were sampled on the same interval: day 0, 2, 4, 6, 8, 10, 12, 14, 17, 20, 23, 26, 30, 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90. The mass of CO₂ evolved from the soil in each incubation chamber was calculated with the gas law using the measured concentration, known volume, and Standard Ambient Temperature and Pressure (SATP). The mass of C from CO₂ (C-CO₂) was normalized by the dry weight of soil in each chamber and cumulative flux curves of were tracked throughout the 90-day incubation experiment. C pools and fluxes were estimated using compartment models fit to cumulative C–CO₂ incubation curves as they describe C mineralization kinetics of soils ¹⁸.

Modeling soil carbon pools using SoilR

To estimate kinetic SOC fractions, Bayesian solutions of two- and three-pool compartment models (Fig. 2) were implemented in R. Models were fit to cumulative C-CO₂ curves using the SoilR package and a generalized model of linear dynamical systems¹⁸:

where I(t) is a time-dependent column vector of inputs to each compartment n, A(t) is a n × n square matrix that contains compartment decomposition rates and transfer coefficients between compartments, and C(t) is a n × 1 vector of C compartment size (pools) at a given time t. Three-pool SOC models can be expressed generically by a system of differential equations of the form:

With initial conditions:

where k_i represents the decay rate in each pool, α_i,j represents the transfer rate of decayed C in pool j that flows to pool i, while γ_i represents the proportion of C in each pool, and I_i represents the input to each compartment. Two-pool models take the same form as Eqs. (2) and (3), except the third row of every matrix, and third column of the A(t) matrix, are removed. In this sense these matrix models can be scaled up or down to include as many compartments as desired, with the caveat that complex models with many compartments and transfers can become untenable.

Equation (1) was thus fit to the cumulative CO₂ curve of each soil sample using a two-step parameter optimization involving initial classical optimization with a Nelder-Mead algorithm³⁴ followed by Bayesian optimization using a Markov Chain Monte Carlo (MCMC) procedure³⁵. Both optimizations were conducted in R³⁶ using a combination of the FME³⁵ and SoilR¹⁸ packages.

Collinearity analysis indicated that two-pool models were the least collinear (i.e., the most likely to result in meaningful solutions) but would still require fixing at least two parameters. From the collinearity analysis, the best two-pool models used pool sizes (γ₁, γ₂) and decay rates (k₁, k₂), but required that the transfer rates (α_2,1, α_1,2) be fixed. To address this issue, we used the transfer rates from a soil fraction-based model developed by Crow et al.³⁷ at one of the 20 sites investigated in this study. The transfer rates, which describe the amount of decayed C that flows through any given transfer path, are difficult to quantify and would require soil density fractionation and fraction-based compartment models, or controlled input and labeling experiments, to quantify, which have not been conducted. Thus, models were informed with best-available data by fixing the forward transfer rate (α_2,1) and backward transfer rate (α_1,2) at 0.85 and 0.15, respectively³⁷.

The resulting kinetic SOCC fraction estimates (i.e., pool sizes and decay rates) are framed here as dependent variables that respond to differences in independent soil physicochemical parameters measured across this agricultural landscape. The slow pool decay rate (k₂) and fast pool decay rate (k₁) from two-pool models are presented in comparison to SOC. To discuss only the most identifiable solutions, the 3-pool models were not presented, and only 2-pool model results were used.

Linear mixed models and model selection

Several statistical models were created to investigate relationships between SOC content, SOC decay rates, and soil physicochemical properties (Table 2) that may promote SOC persistence despite long-term agricultural disturbance. For each C response variable (SOC, k₂, and k₁), three sequential linear mixed models (LMM), with pedon as a random effect, were constructed: (1) using no fixed effect (intercept-only; equivalent to a random effects ANOVA), (2) using only physicochemical properties as fixed effects that showed significant spearman’s ρ to the conditional residuals of the intercept-only LMM, and (3) a final model where root-mineral interaction terms and physicochemical properties with significant spearman’s ρ were further reduced by model selection using the Akaike Information Criterion for small sample sizes (AICc).

Due to computational and theory-based restrictions on model selection, it was necessary to first reduce the possible independent variables (i.e., soil physicochemical properties) to only those that have prediction potential before performing model selection. This was accomplished using the under-defined models, i.e. LMM(1), of each C response variable (SOC, k₂, and k₁) and identifying which physicochemical properties showed spearman correlation to conditional model residuals, like Rasmussen et al.²¹. Soil properties that could explain residual model variance were thus candidates for final model predictors. The conditional residuals were also parsed above and below the median of each climate variable (i.e., rainfall, humidity, and surface temperature) to evaluate climate effects. Correlation significance was assessed using Bonferroni’s corrected α (0.05/308 = 0.000162), where 308 is the number of correlations conducted for each response variable. The correlative analyses allowed visualization of potential relationships between C response variables and physicochemical properties, as well as how climate may mediate these relationships.

Using this information allowed LMM (2) to be built with physicochemical properties that explained some residual model variance. For each response variable, LMM (2) was developed for the sole purpose of creating a global model for model selection by AICc. Interpreting LMM (2) of each response variable is not very productive as many soil metrics are collinear and a model with this many parameters will suffer from model overfitting. Model selection by AICc was thus chosen to identify combinations of predictors that could best-explain C response while avoiding overly complex models and overfitting.

For parsimony, and under the assumption that similar soil properties should control both the decay rates and overall SOC outcomes, a unified list of significant factors across the intercept-only models of SOC, k₂, and k₁ was created as the basis of the third and final model, LMM (3). Root-mineral interaction terms were added to LMM (3) before model selection based on a-priori expectations that roots and mineralogy are major mediators of C input and storage in soils³⁸. Models were selected using AICc in R using the dredge function in the MuMIn package³⁹. Due to the high correlation between many soil properties (e.g., roots and depth, depth and mineralogy, mineralogy and mineral ratios, etc.), the model selection process for SOC also selected between highly correlated predictors to avoid confounding regression coefficient estimates. This was accomplished using logical subsetting in the dredge function of the MuMIn package to eliminate nested models where both correlated predictors occurred.

Climate variables are not included in any LMM. Instead, climate variables are framed as soil forming factors that established soil physicochemical properties and not as direct drivers of SOC storage or decay. For the same reason that collinear soil properties should not be included in a single model, including climate variables that correlate with soil physicochemical properties (e.g., rainfall and clay) will confound regression coefficient estimates. Further, we assumed that direct climate influence on SOC have been constrained by consistent agricultural practices, mono-cropping, and landscape-scale irrigation that have normalized C inputs and losses over 100 years of sugarcane cultivation. We also expect that agricultural disturbance has reduced SOC to minimal levels with remaining C largely persisting due to physicochemical protection mechanisms like aggregation, organo-mineral complexation, and organo-metal ligand bonding. The moisture-adjusted and temperature-controlled incubation also removed direct influence of climate on comparative measurement of SOC dynamics across soil horizons. With these assumptions, we focused on soil physicochemical properties as key controls of SOC and dynamic C response in this system.

Hierarchical clustering to evaluate depth and rainfall gradients

Visualization of relationships within the data between C responses and best predictors by depth and rainfall was accomplished using agglomerative hierarchical clustering. The clustering algorithm used Euclidean distance measures and Ward’s minimum variance method as implemented in the hclust function within the Stats package in R³⁶. Unscaled C responses, final model predictors, and depth data were clustered by minimizing total within-cluster variance, where additional clusters beyond the fifth showed little reduction in variance. The clustering procedure was repeated for surface soils by subsetting C response and final model predictors to only the topsoil, adding associated rainfall data, and clustering again based on minimized within-cluster variance. In the case of surface soils clustered by rainfall, additional clusters past the third showed little reduction in variance. Both clustering procedures (i.e., depth and rainfall) used unscaled data to take advantage of the numerically higher values of depth (0–100 cm) and rainfall (306–1493 mm) compared to other data types. This negated the need to weight or scale the data and resulted in parameters clustered predominantly by the depth or rainfall gradient, respectively. Clear separation of clusters by depth and rainfall were found, which allowed for interpretation of the relationships between C response, soil physicochemical predictors, and depth and rainfall gradients.

Source: Ecology - nature.com