Enhanced carbon storage through management for old-growth characteristics in northern hardwood-conifer forests

Forest management practices emphasizing stand structural complexity are of interest across the northern forest region of the United States because of their potential to enhance carbon storage. Our research is part of a long-term study evaluating silvicultural treatments that promote late-successional forest characteristics in northern hardwood-conifer forests. We are testing the hypothesis that aboveground biomass development (carbon storage) is greater in structural complexity enhancement (SCE) treatments when compared to conventional uneven-aged treatments. Structural complexity enhancement treatments were compared against selection systems (single-tree and group) modified to retain elevated structure. Manipulations and controls were replicated across 2-ha treatment units at two study areas in Vermont, United States. Data on aboveground biomass pools (live trees, standing dead, and downed wood) were collected preand post-treatment, then again a decade later. Species group-specific allometric equations were used to estimate live and standing dead biomass, and downed log biomass was estimated volumetrically. We used the Forest Vegetation Simulator to project “no-treatment” baselines specific to treatment units, allowing measured carbon responses to be normalized against differences in site characteristics affecting tree growth and pre-treatment stand structure. Results indicate that biomass development and carbon storage 10 yr post-treatment were greatest in SCE treatments compared to conventional treatments, with the greatest increases in coarse woody material (CWM) pools. Structural complexity enhancement treatments contained 12.67 Mg/ha carbon in CWM compared to 6.62 Mg/ha in conventional treatments and 8.84 Mg/ha in areas with no treatment. Percentage differences between post-treatment carbon and simulated/projected baseline values indicate that carbon pool values in SCE treatments returned closest to pre-harvest or untreated levels over conventional treatments. Total carbon storage in SCE aboveground pools was 15.90% less than that projected for no-treatment compared to 44.94% less in conventionally treated areas. Results from classification and regression tree models indicated treatment as the strongest predictor of aboveground C storage followed by site-specific variables, suggesting a strong influence of both on carbon pools. Structural enhancement treatments have the potential to increase carbon storage in managed northern hardwoods. They offer an alternative for sustainable management integrating carbon, associated climate change mitigation benefits, and late-successional forest structure and habitat.


INTRODUCTION
Forests globally function as a significant carbon sink, storing~45% of terrestrial carbon (Bonan 2008), yet are a leading source of carbon emissions due to deforestation and forest degradation (Keith et al. 2009). The world's forests have declined an estimated 36% in area (16.4 million km 2 ) over the last 200 yr (Meiyappan and Jain 2012). Although it is widely acknowledged that forest ecosystems of greater maturity and structural complexity maintain high levels of carbon storage (Harmon et al. 1990, Luyssaert et al. 2008, Keeton et al. 2011, Gunn et al. 2014, there is ongoing debate regarding the effects of different silvicultural approaches on carbon storage (Ruddell et al. 2007, Thomas et al. 2007, Nunery andKeeton 2010). Managing for old-growth forest structure has the potential to provide both carbon storage and late-successional habitat co-benefits, and is therefore of particular interest across the northern forest region of eastern North America (D'Amato et al. 2011). In this study, we quantified aboveground biomass and carbon storage in northern hardwood-conifer forests over a 10-yr period following an experimental silvicultural treatment, called structural complexity enhancement (SCE), and determined how this compared both to conventional selection harvesting systems and passive (no-harvest) management. This study also provides an opportunity to validate earlier projections (Keeton 2006) and other model simulations (Nunery and Keeton 2010) using empirical data.
Compared to old-growth and/or late-successional forests, young to mature secondary forests have lower densities of large live and dead trees (Keeton et al. 2007), less vertical and horizontal canopy complexity, and lower diversity of tree species, ages, and sizes (McGee et al. 1999, Franklin et al. 2002, McElhinny et al. 2005, D'Amato et al. 2011). Structurally complex temperate forests are known to maintain higher levels of biological diversity (Lindenmayer et al. 2000), hydrologic regulation (Wirth et al. 2009), and carbon storage (Keeton et al. 2011, Gunn et al. 2014, McGarvey et al. 2015. Over the past decade, several studies have investigated forest management practices specifically designed to promote late-successional/oldgrowth forest characteristics, both in the United States and internationally (Lindenmayer et al. 2000, Bauhus et al. 2009, Ducey et al. 2013, Duveneck et al. 2014. However, less well understood is whether these approaches would also have utility for carbon forestry, and thus, this question is the focus of our paper.
Carbon forestry often employs silvicultural prescriptions intended to improve growth and productivity and/or enhance aboveground carbon storage, the latter typically associated with high levels of stocking in mature trees (Nunery and Keeton 2010). Maintenance of belowground carbon storage is also important, especially given that temperate forests, for example, store 30-50% of their total carbon in soil pools (Pan et al. 2011). However, soil carbon responses to forest management have varied globally (Johnson and Curtis 2011), in some cases showing declines correlated with management intensity (Bucholz et al. 2014) and other times showing little response to management (Nave et al. 2011). Silvicultural treatments that enhance carbon storage while providing other co-benefits, such as latesuccessional biodiversity associated with stand structural complexity, are of interest both domestically and abroad (e.g., Gustafsson et al. 2012, Ducey et al. 2013, Chen et al. 2015. Prior studies (Angers et al. 2005, Bauhus et al. 2009, Puettmann et al. 2015 suggest that silvicultural treatments promoting or accelerating the development of late-successional structure may offer particular potential for this type of multi-functional forestry.
How to best manage forests for carbon storage is actively debated among researchers (Harmon 2001, Fahey et al. 2010, Birdsey and Pan 2015. Some have argued that more intensive management, for example, involving higher cutting frequencies and thus increased growth rates in younger trees and transfer of carbon to wood products, would yield greater levels of net carbon sequestration (Birdsey et al. 2006 andMalmsheimer et al. 2008). This strategy might also produce substitution effects, which are the avoided emissions from using wood products in place of materials requiring more energy, and thus emissions (Eriksson et al. 2007). Others propose reducing harvest frequency and intensity or passive management techniques Keeton 2010, Gunn et al. 2011), which favors development and maintenance of high levels of in situ carbon storage in extant forests (e.g., Krankina andHarmon 1994, Luyssaert et al. 2008). Extended rotations (Gronewald et al. 2010 andSilver et al. 2013) and post-harvest legacy tree retention (Gustafsson et al. 2012, Palik et al. 2014) are examples of this type of approach. While the relative effectiveness of carbon management strategies remains under debate (McKinley et al. 2011), it is clear that shorter rotations and intensified harvesting generally produce less complex stand structures, and if applied broadly enough, reduce the availability of late-successional habitats at the landscape scale (D'Amato et al. 2011, Gronewold et al. 2012, Puettmann et al. 2015. Less intensive management focused on promoting structural complexity may provide an alternative, contributing to climate mitigation through enhanced carbon storage while providing late-successional habitats (see, e.g., Smith et al. 2008, Dove and Keeton 2015, Kern et al. 2016. Carbon sequestration in the world's forests offsets 30% of global CO 2 emissions (Pan et al. 2011), and by some estimates, there is a potential to increase gross terrestrial C uptake by~2 Pg C annually (Birdsey and Pan 2015). Rapidly developing voluntary and compliance carbon markets seek to strengthen this sink capacity by providing a financial incentive for forest carbon projects that generate greenhouse gas emissions offsets (Russell-Roy et al. 2014, Kerchner and. This is encouraging broader adoption of forest management techniques that stock carbon across larger scales (Hoover and Heath 2011). Silvicultural approaches, such as SCE, that integrate carbon, timber, and late-successional biodiversity might have particular utility for generating forest carbon offsets attractive to buyers, for instance, in voluntary or "over-the-counter" markets, interested in co-benefits.
This paper adds to several previous investigations of responses to silvicultural approaches promoting late-successional forest characteristics (Angers et al. 2005, Dyer et al. 2010, Silver et al. 2013, although relatively few studies have explored complementarity between carbon storage and structural complexity objectives. Prior research on the SCE system, more specifically, has focused on harvest effects on stand structure (Keeton 2006), economic tradeoffs (Keeton and Troy 2006), and elements of late-successional biodiversity, including herpetofauna (McKenny et al. 2006), herbaceous plant communities (Smith et al. 2008), and fungal response (Dove and Keeton 2015). Here we report on aboveground carbon pools, explicitly addressing the question of how SCE, as compared to conventional selection systems and passive management, affects carbon storage and fluxes in aboveground pools in northern hardwood-conifer forests. Answering this question will help inform efforts to integrate carbon forestry and old-growth silviculture in both the northern forest region and beyond (e.g., Bauhus et al. 2009, Burrascano et al. 2013). In addition, the carbon responses to low-intensity harvesting systems, like SCE, are uncertain. High levels of retention under SCE might be predicted to negatively influence growth rates in residual trees, despite deliberate creation of small gaps and variable horizontal structure.
We hypothesize that aboveground carbon recovery 10 yr post-treatment will be greater under SCE compared to conventional selection treatments and relative to modeled carbon accumulation potential without treatment. Furthermore, we investigate whether multiple sources of variability influence carbon accumulation outcomes and interact with treatment effects.

Study area
The study areas for this project are located within the Mount Mansfield State Forest (MMSF, 44°30 0 23.03″ N; 72°50 0 11.24″ W) and the University of Vermont's Jericho Research Forest (JRF, 44°26 0 43.70″ N; 72°59 0 44.15″ W) in Vermont, United States (Fig. 1). Both areas fall within the central portion of the Green Mountains, a northern extension of the Appalachian Mountain Range. The MMSF study area spans elevations ranging from 470 to 660 m and is dominated by Peru stony loam soils; the sites at JRF are 200-250 m above sea level (a.s.l.), with soils classified as Adams and Windsor loams sands or sandy loams. Supplementary live tree inventory data from the Forest Ecosystem Research Demonstration Area (FERDA) experiment in New York, United States, were also used to complement existing data in this study. The FERDA experiment contains two study sites, Keese Mill and Visitor Interpretive Center (VIC; 44°25 0 59.6″ N; 74°20 0 36.4″ W), adjacent to Paul Smith's College in Franklin County, New York. Elevations at the FERDA sites ranged from 500 to 540 m a.s.l.; soils are Adams-Colton and Becket-Tunbridge-Skerry complex and are rocky and well drained.
Forests in the Vermont study areas are comprised of mature, 70-to 100-yr-old northern hardwood-conifer species. Dominant overstory species include Acer saccharum (sugar maple), Fagus grandifolia (American beech), Betula alleghaniensis (yellow birch), and Tsuga canadensis (eastern hemlock). There are minor components of Picea rubens (red spruce) at the MMSF study area and Acer rubrum (red maple) and Quercus rubra (red oak) in the dominant canopy at the JRF study area. Over the course of the 20th century, there were four to six recorded management entries (thinning and selection harvesting) in the study areas postestablishment (Hannah 1999), resulting in multi-aged forest structure confirmed through pretreatment tree coring as reported in Keeton (2006).
Dominant overstory species at the FERDA study area include A. saccharum, F. grandifolia, B. alleghaniensis, with minor components of A. rubrum and P. rubens. The FERDA sites were used extensively for agriculture until the early 1900s, after which they were abandoned and reverted to forest cover, and were partially harvested at least once in the mid-20th century.

Silvicultural treatments
This project employs a Before-After-Control-Impact experimental design (Krebs 1999), with structural metrics compared pre-and postharvest and between silvicultural treatments. Pre-treatment data were collected in 2001 and 2002; treatments were introduced to the MMSF and JRF study sites in 2003. This paper reports post-treatment data collected over 10 yr with the last re-measurement in 2013.
The three experimental manipulations included two conventional uneven-aged treatments (singletree and group selection) modified to enhance post-treatment structural retention (see Table 1 for treatment details), and a SCE treatment designed to accelerate the development of latesuccessional forest structure. Treatments were implemented across 2-ha units in a randomized block design, separated by a minimum 50-m buffer (see Fig. 1 for treatment unit layout). Each MMSF and JRF treatment unit contains five randomly placed permanent sample plots that are 0.10 ha in size; plots thus cover 25% of each treatment unit's total area. An important element of the SCE treatment was the target diameter distribution, which was based on a rotated sigmoid form (Goff andWest 1975, O'Hara 1998). In selecting this distribution over a negative exponential or "reverse J" form, the objective was to allocate more growing space and basal area to larger size classes, thereby promoting recruiting of larger trees and a related increase in biomass accumulation over time (see Keeton et al. 2011). The rotated sigmoid distribution was achieved by harvesting to three different q-factors applied to each of three portions of the target diameter distribution (see Table 1 for details). Target residual basal area, in this case a desired future condition, was set at 34 m 2 /ha and maximum diameter at breast height (dbh) to 90 cm, indicative of late-successional structure. Late-successional structure was further enhanced through crown release around larger trees, and silvicultural creation of coarse woody debris, small canopy gaps (0.02 ha mean size), standing dead trees, and tip-up mounds (see Keeton 2006). To ensure harvests were followed as prescribed, including that only stump-marked trees were felled and that all dead trees were retained in all the treatments, a certified forester monitored the harvesting operations daily.
Conventional uneven-aged treatments included single-tree and group selection harvests, which were based on a BDq prescription defining the form of the intended post-harvest diameter distribution. The prescription specified a residual basal areas (B) of 18.4 m 2 /ha, a maximum diameter at breast height (D) of 60 cm, a q factor (the ratio of the number of live stems among each successively larger 5 cm diameter class) of 1.3. These were compared to regional average target residual basal areas of 13.8-16.1 m 2 /ha, 40.6-45.7 cm max D, and q-factors of 1.5-1.7. Single-tree and group selection treatments had the same BDq prescriptions (Table 1), though applied in dispersed or aggregated patterns, respectively. Group selection cutting patches BDq is equal to the residual basal area (B), maximum target diameter (D), and q-factor (q). Adapted from Keeton (2006).
† The q-factor is equal to the ratio of the number of trees in each successively larger size class.
averaged 0.05 ha in size, with nine groups per treatment unit. Groups were well distributed but placed to release desired advanced regeneration; there was light retention of large dead trees and mature beech exhibiting resistance to bark disease (Nectria coccinea var. faginata) within some groups. Supplementary live tree data for conventional (single-tree and group selection) and control treatments were acquired from the FERDA Project. The project was initiated in 1998 and harvested in 2000. We used data from two FERDA sites, Keese Mill and VIC. Forest Ecosystem Research Demonstration Area treatment units are 2 ha with eight permanent plots per unit, each 0.04 in size. The FERDA single-tree and group selection treatments matched MMSF and JRF, with similar target post-harvest residual basal areas (18.4 m 2 /ha), BDq, and selection patch sizes (0.05 ha), and therefore were suitable as replicates. We selected FERDA replicates with pre-treatment basal areas most comparable to the JRF and MMSF study sites. There were additional treatment types tested in the FERDA project, data from which were not used in this study.
The conventional treatments (N = 4 per treatment) were replicated twice each at the MMSF and JRF study sites and twice at each FERDA site. The SCE treatment (N = 4) was replicated two times at MMSF and two times at JRF. There are two unharvested control units at MMSF, two at JRF, and two controls used from the FERDA Project (N = 6).

Field inventory
The field inventory data used in this study focused on measurements needed for aboveground biomass estimations. Within each plot (MMSF and JRF sites), we measured, identified, recorded, and permanently tagged all live and standing dead trees ≥5 cm dbh and ≥1.37 m tall. We recorded decay class (1-9) for all standing dead trees following Sollins et al. (1987) and measured standing dead heights using an Impulse 200 laser range finder (Laser Technology, Englewood, Colorado, USA). Downed log volume by decay class (1-5; Sollins et al. 1987) was estimated following the line-intercept method (Shiver and Borders 1996) for all downed logs ≥1 m in length and ≥10 cm diameter along two 31.62-m transects bisecting each plot at right angles through the center. Diameter at intercept, species, and decay class for each log along the center transects were recorded. For regeneration estimates, we tallied seedlings by species within a 1-m belt along each line-intercept transect. In the FERDA plots, all live and standing dead trees ≥2.54 cm dbh were inventoried and permanently tagged, but downed wood was not inventoried.

Data processing and analysis
Stand structural metrics.-We compared MMSF and JRF field inventory data collected in 2013 with inventory data from 2003 (first year posttreatment) and 2001 (pre-treatment and the year of project initiation) to assess differences in levels of carbon storage pre-and post-treatment and between treatment types. A comparable period pre-and post-treatment was used for the FERDA live tree data. We input all field inventory data (pre-and post-treatment) into the Northeast Ecosystem Management Decision Model (NED-3; Twery and Thomasa 2014) to generate stand structural metrics. These included live, dead, and total tree basal area, stem density, aboveground biomass, live tree quadratic mean diameter, and percent hardwood by basal area (Table 2). Slope and aspect were averaged for each treatment unit, and site index (using sugar maple as the focal species) was determined from pretreatment tree core and height information (MMSF and JRF) and from the Natural Resources Conservation Service Soil Survey (FERDA).
Biomass and carbon quantification.-To quantify live tree and standing dead carbon during each pre-and post-treatment inventory year, we first estimated live and dead tree biomass using the Jenkins et al.'s (2003) group-specific allometric equations embedded in NED-3. Live tree carbon was calculated by dividing the mean biomass for each treatment unit by two. Biomass calculations are the same for both live and dead trees in NED using the Jenkins et al.'s (2003) equations. Consequently, to determine standing dead tree biomass and carbon content, we made deductions to the allometrically derived estimates following the CARB (California Air Resources Board) carbon inventory protocol (Climate Action Reserve 2014). Adjustments reflected the amount of biomass missing (e.g., from breakage, decay, disease) from each dead tree when compared to its living counterpart. Deductions were determined by calculating the difference between the measured standing dead tree height and the pre-treatment inventoried live tree height for the same stem. To estimate former (when live) full height on snapped and dead trees, we used regression equations (R 2 ranged from 0.75 to 0.93) developed from our dataset for species-and site-specific diameter-height relationships. We then applied a density reduction factor following Harmon et al. (2011) correlating with measured decay class to all adjusted standing dead biomass values. Mean standing dead carbon for each unit was calculated using the final adjusted biomass values for individual standing dead trees divided by two. Downed log carbon content was determined following Harmon et al. (2008). Inventoried downed log volumes were adjusted by species specific gravities for each decay class (1-5). Adjusted volumes were then transformed to biomass and adjusted by carbon content by decay class. Species were assigned proportionate to the mean live tree basal area per treatment unit for all unknown species.
Carbon responses to treatments. -Carbon response trends were evaluated first using mean absolute values for each pool (live tree, standing dead, downed log) by treatment. For all carbon response comparisons, single-tree and group selection treatments were grouped into a "conventional" uneven-aged treatment following Keeton (2006). There were no statistically significant differences in stand-level structural outcomes between the selection treatment types, supporting this grouping (see Table 2).
For the second carbon quantification assessing post-harvest carbon storage in each treatment relative to untreated or "baseline" conditions within each unit, we calculated percentage differences between post-treatment and baseline carbon values for all measured pools within treatments. This determined how near to the untreated or "baseline" condition each treatment returned 10 yr post-harvest. We chose the percentage difference metric as a standardized comparison normalizing relative difference between harvested and baseline values across the range of inherent site variability. Percentage differences were calculated following Littlefield and Keeton (2012), modified from Westerling et al. (2006). Percentage differences were calculated as follows: where V T is equal to a post-treatment carbon value, and V B is equal to a baseline carbon value (see next section for "baseline" carbon definitions). Using the above formula, we compared carbon storage in each pool 10 yr post-treatment against baseline values specific to each treatment unit. In this analysis, a zero (0) percentage difference indicates no difference from the baseline conditions; a negative (À) percentage difference indicates below the baseline; and a positive (+) percentage difference indicates above or surpassing the baseline conditions. Therefore, for posttreatment carbon pool comparisons, percentage differences that are closer to or above 0 indicate greater C storage potential relative to a "no-treatment" baseline. Carbon flux for each pool was defined as the amount of C lost or gained over the 10-yr interval post-treatment (MgÁha À1 Áyr À1 ; Harmon 2001. Carbon flux was calculated by determining the difference between mean carbon the year immediately post-treatment and mean carbon 10 yr post-treatment and dividing by 10 (the time span of comparison).
Modeling baseline carbon dynamics under a "notreatment" scenario. -In testing the hypothesis, it was important to control for inherent differences in site quality, initial stocking, species composition, and other factors in order to assess the carbon development potential under a "no-treatment" baseline, which was then compared against the measured carbon responses. In other words, we needed to normalize the measured (i.e., empirical) carbon responses against the inherent tree growth and carbon accumulation potential specific to each unit. Therefore, we used growth and yield modeling to simulate 10 yr of growth based on pre-treatment data and assuming no treatment, disturbance, or management. The northeastern variant of the Forest Vegetation Simulator (NE-FVS) was selected for this purpose because of its wide use in a variety of forest management (Crookston and Dixon 2005) and carbon offset applications (Kerchner and Keeton 2015). Northeastern variant of the Forest Vegetation Simulator is a distance-independent and individual treebased growth and yield model suitable for both even-and uneven-aged stands. Regional validation studies of NE-FVS have shown accurate volume and biomass projections in northern hardwood forests, within 10-15% when simulating forest growth (Yaussy 2000). However, a limitation is that FVS has been shown to have inaccuracies estimating large, live tree growth in northeastern U.S. late-successional and old-growth forests (Gunn et al. 2014, MacLean et al. 2014. In our study, this limitation is acceptable in that the resulting growth projections are conservative, especially for only a 10-yr time interval. Our FVS projections allowed us to test our first hypothesis regarding treatment effects. Stand-level growth simulations in FVS are known to be sensitive to tree regeneration inputs (Ray et al. 2009). Therefore, we evaluated sensitivity in our projections by modeling three different regeneration input scenarios: inventoried pre-harvest regeneration densities, adjusted inventoried regeneration densities (reduced by one order of magnitude), and no regeneration (Table 3). There was a <1-25% range of variability in growth projections between the different regeneration scenarios. We ran a second set of 10-yr FVS simulations using post-treatment data from 2003, modeling each of the three regeneration scenarios. Modeled regeneration scenarios for all the units were compared against measured values for those same units. We chose the no-regeneration scenario for further modeling because it had the closest fit to the measured values (mean = 7% difference). The largest discrepancy between measured and projected values was Jericho Unit 3, where a fine-scale wind disturbance occurred over the sampling interval and further reduced post-treatment basal area.
The no-management baseline for coarse woody material (CWM; standing dead and downed wood) carbon pools was assumed to be equivalent to the pre-treatment values; changes in the CWM pools were also compared against the controls, as was the live tree pool. We did not model CWM development because recruitment into this pool did not change significantly (see Results) in the control units over the 10-yr time interval except for a small amount due to windthrow. This was consistent with previous research showing minimal CWM recruitment from individual tree mortality over this timeframe in eastern deciduous forests (Woodall 2010. Statistical analyses.-To explore the effects of each treatment on carbon response, we tested for statistically significant differences in carbon responses between treatments by pool, and compared the empirical values 10 yr post-treatment to the no-management baseline. For this purpose, we employed one-way ANOVA and post hoc ❖ www.esajournals.org Tukey-Kramer multiple comparisons. Statistical comparisons of means by treatment and pool were made in JMP Pro 11 (SAS Institute Inc. 2013). Shapiro-Wilk tests for normality confirmed normal distribution of data (a = 0.05) and one-way ANOVAs and Tukey-Kramer honestly significant difference post hoc mean comparisons tested for significant differences in carbon pool means pre-and post-treatment. Homogeneity of variance was tested using F tests.
To determine whether variables other than treatment had a significant influence on carbon accumulation, we evaluated the relative influence of multiple independent variables (e.g., treatment type, site productivity, location, and other site characteristics) on the dependent variables (percentage difference carbon per pool). This consisted of classification and regression tree (CART) analyses (Breiman et al. 1984) conducted in S-plus 8.2 (TIBCO Software 2010). Classification and regression tree is a robust nonparametric technique that accommodates both categorical and continuous variables (Littlefield and Keeton 2012). A tree hierarchically ranks the predictive power of multiple independent variables by repeatedly splitting dependent variables into more homogenous groups based on combinations of independent or explanatory variables, which can explain variation within partitioned values of the dependent variable (De'ath and Fabricius 2000). We selected CART as the preferred multivariate method because variance (or deviance) is partitioned into increasingly finer resolution, rather than being assessed across all data points collectively as with certain other methods, such as general linear models. This helps identify secondary predictor variables that may be more influential when paired with a particular combination of variables (and for a given subset of data points), but less so when other variables are operative (or for a different subset of data points). We used a robust set of predictor variables (n/2) representative of site variability (percent hardwood, slope, aspect, location; Table 4). Cost-complexity pruning was used to remove insignificant nodes (a = 0.05).

Carbon responses post-treatment
Our results support the hypothesis that SCE carbon recovery 10 yr post-treatment would be greatest relative to pre-treatment or no-management baseline values when compared to conventional treatments. Comparisons of treatments pre-to postharvest indicate greatest amounts of biomass development (carbon storage; Fig. 2) and greatest carbon fluxes (Table 5) in SCE treatments as compared to  Fig. 3) show the greatest increases occurred in post-harvest SCE carbon relative to pre-harvest levels.
Mean SCE standing dead and downed log carbon post-harvest (2013) were greater in SCE treatment units than in conventional and controls (Fig. 2). Live tree and total C were significantly greater in controls than in conventional units 10 yr post-harvest (P = 0.004) and were also greater in SCE units than in conventional units (see Table 6). Limited significance in C differences between treatments for CWM after 10 yr can be attributed to gradual declines in the SCE units from attrition and decay, as well as some downed log inputs from a few wind-thrown trees confined to the control units at MMSF. Relative to conventional treatments, SCE carbon fluxes were either greater or comparable in live tree, standing dead, and downed log pools, which also supported our hypothesis that C fluxes would be greatest in SCE-treated areas (Table 5). Structural complexity enhancement downed log pools demonstrated the greatest difference in fluxes compared to other treatments, measuring À0.72 MgÁha À1 Áyr À1 compared to À0.33 and À0.05 MgÁha À1 Áyr À1 in conventional and control treatments, respectively (P < 0.05).
Percentage differences between post-treatment and no-management baseline carbon for each pool were lowest in SCE treatments relative to conventional treatments, supporting our hypothesis that SCE treatments will result in C levels closest to untreated or manipulated stand development (Table 6, Fig. 3). Percentage differences were below the baseline (negative) in all treatments in the live tree pool, measuring À17.45% in SCE units and À42.81% in conventional units. Standing dead SCE carbon was again closest to the baseline, with a À65.70% difference compared to À90.20% in conventional units. Downed log carbon under SCE surpassed the baseline with a measured 32.86% increase, yet in conventional units was still below, with a 19.36% decrease.

Effects of site variability on biomass development and carbon storage
It is evident from our CART results that treatment type was most influential on carbon storage across all pools, but that variability in site conditions interacted with treatment in determining carbon response in most situations. Treatment type explained variations in percentage differences in carbon storage at the first and sometimes secondary splits of all trees (Fig. 4). Individual pools demonstrated different responses in carbon storage due to variations in site conditions, demonstrated by relative ranking of secondary predictor variables. Five predictor variables were selected for the final CART models: treatment, aspect, slope, site index, and percent hardwood (percent of hardwood basal area; Table 4).
The live tree carbon CART model (Fig. 4A) primary split (most influential predictive variable) was split between conventional treatments and SCE and control treatments. Carbon storage potential increased moving from conventional treatments to control and SCE treatments. Conventionally treated areas were additionally  influenced by aspect (at the secondary split), with the percentage difference carbon increasing in negative value (less post-treatment carbon relative to baseline carbon) as site orientation moved to the northwest. Percent hardwood and treatment were selected as partitioning points for variance within controls and SCE treatments, with percentage difference increasing in positive value (greater carbon storage) with increase in percent hardwood and from SCE to control treatments.
In both CWM models (Fig. 4B, C), there is a stronger influence of site variation on carbon storage, as indicated by the selection of slope, site index, percent hardwood, and location as predictor variables in the final models. Treatment was selected as the most important predictor of percentage difference between post-harvest and baseline carbon in both models. In the standing dead model (Fig. 4B), the primary node was split between SCE and conventional treatments and controls. Carbon storage potential was greatest in controls and lowest in SCE and conventional treatments. Slope was selected as a secondary predictor for SCE and conventional treatments, with slopes less than 29°having less post-treatment than baseline carbon, or a greater negative percentage difference. As slopes increased, percentage differences in CWM carbon decreased. Location was also ranked as secondary predictor Table 5. Mean annual C flux per pool and treatment over the 10-yr interval post-treatment and significance levels (a = 0.05).  Notes: HSD, honestly significant difference; SCE, structural complexity enhancement; Conv., conventional. The percent differences in values compare the two following Littlefield and Keeton (2012). One-way ANOVA and Tukey-Kramer post hoc analysis results are also presented (a = 0.05). Degrees of freedom = 2. Significant results are in bold (P < 0.05).
variable for standing dead carbon in control treatments, with a greater percentage difference at the MMSF site (higher amounts of standing dead carbon post-harvest than pre-harvest). Percentage difference for downed log carbon was greatest in SCE treatments, with the CART model (Fig. 4C) split between control and conventional treatments and SCE treatments. A secondary predictor of SCE-treated sites was site index, with the percentage difference for carbon increasing with decreasing site productivity. Percent hardwood and treatment were selected as partitioning points for variance among control and conventional treatments.

DISCUSSION
Carbon stocking in aboveground biomass pools in northern hardwood forests has the potential to increase with silvicultural prescriptions that retain elements of stand structure, increase horizontal and vertical complexity, and elevate the availability of CWM. Of the treatments tested in this study, the SCE treatment resulted in aboveground carbon storage levels  closest to untreated controls and modeled nomanagement scenarios. Additionally, after a decade, this treatment maintained and developed greater amounts of carbon storage than the other selection systems tested, likely as a result of elevated post-treatment structural retention and other techniques employed in the SCE treatment (Table 1). We also found site variability to have an important secondary effect on the amount and rate of carbon accumulation in each pool, with carbon storage potential generally increasing with site conditions favoring better growth response to silvicultural treatment, as indicated by our CART models.

Carbon responses to management for old-growth characteristics
Pre-and post-treatment measured carbon outcomes. -Absolute carbon values 10 yr post-treatment were greater in SCE units in all pools relative to conventional treatments. This can be attributed, in part, to a higher post-treatment target residual basal area in SCE units and also to silviculturally created CWM inputs. The absolute carbon values for SCE were comparable to or above published regional values for C stocks and, for some pools, were close to regional averages reported for oldgrowth/late-successional forests. For mature northern hardwoods, the USDA Forest Service (2015) recently reported mean live tree C storage between 60 and 80 Mg/ha, standing dead C between 2 and 4 Mg/ha, and downed log C between 6 and 9 Mg/ha. Other studies specific to northern New England report comparable values (Keeton et al. 2011, Gunn et al. 2014. Bradford et al. (2010), studying mature hardwood forests in northern New Hampshire (maximum age of 120 yr), found 96 Mg/ha of live tree C and 18 Mg/ha C in CWM. Whitman and Hagan (2007) reported higher levels for the same forest type in Maine: 113 Mg/ha live tree C, 10 Mg/ha standing dead C, and 12 Mg/ha downed log C. Aboveground carbon amounts in our experimental units were comparable to regional values, with SCE live tree C measuring 98.22 Mg/ha, standing dead at 3.67 Mg/ha, and downed log at 9.00 Mg/ha. Regional old-growth northern hardwood C stocking has been reported at 116-141 Mg/ha live tree C, 8-22 Mg/ha standing dead C, and 12-18 Mg/ha downed log (Goodburn and Lorimer 1998, Fisk et al. 2002, Whitman and Hagan 2007, Bradford et al. 2010, Keeton et al. 2011, Gunn et al. 2014, McGarvey et al. 2015. Carbon stocking in SCE treatments 10 yr post-treatment was at the upper threshold or above regional mean values, and in some cases approaching regional old-growth stocking levels, indicating the effectiveness of this treatment type in promoting late-successional/old-growth C stocking levels and structure. We found FVS to be a useful model for making short-term forest growth projections and comparing predicted values against empirical results, allowing us to more fully evaluate the effects of each treatment type on C stocking and fluxes relative to site-specific potentials. However, the sensitivity of FVS to regeneration inputs should be noted when using this model for short-term growth projections. Regeneration inputs in shortterm projections seem likely to result in additional variability in outputs because of model behavior, which tends to induce pronounced density-dependent mortality at time steps immediately after user-specified regeneration inputs. The model may therefore function more accurately without the addition of regeneration for short-term (10-yr) forest growth projections.
Management vs. no-management effects on carbon accumulation.-When comparing measured carbon outcomes after treatment with modeled notreatment baselines, SCE percentage differences for all pools were closest to or above the notreatment baseline relative to conventional treatments. This is consistent with the literature predicting accelerated biomass development and recovery of late-successional characteristics following treatments similar to SCE (Bauhus et al. 2009). Management scenarios involving no treatment have consistently shown the greatest total long-term carbon storage, accounting for both in situ forest carbon and the life cycle of wood products (Harmon 2001, Fahey et al. 2010, Nunery and Keeton 2010. However, in our study, the contrast with no management was lowest across all carbon pools under SCE (17.45% less than baseline for live tree C) as compared to the conventional treatments (42.81% less than baseline for live tree C). This finding suggests great potential for low-intensity silvicultural techniques favoring in situ carbon storage in the northern hardwood region, assuming regeneration and other management objectives are met (Gottesman and Keeton 2017), which of course will vary tremendously by site and ownership (Schwenk et al. 2012). The utility of management for in situ C storage must, of course, be considered within the larger context of carbon forestry, which includes approaches emphasizing carbon storage in wood products (FAO 2016) and avoided emissions from substitution of wood products for more energy-intensive materials (Malmsheimer et al. 2008).
The response of downed coarse woody material in this study was particularly promising toward the integration of management for late-successional habitats with carbon storage. Ten years post-treatment, downed log carbon under SCE was significantly higher than the no-treatment baseline and the control units. In addition to providing important habitat for a variety of late-successional organisms (McGee et al. 1999, McKenny et al. 2006, Dove and Keeton 2015 and riparian functions (Keeton et al. 2007, Warren et al. 2009), our results suggest structural complexity approaches have the potential to store significant amounts of carbon in downed woody material, as well.
Carbon flux variations by pool and treatment.-Carbon fluxes were greatest in the live tree and downed log pools following the SCE treatment. These results indicate both a high level of C sequestration (uptake) from accelerated tree growth in response to harvest, and C loss through decay. The latter is likely due to the large input of silviculturally created CWM. Our CWM flux rates are comparable to regional published estimates (Bradford et al. 2010, Gunn et al. 2014). However, we note the difficulty in accurately measuring CWM flux rates due to the combined effects of density, volume, and/or biomass depletion in addition to losses from heterotrophic respiration (Forrester et al. 2015). Additionally, CWM flux is usually greatest within the first 10 yr post-treatment. Live tree and total C flux rates for all treatments were greater than Keeton 2010, Gunn et al. 2014) or comparable to regional estimates (Bradford et al. 2010), with SCE live tree flux measuring higher than conventional treatments. While the conventional treatments also showed elevated levels of tree growth, we found that SCE achieved similar or greater growth responses in overstory trees. This is an important finding relative to the potential for low-intensity treatments of this type both to maintain complex stand structures and to elevate carbon sequestration (see, e.g., Bauhus et al. 2009). While our study does not provide a basis for determining a mechanism for the elevated uptake rates, it is possible this was due to crown release of dominant trees as well as variable canopy openness (or gapiness), both of which were explicit objectives of SCE.

Comparisons of empirical data with modeled forest development projections
This study provided a unique opportunity to compare empirical data with previous simulations and projected outcomes. Using data from one year post-treatment, Keeton (2006) projected that after 50 yr aboveground biomass development in the SCE treatments would be 91.4% of that projected for no treatment. Conventional treatment units were projected to achieve 79.1% of their no-treatment potential on average. Using empirical data from 10 yr post-treatment, carbon accumulation in SCE units is already 84.1% on average of the level simulations project would develop over 50 yr without treatment. This compares to the carbon measured in conventional treatment units, which has achieved 55% of the simulated no-treatment potential for the 50-yr timeframe in Keeton (2006). It is evident from these results that FVS significantly underestimated biomass development in the Keeton (2006) projections. This is consistent with the findings of MacLean et al. (2014) who found that uncalibrated regional variants of FVS tended to under-predict carbon for FIA plots across the northeastern United States. Our findings are contrary to Gunn et al. (2014), however, who found FVS to over-estimate carbon stocks in both latesuccessional and old-growth northeastern forests. Gross carbon under-estimation in FVS may be a result of a lack of late-successional and old-growth data available for model calibration (Gunn et al. 2014). Further research is recommended to reconcile these differences. Finally, total post-treatment aboveground carbon for conventional treatments measured in this study was nearly equal to (<1% difference) comparable treatments projected by Nunery and Keeton (2010). Structural complexity enhancement total carbon measured 10 yr posttreatment was only 7.89% below that of 10-yr projections for no-management scenarios modeled by Nunery and Keeton (2010). These differences in FVS-projected outcomes for northeastern tree growth and C stocking suggest a need to improve model calibration and accuracy, particularly given the wide acceptance of FVS by forest carbon markets (Kerchner and Keeton 2015). Additionally, the effects of natural disturbances, invasive species (e.g., beech bark disease), and climate change (e.g., changes in species distributions) on forest stand development also need to be considered when projecting future forest conditions (Seidl et al. 2014).

Site variability effects on carbon storage potential
Our CART models were consistent in showing treatment type to have the greatest influence on carbon stocking in all aboveground pools measured. Model results also demonstrated the effect of site variability on C in pools. Disparity in carbon storage potential in all pools as a result of differences in site conditions was explained in CART results (Fig. 4), suggesting a relationship between alterations in site conditions and biomass development/carbon storage similar to those described by Nunery and Keeton (2010) and Littlefield and Keeton (2012). It can generally be assumed that C storage potential was directly affected by aspects of site variability, such as percent hardwood, productivity, and slope, which was most clearly evident in the live tree and standing dead CART models (Fig. 4A, B).

Implications for adaptive forest carbon management integrating multiple co-benefits
Multi-functional forest management practices promoting the development of stand structural complexity and associated late-successional habitat characteristics (Keeton 2006, Bauhus et al. 2009) are likely to provide important carbon storage co-benefits based on the results of this study. Disturbance-based management (Seymour et al. 2002, Franklin et al. 2007) promoting legacy tree retention, inputs to coarse woody debris pools, increased vertical and horizontal heterogeneity, and elevated biomass levels are options for maximizing C storage potential (Franklin andPelt 2004, Gustafsson et al. 2012). These provide important co-benefits in terms of habitat function and biodiversity conservation targeted at the full array of temperate forest species, including those associated with late-successional habitats (Lindenmayer et al. 2000, Keith et al. 2009). The re-allocation of diameter distributions to larger size classes supports the growth of large trees, an important element of late-successional forest structure. Recent research (Stephenson et al. 2014) analyzing 403 tropical and temperate trees species indicates that tree growth rate increases continuously with size, as does C sequestration and storage for most trees. Large trees, previously assumed to slow in both productivity and growth rate (Weiner andThomas 2001, Meinzer et al. 2011), function as long-term carbon sinks (Carey et al. 2001). These findings further support the significance of structural retention as a co-benefit to forest carbon storage.
Adaptive silvicultural practices promoting multiple co-benefits, for instance, by integrating carbon with production of harvestable commodities, can contribute to efforts to dampen the intensity of future climate change while maintaining resilient ecosystems (Millar et al. 2007). Prescriptions that enhance in situ forest biomass and thus carbon storage offer one such alternative (Ducey et al. 2013). U.S. forests currently offset approximately 16% of the nation's anthropogenic CO 2 emissions, but this has the potential to decline as a result of land-use conversion and lack of management (EPA 2012, Joyce et al. 2014). While passive or low-intensity management options have been found to yield the greatest carbon storage benefit, assuming no inclusion of substitution effects (Nunery and Keeton 2010) or elevated disturbance risks (Hurteau et al. 2016), we suggest the consideration of SCE to enhance carbon storage. Multiple studies have explored co-benefits provided by management for or retention of elements of stand structural complexity, including residual large living and dead trees, horizontal variability, and downed CWM (Angers et al. 2005, Schwartz et al. 2005, Dyer et al. 2010, Gronewold et al. 2012, Chen et al. 2015. Silvicultural treatments can effectively integrate both carbon and late-successional biodiversity objectives through SCE based on this study and previous research (e.g., Dove and Keeton 2015). Remaining cognizant of the potential for old-growth compositional and structural baselines to shift over time and space with global change-climate impacts on forest growth and disturbance regimes, altered species ranges, and the effects of invasive species-will be important for adaptive management for late-successional functions such as carbon storage.