Impact of fire and harvest on forest ecosystem services in a species‐rich area in the southern Appalachians

To mitigate and adapt to climate change, forest carbon sequestration and diversity of the ecosystem must be included in forest management planning, while satisfying the demand for wood products. The future provisions of ecosystem services under six realistic management scenarios were assessed to achieve that goal. These services were carbon sequestration, types and quantities of roundwood harvested, and different indicators of forest health—biomass of major species, species diversity, and variation of tree age. A spatially explicit forest succession model was combined with statistical analyses to conduct the assessment at the level of both the whole forest landscape and different ecological zones (ecozones) within. An important aspect of this study was to explore the effects of the biophysical heterogeneity of different ecological zones on the outcomes of different management scenarios. The study area was located in an area of the southern Appalachian Mountains in North Carolina with high tree diversity and active forest management activities. Along with a range of management practices, such richness in diversity allowed us to examine the complexity of the interaction between management activities and species competition. The results showed that fire suppression had a greater effect on increasing biomass carbon sequestration than any management scenario that involves harvest and replanting afterward, but at the expense of other indicators of forest health. The effect of fire on species composition was the largest in the xeric parts of the study area. Based on the study results, it was proposed that a low harvest intensity with a mix of fire and fire suppression across the landscape would best balance the need for roundwood products, biomass carbon sequestration, and desirable species composition. This study also demonstrated that the combination of a spatially explicit forest succession model and statistical analyses could be used to provide a robust and quantifiable projection of ecosystem service provisions and possible trade-offs under different manage-

ecosystem must be included in forest management planning, while satisfying the demand for wood products. The future provisions of ecosystem services under six realistic management scenarios were assessed to achieve that goal. These services were carbon sequestration, types and quantities of roundwood harvested, and different indicators of forest health-biomass of major species, species diversity, and variation of tree age. A spatially explicit forest succession model was combined with statistical analyses to conduct the assessment at the level of both the whole forest landscape and different ecological zones (ecozones) within. An important aspect of this study was to explore the effects of the biophysical heterogeneity of different ecological zones on the outcomes of different management scenarios. The study area was located in an area of the southern Appalachian Mountains in North Carolina with high tree diversity and active forest management activities. Along with a range of management practices, such richness in diversity allowed us to examine the complexity of the interaction between management activities and species competition. The results showed that fire suppression had a greater effect on increasing biomass carbon sequestration than any management scenario that involves harvest and replanting afterward, but at the expense of other indicators of forest health. The effect of fire on species composition was the largest in the xeric parts of the study area. Based on the study results, it was proposed that a low harvest intensity with a mix of fire and fire suppression across the landscape would best balance the need for roundwood products, biomass carbon sequestration, and desirable species composition. This study also demonstrated that the combination of a spatially explicit forest succession model and statistical analyses could be used to provide a robust and quantifiable projection of ecosystem service provisions and possible trade-offs under different management scenarios.

INTRODUCTION
The increasing demand for forest roundwood products may be in conflict with the need to mitigate climate change, and to maintain ecosystem diversity, forest carbon sequestration, and species diversity. Thus, there is a pressing need to understand the effects of different types of management on forest dynamics. In particular, the interactions between the ecological processes and management prescriptions are critical because they can result in spatial and temporal differences in the trade-offs between ecosystem services (Seppelt et al. 2011). Some management schemes favor one ecosystem service over others (Schwenk et al. 2012). For example, intensive management such as clear-cutting yields the highest volume of roundwood harvested, but regrowth results in large patches of even-aged trees, which could reduce habitat diversity (Keenan and Kimmins 1993). Management schemes that focus on carbon sequestration may achieve their goal at the expense of other ecosystem services, such as timber production and habitat quality (Seidl et al. 2007, Schwenk et al. 2012). Studies show that some level of timber harvest can enhance biodiversity and increase carbon sequestration (Davis et al. 2009, Fedrowitz et al. 2014. Diverse landscapes with various forest types, topography, and soils provide multiple ecosystem services, and the amount varies greatly depending on locations. Having a robust understanding of such interactions can help predict if a particular policy may succeed or fail to achieve the desired balance (Carpenter et al. 2009).
Systems of coupled human-nature interactions are complex and often involve non-linear processes and spatial variation (Liu et al. 2007). Many previous studies of the trade-offs between roundwood production and forest carbon sequestration have been highly simplified. For example, the spatial dimension, species competition, forest succession, the influences of natural disturbances, and types of roundwood harvested are often not considered (Seidl et al. 2007, Touza et al. 2008, Schwenk et al. 2012, Solarik et al. 2012, Niinim€ aki et al. 2013, M€ akip€ a€ a et al. 2014, Lutz et al. 2016). However, some recent studies (Scheller et al. 2019, Remy et al. 2019, Flanagan et al. 2019) have recognized the importance of these. With new modeling capabilities, it is now possible to analyze the outcomes of forest management in the context of these aspects of forest ecology. In addition to intensity, the location and the timing of roundwood harvest affect both the income of forest owners and carbon accounting.
Following harvest, the time for the carbon to be released to the atmosphere varies depending on the type of harvested roundwood, whether it is sawlog or pulpwood. In general, it can take up to a century for the carbon stored in a sawlog to be completely released to the atmosphere, but only about three years for pulpwood (Smith et al. 2006, Dymond et al. 2012). However, their relative effects on ecosystem carbon sequestration over time remain a major uncertainty (Davis et al. 2009, Nunery and Keeton 2010, McKinley et al. 2011.
Forests that consist of a variety of cover types and stand structures at different successional stages provide a range of habitat niches that support more diverse non-tree species (Farnsworth et al. 2015, Trumbore et al. 2015. Disturbances can create variations in tree species composition and stand structure, as they have important impacts on forest carbon dynamics and vegetation structure (Kurz et al. 2008, Running 2008, Trumbore et al. 2015. Some natural and anthropogenic disturbances have been found to lower carbon stocks, increase emissions, or decrease sequestration rates. They have also been found to affect indicators of forest health, including the forest average ages, the age variability, and species composition. In turn, tree species composition and stand structure affect the richness of non-tree species, such as birds (Goetz et al. 2007), bryophytes, lichens (Neumann and Starlinger 2001), and the herbaceous-layer (Gilliam et al. 1995, Gilliam 2002. The three forest health indicators show how well a forest habitat can support biodiversity (Gao et al. 2014, Trumbore et al. 2015. Therefore, knowing the locations and the impact of disturbances is important because it can help allocate different treatment plans in different locations to improve forest health (Thompson et al. 2011, Martin et al. 2015, Creutzburg et al. 2017, Krofcheck et al. 2017, 2018. It is particularly interesting to conduct the analysis of ecosystem service trade-offs in a species-rich area, because tree species diversity promotes productivity in temperate forests mostly through strong complementary effect between species (Morin et al. 2011). To conduct such analysis in a species-rich forest, a forest would be subdivided into smaller ecological zones as the analysis unit. An ecological zone (ecozone) is where climate and soil are assumed to be homogenous (Simon et al. 2005, Scheller et al. 2007. It is valuable for managers to have research results in finer ecozone scale in addition to the landscape scale (USDA Forest Service at North Carolina 2011). To illustrate this concept, an area with active management and high diversity was chosen. It was located in the southern Appalachian Mountains of the United States: It included 12 different vegetation communities as defined by Simon et al. (2005), 20 common tree species, and a range of management practices within an area small enough to model with operational-level details.
The objectives of this study were to further our understanding of the impacts of fire, as well as various harvest and replanting prescriptions, on the following: (1) carbon sequestration in the standing forest biomass; (2) the amount of roundwood (sawlog and pulpwood) harvested; and (3) forest species diversity, composition, and age structure, which are indicators of forest health. In this study, it was hypothesized that no fire and intensive harvest would maximize the amount of roundwood produced at the expense of biomass carbon sequestration and the quality of the forest habitat to support biodiversity. Furthermore, when the harvest and replanting intensity decreased, the volume of total harvested roundwood would decrease, but the forest was expected to have a higher proportion of ecologically desirable species. In addition, the inclusion of fire was expected to decrease the average forest age while opening up the canopy to allow more early-successional species to grow.

Study area
The study area is located in the Grandfather Ranger District (GRD) within the Pisgah National Forest. It is about 777 km 2 in size and is located in the eastern edge of the Southern Blue Ridge Ecological Province, in North Carolina in the southeastern USA (Fig. 1). The elevation ranges from 314 to 1810 m above sea level, changing sharply within short distances (USDA Forest Service at North Carolina 2011). The steep topography, in combination with the soil type and climate, results in a diverse vegetation composition in the area (Pittillo et al. 1998). The GRD has 12 ecozones ( Fig. 1) that belong into three ecological groups: "Northern hardwood forest," "Mesic forest," and "Xeric forest" (Simon et al. 2005). The GRD is rich in plant and animal diversity and is facing pressure for increased timber harvest (Fox et al. 2007).
According to the managing agencies of the GRD, due to historical management methods, five ecozones have "departed from desirable conditions" (USDA Forest Service at North Carolina 2011). Those ecozones include "Pine-Oak Heath," "White Pine-Oak Heath," "Shortleaf Pine-Oak," "Rich Cove," and "Acidic Cove." All of them lack quality early-successional habitat, diverse tree age structure, open canopy system, and high-quality soil.
Fire had been the most frequent natural disturbance in the southern Appalachian region before the U.S. Forest Service started actively suppressing fires in the early 1900s, but fire suppression is being re-examined more recently by land managers in the area because of the unintended consequences of reducing the diversity and natural variability across the landscape (Waldrop et al. 2016, Lafon et al. 2017. The exact impact of fire on vegetation, with climate change, in the southern Appalachian region is a current active field of research: for example, Lafon et al. (2007), Flatley et al. (2011), Flatley et al. (2013, and Waldrop et al. (2016). An estimation of one of the possible natural fire scenarios was used in the model simulation of this study.
Field data.-Field plot data were used for model calibration and validation. The field data were obtained from the USDA Forest Service (USFS) Forest Inventory Analysis (FIA) data version 6.0.1 (USDA Forest Service 2016), available from 1974 to 2015. The FIA consists of assessments of field plots of 0.4 ha that include information on tree age, species, site condition, forest type, and estimated biomass for individual trees. The sampling intensity of the FIA plots is, on average, one plot for every 2428 ha, which was designed to provide an unbiased representation of forest types. Plots are censused every five to seven years (Riemann et al. 2010). The biomass data for FIA plots are estimated using allometric equations provided by Jenkins et al. (2003). This study used the plot data of live trees located in areas classified as forestland, from 2002 to 2015.
In order to preserve confidentiality for private forest data, FIA varies the locations of the plot centers randomly to within a radius of 1.6 km from the exact location. Some plots located in private forests are swapped with plots of similar conditions, but within the same county. Since the plots are swapped with nearby plots of similar forest conditions, the data remain valid for the forest type measured. Therefore, not having the precise location information could only marginally affect the use of the FIA data in this study.
Species.-Twenty tree species were chosen for this analysis (Table 1). They compose the largest portion of biomass, out of the 88 species reported by the FIA data in the southern Appalachian region in North Carolina. The selected species occur at different successional stages, with some of them being resistant to fire and some of them being shade-loving. To understand better the impact of the management on forest health and the income from the forest, as well as to better model the species that would be planted and harvested in different management scenarios, two rankings were assigned to all the species: the ecological rank, based on the ecological role of each species; and the commercial rank, based on the type of wood sold and products made. The ecological rank of the species was based on the potential habitat and food that a species can ❖ www.esajournals.org provide to non-tree species; the commercial rank of the species was based on the potential commercial use of each species harvested and the average real prices of different types of roundwood in the recent five years (Appendix S1: Section 5). Both ranks were provided by local experts working in the National Forests in North Carolina, including forest silviculturists, land managers, and forest ecologists (J. A. Rodrigue, personal communication, Oct 2015). Based on the ranks, species were grouped into 4 classes based on ecological desirability (Table 1): from Class A, the most desirable, to Class D, the least. Class A consisted of pitch pine, white oak, and eastern hemlock, which are late-successional, climax species. For example, eastern hemlock is a keystone species in riparian forests and acidic (Brantley et al. 2013). Class B consisted of some mast producers and climax species. Class C consisted mostly of early-successional species. The three species in Class D, yellow poplar, eastern white pine, and red maple, are aggressive species with capability to displace other species. Species were also classified into 5 classes based on the commercial value: from Class I, the most economically valuable, to Class V, the least. Yellow poplar and eastern white pine belonged to Commercial Class I, but were also two of the three species in Ecological Class D.

Management targets and scenarios
Management prescriptions and goals can differ between forest owners. Harvest rates and replantation currently vary owing to the different management objectives of different spatial units (National Forests in North Carolina and USFS 2016). The forest in the GRD consists of private, public, and Congressional Designated Roadless areas (Appendix S1: Fig. S3). Management practices in private forests were assumed to aim at maximizing the economic return from selling harvested roundwood, but the harvest prescriptions varied depending on the protection status of the land unit as classified by the North Carolina GAP study (McKerrow et al. 2006); public forests were subdivided into smaller patches called the "Management Areas" (MAs) that took into account the specific ecological and logistical  Notes: "Comm. Rank" and "Eco. Rank" stand for Commercial Rank and Ecological Rank, respectively. The rankings were provided by local experts based on resulting harvested timber products and the contribution to its community of each species. Smaller rank index indicates greater commercial or ecological importance. Based on the ranks, the species were grouped into four ecological classes (Eco. Class) and 6 commercial classes (Comm. Class). The three species with the lowest ecological rank are considered as undesirable aggressive species that are capable to displace other species under current forest disturbance patterns.
concerns, and were managed accordingly. The MAs of the public forests were designated by the National Forests in North Carolina and USFS (2016). No harvest takes place in the Congressional Designated Roadless areas, where the U.S. federal law prohibits any forest management activity including logging or road construction.
Six management scenarios were simulated. The scenarios consist of two harvest and planting prescriptions and no prescription, each with and without fire: (1) no harvest and planting ("reference"); (2) clear-cutting and planting economically valuable species ("Aggressive"); and (3) species selection harvest and planting ecologically valuable species ("Moderate"). The two prescriptions, the "Aggressive" and the "Moderate," each consisted of different degrees of harvest and planting that varied depending on the landownership (Table 2). Specific management prescriptions were applied to the forest stands depending on the MAs. The harvest intensities of both prescriptions were similar to those the land managers were planning to implement (J. A. Rodrigue, personal communication, Jan 2016). The spatial information of the stands in the GRD was obtained from the Continuous Inventory of Stand Conditions (CISC) data based on the National Forest Stands for the Southern Appalachian Assessment (SAA) Study Area (SAMAB, 1996). Although the Aggressive prescription consisted of various levels of clear-cutting, its harvest rates in different MAs, calculated based on the estimate from the historical forest disturbance maps in between 1990 and 2005 of the GRD, were still lower than those of the historical clear-cutting in the Appalachian regions (Yarnell 1998). The historical disturbance maps in the GRD were obtained from the forest disturbance history maps derived from the Vegetation Change Tracker (VCT) algorithm (Huang et al. 2010). The major disturbances in the area have been fire, harvest, and insect outbreak (USDA Forest Service at North Carolina 2011), but the disturbances detected at 30 m resolution were mostly fire and harvest, which were the disturbance types simulated in this study. The "Moderate" prescription simulated the practice of the U.S. Forest Services in the National Forest, which can be viewed as restoration management regimes.
The MAs for private forests were designated based on the protection status of the land unit. The harvest intensity of the private forest located in protected areas was the same as that of the public forest with the same protection status. Prescription details in each MA and protection status are provided in Appendix S1: Section 6. Harvest events occurred every five years in the locations within each MA where the criteria of harvest detailed in the prescription were satisfied. The current complete harvest rotation of the southern Appalachian area of 80 yr was adopted in this study (Fox et al. 2007).
A no-fire scenario was simulated that resembled the condition under complete fire suppression management currently undertaken by the USFS (Flatley et al. 2013). For comparison, a natural fire scenario was also simulated. It was based on an estimation of historical fire (Xi et al. Note: "Roadless" and "Ecology" are public forest areas with special designated protected status where minimal or no harvest activities can take place, compared to other publicly owned forest. 2009), where about 0.4% to 2.6% of each ecozone was burned in each fire event (see Appendix S1: Section 7 for details).

Ecological model
Forest change over 100 yr was simulated using the forest succession model, Landis-II (Scheller et al. 2007) at 150 m resolution (2.25 ha grid). Landis-II is a mechanistic, stochastic, forest succession model that can be used to simulate forest response to different disturbances and provide geospatial results that are based on the interactions of individual trees and their neighbors including species productivity, mortality, reproduction, natural disturbances, and management activities Mladenoff 2004, Scheller et al. 2007). The stochastic components of the model that were used in this study came from fire ignition and seed dispersal (Scheller and Mladenoff 2004). Since it is a spatially explicit model, Landis-II can simulate the spatial aspect of inter-and intraspecies competition, including processes such as seed dispersal in the landscape and succession processes. The Forest Carbon Succession extension version 2.0 (Dymond et al. 2016) was used to account for the carbon and biomass dynamics in the forest during the succession process. It also provides information on the carbon in the trunk and branches of the tree that goes to the harvested product pool (Dymond et al. 2016).
In addition to the core model, two extensions of the model were used to simulate the impact of the harvest and fire disturbances: the Base Harvest and the Base Fire. The Base Harvest version 3.0 allows the simulation of different harvest methods, areas harvested, and rotation length (Gustafson et al. 2000). Wildfires were generated based on the physical characteristics of each ecozone in the Base Fire extension version 3.0.3 (He and Mladenoff 1999). Parameterization of the fire information, such as the probability of ignition, and the fire cycle of each ecozone were based on the estimation from the studies in the southern Appalachians before human settlement by Lafon et al. (2017), as well as Wade et al. (2000), Xi et al. (2009), Fesenmyer andChristensen (2010), and Flatley et al. (2013). Fire cycle in each ecozone for each management scenario is provided in Appendix S1: Section 7.

Model
parameterization, calibration, and uncertainties.-Landis-II models forests as a grid of interacting cells (Scheller et al. 2007), with each cell assigned to an ecozone (Fig. 1). The initial conditions of the species and age cohort composition in each cell were calculated based on the FIA data in between 2002 and 2010. Then, the model was simulated without any disturbance for 26 yr. Output from the model was in units of mass of carbon and was converted to biomass by assuming carbon is half of total biomass (Smith et al. 2006). The value of the aboveground live trees biomass in each cell changes with time in the model simulation based on species, the impact of disturbances, growth, and age-dependent mortality. Details of the initial values of climate, soil, and species are provided in Appendix S1: Section 1. The values of the live tree aboveground biomass density chosen from the FIA database and used for model calibration and validation were those located in forestland where there had not been any disturbance between 1990 and 2015: Those between 1990 and 2010 were used for calibration and those between 2011 and 2015 for validation of the model results. Calibration was done at the plot and species levels, while validation was done at the plot level. Details on the FIA data used for calibration and validation are provided in Appendix S1: Section 1.13.
To capture the stochastic components of the Landis-II model, each management scenario was simulated five times to estimate the betweenrun variability (Thompson et al. 2011). The standard deviation of the mean of the results of the between-run variability for each management scenario was used to determine the significance of the difference between management scenarios. The contributions of the uncertainties in input parameters were analyzed using sensitivity analysis. The inherent uncertainties in the model projection results due to the uncertainties of the input parameters, the calibration, the accuracy of the model projection, etc., were assumed to be the same for each scenario. In all the six management scenarios that were simulated in this study, the two scenarios where no harvest and planting occurred served as reference scenarios and were used to compare the impact of fire and harvest and planting prescriptions.

Analysis of the simulation results
Carbon sequestration was simulated for the aboveground standing biomass, defined by the net change of carbon in the forest biomass between simulation years 1 and 100. The net change was divided by 100 to obtain the average annual carbon sequestration. The volume of roundwood harvested was estimated by converting the annual carbon flux to the harvested product pool into biomass, then into wood volume. An average specific gravity of both hardwood and pulpwood in the southeastern United States, 0.49 t/m 3 (Smith et al. 2006), was used to convert from biomass to volume units, in which commercial roundwood production is measured. Harvested roundwood was divided into sawlog and pulpwood harvested. For species that belong to Commercial Classes I to III, the harvested trunk was considered as sawlog while the branch portion was considered as pulpwood, except for harvested Virginia Pine, and all species that belong to Commercial Classes IV and V were considered to be all harvested pulpwood. Harvest density is a useful indicator for the productivity of the land and was calculated by dividing the harvested volume by the area of harvest.
The results could be used to estimate the effect of the simulated forest on the diversity of the landscape if all organisms were present. Four indicators were established to quantify the quality of the forest habitat to sustain biodiversity: (1) the relative contribution of biomass in each ecological class to the total biomass; (2) the Shannon diversity index based on the biomass of each species; (3) the average age; and (4) the standard deviation of age (Gao et al. 2014, Farnsworth et al. 2015, Trumbore et al. 2015. Older forests and those having a wider range of tree ages can provide more distinct niches, and hence greater diversity (Spies 2004, Gao et al. 2014). To measure for the higher probability of niche diversity, the distribution of the species biomass and the average and the standard deviation of ages of all species for the entire GRD and for the individual ecozones at the end of the simulation year were also assessed. Biomass of each tree species in each management scenario was plotted to show the distribution of the species biomass. Slope of the biomass distribution shows the species equitability: The steeper the slope, the more unevenly distributed the species. The Shannon diversity index (Shannon 1948) was used to assess the tree species abundance and richness after 100 yr. Landis-II does not track individual trees, but models the species present, their biomass, and also the number of cohorts of different ages (Scheller et al. 2007), so the proportion of the biomass of each species was used in the calculation of a modified Shannon diversity index (H') as follows: where p i is the proportion of biomass of species i. The use of relative biomass in the Shannon diversity index rather than the relative species abundance results in an index of diversity that is more sensitive to the relative biomass (Dickman 1968), and energy distribution among species (Wilhm 1968), which is more appropriate for this study. For example, it corrects the discrepancy which would arise where a species with high number of individuals but only small relative biomass, such species would otherwise dominate the diversity index values. However, the Shannon diversity index does not account for the ecological values of species (Hurlbert 1971) nor the species age structure. But both the ecological values of the dominant species and the distribution of the species age structure are also important indicators of the likely ability of a forest habitat to sustain biodiversity (Noss 1999, Tews et al. 2004, Leinster and Cobbold 2012, which were also analyzed in this study.
To compare the Shannon diversity index between each pair of management scenarios, bootstrapping was used to generate 1000 samples of the Shannon diversity index for each scenario. The biomass of each species in each management scenario for all of the five simulations was used to calculate the Shannon diversity index samples used in the bootstrapping. Pairwise t-tests were used to examine if the differences of the indices between each pair of management scenarios were significant (Hutcheson 1970), where the P-values were adjusted by the Benjamini-Hochberg method (Benjamini and Hochberg 1995) with the false discovery rate of 0.01. For the differences that were significant, the effective numbers of species, calculated by taking the exponent of the Shannon diversity index, were calculated, so that their differences could be ❖ www.esajournals.org compared on a linear scale (MacArthur 1965, Jost 2006, Leinster and Cobbold 2012. Effect of hemlock woolly adelgid infestation.-The eastern hemlock is threatened by an invasive insect, the hemlock woolly adelgid (HWA) (Elliott and Vose 2011). The FIA biomass data for the eastern hemlock may not reflect the overall damage by HWA, because of the frequency of revisiting the same plots (5-7 yr) and the possible selection bias because the data report the biomass of only healthy trees and not those killed by HWA. Species composition of the area changed drastically in the complete absence of the eastern hemlock as a result of HWA infestation in the GFD, as shown in a simulation study (Birt et al. 2014). In order to understand the impact of the reduction of eastern hemlock biomass due to HWA infestation, sensitivity analysis on the management scenario that had the highest proportion of resulting eastern hemlock biomass was conducted, as that scenario would be the most sensitive to any change on the eastern hemlock parameters. The eastern hemlock parameters were altered in three ways to test the resulting species and age composition: (1) no planting of eastern hemlock; (2) reducing the maximum biomass value of eastern hemlock by 70%, based on observed eastern hemlock biomass loss from the HWA infestation of 1987 in southern New England (Small et al. 2005); and (3) both no planting and maximum biomass reduction.

Statistical analyses
Besides directly comparing the Landis-II results, a regression analysis that controlled the impacts of the area of each ecozone and heterogeneity between ecozones was performed to assess the impact of each management scenario on the four indicators that were used to quantify the quality of the forest habitat. Those four indicators were as follows: (1) biomass of each ecological class; (2) the Shannon diversity index; (3) the average age; and (4) the standard deviation of age. Regression analysis was performed because it can provide quantification of the effects of various management scenarios. Moreover, since species biomass proportions are compositional data, where each of the proportion is nonnegative and summed up to one, using ordinary least squares (OLS) regression can give misleading results (Aitchison 1982). Therefore, a Dirichlet component regression approach was used in addition to the OLS regression (Aitchison 2003), because it takes into account the compositional nature of the data. In both of the regressions, the ecological classes of species, the Shannon diversity index, age, and its standard deviation were the dependent variables; fire, the harvest and planting prescriptions, and the areas of each ecozone were the independent variables. Natural logarithm transformation was taken for the dependent variables so that the data would satisfy the normal distribution assumption for the regression analysis. Moreover, in an effort to check the model for errors with OLS assumptions, tests were run for high multicollinearity in the regressors, looking for heteroskedasticity and spatial autocorrelation of the error term, using the variance inflation factors (VIF), White's heteroskedasticity-robust standard errors, and Moran's I.

Calibration, validation, and sensitivity analyses
For species-level calibration of Landis-II with the FIA data, the adjusted R-squared of the regression line was 0.65 (P < 0.01) (Fig. 2). Statistical analyses of both calibration and validation at the plot level showed that the model simulations were statistically similar to the FIA data: The simulated average aboveground biomass density was 162.5 t/ha, while that of the FIA was 171.7 t/ha. Since the Fisher F-test showed that the variances of the two distributions were not homogeneous (P = 0.003), a Welch two-sample ttest was used (Ruxton 2006). The results showed that the two groups were similar (P = 0.57). Details of various testing methods and comparison are available in Appendix S1: Section 1.13.
Two variables, the maximum ANPP and the maximum biomass of each species, were varied by 10% to test the sensitivity of the model. Both variables have been shown to be the most influential in Landis-II (Thompson et al. 2011, Simons-Legaard et al. 2015. Overall, the aboveground biomass density showed <10% change with 10% change with the two tested parameters. Details of the sensitivity analysis of this study are provided in Appendix S1: Section 4. Heteroskedasticity and spatial autocorrelation were not found in the data, so the results of the regression analyses could be interpreted without any further transformation in the data. The results of the various tests conducted showed that there was no high multicollinearity in the regressors and no heteroskedasticity and spatial autocorrelation of the error term. The test results showed that the methods used in the statistical analyses were valid. Details of the tests performed and the results are provided in Appendix S1: Section 10.

Effects of management on the long-term ecosystem services
Biomass carbon sequestration and roundwood harvested.-The simulation of the entire GRD showed that fire control was the most important factor influencing carbon sequestration in the aboveground biomass (Table 3). When there was no fire, there was a net annual carbon gain in the standing biomass. Harvest and planting activities resulted in less carbon sequestered than those for the reference scenarios with the same fire suppression status. However, when harvest took place, some of the carbon in the aboveground biomass was transferred to the harvested sawlog, about 33% of which serves as a longterm carbon pool for storing carbon for up to 100 yr in the end-use products such as houses (Smith et al. 2006). The difference between the scenarios with the Aggressive prescription and the reference scenarios in annual live biomass  Table 1. Orange line represents the 1:1 relationship and the red line represents the linear regression of the model and the FIA data. Slope of the regression line is 0.90, and the adjusted r-squared was 0.65, with P-value < 0.01. The outer pair of the red dash lines is the 95% prediction interval, while the inner pair is the 95% confidence interval.
carbon sequestration was about 13% less, while that for the Moderate prescription was about 5% less.
Although harvested sawlog provides more benefits in terms of serving as a carbon storage pool and generating profit, both harvested sawlog and pulpwood provide income. The percentage of harvest area was specified in the model; however, the actual volume of harvest depended on the number of species that satisfied the specified harvesting criteria, such as the requirements on species and age. For the scenarios with Aggressive prescription, there was no significant trend in the harvested density of pulpwood, and the decreasing trend in the harvested density of sawlog was not significant when there was fire. However, when there was no fire, the decreasing trend in the harvested density of sawlog was very significant with a rate of 0.02 m 3 Áha À1 Áyr À1 (P < 0.0001; adjusted r 2 = 0.91). For both of the scenarios with Moderate prescription, there was no significant trend in the harvested density for sawlog, but that for pulpwood increased significantly at a rate of 0.02-0.03 m 3 Áha À1 Áyr À1 (P < 0.001) (Fig. 3).
Species composition and diversity. -The Shannon diversity index in the initial condition was 2.45. The Shannon diversity indices reflect the evenness of the tree species biomass distribution and the number of species. Fire resulted in more diverse forest but with less ecologically preferable species. Dominant species were killed in fire, equivalent to resetting the forest to an earlier successional stage, and hence resulted in a more even distribution of the biomass among species (Fig. 4). The values of the diversity index after 100 yr decreased for the No Disturbance, and the Moderate with no fire scenarios, but increased for all the scenarios with fire or with Aggressive prescription (Table 3). Fire significantly increased diversity by about 3.2% compared with scenarios with no fire, as shown by the regression analysis (Appendix S1: Table S10).
The species distribution of the landscape in the initial condition was dominated by chestnut oak and other species of lower ecological ranks (Appendix S1: Fig. S1). After 100 yr, in the scenario where there was no fire and no planting and harvest, the landscape was dominated by American beech and eastern hemlock. However, in the scenarios with fire, especially in Xeric forest, no one species dominated, with the forest consisting mostly of early-successional species in Ecological Class C. This change was at the expense of ecologically preferred Class A and ecologically least preferred Class D. In other words, fire only affected the proportion of earlysuccessional species (Ecological Class C), as shown in Appendix S1: Table S8. Thus, it can be concluded that the ratio of Ecological Classes A and D appeared independent of fire status and that Ecological Classes A and D were replaced by Ecological Class C in approximately equal proportion in case of fire. It was even more evident when compared to scenarios with no fire.
The Aggressive scenario had significant but different impacts depending on whether fire was present. It significantly decreased the Shannon diversity index by 53% with fire and significantly increased it by 2.6% with no fire (see Appendix S1: Section 9 for regression analysis details). The biomass of the dominant species decreased with fire ( Fig. 4), increasing in the evenness of the species biomass distribution, as indicated by the higher

Notes:
The reported values are the average AE standard deviation of the 5 model runs. Negative numbers indicate a loss. "Fire" indicates the presence of fire; "Avg. age" is the average tree age in year; "SD age" is the standard deviation of age in year.
value of the Shannon diversity index for the scenarios with fire. In those scenarios, any planting or harvesting decreased the evenness of the species biomass distribution, because they favored the newly planted species over the others. However, for the scenarios where fire was absent, later successional species became dominant. Any management activities that did not favor the dominant species, such as planting, increased the biomass of other species and the evenness of the species biomass distribution (Fig. 4).
When the input parameters related to eastern hemlock were altered, the species composition of the landscape changed drastically. The most significant change was when the allowed maximum biomass of eastern hemlock was reduced by 70% and no planting: The eastern hemlock was no longer the dominant species in the landscape being replaced by deciduous species. One such species was the ecologically more desirable white oak (QUAL), which was replaced by the less ecologically desirable, short-lived scarlet oak (Fig. 5). In other cases, the proportion of the biomass of eastern white pine (PIST), the least ecologically desirable, increased significantly. The regression analysis showed that when the allowed maximum biomass of eastern hemlock decreased, the Shannon diversity index increased significantly (about 1.7%), but the proportion of species belonging to Ecological Class A decreased significantly (about 62%) and the early-successional, Ecological Class C, increased significantly (about 42%) (see Appendix S1: Table S11).
Average and the variation of age.-Ecozones with even-age tree distributions were especially sensitive to fire and the Aggressive prescription, especially the Xeric forest. Fire significantly reduced the average age of trees by about 38% and its variation by about 18%. With fire, the Aggressive prescription further reduced the standard deviation of tree age significantly by about 44% (see Appendix S1: Section 9 for regression analysis details). Results showed that the Xeric forest ecological group, ecozones 5 ("Pine-Oak Heath"), 6 ("White Pine-Oak Heath"), and 11 ("Shortleaf Pine"), all consisting of large blocks of even-aged forest (USDA Forest Service at North Carolina 2011), were especially sensitive to fire-the average and standard deviation of age were much lower with fire. In the case of no fire, the ages of the forest were particularly sensitive to the Aggressive prescription (Fig. 6).

Implications for conservation policy
The results showed that the aboveground live biomass carbon transferred to the harvested sawlog was a significant carbon sink. Harvested sawlog production serves as an indicator of the longterm carbon sequestration potential of the forest under different management scenarios. Carbon loss in the live biomass was transferred to the harvested sawlog pool, of which a small portion can remain for between 60 and 100 yr (Dymond et al. 2012, Hoover et al. 2014). In the case when there was no fire, the forest served as a net carbon sink. In this study, when there was no fire, the amount of carbon that was transferred to the long-term carbon pool was not that large, of about 6%, and it would eventually be released to the atmosphere at the end of its useful life. However, in the scenarios where fire occurred, the whole forest landscape became a carbon source. In this study, the amount of carbon that was transferred to the long-term carbon pool ranged from about 13% in the Moderate prescription to about 52% in the Aggressive prescription. Compared to the scenario where no harvest had taken place and the total carbon loss in the forest aboveground live biomass was released to the atmosphere, the harvested sawlog could at least serve as a long-term carbon pool that defers the release of carbon to the atmosphere for some time. This is an important component in more accurate carbon accounting. Although intensive roundwood harvest would bring more income to forest landowners and provide a larger potential long-term storage of carbon, the significant decreasing trend of the volume of the harvested sawlog per hectare in at least a prescription similar to the Aggressive prescription showed that such harvest intensity may not be sustainable for a long term. The annual average harvest volume per hectare of the whole area ranged between 2.60 and 3.45 m 3 Áha À1 Áyr À1 for pulpwood and between 3.62 and 4.95 m 3 Áha À1 Áyr À1 for hardwood (Table 3). The results are similar to the finding in Ling et al. (2016), where harvested roundwood volume per hectare of forestland in the same region between 1986 and 2009 was about 0.08 to 7.35 m 3 Áha À1 Áyr À1 . Although the harvest volume per hectare is similar for both prescriptions, the harvested areas for the Aggressive prescriptions were greater than those in the Moderate prescriptions. Also, the results in this study, as well as other field studies conducted in the area, showed that the density and the basal area of trees decreased under repeated intensive harvest (Leopold et al. 1985, Elliott andSwank 1994).
Conserving eastern hemlock is very important in preserving the species composition of a forest similar to that in the GRD. The results of the simulation of the effect of the hemlock woolly adelgid (HWA) infestation found that with HWA infestation, eastern hemlock would be replaced by less desirable, deciduous species (Small et al. 2005, Brantley et al. 2013, Birt et al. 2014. Such drastic changes in species composition could exert a strong negative impact on the hydrological cycle (Brantley et al. 2013), on non-tree species diversity (Small et al. 2005), and on the carbon and nitrogen cycles (Nuckolls et al. 2009). Furthermore, the replacement tree species that became dominant, American beech and scarlet Fig. 5. Species biomass distribution after 100 yr of simulation for the three cases of altering the parameter for eastern hemlock in the management scenario where there was no fire and the Aggressive prescription. Each management scenario was simulated for five times to account for the between-run variability of the model. Error bars represent three standard deviations from the mean biomass of each species calculated from the five simulations. Green bars indicate the three most ecologically preferred species and the red the least ecologically preferred species. oak, are both deciduous. Not only would that affect the growth of other understory vegetation, but also it could affect the temperature of runoff the streams, affecting the insect population and aquatic ecosystem (Snyder et al. 2002, Ross et al. 2003. In addition, the results showed that moderate harvest and letting fire occur, especially in ecozones that belong to the "Northern hardwood" ecological group, would provide a forest habitat with increased diversity in tree species and ages. As a result, that can be expected to provide more niches for a greater variety of species. However, the amount of roundwood harvested and biomass carbon would be significantly lower. The lack of harvest and planting lowered the diversity of the forest, but the Aggressive prescription that included clear-cutting also decreased the diversity, because both management methods resulted in uniform stand ages of the forest. The impacts of those changes on the diversity of herbaceous species would vary depending on the ecozone (Gilliam et al. 1995, Gilliam 2002. These trade-offs have been reported in a species-rich area in Oregon (Creutzburg et al. 2017).
Xeric forest in the area was sensitive to the fire regime, as recognized by the Nantahala and Pisgah National Forests Land Resource Management Plan (National Forests in North Carolina and USFS 2016). This study results agree with their finding: The three Xeric forest ecozones Pine-Oak Heath (ecozone 5), Shortleaf Pine-Oak (ecozone 6), and Shortleaf Pine (ecozone 11) were especially sensitive to fire (Fig. 6). Fire has also been reported to promote high levels of Table Mountain pine and pitch pine (Lafon et al. 2007), but also decreased both the average and standard deviation of the tree age in those ecozones decreased by about 50% (see Appendix S1: Table S8).

Limitations
A simplified fire simulation approach for management scenarios analysis was used in order to avoid introducing unnecessary model complexity that would be beyond the scope of this paper. The parameters used in the fire simulation were estimated based on existing literature that may not be accurate, and the Landis-II Base Fire extension is based on simple assumptions such as the relationship between fire intensity and Fig. 6. Age distribution for each ecozone for the initial condition and scenarios with different fire suppression status: (a) fire present; (b) no fire. Each number represents an ecozone: "Northern Hardwood" (1), "Spruce-fir" (7), and "High-Elevation Red Oak" (8) belong to the "Northern hardwood forest" group; "Mesic Oak-Hickory" (2), "Acidic Cove" (4), "Rich Cove" (10), and "Dry Oak-Hickory" (12) belong to the "Mesic forest" group; and "Oak Heath" (3), "Pine-Oak Heath" (5), "White Pine-Oak Heath" (6), and "Shortleaf Pine-Oak" (11) belong to the "Xeric forest" group. Specifically, the ecozones belonging to the "Northern hardwood forest" group are located in high-elevation environments, while the ecozones that belong to the other two ecological groups are located in the low-elevation environments (Simon et al. 2005). time since last fire, and does not account for fuel types and their spatial pattern (Sturtevant et al. 2009). This study simulated scenarios with fire and without fire. However, some fire does currently occur in the GRD. Although only one percent of the prescribed burns escaped, the Forest Service does undertake some low-to mediumdensity burning to thin out leaves and woody debris (USDA Forest Service 2018). Future studies using the Dynamic Fire extension of Landis-II model, which accounts for fuel types, topography, and weather, may give a better quantification of the impact of climate change, but additional parameters will be required.
The impact of climate change was not modeled. In this study, temperatures and atmospheric CO 2 concentration were assumed to fluctuate within the range of the past 50 yr. Any increase in the values of these may increase tree growth, and possibly fire size and frequency (McKenzie et al. 2004, Westerling et al. 2006, Littell et al. 2010). This will affect not only the species composition (Remy et al. 2019) and the feedback between vegetation and fire ), but also the economic values of roundwood and hence the choice of species to harvest and plant. Climate change affects the supply of the types of roundwood, and hence its pricing. A detail analysis is needed to accurately model the impact of climate change on the harvested species pattern. Our study presented the possible impact of the scenarios without climate change, which can be served as a foundation to future research that investigates the effect of climate change on fire, species availability, and also harvested roundwood.
Thinning was not simulated in this study, although it is a very common practice that allows the trunk-wood of desirable trees to increase faster (Keyser and Brown 2014). Accounting for the carbon flux to the harvested wood product pool from thinning and the impact of climate change coupled with different management regimes in a species-rich area could be included in a future study.
In the model projection of 100 yr, there may be some uncertainties in the model that were not accounted for during the process of data collection, calibration, and projection (Sklar and Hunsaker 2001). Although sensitivity analyses, calibration, and validation were done to account for uncertainties, the statistical analyses in the study compared the results in each management scenario only with those in its corresponding reference scenario.
One of the limitations of this study is that we used a number of different indicators with different units. It can be more direct to evaluate ecosystem services after converting each into the same unit, but it often involves additional complex analyses or is infeasible. In such cases, reporting the outcomes in biophysical terms is often preferable (Guerry et al. 2015). Statistical analyses that account for the spatial heterogeneity allowed the evaluation of the impact of each management scenario with useful insight for comparison.

CONCLUSIONS
For our study of a region in the southern Appalachians, the results showed that the presence of both eastern hemlock and fire has an important influence on the ecosystem. Despite all the limitations, the simplified fire model could provide a quantifiable impact assessment on fire restoration, achieving the goal to let the forest managers understand the impact of letting some fire return to in the landscape. The results of this study distinguished the different types of harvested roundwood, which, though important for both carbon accounting and consideration of the potential income of forest owners, often is not included in forest succession dynamics or ecosystem service trade-off studies. The results also suggested some level of harvesting and burning would be beneficial to the overall long-term forest health while satisfying the demand on roundwood. The methodologies used in this study could be applied to different places to assess the longterm ecosystem service trade-offs and to assist management decisions.
The restoration project of the GRD planned to restore natural fire, especially in certain ecozones where the forest health had been negatively impacted by years of fire suppression (USDA Forest Service at North Carolina 2011). The results showed that for ecozones that suffered from lack of quality early-successional habitat and even-aged trees, such as the Pine-Oak Heath (ecozone 5), restoring fire would have a big impact of the average and standard deviation of forest age (Fig. 6). The result of this study also showed that some level of sawlog harvest helped delay some carbon that would have otherwise been released to the atmosphere immediately in the case of fire.