DOI: 10.5061/dryad.0k6djhb5k
title: Empirical data and model simulations of the effect of repeated hurricanes on soil carbon dynamics in a humid tropical forest

description: 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 d13C 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.
description: <em>Site description</em> 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<sup>-1</sup> 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
 &gt; 4 cm diameter at breast height (Weaver & Murphy, 1990),
 and dominated by tabonuco (<em>Dacryodes excelsa</em>
 Vahl<em>)</em>. 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.
 <em>Hurricane occurrence </em>
 <strong>Figure 1: Timeline of major hurricanes that have
 affected Luquillo Experimental Forest between sampling dates.
 </strong> 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). 
 <em>Sampling</em> 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 <em>et al</em>.
 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. <em>Soil Density
 fractionation</em> 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<sup>-3</sup> 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%. <em>Soil C and N and
 δ<sup>13</sup>C</em> 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
 δ<sup>13</sup>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 (&lt; 2 mm) and
 dividing by depth and the bulk density as measured in 1988 (Silver et al.,
 1994; Throop et al., 2012).
 <em>Radiocarbon</em> Homogenized
 soil samples were combusted to CO<sub>2</sub> 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<sub>2 </sub>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
 Δ<sup>14</sup>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 Δ<sup>14</sup>C precision of 3.2 %. We report the
 corrected Δ<sup>14</sup>C value and
 ΔΔ<sup>14</sup>C, which is calculated as
 Δ<sup>14</sup>C of the sample minus
 Δ<sup>14</sup>C of the atmosphere, to account for rapidly
 changing atmospheric Δ<sup>14</sup>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<sup>14</sup>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<sup>14</sup>C values, or
 DD<sup>14</sup>C, because atmospheric
 D<sup>14</sup>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<sup>14</sup>C is nonlinear, and thus the
 DD<sup>14</sup>C in 2018 soil will be less sensitive to
 short-term shifts in D<sup>14</sup>C inputs than the samples
 from 1988. <em>Carbon age and transit time
 modeling</em> 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,
 <em> </em> <em>           
                         </em>(2)  dC(t)/dt = [Litter Inputs; 0; 0] +
 [-<em>k</em><sub>FLF</sub>, 0, 0 ;
 a<sub>21</sub>, -<em>k</em><sub>OLF</sub>, 0; 0, a<sub>32</sub>, -<sup>k</sup><sub>HRF</sub>] * [C<sub>FLF</sub>; C<sub>OLF</sub>; C<sub>HF</sub>] where <em>k</em><strong> </strong>is the first-order decay constant for each pool, <em>a</em> is the C transfer rate between pools (<em>i.e. a<sub>21</sub> </em>is the transfer from FLF (pool 1) to OLF (pool 2) and <em>a<sub>32</sub></em> is the transfer from OLF (pool 2) to HF (pool 3)), and <em>C </em>is the C stock of each pool.<strong> </strong>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). <strong>Figure 2: Hypothesized flow of C in soils. </strong> Free light fraction (FLF) C (pink) is either decomposed (at cycling rate -<em>k<sub>FLF </sub>* FLF</em>) or transferred to the occluded light fraction pool (OLF, blue) with the transfer proportion defined by <em>a<sub>21</sub></em>. Carbon transfer between the OLF and heavy fraction (HF, purple) is defined by transfer coefficient <em>a<sub>32</sub></em>, and is respired from these pools at cycling rates -<em>k<sub>OLF</sub>* OLF</em> and <em>-k<sub>HF</sub>* HF</em>, respectively. Figure adapted from Sierra et al. (2012). Soil D<sup>14</sup>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<sup>14</sup>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<sup>-2</sup> 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<sup>-2  </sup>(Hall et al., 2015), suggesting that total root C inputs of approximately 450 g C m<sup>-2 </sup>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<sup>-2­</sup>), and 5–8% of biomass is in the 35–60 cm depth (~40 g C m<sup>-2­</sup>) (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<sup>-2­</sup> to the surface in 1989 (Hurricane Hugo) and 1998 (Hurricane Georges) (Scatena et al., 1993; Silver et al., 1996) and 1611 g C m<sup>-2  </sup>in 2017 (Hurricanes Irma and Maria) (Liu et al. 2018a). The same soil D<sup>14</sup>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 (<em>-k­­<sub>FLF</sub>,</em><sub> </sub><em>-k­­<sub>OLF</sub>,<sub> </sub>-k­­<sub>HF</sub>,<sub> </sub>a<sub>21</sub>, a<sub>32</sub></em>) were constrained using a cost function to accept or reject potential parameter sets over 1000 iterations, based on observed D<sup>14</sup>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 <em>modMCMC() </em>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 <em>systemAge() </em>and <em>transitTime()</em> 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 <em>systemAge() </em>and <em>transitTime()</em> 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. <em>Statistics</em> 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<sup>13</sup>C, DD<sup>14</sup>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<sup>13</sup>C and DD<sup>14</sup>C, FLF DD<sup>14</sup>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 <em>P</em> &lt; 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. 
description: # Empirical data and model simulations of the effect of repeated
 hurricanes on soil carbon dynamics in a humid tropical forest --- ##
 Description of the Data and file structure ### All Bisley Data.csv This is
 a CSV file containing empirical data from 1988 and 2018 sampling dates in
 the Bisley Watershed of Luquillo Experimental Forest. * Year: sampling
 year * Topography: One of four positions, either Ridge, Slope, Upland
 Valley or Riparian Valley * Depth: Midpoint in cm of depth horizon (e.g.
 -5 indicates 0 - 10 cm) * Top of horizon: (cm) * Bottom of horizon: (cm) *
 Fraction: Density fraction - HF (heavy fraction), OLF (occluded light
 fraction), or FLF (free light fraction) * Sample: Full sample name
 including all descriptors * d13C: stable isotope delta 13 C in per mile *
 fraction modern: Fraction of modern carbon as measured by accelerator
 spectrometer radiocarbon data * plus_minus: analytical standard error of
 fraction modern * Del14C: Capital Delta 14C, or radiocarbon signature (per
 mille) * plus_minus: analytical standard error of Del14C * 14C age: age in
 years * plus_minus: analytical standard error of 14C age * Percent N:
 Nitrogen concentration of the sample * s.e. per N: analytical standard
 error of %N * Percent C: Carbon concentration of the sample * s.e. per C:
 analytical standard error of %C * C (g C/m^2): Carbon stock as calculated
 by bulk density given in Silver et al. 1994 * DepthID: depth horizon in
 cm, either 0- 10, 10- 35, or 35 - 60 * Core: Core ID within Bisley
 watershed sampling grid, horizontal and vertical identifiers ##### NOTES
 on blank data: 1\) Very low mass of sample precluded isotopic analysis for
 the following samples: * FLF of 1988 sample of 360, 600 Ridge soil (35-60
 cm) * FLF of 1988 sample of 200, 600 Valley Upland soil (35-60 cm) * FLF
 of 2018 sample of 200, 560 Slope soil (35- 60 cm) * for this sample,
 practically no FLF is present - not enough for C/N analysis either 2\)
 Where s.e. of C and N data is blank, these samples did not have enough
 mass for more than two replicates to run on the elemental analyzer 3\)
 plus_minus of D14C is not given for &gt;Modern samples ### allrpdfs.csv
 This file gives probability density output for 200 out of the 10,000 MCMC
 runs. * Column B: iter = iteration, or ID for mcmc run number * Column C:
 age = age of pool, 0 to 200 years * Column D: rFLFdens = probability
 density for free light fraction at each age in column C * Column E:
 rOLFdens = probability density for occluded light fraction at each age in
 column C * Column F: rFHFdens = probability density for heavy fraction at
 each age in column C * Column G: rsadens = probability density for bulk
 soil at each age in column C * Column H: rttdens = probability density for
 bulk soil transit times at each transit time year in column C * Column I:
 Type = mcmc (Markov chain Monte Carlo) * Column J: Scenario = either
 constant inputs or hurricane inputs (i.e. with a large pulse of inputs at
 each hurricane, followed by reduced inputs for 6 years) * Column K: Depth
 = soil horizon (cm)