Significant Increase in Ecosystem C Can Be Achieved with Sustainable Forest Management in Subtropical Plantation Forests

Subtropical planted forests are rapidly expanding. They are traditionally managed for intensive, short-term goals that often lead to long-term yield decline and reduced carbon sequestration capacity. Here we show how it is possible to increase and sustain carbon stored in subtropical forest plantations if management is switched towards more sustainable forestry. We first conducted a literature review to explore possible management factors that contribute to the potentials in ecosystem C in tropical and subtropical plantations. We found that broadleaves plantations have significantly higher ecosystem C than conifer plantations. In addition, ecosystem C increases with plantation age, and reaches a peak with intermediate stand densities of 1500–2500 trees ha−1. We then used the FORECAST model to simulate the regional implications of switching from traditional to sustainable management regimes, using Chinese fir (Cunninghamia lanceolata) plantations in subtropical China as a study case. We randomly simulated 200 traditional short-rotation pure stands and 200 sustainably-managed mixed Chinese fir – Phoebe bournei plantations, for 120 years. Our results showed that mixed, sustainably-managed plantations have on average 67.5% more ecosystem C than traditional pure conifer plantations. If all pure plantations were gradually transformed into mixed plantations during the next 10 years, carbon stocks could rise in 2050 by 260.22 TgC in east-central China. Assuming similar differences for temperate and boreal plantations, if sustainable forestry practices were applied to all new forest plantation types in China, stored carbon could increase by 1,482.80 TgC in 2050. Such an increase would be equivalent to a yearly sequestration rate of 40.08 TgC yr−1, offsetting 1.9% of China’s annual emissions in 2010. More importantly, this C increase can be sustained in the long term through the maintenance of higher amounts of soil organic carbon and the production of timber products with longer life spans.


FORECAST model description
The ecosystem management simulation model FORECAST (Kimmins and others 1999) has been used as a long-term management evaluation tool in several types of forest ecosystem (e.g., Morris and others 1997;others 2000, 2003;Seely and others 2002;Welham and others 2002), including tropical and sub-tropical plantations (Bi and others 2007;Blanco and González 2010).
Evaluation exercises have demonstrated the reliability of this model (Blanco and others 2007;Seely and others 2008;Blanco and González 2010). FORECAST was specifically designed to examine the impacts of different management strategies or natural disturbance regimes on long-term site productivity. The projection of stand growth and ecosystem dynamics is based on a representation of the rates of key ecological processes regulating the availability of, and competition for, light and nutrient resources ( Figure S1). The rates of these processes are calculated from a combination of historical bioassay data (biomass accumulation in component pools, stand density, etc.) and measures of certain ecosystem variables (e.g. decomposition rates, photosynthetic saturation curves) by relating 'biologically active' biomass components (foliage and small roots) to calculations of nutrient uptake, the capture of light energy, and net primary production. Using this 'internal calibration' or hybrid approach, the model generates a suite of growth properties for each tree and plant species to be represented. These growth properties are subsequently used to model growth as a function of resource availability and competition (Kimmins and others 1999). They include (but are not limited to): 1) Photosynthetic efficiency per unit foliage biomass based on relationships between foliage biomass, simulated self-shading, and net primary productivity after accounting for litterfall and mortality; 2) Nutrient uptake requirements based on rates of biomass accumulation and literature-or field-based measures of nutrient concentrations in different biomass components on different site qualities; 3) Light-related measures of tree and branch mortality derived from stand density input data in combination with simulated light profiles. Light levels at which foliage and tree mortality occur are estimated for each species.
Soil fertility in FORECAST is represented based on empirical input data describing decomposition (mass loss) rates and changes in chemistry as decomposition proceeds. These data allow for the calculation of nutrient release from litter and humus ( Figure S1). Nutrient uptake demands of different species on sites of different fertility are based on observed biomass accumulation rates and tissue nutrient concentrations on these sites, allowing for internal cycling of nutrients. The calculated uptake demand by the observed growth rates on sites of different productivity permits a definition of nutritional site quality. This assumes that moisture is not the major limiting factor, or that, if it is limiting, it acts dominantly through soil processes that determine nutrient availability. In the humid climates that characterise the Chinese fir region this assumption is felt to be reasonable. Figure S1. A schematic representation of the ecosystem compartments and transfer pathways represented in FORECAST (adapted from Kimmins et al. 1999).
Carbon allocation in response to soil fertility and tree/plant nutrition is based on empirical biomass ratios and biomass turnover rates (e.g., number of years of leaf retention for evergreens) for sites of different fertility (e.g., different site nutritional quality), and on literature or locally-obtained values for variation in fine root turnover along fertility gradients. FORECAST performs many of its calculations at the stand level but includes a submodel that disaggregates stand-level productivity into the growth of individual stems with user-inputted information on stem size distributions at different stand ages. Top height and diameter at breast height (DBH) are calculated for each stem and used in a taper function to calculate total and individual gross and merchantable volumes.

Model application
FORECAST has four stages in its use: 1) data assembly, input and validation; 2) establishing the ecosystem condition for the beginning of a simulation run (by simulating the known or assumed history of the site); 3) defining a management and/or natural disturbance regime; 4) simulating this regime and analyzing model output. The first two stages represent model calibration. Calibration data are assembled that describe the accumulation of biomass (above and below-ground components) in trees and minor vegetation for chronosequences of stands developed on sites that vary in nutritional quality. Tree biomass and stand self-thinning rate data are often generated from the height, DBH and stand density output of traditional empirical growth and yield models in conjunction with speciesspecific component biomass allometric equations. To calibrate the nutritional aspects of the model, data describing the concentration of nutrients in the various biomass components are required.
FORECAST also requires data on the degree of shading produced by different quantities of foliage and the response of foliage to different light levels (this information is derived from literature values, field measurements, or simulation models). A comparable but simpler set of data on minor vegetation must be provided if the user wishes to represent this important ecosystem component (e.g., Royo and Carson 2006). Data are obtained from the literature or field measurements. Lastly, data describing the rates of decomposition of various litter types and soil organic matter are required for the model to simulate nutrient cycling. A detailed description of the input data requirements can be found in Kimmins and others (1999), or from the corresponding author.
With the calibration data obtained from different sources, the model calculates the annual rates of different ecological processes (tree growth, litterfall production, mortality, etc.) based on the historical data on tree growth and density provided by the user. Therefore, for each plant species for which historical data are provided, the total net primary production (TNPP) that occurred for each annual time step (t) is calculated with Eq. 3. TNPP t = Δbiomass t + litterfall t + mortality t where Δbiomass t = the sum of the change in mass of all the biomass components of the particular species in time step t; litterfall t = the sum of the mass of all ephemeral tissues that are lost in time step t (e.g., leaf, branch, bark and reproductive litterfall, and root death), and mortality t = the mass of individual plants that die in time step t. Change in biomass (Δbiomass t ) in each time step is derived from a series of age-biomass curves created with empirical data. Litterfall is calculated using userdefined values based on empirical litterfall rates. Mortality is derived from a series of age-stand density curves created with empirical data (for a detailed description on mortality simulation in FORECAST, see Kimmins et al. 1999). Mortality is calibrated through two different parameters: curves of historical stand density for different ages and the proportion of mortality that is due to nonintraspecific competition factors.
The model also estimates the shade-corrected foliage N content (SCFN), which represents the amount of fully illuminated foliage N that was required to produce the calculated historical TNPP. To estimate foliage shading, FORECAST simulates canopy foliage biomass as a "blanket" that covers the stand and that is divided in several layers of 0.25 m height, each of them increasingly darker from the top to the bottom of the canopy. The light absorbed by each layer is calculated based on the foliage biomass present in each time step and a user-defined empirical curve of foliage mass-proportion of full light. Once an estimation of self-shading has been completed for a particular time step using the method described above, FORECAST calculates a foliar N content adjusted for the effects of selfshading (Eq. 4 and 5).
(4) FN t,i = foliage biomass t,i x foliar N concentration where FN t,i = mass of foliage nitrogen in the ith quarter-meter height increment in the live canopy at time t, PLSC i = photosynthetic light saturation curve value for the associated light level in the ith quarter-meter height increment in the live canopy, n = number of quarter-meter height increments in the live canopy at time t. The mean photosynthetic rate of the foliage in canopy level i is calculated by combining simulated light intensities in canopy level i with input data that define photosynthetic light saturation curves for the foliage type in question. Finally, the driving function curve for potential growth of a given species in FORECAST is the shade-corrected foliar nitrogen efficiency (SCFNE) calculated for each annual time step (t) with Eq. 6: When data describing the growth of a species on more than one site quality (i.e. nutrient availability) are provided, SCNFE function curves will be generated during the calibration stage for each site quality. To calculate the nutritional aspects of tree and plant growth, FORECAST requires data on nutrient concentration in each different tree organ. Nutrient dynamics in this study were restricted to nitrogen (the most limiting nutrient at this region (Wang et al. 2013, Bi et al. 2007, Blanco et al. 2012, Wei et al. 2012).
The combination of light and nutrient limitation is usually not enough to explain complex ecological patterns through models, and also including understory vegetation in the simulations is recommended (Kimmins et al. 2008). Therefore, a comparable but simpler (e.g. no data on bark, wood, mortality, etc.) set of data for understory vegetation must be provided to represent this ecosystem component. Lastly, data describing decomposition rates for various litter and humus types are required to simulate nutrient cycling. Decomposition rates are defined by the user (using values from empirical studies) and are affected by site quality, which in turn is defined depending on nutrient and water availability. Snags and logs are tracked by placing them into different categories depending on their original sizes (with slower decomposition rates for snags and for stems with larger sizes).
The second stage of calibration requires running the model in "set-up" mode to establish initial site conditions. In this stage, the model is run with nutrient feedback turned off to allow it to accumulate vegetation, litter and soil organic matter representative of the site(s) to be modeled, and which reflects the historical patterns of accumulation. This is typically achieved by simulating the known or estimated natural disturbance and/or management history of the site (see Seely and others 2002; or Blanco and others 2007 for a detailed description of this process).
After calibrating, estimating the historical ecological rates, and creating the initial conditions, the model is ready to simulate each particular scenario. During the simulation stage, for each annual time step, the annual potential growth (APG) of vegetation is driven by the photosynthetic production of the foliage biomass (Eq. 7). The productive capacity of a given quantity of foliage biomass (photosynthetic rate) is assumed to be dependent on foliage nitrogen content corrected for shading created by the canopy of the simulated site (SCFN t *). SCFN t * is different from the SCFN t that was previously calculated during the internal calibration stage. During the simulation stage the canopy simulated corresponds to the site defined by the user for that particular scenario, which can be different from the empirical canopy data used (i.e. different stand density) during the calibration stage, and therefore SCFN t * is particular for each simulation.
where: APG (t+1) = annual potential growth for a given species in the next time step. During the simulation stage, the model interpolates between the different curves of SCFNE calculated before to find the site quality of the simulated site. Nutrient uptake requirements to support APG are calculated based on rates of biomass growth and data on nutrient concentration in the different biomass components. Nutrient availability is calculated based on empirical data describing litter and humus decomposition rates, changes in chemistry as decomposition proceeds, and the size of nutrient pools in the mineral soil and humus (cation exchange capacity (CEC) and anion exchange capacity (AEC), respectively). If the availability of nutrients for each time step is less than required to support APG, vegetation growth is limited by nutrients and the realized annual growth is lower than APG.
Nitrogen cycling in FORECAST is based on a mass balance approach ( Figure 4) where N can exist in three distinct pools: 1) the plant biomass pool; 2) the available soil nutrient pool, and 3) the soil organic matter/forest floor pool. Inputs and outputs of N to the ecosystem are simulated in a fourstage process for each annual time step. The "available N" pool in FORECAST can be assimilated to represent the interchangeable N present in the soil during one year as NH 4 + , NO 3 or labile organic N fractions with turnover rates shorter than one year. N deposition and N fixed by bryophytes and other microorganisms are simulated as constant annual N fluxes that directly reach the soil solution and are incorporated into the available N pool. Annual values of available N are calculated by simulating consecutively the different inputs and outputs of the biogeochemical cycle: deposition, fertilization, seepage, leaching, mineralization, immobilization ( Figure 4). The simulation of each of these fluxes in FORECAST has been described in detail before (Kimmins et al. 1999, Blanco et al. 2012. The definition of site fertility based on N availability assumes that soil moisture is not limiting in these sites (Blanco et al. 2012, Wei et al. 2012, Wang et al. 2013). However, soil moisture is still implicitly affecting the simulation by the use of the parameter "maximum foliage per tree" which is directly correlated with soil moisture availability (Kimmins et al. 1999).

Figure S2.
Estimation of available N in FORECAST in each annual time step.
Step 1: geochemical inputs were calculated, with all the forms of N lumped together.
Step 3: Plants uptake the available N.
Step 4: Soil N remaining for next time step is calculated by subtracting the remaining N from the soil CEC (for ammonium) or AEC (for nitrate). The N excess was assumed to be lost via leaching.
Carbon and nitrogen cycles are linked through the use of the foliar N efficiency as the driving function of the model (amount of biomass generated in a year per kg of foliar N). Therefore, a limitation in N uptake will result in a reduction of foliar N, reducing biomass produced by the trees.
Nutrient uptake demands on sites of different N fertility are based on observed biomass accumulation rates and tissue nutrient concentrations on these sites, allowing for internal cycling of nutrients.

FORECAST model evaluation for Chinese fir and Phoebe bournei plantations
Published field data on several chronosequences at different site qualities were used to evaluate FORECAST performance for Chinese fir plantations in SE China (Tian 2003, Rong et al. 2008). Data to evaluate FORECAST performance for Phoebe bournei plantations were obtained from literature (Ma et al. 2008, Peng 2003, Peng 2008a, Peng 2008b, Wu 2009, Liao et al. 1989, Chen et al. 2007, Cai 2009, Liang et al. 2009Li 2003, Long et al. 2011, sun 2008, tong 2010, Wei an Ma 2006, Zhang and Wu 2007. For our purposes we used the data described for a good site (27 m dominant height at stand age 50 years) and in areas with low levels of N deposition (5 kg ha -1 y -1 ), assuming that this is the historical level of N deposition in which the forests described by Tian (2003) and Rong et al. (2008) grew (Wei et al 2012) To assess the performance of FORECAST relative to field observations, data pairs of observed vs. predicted were subjected to graphical comparisons, assessments of average and absolute biases, and measures of goodness-of-fit (Blanco et al. 2007). A linear regression of predicted vs. observed values was fitted to calculate the coefficient of determination (r 2 ). In addition, two different indices were calculated.
The first performance index was Theil's inequality coefficient (Theil 1966): where D i is the difference between Observed i and Predicted i, and n is the number of data pairs. U can assume values of 0 and greater. If U = 0 then the model produces perfect predictions. If U = 1 the model produces predictions of system behaviour that are not better than assuming the system does not change. If U > 1, then the predictive power of the model is worse than the no-change prediction. The second index was modelling efficiency (ME) (Vanclay an Skovsgaard 1997): This statistic provides a simple index of performance on a relative scale, where ME = 1 indicates a perfect fit, ME = 0 reveals that the model is no better than a simple average, while negative values indicate poor performance. Finally, the critical error e* was calculated for two different confidence levels 64 . This error can be interpreted as the smallest error level, in absolute terms, which will lead to the acceptance of the null hypothesis (i.e. that the model is within e* units of the true value). If e* is lower than the accuracy level defined by the model user (the minimum acceptable difference between observed and modelled values), then the model is accepted as suitable for the model user's needs. All statistical analyses were carried out using JMP version 5.0.1.2 from SAS Institute.
All the indices of model performance indicated that FORECAST produced acceptable predictions (Table S1), although the results were better for Chinese fir than for Phoebe bournei. Model predictions were better for the first half of the rotation (until year 25), after which there was a slight tendency to underestimate DBH, aboveground biomass and forest litter mass ( Figure S4). Critical error values were low, less than 10% of the maximum value of all variables except forest floor litter mass, for which the values were 20% and 16% of the maximum mass for confidence levels of 95% and 80%, respectively.   Low average biases and absolute differences, together with high Pearson's r values, indicate acceptable agreement between observed and predicted values. However, it has been argued that r is not the most reliable measure of model performance (Power 1993) because it is not related to the "perfect fit" line (the line in which observed equals predicted). As a consequence, this coefficient is more about a model's capacity to use the calibration data set to reduce differences between observed and predicted values rather than a measure of the accuracy of a model's predictions. A different measure of model performance is given by Theil's U coefficient, whose values were always lower than 1, indicating that the model always performed better than a general average value such as that provided by traditional growth and yield tables. Modelling efficiency, recommended as a more adequate measure of model performance (Power 1993, Mayer and Butler 1993, Smith et al. 1997) was close to 1 for all variables, indicating acceptable agreement between observed and simulated values. forests (Seely et al. 2008). Therefore, FORECAST appears to be a valuable tool for studying ecological processes in forest ecosystems in situations when robust predictions are needed.

References used for the review analysis
The following tables provide the details of all the studies used to estimate carbon pools in subtropical plantations.