Increasing hurricane frequency and intensity with climate change is likely to affect soil organic carbon (C) stocks in tropical forests. We examined the cycling of C between soil pools and with depth at the Luquillo Experimental Forest in Puerto Rico in soils over a 30-year period that spanned repeated hurricanes. We used a non-linear matrix model of soil C pools and fluxes (“soilR”) and constrained the parameters with soil and litter survey data. Soil chemistry and stable and radiocarbon isotopes were measured from three soil depths across a topographic gradient in 1988 and 2018. Our results suggest that pulses and subsequent reduction of inputs caused by severe hurricanes in 1989, 1998, and two in 2017 led to faster mean transit times and younger mean ages of soil C in the particulate, occluded, and mineral-associated soil organic matter pools at 0–10 cm and 35–60 cm depths relative to a modeled control soil with constant inputs over the thirty years. Between 1988 and 2018, the occluded C stock increased, and d¹³C in all pools decreased, while changes in particulate and mineral-associated C were undetectable. The differences between 1988 and 2018 suggest that hurricane disturbance results in a dilution of the occluded light C pool with an influx of young, debris-deposited C, and possible microbial scavenging of old and young C in the particulate and mineral-associated pools. These effects led to a younger total soil C pool with faster mean transit times. Our results suggest that increasing frequency of intense hurricanes will speed up rates of C cycling in tropical forests, and eventually lead to net losses of C from tropical forest soils.

Site description

Soils were sampled from the Bisley Experimental Watershed of the LEF, Puerto Rico (18.3157 deg. N, 65.7487 deg W), a Long-Term Ecological Research and Critical Zone Observatory and Network site (https://luq.lter.network). The mean maximum daily temperature at Bisley was 27 ºC between 1993 and 2010 (Gonzales, 2020), with little seasonality. The mean annual precipitation at Bisley was 3883 (± 864 s.d.) mm y^-1 from 1988 through 2014 (González, 2017; Murphy et al., 2017). Rainfall occurs all year, though January through April experience slightly less precipitation than other months (Heartsill-Scalley et al., 2007). The site is a humid tropical forest with a diverse tree community of approximately 170 species > 4 cm diameter at breast height (Weaver & Murphy, 1990), and dominated by tabonuco (Dacryodes excelsa Vahl). Elevation of Bisley spans from 261 m a.s.l. at the base to 450 m a.s.l. on the ridges (Scatena, 1989).

Soils in Bisley are derived from volcaniclastic sediments of andesitic parent material (Scatena, 1989).  Ridge soils are classified as Ultisols (Typic Haplohumults), while slope soils are classified as Oxisols (inceptic and Aquic Hapludox), and valley soils are classified as Inceptisols (Typic Epiaquaepts) (Hall et al., 2015; McDowell et al., 2012; Scatena, 1989). Detailed site descriptions can be found in Scatena (1989), Heartsill-Scalley et al (2010), and McDowell et al (2012). Here we refer to soil organic C (SOC) and soil C interchangeably because there is no detectable inorganic C in these soils.

Hurricane occurrence 

Figure 1: Timeline of major hurricanes that have affected Luquillo Experimental Forest between sampling dates.

Nine major hurricanes (category 3 or higher) have impacted Puerto Rico between 1851 and 2019 (López-Marrero et al., 2019), and five of these hurricanes have impacted the LEF. Until 1998, hurricanes had historically directly impacted the LEF approximately every 60 years (Scatena & Larsen, 1991). Before the initial sampling campaign of this study, Hurricane San Ciprián in 1932 was the most recent storm to cause major disturbance to the LEF (Scatena & Larsen, 1991).  However, since sampling in 1988, four major hurricanes have impacted the forest (Figure 1). Hurricane Hugo (Category 3-4) in 1989, Hurricane Georges (Category 3) in 1998, and Hurricanes Irma and Maria (Categories 5 and 4, respectively) within two weeks in 2017. The trajectory and windspeeds of all these hurricanes caused widespread defoliation. Litterfall historically takes over five years to return to pre-hurricane levels (Scatena et al., 1996). 

Sampling

Sample collection occurred in 1988 and again in 2018. In both years, samples were collected from three depths: 0–10 cm (the A horizon), 10–35 cm (all of the B1 horizon and part of B2), and 35–60 cm (B2 to C) using an 8 cm diameter soil auger. Soils in this study were sampled at three separate sites at least 40 m from one another for each of three topographic locations, ridge, slope, and upland valley. Two separate cores were taken from a fourth topographic location in the riparian valley, that characterized a smaller proportion of the area of these watersheds (Scatena & Lugo, 1995). Riparian valley sites were ephemeral streambeds with a high boulder presence that limited sampling to less than 25 cm depth in one case. Sampling sites from 1988 were marked with flags, and samples from 2018 were collected from within 15 m of the same locations as the replicates from 1988, for consistency.

Samples collected in 1988 were analyzed for bulk density, pH, soil moisture, and a suite of soil chemical properties (see Silver et al. 1994). Samples were then air-dried and stored in closed Ziploc bags within paper bags in a storage facility in Richmond, CA, USA before density fractionation in 2018. Fresh samples collected in 2018 were also characterized for pH, soil moisture, and soil chemistry. Approximately 3 g subsamples from each fresh sample in 2018 were immediately extracted with 45 mL of 0.2 M sodium citrate/0.5 M ascorbate solution, shaken for 16 hours, then centrifuged and the supernatant decanted to measure concentrations of poorly crystalline iron (Fe) oxides. Within two days of being double-bagged in Ziploc bags, fresh samples were further subsampled and analyzed for pH in a 1:1 soil-to-water slurry (Thomas, 1996) and for gravimetric soil moisture by oven-drying ~10 g subsamples at 105 ºC until a constant weight. Soil samples were air-dried before further processing and analysis. Air-dried soils from both sampling years were sieved to 2 mm and large roots were sorted out.

Soil Density fractionation

Soil was fractionated by density following the method of Swanston et al. (2005), as modified by Marin-Spiotta et al., (2009). Approximately 20 g of air-dried soil was added to centrifuge tubes. Sodium polytungstate (SPT, Na6 [H2W12O40] TC-Tungsten Compounds, Bavaria, Germany) in solution of density 1.85 g cm^-3 was added to centrifuge tubes and agitated before centrifuging. The density of the SPT followed previous studies from this and nearby sites to allow direct comparison (Gutiérrez del Arroyo & Silver, 2018; Hall et al., 2015). Particulate organic matter floating at the surface after centrifugation, the free light fraction (FLF), was aspirated and then rinsed with 100 ml of deionized water 5 times on a 0.8 µm pore polycarbonate filter (Whatman Nuclepore Track Etch Membrane, Darmstadt, Germany). Rinsed FLF was oven-dried at 65 ºC until weight had stabilized. The remainder of the sample was combined with 70 ml of additional SPT and mixed using an electric benchtop mixer (G3U05R, Lightning, New York, NY, USA) at 1700 rpm for 1 min and sonicated in an ice bath for 3 min at 70% pulse (Branson 450 Sonifier, Danbury, CT, USA). Sonication is intended to disrupt soil structure and liberate organic matter that has been occluded in aggregates. The sonicated slurry was centrifuged again, and the light fraction at the surface, the occluded light fraction (OLF), was aspirated, rinsed, and dried using the same method as for the FLF. The remaining soil pellet was considered the heavy fraction (HF), or mineral-associated organic matter fraction. The HF was rinsed by thoroughly mixing with 150 ml of deionized water in the centrifuge tube, centrifuging, and removing the supernatant repeatedly until the fraction had been rinsed 5 times. The rinsed HF was oven-dried at 105 ºC until weight stabilized. The average mass recovery was 98%.

Soil C and N and δ¹³C

Dried bulk and HF soils were homogenized separately using a Spex Ball mill (SPEX Sample Prep Mixer Mill 8000D, Metuchen, NJ). The FLF and OLF were homogenized separately by hand using a mortar and pestle. All homogenized samples were then analyzed at U. C. Berkeley for C and N concentrations on the CE Elantech elemental analyzer (Lakewood, NJ) and for δ¹³C in the Stable Isotope Laboratory at UC Berkeley, using a CHNOS Elemental Analyzer interfaced to an IsoPrime 100 mass spectrometer (Cheadle Hulme, UK), with a long-term external precision of 0.10 %.  Soil C stocks were calculated by multiplying the C concentrations (%) by the oven-dry mass of bulk soil (< 2 mm) and dividing by depth and the bulk density as measured in 1988 (Silver et al., 1994; Throop et al., 2012).

Radiocarbon

Homogenized soil samples were combusted to CO₂ in sealed glass tubes along with silver (Ag) and copper oxide (CuO) at the Center for Accelerator Mass Spectrometry at Lawrence Livermore National Lab. The CO₂ was then graphitized on Fe powder under pressurized hydrogen gas (Vogel et al., 1984). Graphite was pressed into aluminum targets and run on the Compact Accelerator Mass Spectrometer for radiocarbon analysis (Broek et al., 2021). Radiocarbon is reported in Δ¹⁴C, following Stuiver & Polach (1977), and calculated based on the fraction of modern isotope composition, corrected for the year of sampling, and corrected for mass-dependent fractionation with observed δ13C values of the sample. The compact AMS had an average Δ¹⁴C precision of 3.2 %. We report the corrected Δ¹⁴C value and ΔΔ¹⁴C, which is calculated as Δ¹⁴C of the sample minus Δ¹⁴C of the atmosphere, to account for rapidly changing atmospheric Δ¹⁴C during the study period. Atmospheric radiocarbon has been decaying nonlinearly since the peak of weapons testing in the 1950s. Radiocarbon signatures in the soil are strongly influenced by the atmospheric D¹⁴C signature, making them useful for modeling soil C age and transit time, especially since the 1950s. To compare the contribution of modern C between 1988 and 2018, it is useful to take the difference between soil and atmospheric D¹⁴C values, or DD¹⁴C, because atmospheric D¹⁴C declined between 1988 (98 %) and 2018 (4.4 %) in Northern Hemisphere Zone 2 (Hua et al., 2013). We note that the decline in atmospheric D¹⁴C is nonlinear, and thus the DD¹⁴C in 2018 soil will be less sensitive to short-term shifts in D¹⁴C inputs than the samples from 1988.

Carbon age and transit time modeling

Transit times and ages of C were modeled with the package “SoilR” (Sierra et al., 2012, 2014) in R, version 4.0.2. The change in C density fractions over time, termed C flow, was modeled using a 3-pool structure with a series flow matrix, under the simplifying assumption that C flows from the litter pool to the FLF, where it is sequentially transferred into the OLF and HF pools (Figure 2). The model structure is depicted in basic form in equation 1,

                                                                            (1)  dC(t)/dt = Inputs - k*C

 in matrix form with explicit pools in equation 2,

                                      (2)  dC(t)/dt = [Litter Inputs; 0; 0] + [-k_FLF, 0, 0 ; a₂₁, -k_OLF, 0; 0, a₃₂, -^k_HRF] * [C_FLF; C_OLF; C_HF]

where k is the first-order decay constant for each pool, a is the C transfer rate between pools (i.e. a₂₁ is the transfer from FLF (pool 1) to OLF (pool 2) and a₃₂ is the transfer from OLF (pool 2) to HF (pool 3)), and C is the C stock of each pool. The transitTime and systemAge functions within the “soilR” package use this model structure to solve for the distribution of ages (time since entry) of each pool, and the distribution of transit times (times between entry and exit from the bulk soil) (Sierra et al 2016). Distributions of age and transit time were time-independent and did not assume a specific distribution (Sierra et al., 2014, 2017).

Figure 2: Hypothesized flow of C in soils.

Free light fraction (FLF) C (pink) is either decomposed (at cycling rate -k_FLF * FLF) or transferred to the occluded light fraction pool (OLF, blue) with the transfer proportion defined by a₂₁. Carbon transfer between the OLF and heavy fraction (HF, purple) is defined by transfer coefficient a₃₂, and is respired from these pools at cycling rates -k_OLF* OLF and -k_HF* HF, respectively. Figure adapted from Sierra et al. (2012).

Soil D¹⁴C and C stock mean and standard deviations from each time point, depth, and fraction were used to constrain the matrix model describing the movement of C through three soil pools and losses of C from each pool. Topography was not a strong predictor of patterns in D¹⁴C, C stocks, or C fractions, so samples from all topographies were aggregated for model simulations. The model used mean observed C content in each pool for each depth in 1988 as initial conditions for SOC stocks. Above and belowground litter inputs at 0–10 cm were assumed to be 900 g C m^-2 in non-hurricane or hurricane recovery years, based on observations from the same site (Liu et al., 2018; Scatena et al., 1996; Silver et al., 1996; Vogt et al., 1996).

Inputs to the 10–35 cm and 35–60 cm depths were estimated using observations of live fine roots on the surface and typical root distribution in the forest (Silver & Vogt, 1993). Total root input is approximately threefold the input of fine roots alone (McCormack et al., 2015; Yaffar & Norby, 2020), and live fine roots in the 0–10 cm depth had a mean biomass of 80 - 250 g C m^-2  (Hall et al., 2015), suggesting that total root C inputs of approximately 450 g C m^-2 to the surface would be well within the expected range. Root inputs below 0–10 cm were estimated assuming that inputs follow the typical distribution of root biomass in Puerto Rican tropical forests, with 60–70% of root biomass in 0–10 cm, an additional 20-30% of biomass in 10–35 cm (~135 g C m^-2­), and 5–8% of biomass is in the 35–60 cm depth (~40 g C m^-2­) (Silver & Vogt, 1993; Yaffar & Norby, 2020).

The model was parameterized under two scenarios for each depth: 1) constant inputs, assuming a steady-state undisturbed forest, and 2) hurricane inputs, which simulated the input fluxes from defoliation during the three major hurricanes, followed by a subsequent reduction in litter inputs and then litterfall increasing linearly to pre-hurricane inputs over 6 years (Scatena et al., 1996; Silver et al., 1996; Vogt et al., 1996). Hurricane inputs were imposed as an additional pulse of litter inputs to each depth interval, declining with depth.  The 0–10 cm interval received 100% of the surface input pulse, the 10–35 cm depth received a pulse of root inputs equivalent to 30% of the surface input pulse, and the 35–60 cm depth received root inputs equal to 10% of the surface input pulse. Surface litter pulses under hurricanes were specified according to measured litterfall values and were 42.5 g C m^-2­ to the surface in 1989 (Hurricane Hugo) and 1998 (Hurricane Georges) (Scatena et al., 1993; Silver et al., 1996) and 1611 g C m^-2  in 2017 (Hurricanes Irma and Maria) (Liu et al. 2018a). The same soil D¹⁴C and C stock observations were used to constrain the model under each scenario, with only the input regime varying.

Parameters of the transfer matrix (-k­­_FLF, -k­­_OLF, -k­­_HF, a₂₁, a₃₂) were constrained using a cost function to accept or reject potential parameter sets over 1000 iterations, based on observed D¹⁴C and C stock means and standard errors from both time points (1988 and 2018). A Markov chain Monte Carlo (MCMC) simulation initialized with cost-optimized parameters was run to assimilate observed data and optimize parameter choices to the observations using function modMCMC() from R package “FME” (Sierra et al., 2014; Soetaert & Petzoldt, 2010). The MCMC was iterated over at least 20,000 simulations or until parameter solutions converged according to the trace, which was over 100,000 iterations at the 35–60 cm depth. The first half of the iterations was considered the burn-in period before the chain started to converge near an equilibrium, and these iterations were discarded in calculations of optimal parameters. The model output for the surface soils of the HF pool was validated using published radiocarbon values from the mineral-associated fraction (the only fraction analyzed) of samples from the site taken in 2012 (Hall et al., 2015). 

Bulk and pool soil C age and transit time density distributions and mean values were calculated using the systemAge() and transitTime() functions from the “SoilR” package. Mean density distributions were calculated using the mean parameter set given from the MCMC analysis. Standard deviation from the mean was calculated using the systemAge() and transitTime() functions on 200 sets of five parameters selected randomly within one standard deviation of the mean of each parameter given as output from the MCMC. Lower and upper limits of SOC ages and transit times were calculated using the upper and lower ranges of these iterations.

Statistics

Statistics were run in R, version 4.0.2 (R Core Team, 2020). The statistical model selection followed the recommendations of Zuur et al (2009). Statistical models were chosen using a linear mixed effects model in package “lme4”, with random slopes accounting for the influence each core, or sampling site, had on the response variable values as they varied with depth. This random effect of the core site on the depth effect was evaluated using a restricted maximum likelihood approach and was included in the initial evaluation of all model comparisons. Linear mixed effect models included year, topographic position, depth, and interactions as fixed factors, and the depth effect of each core as a random factor for each of the response variables: C concentration, N concentration, d¹³C, DD¹⁴C. In evaluations of some response variables with AIC and BIC criteria, the random effect no longer enhanced the model, and model comparison proceeded using ANOVAs of linear models without random effects. Topographic effects on C concentrations are discussed in the supplemental information. Model assumptions were evaluated using the check_model function in R package “performance”, to check for multicollinearity, normality of residuals, homoscedasticity, homogeneity of variance, influential observations, and normality of random effects.

In the cases when random effects were significant (bulk soil d¹³C and DD¹⁴C, FLF DD¹⁴C and HF C and N concentrations), fixed effects were chosen using ANOVA of subsequent models using maximum likelihood estimation, with the random effects held constant. Once fixed effects were established, the model was re-fitted using a restricted maximum likelihood approach to report model estimates, and an ANOVA was run to determine the significance of the response variable. In all cases, P-values were estimated using Tukey’s honest significant post-hoc test to assess significant differences between variables, in the package “agricolae” in R, and contrasts and standard errors of contrasts were estimated using lsmeans() function in package “lsmeans” in R. Values of P < 0.10 were reported as significant unless otherwise specified. The topographic position was not a significant predictor for most variables, so results are reported as means aggregated across positions.