Ecological Archives A023-096-A1

April M. Melvin, Jeremy W. Lichstein, Christine L. Goodale. 2013. Forest liming increases forest floor carbon and nitrogen stocks in a mixed hardwood forest. Ecological Applications 23:1962–1975. http://dx.doi.org/10.1890/13-0274.1

Appendix A. Modeling the dynamics of the Oe and Oa soil carbon stocks.

Because we did not measure all C stocks and fluxes needed to parameterize the model (Eqs. 1–2), we combined our data with values reported by Fahey et al. (2005) (hereafter, "Fahey"). Fahey report a complete soil C budget for forest floor and mineral layers for the Hubbard Brook Experimental Forest, New Hampshire, USA. We used the ratios of Oe:Oa stocks and fluxes observed in our control plots to partition Fahey’s forest floor fluxes into separate Oe and Oa layers. We then solved the system for equilibrium as the starting point for simulating 19 years of post-liming C dynamics (corresponding to the 19 years between the liming treatment at our study site and the time of our sampling). Starting our simulations at equilibrium eliminated liming-independent transient dynamics and allowed us to isolate the model's response to liming.

The underlying assumption of our modeling approach was that the ratios observed in our study (ratios of Oe:Oa stocks and fluxes, and ratios of limed:control stocks and fluxes within the Oe and Oa layers) were applicable to the Fahey system. To illustrate the implications of this assumption, let F and S be a flux and stock in the Fahey system, and let βF and βS be limed:control ratios of the corresponding flux and stock in our study. We assume that if liming were applied to the Fahey system, then the post-liming flux would be βFF throughout the 19-year simulation, and the predicted stock after 19 years would be βSS. Alternatively, we could have adjusted the Fahey fluxes so that they yielded an equilibrium condition equal to the mean stocks in our control plots, and then compared the model predictions (19 years post-liming) directly to the observed stocks in our limed plots. This alternative approach would have yielded similar qualitative results to those we report, but would have been more complicated to implement.

Estimating control parameter values

Here we describe the methods used to estimate control (non-limed) parameter values in Eqs. 1–2, which were used to define the equilibrium starting values for simulating the response to liming. All Fahey values come from their Fig. 4. All values from our study come from our control (non-limed) plots. Lit and CWD were taken directly from Fahey. Roote and Roota were estimated by multiplying Fahey’s root mortality flux into the forest floor by the fraction of forest floor fine roots in our control plots located in the Oe and Oa layers, respectively. Rhizoe and Rhizoa were estimated in a similar way; i.e., we multiplied Fahey’s rhizosphere flux into the forest floor layer by the Oe and Oa fine-root ratios from our control plots.

To estimate Oe and Oa heterotrophic respiration rates (re and ra; yr−1), we combined Fahey’s forest floor respiration flux (Fffr; t C ha−1 yr−1) and forest floor organic matter (Cff; t C ha−1) with the following observations from our study: the fractions of forest floor organic matter in the Oe and Oa layers (fOMe and fOMa = 1 − fOMe, respectively), and the Oa:Oe ratio of soil basal respiration rates (ρae-r). We assume that ra = re×ρae-r, so we need only estimate re. The forest floor respiration flux is

Fffr = re×Cff×fOMe + ra×Cff×fOMa = re×Cff×fOMe + re×ρae-r×Cff×fOMa

which can be rearranged to yield

 re = Fffr/(Cff×fOMe + ρae-r×Cff×fOMa).

To estimate the transfer rate from the Oe to the Oa layer (τea; yr−1), we first note that the turnover rate of Oe carbon (αe; yr−1) is equal to the sum of re and τea. The half-life of structural litter is about one year (Bolker et al. 1998), so αe ≈ ln(2), where “ln” is the natural logarithm. This yields the estimate τea = ln(2) – re.

Modeling effects of liming on C stocks

To assess if liming effects on C fluxes inferred from our study were sufficient to account for observed effects on C stocks 19 years after liming, we simulated Eqs. 1–2 with a 0.1-year time-step for 19 years using parameter values that reflect our observed liming effects. Specifically, we multiplied Lit by the limed:control ratio of litterfall observed in our study; we multiplied (Roote + Rhizoe) by the limed:control ratio of Oe fine-root mass observed in our study; we multiplied (Roota + Rhizoa) by the limed:control ratio of Oa fine-root mass observed in our study; we multiplied (re + τea) by the limed:control ratio of Oe basal soil respiration rate observed in our study; and we multiplied (ra + τam) by the limed:control ratio of Oa basal soil respiration rate observed in our study. The key simplifying assumptions required to infer the above flux ratios from our observations are: (1) no liming effects on root lifespan (so that root-biomass ratios and root-productivity ratios are equal); (2) equal liming effects on root production and rhizosphere transfers; and (3) equal liming effects on soil respiration and transfer rates between soil layers (i.e., a single decomposition-rate effect per soil layer).

We assumed that liming effects on fluxes were constant throughout the 19-year post-liming period. Alternative scenarios (not presented), in which liming effects on fluxes were greatest immediately after liming (year 0) and then declined to the inferred year-19 effects, yielded closer correspondence between predicted and observed C stocks than depicted in Fig. 5. Thus, there may be greater consistency between observed flux and stock effects than implied by the modeling results in Fig. 5.

To specify the initial condition for simulating liming effects, we solved the control (non-limed) system (see previous section) for its equilibrium to eliminate any transient effects unrelated to liming. The stocks and fluxes reported by Fahey reflect a system close to equilibrium. To solve for the exact equilibrium, we used the control parameter values, set Eqs. 1–2 equal to zero, and solved for Ce and Ca.

Sensitivity analysis

Numerous simplifying assumptions were required to estimate the control and limed parameter values. To quantify sensitivity of model predictions to parameter values, we performed a series of simulations in which we either increased or decreased each parameter value by 10%. This sensitivity analysis indicated that model predictions were most sensitive to liming effects on fluxes (i.e., the limed:control ratios described above), and only mildly sensitive to control parameter values. Therefore, we ignored uncertainty in control parameter values and focused on uncertainty in the limed:control flux ratios.

Uncertainty analysis

We used Monte Carlo methods to quantify uncertainty in model predictions due to uncertainty in the limed:control flux ratios described above (see Modeling effects of liming on C stocks). We generated 10,000 random values for each limed:control flux ratio by generating 10,000 random limed and control values and calculating their ratios. Each limed and control value was drawn from a normal distribution whose mean and standard deviation were equal to the treatment mean and standard error, respectively, from the mixed model described in Methods. Thus, we performed 1000 simulations, each of which predicted how Oe and Oa C stocks changed over a 19-year post-liming period. All simulations had the identical initial condition, and differed only in the limed:control flux ratios. For each time step, we quantified the mean, 2.5 percentile, and 97.5 percentile of the model predictions, which are reported as the black curve and gray shaded region in Fig. 5.

We also used Monte Carlo methods to quantify uncertainty in the observed liming effect on C stocks at 19 years, and in the difference between the predicted and observed C stocks. To calculate observed Oe and Oa C stocks appropriate for comparison to our model predictions at 19 years, we multiplied the limed:control C stock ratios observed in our study by the initial (year 0) values of Ce and Ca used in the simulations (which reflect the equilibrium, non-limed state of the Fahey system). As above, we generated 10,000 random draws of Oe and Oa C stocks from their limed and control normal distributions to yield 10,000 ratios. These ratios, multiplied by the initial values of Ce and Ca, are reported as the point and error bars at year 19 in Fig. 5. We randomly paired each of these 10,000 observed C stock values with one of the 10,000 year-19 predictions to yield a sample of 10,000 differences between predicted and observed C stocks. The proportion of these differences that are positive, multiplied by two, yields a two-tailed P value (reported in Fig. 5) for rejecting the null hypothesis that the observed liming effects on fluxes are sufficient to explain the observed liming effects on stocks.

Literature Cited

Bolker, B., S. Pacala, and W. Parton. 1998. Linear analysis of soil decomposition: insights from the CENTURY model. Ecological Applications 8:425–439.

Fahey, T. J., G. L. Tierney, R. D. Fitzhugh, G. F. Wilson, and T. G. Siccama. 2005. Soil respiration and soil carbon balance in a northern hardwood forest ecosystem. Canadian Journal of Forest Research 35:244–253.


[Back to A023-096]