A Flexible Hybrid Model of Life Cycle Carbon Balance for Loblolly Pine (pinus Taeda L.) Management Systems

In this study we analyzed the effects of silvicultural treatments on carbon (C) budgets in Pinus taeda L. (loblolly pine) plantations in the southeastern United States. We developed a hybrid model that integrated a widely used growth and yield model for loblolly pine with published allometric and biometric equations to simulate in situ C pools. The model used current values of forest product conversion efficiencies and forest product decay rates to calculate ex situ C pools. Using the model to evaluate the effects of silvicultural management systems on C sequestration over a 200 year simulation period, we concluded that site productivity (site quality), which can be altered by silviculture and genetic improvement, was the major factor controlling stand C density. On low productivity sites, average net C stocks were about 35% lower than in stands with the default average site quality; in contrast, on high quality sites, C stocks were about 38% greater than average productivity stands. If woody products were incorporated into the accounting, thinning was C positive because of the larger positive effects on ex situ C storage, rather than smaller reductions on in situ C storage. The use of biological rotation age (18 years) was not suitable for C sequestration, and extended rotation ages were found to increase stand C stock density. Stands with an 18-year-rotation length had 7% lower net C density than stands with a 22-year-rotation length; stands with a 35-year-rotation length had only 4% more C than stands harvested at age 22 years. The C sequestered in woody products was an important pool of C storage, accounting for ~34% of the average net C stock. Changes in decomposition rate, associated with possible environmental changes 750 resulting from global climate change, affected C storage capacity of the forest. When decay rate was reduced to 10% or increased to 20%, the C stock in the dead pool (forest floor and coarse woody debris) was reduced about 11.8 MgC·ha −1 or increased about 13.3 MgC·ha −1 , respectively, compared to the average decay rate of 15%. The C emissions due to silvicultural and harvest activities were small (~1.6% of the gross C stock) compared to the magnitude of total stand C stock. The C model, based on empirical and biological relationships, appears appropriate for use in regional C stock assessments for loblolly pine plantation ecosystems in the southern U.S.


Introduction
Carbon dioxide (CO 2 ) emissions from fossil fuels have increased on average at a rate of 1.9% per year over the last decades [1], with atmospheric CO 2 concentration rising at approximately 1.4 parts per million by volume per year [2].As Sundquist et al. [3] pointed-out, mitigation of atmospheric CO 2 requires an approach that combines reductions in CO 2 emissions with increasing CO 2 storage.One of the most effective mechanisms for offsetting carbon (C) emissions is the fixation of atmospheric CO 2 into plant tissues [3,4].In the United States (U.S.), forests represent over 90% of the terrestrial C sink, equivalent to 12 to 16% of the U.S. greenhouse gas (GHG) emissions [2], and southern pines represent around 36% of the terrestrial C stock in the conterminous U.S. [5].The applications of sustainable forest management systems can increase C sequestration potential of managed forests [6].For example, forests in the southeast and south-central U.S. could potentially capture CO 2 equivalent to about 23% of GHG emissions of that region [7].
Managed southern pine forests have played a large role in C sequestration in the U.S. [8].The land area of southern pine plantations has increased steadily over the past half century [6], with plantations today occupying more than 13 million ha [9].Intensive management of southern pine plantations, using competition control, fertilization, and superior genotypes, can now increase productivity four-fold, compared to mid-1950s plantations [9].Rotation length is known to affect C stored in forest stands [10,11] and the use of extended rotations has been proposed as an effective way to increase C sequestration for different forest types such as scots pine (Pinus sylvestris L.) [12], Douglas-fir (Pseudotsuga menziesi (Mirb.)Franco) [13] or slash pine (Pinus elliottii Engelm.var.elliottii) [14].
In managed forests, C stocks can be divided into two major pools: in situ C in standing biomass (above and below ground) and soil organic matter, and ex situ C sequestered in products created from harvested wood [15].Sustainable forest management has the potential to greatly influence both in situ and ex situ C pools [6].
In this study, we assessed the effects of silvicultural treatments on forest C stocks per unit area in Pinus taeda L. (loblolly pine) plantations established in the southeastern U.S. Coastal Plain and Piedmont physiographic regions.Following the approach of Gonzalez-Benecke et al. [14], we developed a hybrid model that accounted for both in situ and ex situ C stocks.In situ C pool dynamics were determined using growth and yield models for loblolly pine [16,17], combined with allometric and biometric equations; ex situ C pool dynamics were determined using current values of industrial conversion efficiencies and product-specific decay rates [14].The model, which also accounted for C costs of silvicultural operations, did not include changes in soil C. We used the model to investigate net C balance, expecting that: (i) extended rotation lengths would increase time-averaged C stocks; (ii) regimes that incorporated longer rotations and thinning would maximize C accumulation in ex situ C pools; (iii) increasing stand productivity through intensive silviculture would increase net C storage; and (iv) under similar site and silvicultural treatments, loblolly pine would accumulate more C than slash pine plantations.

Models
Allometric and biometric equations were combined with growth and yield models to estimate C stocks and dynamics for loblolly pine plantations in the southeast U.S. We used published loblolly pine growth and yield equations that could be broadly applied to Lower Coastal Plain, Upper Coastal Plain and Piedmont regions [16][17][18].These equations predicted stand growth in basal area (BA, m 2 •ha −1 ), dominant height (H d , m), number of surviving trees per hectare (N ha , trees•ha −1 ), quadratic mean diameter (QMD, the diameter of the trees of mean BA, cm) and total stem volume (V, m 3 •ha −1 ) , using as inputs number of trees planted per hectare (N ha ), site index (SI, the height (m) reached by the stand's dominant trees at a reference age of 25 years), type of site preparation (no soil preparation or bedding), weed control treatment (age and method of application, broadcast or banded), fertilization treatment (age and amount of N and P applied), and thinning timing and intensity (as a removal percentage from the surviving trees before thinning).
At each age, allometric equations (Table 1) [19][20][21] were used to estimate coarse root, stem, crown and total aboveground stand biomass from QMD and N ha simulated by the growth and yield model.Biomass estimations were corrected using log-bias correction following back-transformations [22].Fine root biomass was estimated from fine root/needle mass ratios reported for stands younger than 5 years [23], and for stands 5 years and older [21].Loblolly pine projected leaf area index (LAI, the ratio of leaf surface area supported by a plant to its corresponding horizontal projection on the ground, m 2 •m −2 ) and litterfall (B L ; Mg•ha −1 •year −1 ) were estimated from the model reported by Gonzalez-Benecke et al. [24].Forest floor biomass accumulation (B F ; Mg•ha −1 ) was determined as in Gonzalez-Benecke et al. [24] as the sum of yearly litterfall inputs corrected for decay loss using an equation to estimate decay rate of the forest floor [25].This equation used site coordinates (latitude and longitude, in decimal units) as inputs to estimate decay rate and was in good agreement with decay rates reported for southern pines in the southeastern U.S. [26][27][28].Understory biomass accumulation (B U ; Mg•ha −1 ) was estimated from published equations [14] that used a 17-year chronosequence of understory and litterfall biomass [29].Harrison and Borders 1996 [16] Note: B F is foliage biomass (kg•tree −1 ); B B is branch biomass (kg•tree −1 ); B C is crown biomass (kg•tree −1 ); B S is stem plus bark biomass (kg•tree −1 ); B AG is aboveground biomass (kg•tree −1 ); B TR is tap root biomass (kg•tree −1 ); B CR is coarse root biomass (kg•tree −1 ); B FR is fine root biomass at land area basis (Mg•ha −1 ); B F is foliage biomass at land area basis (Mg•ha −1 ); LAI is mean annual projected leaf area index (m  Standing dead trees, estimated from mortality equations [19], were incorporated into the dead component of total biomass.A large fraction of the stand mortality occurred in diameter classes below the median, due to the effects of resource competition of suppressed and weakened trees [4,30].Using diameter distribution models reported for unthinned and thinned plantations [16], the diameters at the 25th percentile (cm) were determined for each age.This diameter class was assumed to better represent the observed diameter class of dying trees [30].Biomass of dying trees was computed in the same way as standing biomass, but diameters at the 25th percentile from the previous year were used instead of QMD to estimate individual tree biomass.
At the time of thinning, reductions in loblolly pine LAI were set to be proportional to reductions in BA due to thinning and, therefore, litterfall, forest floor and understory biomass were affected due to their LAI-dependence [24].At thinning and final harvest (clear-cut), logging slash (root and crown biomass plus stem residues) from harvested tress were also included into flux calculations and incorporated into the dead biomass pool.Stem residues were obtained by assuming a harvest efficiency of 88% of V [31][32][33][34].It was assumed that C storage in soil was not affected by forest management in southern pines plantations [7,[35][36][37][38], therefore the model did not include changes in soil C. Carbon mass (MgC•ha −1 ) was calculated by using an average C content of 50% for loblolly pine and understory biomass components [6,39,40].

Model Validation
Validation of model results was carried out by comparing model outputs against data reported in the peer-reviewed literature for: (i) total above-plus belowground C accumulation in loblolly pine biomass (T C , MgC•ha −1 ); (ii) aboveground C accumulation in live loblolly pine biomass (AG C , MgC•ha −1 ); (iii) belowground C accumulation in live loblolly pine biomass (BG C , MgC•ha −1 ); and (iv) net ecosystem production (NEP; MgC•ha −1 •year −1 ).Results reported in 15 publications including both above-and belowground loblolly pine biomass, for stands between 3 and 48 years old, were considered in our approach.On each validation study, the authors estimated biomass directly using inventory and local biomass functions.For each literature value, the reported initial stand characteristics such as planting density, SI and management activities (site preparation, weed control, fertilization and thinning) were used as model inputs.Loblolly pine's T C , AG C and BG C were validated in two ways, depending on the availability of information: (i) whole-model validation, using reports that included planting density, SI and management activities [25,[41][42][43][44][45][46][47][48]; or (ii) allometry validation, using reports where SI was not reported [27,[49][50][51][52][53].The former enabled us to run the model from planting year to age when biomass was determined.For allometry validation from papers which did not report initial stand characteristics, initial model inputs were adjusted to achieve the final stand characteristics of each report, such as H d , DBH, BA, N ha .The latter approach can be considered a validation of the biomass equations used by the model.Even though two of the studies used for validation [21,23] utilized some of the equations included in the model (see Table 1), we contend that the independent estimation of N ha , DBH and all other C pool components validate the inclusion of those cases.The values of NEP reported by Oren et al. [54] for 16 to 21-year old stands were compared with our model estimates.A summary table of all data used for C accumulation validations is reported in Appendix 1. Locations of validation sites are also shown in Figure 1.Validation of litterfall and forest floor outputs were developed/tested in related research by Gonzalez-Benecke, et al. 2011 [24].Validation simulations were performed with the growth and yield functions appropriate for the physiographic region of each validation study.Sites with black-filled symbols were used for the allometric validation, sites with white symbols were used for the whole model validation, and sites with black and white circles were used for the NEP validation.

Ex Situ Wood Products Pools
Depending on stem DBH and merchantable diameter, harvested roundwood (from thinnings or clear-cuts) was assigned to three main product classes: sawtimber (ST), chip-and-saw (CNS) and pulpwood (PW), using the model proposed by Harrison and Borders [16].Merchantable outside-bark stem volume was calculated at each stand age and product volume was transformed to biomass (Mg• ha −1 ) by multiplying the volume by the average specific gravity (SG) at a given age from equations derived from the model [7].The exponential model that correlates SG and age is shown in Table 1.A C content of 50% was used to calculate C mass of each product type.

Silvicultural Management Scenarios
The effects of silvicultural treatments and rotation lengths on C sequestration were analyzed by simulating C flux under four different scenarios for standard conditions of loblolly pine plantations established in the southeastern US lower Coastal Plain.Initial variables were set as follow: base SI = 22 m, N ha = 1500 trees•ha −1 , bedding, weed control at planting and at age 1, and NP fertilization at age 5 (135 Kg•ha −1 N + 28 Kg•ha −1 P), 10 (225 Kg•ha −1 N + 28 Kg•ha −1 P) and 15 years (225 Kg•ha −1 N + 28 Kg•ha −1 P).After fertilization and weed control treatments, the model estimated that the effective SI was increased to 24.3 m.Initial C accumulated from the previous rotation in coarse root debris (~21.6 MgC•ha −1 ; [41,42,49]) and the forest floor and aboveground coarse woody debris (~33.83MgC•ha −1 ; [61]) was assumed to be 55.4 MgC•ha −1 .This value was similar to the amount of 60.0 MgC•ha −1 reported for a slash pine plantation following a clearcut harvest [62].Thinning intensity removal was set at 33% of the living trees.Details of the four management scenarios used in the modeling effort are shown in Table 3. Silvicultural-scenario simulations were performed with the growth and yield functions for the lower Coastal Plain physiographic region.

Carbon Emissions of Transportation and Silvicultural Activities
Emissions of C by silvicultural activities were determined from Markewitz [63] and White et al. [64].Carbon release in transportation of raw material from the forest to the mill was estimated according to White et al. [64], assuming an average haul distance of 100 km from the forest to the mill, load per logging truck of 24 m 3 , and fuel economy of the diesel logging truck of 2.6 km l −1 .Details of C emissions are presented in Table 4. Average for 24 m 3 load capacity 0.0026 (2)   Note: (1) Carbon use in fertilization includes production, packing, transportation and application [63]; (2) Carbon use for transportation is expressed in MgC used per m 3 transported [64].

Sensitivity Analyses
Sensitivity analyses were carried out to determine the effects of changes in key parameter estimates on total C balance.The effects of site quality were assessed by evaluating the model under contrasting SI values of 15 and 30 m, which corresponded to the full range of site qualities observed in loblolly pine plantations in the southeastern US.After fertilization and weed control treatments, the model estimated that the effective SI increased from 15 to 18.2 m, and from 30 to 32.3 m.Initial stand density effects were evaluated by running the model under contrasting planting densities of 750 and 2250 trees•ha −1 .Forest floor accumulation and total net C stock were evaluated by using decomposition rates of 10 and 20%; values beyond the extremes of 12 and 17% reported for this region [26,27,65].In addition to the silvicultural management scenarios tested, rotation length effects were assessed by evaluating the model under the PP scenario for 18 and 35 years.Average product life span was evaluated by changing the proportion of products in different life span classes.In the case of ST and CNS, the proportion of products in the long life classes (50 years) were changed by 25% (step up and down), distributing the residual proportion in equal parts to the rest of the life span classes.Sensitivity analyses to industrial conversion efficiencies were not considered due to their low impact on ex situ C stocks [14].

Thinning Effects on Loblolly Pine C Stocks
The effects of thinning (as a percentage of living trees removed) were assessed by evaluating the model under different combinations of thinning age (8, 12 and 16 years) and intensity (20, 40 and 60% of living trees removal).Stand and management characteristics were set as those described previously for ST1 (i.e., base SI = 22 m, N = 1500 trees•ha −1 , bedding, weed control at planting and at age 1, and NP fertilization at age 5-10-15 years).

Comparison Between Loblolly and Slash Pine C Stocks
We compared estimates of loblolly pine C stocks with those of slash pine using a previously reported model [14] that was updated with relationships used to estimate LAI, litterfall and forest floor accumulation [24].Initial parameter estimates such as site index, planting density, site preparation, weed control, fertilization and rotation length were set equal for both species.An analysis for unthinned stands with rotations length of 22 and 35 years was performed.For both species, two NP fertilizations were set at ages 5 and at age 10 years (120 kg•ha −1 diamonium phosphate + 210 kg•ha −1 urea, and 120 kg•ha −1 diamonium phosphate + 430 kg•ha −1 urea, respectively).

Model
Net C stock was defined as: Net C stock = Total C in situ (C stored in living loblolly pine trees + understory + forest floor + coarse woody debris + standing dead trees) + Total C ex situ (C stored in wood products ST + CNS + PW) − Total C cost (silvicultural activities, including transportation of supplies).For all simulations, estimates of average C stock were reported as the average of all yearly values from the first ~200 years of management, stopping the simulation at the end of the rotation closest to the 200 year endpoint (not stopping the simulations midway into a rotation).For the scenarios with rotation length of 18, 22 and 35 years, the number of rotations simulated was 11, 9 and 6, respectively.This simulation length was chosen to be sufficiently long to approach steady state values for ex situ pools, while remaining within plausible bounds for consideration of future forest management scenarios.While it is likely that "true" steady state would be achieved for simulations of 500-1000 years, these time periods are far beyond the time horizon over which forest management planning occurs.

Statistical Analysis
Four measures of accuracy were used to evaluate the "goodness of fit" between observed and predicted (simulated) values for each variable from the dataset obtained in the model validation [66][67][68][69][70]: (i) Mean absolute error (MAE); (ii) Root mean square error (RMSE); (iii) Mean bias error (Bias); and (iv) Pearson product-moment correlation coefficient (R).

Model Validation
There was no statistical difference between average NEP modeled (NEP EST ) and measured NEP (NEP EC ) from ages 15 to 21 years (P = 0.77), averaging 4.46 and 4.63 MgC•ha −1 •year −1 for NEP EST and NEP EC , respectively (Table 5; Figure 2).Predicted values were in the mid-range of variation of observed values.Estimations of C stock in aboveground (AG C ), belowground (BG C ) and total (T C ) live pine were well correlated with values reported for stands of different ages and productivity (R 2 > 0.94; Figure 3a,  b and c).In all cases, the intercept and slope of these relationships were not different from zero and one, respectively (P > 0.27).5).For the whole-model validation of C stock estimations, MAE and RMSE ranged between 15.3 to 16.7%, and 12.8 to 14.3% of the observed C stock values, respectively.On the other hand, MAE and RMSE for the validation of the allometric equations for biomass estimation (that was carried out with sites were SI was not reported), ranged between 11.1 to 23.2%, and 8.6 to 21.4% of the observed C stock values, respectively.In all cases, BG C estimations presented the larger differences between the observed and predicted values.The Bias of the validation of allometry ranged between 1.5% under-estimations for AG C to 6.5% over-estimations for BG C , whereas Bias for the whole-model validation ranged between 0.1% under-estimations for AG C to 8.5% over-estimations for BG C , with no clear tendency to under-or over-estimate the results.Estimated and observed values of C stock were highly correlated, with R values greater than 0.94.In the case of NEP, the model was less accurate than the C stock estimates, even though the mean was not different between the observed and predicted values.MAE and RMSE were 28.1 and 30.5 % of the observed values, respectively.The Bias showed a 3.7% underestimation.

Silvicultural Management Effects on C Sequestration
One-thinned or two-thinned stands with longer rotation ages (35 years) stored only 4% or 1% more C, respectively, than unthinned stands harvested at age 22 years.Under conditions used in the simulations, average net C stock, which corresponded to total C in situ (loblolly pine + understory + forest floor + coarse woody debris + standing dead trees) plus total C ex situ (C in woody products ST + CNS + PW) minus total C cost of silvicultural activities (including transport), averaged 176,179,183 and 178 MgC•ha −1 for PP, ST1, ST2 and ST3, respectively, for a 200 year simulation period (Figure 4).In situ C stock accounted for between 62% and 71% of the gross C sequestration (not including silvicultural C costs) across silvicultural regimes.The magnitude of emissions associated with silvicultural activities (including transportation) was between 1.7% and 1.8% of the gross C stock.The relative impact on C sequestration for the different woody products depended on the silvicultural regimes; ST accounted for 40%, 60%, 75% and 74% of total C ex situ for PP, ST1, ST2 and ST3, respectively.In contrast, CNS followed an opposite trend, accounting for 56%, 37%, 23% and 23% for the PP, ST1, ST2 and ST3 silvicultural regimes, respectively.Across different silvicultural management systems, the forest floor + dead trees + coarse woody debris (FF D ) components averaged 42 MgC•ha −1 (between 31% to 39% of total in situ C stock) and the understory averaged 1.0 MgC•ha −1 (less than 1% of total in situ C stock).From a total silviculture C cost perspective, fertilization accounted for 0.92 MgC•ha −1 (three fertilizations), representing between 28% and 30% of the total silvicultural C emissions.Harvest and transportation of woody products accounted for between 1.42 and 1.84 MgC•ha −1 , representing about 47 to 55 % of the total silvicultural C emissions.
In general, after ~200 years, C flux in the woody products converged to stable values, reaching quasi-equilibrium minimum and maximum values.At their respective rotation ages, in situ C stocks were 166, 161, 169 and 158 MgC•ha −1 for the PP, ST1, ST2 and ST3 scenarios, respectively (Figure 5).From that total, the C stock in living loblolly pine and the understory was 128, 126, 131 and 124 MgC•ha −1 for the same silvicultural regimes, respectively (data not shown).Total woody products C stock increased each rotation from 80, 77, 95 and 90 MgC•ha −1 during the first rotation, up to 142, 150, 153 and 150 MgC•ha −1 at the end of the rotation closest to the 200 year endpoint, for the PP, ST1, ST2 and ST3 scenarios, respectively (Figure 5).This result implies that in the long term, for extended rotations, C storage fluxes in woody products were similar than the amount of C stored in living loblolly pine.Differences in tree size (diameter and height) and number of trees remaining due to different thinning and rotation age scenarios created different woody products pools that had different life spans.While PW represented 4.2%, 3.4%, 2.3% and 2.6% of the total C extracted in products at harvest, ST accounted for 40%, 60%, 75% and 74% of that C for the PP, ST1, ST2 and ST3 scenarios, respectively.In general, C stored in products derived from PW (e.g., paper, packing material, office supplies, etc.) had a negligible effect on net C sequestration; between harvest events (thinning or clear cutting) the amount of C stored in products derived from pulpwood diminished towards zero, while C stored in CNS and ST increased between harvests (Figure 6).

Sensitivity Analysis
With all other parameters held constant, site quality (or potential productivity) reflected in SI of the stand, was the major factor controlling C sequestration (Table 6).For example, on low productivity sites (e.g., base SI = 15 m), average net C stocks were about 28% lower than with the default site quality (base SI = 22 m).In contrast, for high quality sites (e.g., base SI = 30 m), C stocks across silvicultural regimes averaged about 38% greater than base SI = 22 (Table 6).When SI was set equal to 15 m, ex situ C stocks were largely reduced for sawtimber-oriented silvicultural regimes.In situ C stocks were reduced 21.6, 20, 25.6 and 24 MgC•ha −1 , and ex situ C stocks were reduced 28, 30, 26.4 and 26.8 MgC•ha −1 , for PP, ST1, ST2 and ST3, respectively.By comparison, for base SI = 30 m, in situ C stocks increased 25, 23, 30 and 28 MgC•ha −1 , and ex situ C stock augmented 44.6, 44.7, 39.3 and 40 MgC•ha −1 , for the PP, ST1, ST2 and ST3 management regimes, respectively.The silvicultural treatments included in our analysis, that improved SI by 3.2, 2.3 and 2.3 m, for stands with base Si of 15, 22 and 30 m, respectively, produced an increase in net C stock of 12.5, 12.4 and 12.8 MgC•ha −1 , respectively, across management regimes (Table 6).The effects of planting density were largely reflected in in situ rather than ex situ C pools.Reducing the initial planting density decreased net C stocks up to 10%, while increasing the planting density enhanced net C stocks by about 6% among silvicultural regimes.By lowering planting density from 1500 trees•ha −1 to 750 trees•ha −1 , the average net C stocks decreased 16.6,17.3,19and 18.9 MgC•ha −1 for the PP, ST1, ST2 and ST3 management systems, respectively.This reduction was explained principally by a decrease in the in situ C stocks of 11.7, 10.5, 10.7 and 10.3 MgC•ha −1 for the same silvicultural regimes.The effects on ex situ C stocks were substantially smaller, producing gains of 2.4 MgC•ha −1 for the PP, but reductions of 0.1, 2.2 and 3.1 MgC•ha −1 for the ST1, ST2 and ST3 management systems, respectively.By increasing the initial planting density from 1500 trees•ha −1 to 2250 trees•ha −1 , average net C stocks were increased 10.7, 9.5, 11.2 and 11 MgC•ha −1 compared to the default initial planting density for the PP, ST1, ST2 and ST3 treatments, respectively.For the high initial planting density, in situ C stocks were increased 8.1, 7.2, 6.8 and 6.6 MgC•ha −1 for the PP, ST1, ST2 and ST3 silvicultural management systems, respectively.The C storage in woody products was reduced by 2.7, 2.4 MgC•ha −1 for the PP and ST1 management regimes, but was increased by 0.2 and 0.7 MgC•ha −1 for the ST2 and ST3 scenarios, respectively.
Shortening the rotation length to 18 years under the PP scenario (unthinned), which corresponded to the biological rotation age for that silvicultural scenario (i.e., the time when mean annual increment equals periodic annual increment), reduced the average net C stock to 164.2 MgC•ha −1 .This 7% decrease, when compared with the default 22-year-rotation, was caused mainly by reductions in the in situ rather than ex situ C stocks (10.8 and 1.4 MgC•ha −1 reductions, respectively).When the rotation age was extended to 35 years, average net C stocks increased to 183.2 MgC•ha −1 .
When the C decay rate was reduced to 10% or increased to 20%, average net C stocks increased or decreased between 7.5 and 6.6%, respectively, across silvicultural regimes (Table 6).Assuming a decay rate of 10%, C storage in the FF D increased to about 13.3 MgC•ha −1 , across silvicultural regimes, compared to a decay rate of 15%.When the decay rate was increased up to 20%, the C stocks in the FF D were reduced to about 11.8 MgC•ha −1 .
Variations in average life span of ST and CNS woody products affected net C stock.On the other hand, paper product life span had a smaller effect on net C storage (Table 6).When the average life span of ST was increased by increasing the product proportion in the long-lived class (half-life 50 years, Table 2) from 50% to 75%, average net C stocks were increased 7.3, 12.6, 11.9 and 12.6 MgC•ha −1 for the PP, ST1, ST2 and ST3 management regimes, respectively; this represented an increment of between 4.1 and 7.1% (Table 6).On the other hand, when the ST proportion of long life span products was stepped-down to 25%, reductions in C increments of the same magnitude were observed.The impact of the CNS half-life on net C stocks was more important in the PP than the sawtimber-oriented regimes, changing average net C stocks by 12.2, 9.2, 4.3 and 4.6 MgC•ha −1 for the PP, ST1, ST2 and ST3 management scenarios, respectively (Table 6).When all PW products were set to have a life span of one year, average net C stocks was reduced between 0.42 and 0.86 MgC•ha −1 .Conversely, when the PW products life span was increased by assuming that 67% of PW products would have a life span of four years, the average net C stock was increased between 0.45 and 0.89 MgC•ha −1 .

Thinning Effects on Loblolly Pine C Sequestration
In general, for a base SI = 22 m stand with 1500 trees•ha −1 and managed under a 22-year rotation, thinning had a positive on net C stock.The negative effect of thinning on the in situ C stock was largely counteracted by increasing the ex situ C stock, producing an increase in average net C stock (Figure 7).
For any given thinning intensity, there was a small effect due to the age of thinning on the in situ C stock (Figure 7a).For example, using a thinning intensity of 20% at ages 8, 12 and 16 years reduced the in situ C stock by 3.1, 3.4 and 3.1 MgC•ha −1 , respectively (Figure 7a).Increments in thinning intensity produced a quasi-constant decline in in situ C stock, independent of thinning age.For instance, in situ C stocks were reduced 3.2, 7.9 and 13 MgC•ha −1 , when thinning intensities were set at 20, 40, and 60% removal of living trees, respectively (Figure 7b).Ex situ C storage had an opposite response to thinning; the more intensive the thinning regime, the more gain in woody products C storage.Ex situ C storage was slightly increased by thinning age for any given thinning intensity.For example, when thinning age was set at 8,12 and 16 years, with a thinning intensity of 20%, the ex situ C stock was increased by 5.6, 5.2 and 4.7 MgC•ha −1 , respectively; when thinning intensity was set to 60%, ex situ C stocks were increased 19.4,19.9 and 17.7 MgC•ha −1 , respectively (Figure 7c).Increments in thinning intensity produced a quasi-constant increment in ex situ C stock independent of thinning age.Ex situ C stocks were increased 5.2, 11.5 and 19 MgC•ha −1 , when thinning intensities were set at 20, 40 and 60% removal of living trees, respectively (Figure 7d).Due to the larger positive effects on ex situ C storage, rather than smaller negative effects on in situ C storage, net C stocks were slightly increased by thinning (Figure 7e and f)., d, f) on In Situ, Ex Situ and average net C stocks for loblolly pine plantations over a ~200 year simulation period.

Comparison in C Sequestration between Loblolly and Slash Pine Stands
Under the default parameters used for simulations in unthinned loblolly pine stands harvested at age 22 years, loblolly and slash pine net C stocks were generally similar, albeit on an absolute basis loblolly pine accumulated 1.3 Mg •ha −1 more C than slash pine (Figure 8).This 1% difference in favor of loblolly pine was explained mainly by increases in stored C in the living tree biomass (2.8 MgC•ha −1 ) and FF D (3.8 MgC•ha −1 ) pools.It should be noted, however, that a 4.6 MgC•ha −1 decrease in the loblolly pine woody products pool partially counteracted that gain.Understory biomass was 0.7 MgC•ha −1 larger in slash pine than in loblolly pine stands.When the rotation length was increased to 35 years without thinning, slash pine C stocks exceeded loblolly pine, averaging 11.9 MgC•ha −1 more net C stock.This difference was explained mainly through increases in net C stored in the living Stand Age (years)

Discussion
Loblolly pine plays an important role in mitigation of CO 2 emissions due to its high productivity and extensive planting throughout the southeastern U.S. Accurate determinations of C stocks and the understanding of factors controlling C dynamics in loblolly pine plantations are essential for C offset projects and the development of sustainable management systems.To validate the model predictions, we utilized 15 peer-reviewed publications that reported estimates of both above and belowground biomass accumulations in loblolly pine stands.Several other publications were also found that reported only aboveground biomass estimates ( [71][72][73][74][75], among many others); however, we decided to validate the model using only those published studies that included both above-and below-ground estimates measured at the same time for a given site.Further, the model was tested on stands of contrasting site quality (SI ranges between 15.5 to 32 m) that ranged in age from 3 to 48 years, across different physiographic regions that cover the planting range of the species (Figure 1; Appendix 1).The strong agreement between observed and predicted values supported the robustness of the model and its usefulness for assessing the effects of forest management activities on C sequestration (in situ, and ex situ) for the most important commercial tree species in the southeastern US.
The dominant factor controlling C sequestration in loblolly pine plantations was site quality as it interacted with different management scenarios.Similar responses have been reported for slash pine [14], Pinus radiata D. Don and Pinus pinaster Ait.[76].Across the silvicultural regimes evaluated there were similar changes in in situ C stocks, but ex situ C storage was largely affected by site quality changes and longer rotation lengths oriented toward sawtimber production.Specifically, our results indicated that: (i) long rotation oriented-sawtimber silvicultural regimes were not as effective for C storage on low quality sites, but they were more effective on high quality sites; and (ii) ex situ carbon pools were important considerations when evaluating the C offset potential of different management systems.Over the past 50 years, the productivity of planted loblolly pine stands has tripled [9,77], implying an important effect of modern forest management systems (interaction between genetic improvement, seedling culture, and nutrient and competition management) not only on volume yield, but also on total C storage.The rate of development and implementation of technology that increases production rates will determine the contribution of future plantation C storage.The effects of site quality changes in C storage were exacerbated by the effects of ex situ C sequestration because increases in site quality not only increased total standing volume (or stand biomass), but also increased the proportion of trees making valuable product grades that had a long life span.

Species -Rotation Length
Increasing initial planting density in the range tested in this study had a positive effect on net C storage, and the effects of planting density on C storage were most apparent in the in situ C pool, affecting both living tree biomass and FF D biomass accumulation.Even though raising the planting density increased the proportion of fixed C used in stem production in loblolly pine [78], this effect was not reflected in the ex situ C pool.As planting densities increased, there was a tendency to decrease sawtimber products yields, affecting the average ex situ C pools; however, the increase in forest floor, coarse woody debris and total living tree C storage largely counteracted that negative effect.Morton [79] concluded that utilizing higher planting densities maximized C sequestration and with a market value for C credits, land expectation values were maximized by utilizing higher planting densities as well.
Increasing the rotation length increased C stock in loblolly pine stands.We estimated similar net increments of 6.9 MgC•ha −1 when rotation lengths were increased from 22 to 35 years for unthinned and thinned stands.Most of this increment was accounted for in the in situ C storage.In contrast, when rotation lengths were shortened to 18 years, net C stock was reduced.Other reports for different conifer species [12][13][14] indicated similar effects of rotation length on C storage i.e., extended rotations increased C sequestration in conifer forest plantations.Longer harvesting cycles represent one of the major management strategies used to increase forest C density [80].Nevertheless, the inclusion of biomass harvest for fossil fuel offset might change our conclusions, especially when shorter rotations include provisions for a technology bump at the end of each rotation.Further research is needed in this area, and this model is a tool to address these types of questions.
In our study, under a low decomposition rate of 10%, the C stock in the FF D increased about 13.3 MgC•ha −1 ; under a higher decomposition rate of 20%, C stocks in the FF D decreased about 11.8 MgC•ha −1 .As changes in soil temperature and moisture regimes can affect litter decomposition rates [81], it is probable that future environmental changes would be expected to alter litter decomposition rates in loblolly pine stands, and thereby affect the C storage capacity of the forest.
When stand C density was compared between loblolly and slash pine stands under similar levels of site quality (base SI = 20 m) and silvicultural inputs over a 22-year rotation length, living pine C stocks of loblolly and slash pine were generally similar (loblolly pine was 1% greater than slash pine).Under longer rotations (i.e., 35 years), the average net C stocks of both species remained similar; however, slash pine was about 6% greater than loblolly pine.For thinned stands with similar SI as our simulations, loblolly pine at age 31 years had about 8% more standing volume than slash pine [82].For unthinned stands at age 28 years, the differences between loblolly and slash pine standing volumes were about 1% [83].The first study to our knowledge that compared C stocks between southern pine species in mature stands was from Vogel et al. [84].The authors reported that at age 26 years, living tree C stocks of loblolly pine were larger than slash pine when nutritional limitations were eliminated through fertilizer additions.Under nonfertilized conditions, even with sustained elimination of understory vegetation, living tree C stocks were larger for slash pine than for loblolly pine.As nutritional demands and the responses to fertilization for loblolly pine tend to be larger than slash pine [77], differences in nutrient requirements and nutrient use efficiency between the two species should be taken in account when developing sustainable and ecological forestry regimes.In our analysis, the fertilization regime included two applications, which may not be sufficient to support the demands of loblolly pine, especially under longer rotations scenarios.
For the ~200 year simulation period, thinning decreased in situ C stocks, but increased ex situ C stocks more, resulting in slight increases in net C stocks.Most of the studies that have addressed the impacts of thinning on C budgets in pine ecosystems have reported the responses only on living pine biomass [76,[85][86][87][88]; few studies have reported the impacts of thinning on total in situ C [89][90][91][92].All the previously cited studies concluded that there was a reduction in pine stand or total in situ C after thinning.The negative effects of thinning were also reported for NEP fluxes for Pinus ponderosa Dougl.forests [90,91].In both studies, the authors reported a decrease in NEP for thinned stands.To our knowledge, the current study represents one of the few that incorporated the ex situ C pool into the analysis of thinning effects on carbon sequestration of forest plantations.Garcia-Gonzalo [93], in a similar analysis that included ex situ C pools for mixed coniferous stands in Finland, reported a net reduction between 25 and 33 MgC•ha −1 in trees and a net increase between 30 and 45 MgC•ha −1 in harvested timber.Even though the wood extracted in thinning was primarily pulpwood, which impacted ex situ C sequestration, increased growth of residual trees due to thinning promoted the production of larger tree size classes at final harvest.These long-lived products increased the ex situ C pool, overcompensating for the reduction in in situ C associated with thinning.When ex situ pools were considered, the possible economic benefits of thinning were not in opposition to maintaining or increasing net C stock.This, and other evidence, suggests that there is no simple inverse relationship between the amount of timber harvested from a forest and the amount of C stored [94].
It should be noted that estimates of the effects of scenarios on net C stocks are somewhat sensitive to the length of simulation.This is because ex situ C stocks take centuries to reach equilibrium levels (200-400 years for the scenarios we tested, data not shown).Therefore, as longer and longer simulation periods are compared, the effects of differences in in situ C stocks on net C stocks tend to be dampened by converging magnitudes of ex situ C stocks.This makes it important to carefully consider the inference space desired for a particular simulation.For relatively short-term management decisions, it may be appropriate to run simulations for only one or two rotations.For decision-making associated with C offset accounting contracts, 50-150 years may be appropriate.As mentioned previously, we chose 200-year simulations as a length of time long enough to allow development of near-equilibrium ex situ stock levels, but not so long that the time period moved out of the realm of realistic forest management planning.
Under the Intergovernmental Panel on Climate Change [95], reporting of C stock in wood products is not mandatory, but the enhancement of that C pool could provide important GHG emission offsets.Results from this study suggested that extending the half-life of PW products had only a marginal effect on C stock.Similar results were reported in other studies [14,56,58].The ex situ C pools could be influenced by both the final utilization of particular products, and also by substituting wood for more C intensive materials.If waste wood and forest biomass residues were used as substitutes for fossil fuels [96], or if long lasting wood products took the take the place of more C intensive materials like concrete or steel [97], then the mitigation impacts of ex situ C stocks could be even larger.

Conclusions
In this article we provide a model for quantifying C sequestration for loblolly pine plantations under varying management conditions in the southeastern U.S. Coastal Plain and Piedmont regions.The model performed accurately when it was tested against reported C measurements over a wide range of stand ages and site qualities.Using the model to evaluate the effects of silvicultural management systems on C sequestration over a 200 year simulation period, we conclude that: (i) site quality (productivity), that can be altered by silviculture and genetic improvement, represents the major factor controlling average net C stock in loblolly pine plantations; (ii) if woody products were incorporated into the accounting, thinning tends to be C positive because of the larger positive effects on ex situ C storage; (iii) shorter rotations (biological rotation age) were not as effective for C sequestration as extended rotations that increased average net C stock; (iv) C sequestered in woody products accounted for ~34% of the net C stock; (v) changes in decomposition rates associated with possible global climate change could affect C storage capacity of the forest; and (vi) emissions due to silvicultural and harvest activities were small compared to the magnitude of the total stand C stock.

Figure 1 .
Figure 1.Location of data validation sources in the southeastern U.S. Shaded area indicates the natural range of loblolly pine.Sites with black-filled symbols were used for the allometric validation, sites with white symbols were used for the whole model validation, and sites with black and white circles were used for the NEP validation.
−1 •year −1 ); AG C is loblolly pine aboveground C accumulation (MgC•ha −1 ); BG C is loblolly pine belowground C accumulation (MgC•ha −1 ); T C is loblolly pine total C accumulation (MgC•ha −1 );is the mean observed value; is the mean predicted value; n is the number of observations; MAE is the mean absolute error; RMSE is the root of mean square error; Bias is the absolute bias estimator; R is the Pearson's correlation coefficient.

Figure 2 .
Figure 2. Time series of simulated versus observed NEP for a loblolly pine plantation located in the transition between the upper Coastal Plain and Piedmont regions in North Carolina.

Figure 3 .
Figure 3. Validation of allometric equations (black filled symbols) and whole model (open symbols) outputs: (a) Simulated versus measured total C stock (T C ) in living loblolly pine; (b) Simulated versus measured aboveground (AG C ) living loblolly pine C stock; (c) Simulated versus measured belowground (BG C ) living loblolly pine C stock.

Figure 4 . 8
Figure 4. Average carbon stock for loblolly pine plantations for a ~200 years simulation period under different silvicultural scenarios.

Figure 5 .
Figure 5. Annual carbon stocks for loblolly pine plantations under different silvicultural scenarios for a 200-year simulation period.

Figure 6 .
Figure 6.Annual carbon stocks in woody products for loblolly pine plantations under different silvicultural scenarios for a 200-year simulation period.

Figure 7 .
Figure 7. Effects of stand age at thinning (a, c, e) and thinning intensity (b, d, f) on In Situ, Ex Situ and average net C stocks for loblolly pine plantations over a ~200 year simulation period.
•ha −1 ) and ex situ biomass (5.3 MgC•ha −1 ) pools.Understory C accumulation was 0.9 MgC•ha −1 larger in slash pine than in loblolly pine stands, and the FF D was 2.1 MgC•ha −1 larger in loblolly than in slash pine stands.

Figure 8 .
Figure 8.Average carbon stock for unthinned loblolly (LOB) and slash (SLA) pine plantations for a ~200 year simulation period under two different rotation lengths (22 and 35 years).

Table 1 .
Allometric equations used to estimate loblolly pine aboveground (needles, branches and stem plus bark) and belowground (tap, coarse and fine roots), forest floor and understory biomass. 2

Table 2 .
Proportions of wood products by life-span category.
Note: ST is Sawtimber; CNS is Chip and Saw; PW is Pulpwood; Values in parenthesis indicate average life span for class (years).

Table 3 .
Stand management scenarios.Th 10 is commercial thinning of 33% of living trees at age 10 yr; Th 17 is commercial thinning of 33% of living trees at age 17 yr.All scenarios included bedding, weed control at planting and at age 1, and NP fertilization at age 5, 10 and 15 yr.

Table 4 .
Carbon used in silvicultural activities and product transportation to mill gate.

Table 5 .
Summary of model evaluation statistics.

Table 6 .
Sensitivity analysis of average net carbon stock for selected parameters under different silvicultural scenarios over a ~200 year simulation period.
Note: Average net carbon stock (MgC•ha −1 ) is the average of a ~200 year simulation period and Δ% is the percentage deviation from default parameter values used (Effective site index = 24.3m; Planting density = 1500 trees•ha −1 ; Rotation length = 22 years; Decay rate = 15%; ST in long life class = 50%; CNS in long life class = 25%; PW in medium-short life class = 33%).