Carbon stocks and accumulation rates in Pacific Northwest forests: role of stand age, plant community, and productivity

Forest ecosystems are removing significant amounts of carbon from the atmosphere. Both abiotic resource availability and biotic interactions during forest succession affect C accumulation rates and maximum C stocks. However, the timing and controls on the peak and decline in C accumulation rates as stands age, trees increase in size, and canopy gaps become prevalent are not well-understood. Our study examines measured change in live and dead woody C pools from 8767 inventory plots on 9.1 million ha of Pacific Northwest National Forest lands to determine how the balance of tree growth, mortality, and dead wood decomposition varied by stand age, plant community type, and site productivity; and to compare the contribution of different tree sizes to C accumulation. Maximum non-mineral soil C for old-growth stands varied significantly by productivity class within plant community types, but on average stands accumulated 75% of maximum stocks by age 127 ± 35 yr. We did not see a decline in net primary production of wood (NPPw) with age in moderate and low-productivity classes, but found a 33% reduction in high-productivity classes. Mortality increased with stand age such that net change in live tree biomass, and change in total woody C, was not significantly different from zero in old-growth stands over age 400 (0.15 ± 0.64 Mg C·ha−1·yr−1 for woody C). However, significant though modest C accumulation was found in forests 200–400 yr old (0.34–0.70 Mg C·ha−1·yr−1, depending on age class). Mortality of trees >100 cm diameter exceeded or equaled NPPw, but trees were growing into the larger sizes at a high-enough rate that a net increase in large tree C was seen across the region. Although large trees accumulated C at a faster rate than small trees on an individual basis, their contribution to C accumulation rates was smaller on an area basis, and their importance relative to small trees declined in older stands compared to younger stands. In contrast to recent syntheses, our results suggest that old-growth and large trees are important C stocks, but they play a minor role in additional C accumulation.


Introduction
Forest ecosystems play a major role in the global carbon cycle because they can attain high levels of carbon storage, and can gain or lose carbon relatively rapidly (McKinley et al. 2011). These carbon (C) stocks are present in aboveand below-ground live, dead, and highly decayed vegetation components that cycle with the atmosphere on time scales of days to centuries.
Carbon stocks and accumulation rates in Pacific Northwest forests: role of stand age, plant community, and productivity v www.esajournals.org

GRAY ET AL.
Although approximately one-quarter of the C in the atmosphere is fixed by photosynthesis in a year, relatively little ends up in carbon pools with long residence times such as trees or soil (Falkowski et al. 2000). Understanding the magnitude and drivers of C flux between forests and the atmosphere has been a focus of research given concerns about the effects of rising levels of atmospheric carbon dioxide on climate change (IPCC Core Writing Team 2007). Much of the terrestrial annual carbon accumulation in recent decades is believed to be occurring in northern hemisphere forests, but uncertainty remains about where and by how much (Pan et al. 2011, Hayes et al. 2012).
In addition to combustion and harvest, the rate at which different forest types store and release C is determined by available resources and environmental conditions, which control net primary production (NPP) and heterotrophic respiration (Rh). The species composition of different forests may also determine the maximum stocks attainable on a site due to species' differences in longevity, growth, and decomposition rates Franklin 1979, Harmon et al. 1986). Predicting relevant sitespecific resource levels given wide variation in climate, topography, and soils is challenging (e.g., Coops et al. 2012). In practice, forest managers often use empirical relationships with site index and community classifications, which integrate conditions at individual locations, to predict forest productivity and evaluate alternative management actions (e.g., Hemstrom et al. 1987, Barnes et al. 1998, Hoover and Rebain 2011. This approach is a useful framework for assessing C flux. The net rate of C accumulation also changes with forest age and successional stage, depending on relative rates of NPP vs. Rh. A net loss of C usually occurs in the initial period of postdisturbance secondary succession, dominated by the decomposition of remnants of the previous stand (Janisch and Harmon 2002). As new forest vegetation dominates the site, C accumulation tends to peak early in stand development as NPP reaches a maximum soon after tree canopy closure, and then declines as stands age due to declines in NPP and increases in decomposition of dead matter produced by the new stand (Rh) (Ryan et al. 2004, McKinley et al. 2011. The hypothesized causes of decline in NPP in older stands include: increased autotrophic respiration as live biomass accumulates (but see Ryan et al. 2004); changes in allocation of carbon to less quantifiable below-ground structures, or inability of subordinate trees in multilayered old stands to use resources as efficiently as overstory trees (Binkley et al. 2002); or inability of large trees to obtain and/or distribute water and nutrients to the foliage in large, tall crowns (McDowell et al. 2002). Despite their low net C accumulation rates, old-growth forests provide useful estimates of maximum C storage and reference points for evaluating the effects of human and natural disturbances on regional C stocks (Smithwick et al. 2002).
Old-growth forests store large amounts of C per unit area, but a recent study suggests that substantial rates of C accumulation may be more common than previously thought (Luyssaert et al. 2008). However, mortality of large trees in old-growth stands may be increasing in the western United States (van Mantgem et al. 2009), with the potential to reduce NPP. Tree mortality resulting in canopy gaps in older forests may either lead to lower rates of C accumulation as space is not fully occupied by photosynthesizing plants (Coomes et al. 2012), or to greater rates because live biomass and NPP recovers quickly, while dead biomass decomposes slowly (Luyssaert et al. 2008). Decomposition rates of standing dead trees and downed wood differ, so the timing of gap events and snag fall rates affect stand-level Rh. Several of these hypotheses about NPP and stand age focus on absolute tree size as well as within-stand relative position (i.e., dominant or subordinate). Despite small increments in diameter, individual large trees accumulate more mass annually than small trees (Stephenson et al. 2014). However, large trees are typically less numerous than small trees and take up more space within a stand, so C accumulation on an area basis may not be that different.
Most of the empirical research on forest C stocks and accumulation rates is based on chronosequences, experiments, long-term measurements of a few plots, or one-time inventories. Many of these studies are conducted in archetypal stands selected to minimize confounding factors, which usually means they are fully stocked and have not been disturbed recently (Botkin and Simpson 1990). The focus of forestry research has often been on the more productive v www.esajournals.org GRAY ET AL.
(i.e., economically valuable) vegetation types in a region, and many of the experimental studies on NPP and succession have been done in relatively simple or single-species stands (e.g., Ryan et al. 2004). As a result, it is not clear how well the patterns and processes found by those studies apply to the larger population of forests in a region that are rarely fully stocked, occur on a wide variety of edaphic conditions, and have likely experienced a number of management, biotic, and abiotic mortality events during their development.
In addition, few landscape and regional assessments of C flux have been based on measured, as opposed to modeled, changes in live and dead woody C pools.
Our objectives in this study were to assess the role of stand age, plant community type, and productivity on forest C stocks (excluding organic C in mineral soil) and the net changes in woody C over a diverse range of forest conditions. We analyzed a systematic inventory of National Forests in the Pacific Northwest, United States, with repeat measurements of most aboveground C pools. Previous work with these data summarized regional patterns of C stocks and change aggregated to landscape scales (Gray and Whittier 2014), while this study is focused on patterns and processes related to stand-level attributes. Specific null hypotheses include: (1) The rate at which stands accumulate C stocks, and the balance of tree growth, mortality, and decomposition, do not vary by plant community type and site productivity; (2) Most carbon accumulation occurs in young stands, and net accumulation in oldgrowth stands is essentially zero; and (3) Small trees contribute as much to the change in live C stocks within stands as do large trees. We focus on the dynamics of relatively undisturbed stands (i.e., no logging or fire between measurements), but also provide some comparisons with all stands in the region combined to determine the net effects of disturbance.

Study area
We assessed C stocks and their change on the 9.1 million ha of forested federal land administered by the Pacific Northwest (PNW) Region of the National Forest System (NFS). These lands are found primarily in the states of Oregon and Washington as well as parts of California and Idaho, USA, between 41.8 ° and 49.0° N latitude and 116.3 ° and 124.7° W longitude. NFS lands in this region occur in a great variety of conditions, with annual precipitation ranging from 25 to over 350 cm, mean annual temperature from −1 to 12 °C, and elevations from 0 to 3300 m above sea level (Franklin and Dyrness 1973). We grouped individual plots into 10 plant association zones (PAZs ; Table 1) designated by the climax tree species as classified by field crews using local NFS guides (Hall 1998), and assessed geographic outliers by examining tree records and overlaying plot locations with a model of PAZs (Henderson 2009). Twenty-three percent of NFS forest land was "reserved" (i.e., where management for production of wood products was precluded; 85% of this was designated Wilderness) with the majority found in the ABLA, TSMEpark, and ABAM zones (Abies lasiocarpa, Tsuga mertensiana and alpine parklands, and Abies amabilis zones, respectively; Table 1).
The data and compilation methods we used for this study are similar to those used in Gray and Whittier (2014). We summarize the methods that were the same and provide full details on calculations and variables not used in that paper.

Field data
The primary data used in this study were collected by the US Forest Service for a strategic inventory of vegetation conditions on all NFS lands in the PNW Region using a probabilitybased sample design (Max et al. 1996). The sample consisted of a systematic square grid at a 5.47 km spacing across all lands, and a denser grid at a 2.74 km spacing outside of designated Wilderness areas, providing a sample density of one plot per 3000 and 750 ha, respectively. Plots were installed using the Current Vegetation Survey (CVS) design (Max et al. 1996(Max et al. ) between 1993(Max et al. and 1997 and remeasured between 1997 and 2007 ("time 2") in four spatially-and temporally balanced panels. The CVS plot remeasurement period ranged from 1 to 14 yr with a mean of 7.1 yr. To avoid high sample errors associated with estimating annual rates of change from short remeasurement periods on small numbers of plots, we only used plots from the last three panels, which were remeasured more than 2 yr after installation. The same grid of points was also measured with the nationally standardized Forest Inventory and Analysis (FIA) design starting in 2001 (USDA Forest Service 2006); we applied the FIA land classification distinguishing forest from non-forest, and FIA measurements of organic soil horizons (i.e., duff and litter depths measured at eight points per plot), to the data in this study.
The CVS plot design consisted of a cluster of five points within a 1-ha circle, with nested circular plots of different sizes used to measure live and standing dead trees of different sizes. Crews sampled downed wood with the line-intercept method and estimated line-intersect cover of shrub and forb vegetation on a 15.6 m transect at each point. Changes in understory biomass appeared unreliable due to changes in vegetation protocol (Gray and Whittier 2014); we therefore do not report estimates of change in understory vegetation C.

Data classification and compilation
There were 8767 plots within NFS lands that had forested conditions measured 3 or more years apart. Forest land is defined as land areas ≥0.4 ha that support, or previously supported, ≥10% canopy cover of trees and were not primarily managed for a non-forest land use. Systematically placed plots often straddle different conditions; we only used portions of plots classified as forest in our analysis. We assigned values for stand age, site index, and forest type based on the FIA compilation of the FIA sample of the same plot (e.g., Donnegan et al. 2008) where available, or from an FIA compilation of the first CVS plot measurement (Waddell and Hiserote 2005). We estimated site productivity in terms of production of wood at culmination of mean annual increment (MAI, m 3 ·ha −1 ·yr −1 ) from the site index tree measurements (Hanson et al. 2002). We identified the nature and year of human and natural disturbances on each plot from a combination of field records and spatial data layers. Most analyses in this paper were based on a set of "undisturbed" plots that consisted of plots where tree cutting or fire had not occurred between measurements (N = 7647).
We used a combination of CVS status codes, estimated growth rates, and disturbance information to validate changes in individual tree status, calculate growth rates, and estimate diameters and heights of trees when they died (N = 1 008 943 tree records). Estimates of aboveand below-ground live tree and standing dead tree woody C used the procedures documented in Woodall et al. (2011) and tree foliage and coarse-root ratios from Jenkins et al. (2003). Above-and below-ground biomass estimates v www.esajournals.org GRAY ET AL.
for standing dead trees were reduced to account for decay and proportional bark and branch loss from Harmon et al. (2011). A standard "trees per hectare" expansion factor derived from the appropriate fixed-area plot size was used to convert individual live tree and standing dead tree C to an area basis (Mg/ha), and to calculate ratios of means for selected tree attributes (e.g., growth per tree or per unit basal area). We calculated C in downed dead woody material (DWM) from line-intercept diameter (for pieces >7.6 cm diameter) and counts by diameter class (for smaller pieces) using the procedures in Woodall and Monleon (2008) and densityreduction constants by decay class from Harmon et al. (2011). Forest-floor C was calculated from the annual FIA duff and litter depth measurements using forest type-specific bulk densities (Woodall and Monleon 2008). Additional details on C calculations are found in Gray and Whittier (2014).
We calculated three components of change for live tree C. "NPP w " consisted of the woody mass increment of trees alive at both measurements, the increment of trees up to the estimated size and year of mortality or harvest, and the size of trees that grew over the minimum 2.5 cm DBH threshold during the remeasurement interval (equivalent to the forestry term "gross growth"). Differences in estimates of NPP w among stands should reflect differences in NPP, though we did not attempt to estimate turnover in foliage and fine roots to calculate NPP. The "Mortality" and harvest components of change were the estimated mass of trees when they died. The net change in live trees (∆Live), is NPP w minus mortality (equivalent to the forestry term "net growth"). The net change in snags (∆Snags) and downed woody material (∆DWM) were based on differences in C between time 1 and 2 measurements. We refer to the net change in all pools (∆Live + ∆Snags + ∆DWM) as net change in woody carbon "Net∆C w ". To put all plots on the same footing, components of change were annualized (i.e., Mg C·ha −1 ·yr −1 ) by dividing by the plot measurement interval. Most C results are presented on a per unit area basis (Mg C/ha), though occasionally we refer to estimates of total C for a stand type (Tg C). Most available volume and biomass equations and decay-reduction factors are based on limited data with unknown bias when applied to regional analyses (Temes-gen et al. 2015); we did not attempt to incorporate these additional errors into our sample error estimates.

Statistical analyses
Statistical analyses of survey estimates differ from approaches commonly applied to designed studies (e.g., ANOVA). We calculated means and variances for all values using the equations for double-sampling for stratification based on sampled condition classes (Cochran 1977, Scott et al. 2005. Essentially, stratum weights for each plot estimate its proportional contribution to the overall population, taking into account differences in sampling intensity due to design and inaccessible plots, and the variation in remotely sensed vegetation attributes, to improve the precision of plot-based estimates. Strata were defined using national forest boundaries, Wilderness boundaries, and classified Landsat satellite imagery (Dunham et al. 2002, Homer et al. 2004. The ratios of the number of pixels in each stratum to the known area of sampled NFS lands in each estimation unit were used to assign population weights to each plot (e.g., MacLean 1972), and ratio estimates (e.g., Mg/ ha) and their variances were calculated using the ratio of means estimator (Scott et al. 2005).
Plots are merely samples of forest stands and in most cases do not describe the mean and variance of conditions within a single stand with much precision. Therefore, we grouped plots with similar stand attributes (e.g., PAZ, MAI, stand age) and estimated mean C values using ratio of means (e.g., total C in a group of plots divided by the total area sampled by that group). The statistical significance of estimated means and variances was assessed with the Type I error for the Z-statistic (Zar 1984). We report most error estimates as the 95% confidence interval around the mean. We also calculated the Holm-Bonferroni adjustment to estimate the "family-wise" error rate for multiple comparisons (Ludbrook 2000) to identify cases where a value of interest that was not within a confidence interval was not significant (identified as P′).
We grouped plots within each PAZ into classes with similar estimated productivity (i.e., MAI) and stand age for analysis and display, based primarily on available sample size for most of the groups. Three MAI classes described as low, v www.esajournals.org GRAY ET AL. medium, and high grouped plots with MAIs of <3.5, 3.5-8.4, and >8.4 m 3 ·ha −1 ·yr −1 , respectively. Stand age classes were 0-20, 20-40, 40-60, 60-80, 80-100, 100-125, 125-150, 150-175, 175-200, 200-250, 250-300, 300-400, and >400 and identified by the mid-point ages of each class: 10,30,50,70,90,113,138,163,188,225,275,350,450, respectively (age intervals increase with age because precision of stand age estimates and sample sizes decline with age). Most analyses present empirically based estimates by categories of interest to assess differences in C across a diverse region. Alternative analyses applying models to continuous variables would require introducing a number of assumptions, including the shape of the relationships between dependent and independent variables, and the nature of the interactions among independent variables. In addition, many of the responses (e.g., NPP and changes in dead wood with age) do not seem to conform to simple linear or nonlinear curves. We provide graphs of many of the responses with continuous variables in the Appendix for context.
To estimate maximum C stocks for different forest conditions and the age at which stands attain 75% of the maximum, we regressed total C (not including mineral soil) on stand age for each PAZ × MAI class group. We used a cumulative two-parameter Weibull model used extensively to model stand-level growth in terms of basal area and volume (Curtis 1967, Hanus et al. 1999): allCdens = exp(a0 + a1 × (stdage a2 )) where allCdens was total C (Mg C/ha), stdage was stand age (yr), and a0, a1, and a2 were parameters to be estimated. Modeling was done using proc NLIN in SAS (SAS Institute 2008) and data sets were restricted to stand age 300 (150 for the PICO-Pinus contorta-PAZ) to avoid problems with extrapolation of models into regions of sparse data. Model predictions were compared to splines to evaluate potential problems with model fits.
To assess the importance of different sizes of live trees to C accumulation, we used absolute diameter classes to assess change across the region, as well as relative diameter classes to assess changes in growth and dominance in stands of different ages and on different PAZs, because absolute tree sizes varied greatly among forest types. Five relative diameter classes (RDCs) of equal basal area (m 3 /ha) were defined on each plot, populated with trees of increasing diameter. Growth and mortality rates of absolute and relative diameter classes were calculated for each class, as well as recruitment from one size class into the next through diameter growth. Total C changes for absolute diameter classes (Mg C·ha −1 ·yr −1 ) were calculated by multiplying individual tree measurements by their trees-perhectare expansion, summing by diameter class per plot, multiplying by plot stratum weights to estimate totals by class and estimating means by dividing with total (weighted) area sampled. Estimating C change for relative diameter classes was similar, but compared different metrics of occupancy because tree density is commonly characterized in a variety of ways (e.g., numbers of trees, basal area, crown area, or biomass). For example, instead of dividing with area sampled as shown above, sums by RDC were divided with total (weighted) number of trees to get per tree change, or with total (weighted) basal area to get per unit basal area change.

Results
The ten PAZs are presented in the results sorted from low mean C stocks (Mg C/ha) to high, as sampled and estimated in this study ( Table 1). The lowest-C PAZs tended to be the driest and/or highest in elevation. Douglas-fir (Pseudotsuga menziesii (Mirbel) Franco) was an important or dominant component of many of the PAZs (species composition and selected attributes by PAZ are shown in Appendix Figs. A1 and A2).

Carbon stocks and chronosequence
The maximum mean C stocks (live and dead woody pools, tree foliage, understory vegetation, and forest floor combined; not including mineral soil) and the apparent rate at which it was reached varied by PAZ and productivity (MAI) class. The Weibull model results identified significant differences in the maximum C attainable by PAZ × MAI class (Fig. 1). Trends of mean C with stand age were not always smooth and were affected by sample size for some PAZ × MAI class combinations (see Appendix  Fig. A3). The modeled curves of C up to stand age 300 (150 for PICO zone) were able to smooth the empirical irregularities comparable to fitting splines and identified differences in the v www.esajournals.org GRAY ET AL. maximum attainable C over time (see Appendix Fig. A4). Maximum C was significantly greater in more productive MAI classes than in less productive MAI classes within the ABAM, TSHEPISI (Tsuga heterophylla and Picea sitchensis), and ABCOGR (Abies concolor and A. grandis) PAZs (Z = 1.99, 4.64, and 5.09; Holm-Bonferroni adjusted P′ = 0.023, <0.001, and <0.001, respectively), and varied among PAZs. The apparent rate of C accumulation (i.e., the steepness of the curve) also differed among PAZs, with TSMEpark and ABCOGR showing the oldest stand age to attain 75% of maximum C, and PIPO (Pinus ponderosa) and JUOC (Juniperus occidentalis) zones the youngest (Z = 5.12, P′ < 0.001; Table 2). The mean stand age required to reach the 75% level across PAZs was 127 ± 35 yr (95% confidence interval).
The proportion of total C in different pools varied significantly among PAZs and was likely affected by inherent differences in vegetation types as well as differences in recent and past disturbance and logging. Most of the non-mineral soil forest C was in the live tree pool, which showed the greatest differences in C among PAZs as indicated by their 95% confidence intervals ( Fig. 2A). The largest C stocks were found in the TSHEPISI and ABAM zones, and were more than twice as high as the those of the five lowest-density zones (Z = 4.49, P′ < 0.001). The ABAM zone tended to have older stands (median = 137 yr), while the TSHEPISI zone tended have stands with higher estimated MAI (median = 8.9 m 3 ·ha −1 ·yr −1 ), than other PAZs (medians = 94 yr and 4.6 m 3 ·ha −1 ·yr −1 ; Fig. 1. Predicted maximum non-mineral soil carbon at stand age 300 (150 for PICO) by plant association zone (PAZ; see Table 1) and MAI class. Bars are 95% confidence intervals around the mean. Note: High-MAI classes excluded for PIPO, PSME, and ABLA due to low N. PAZs sorted from low age to high. see Appendix Fig. A2). The proportion of total stand C in the non-live tree pools differed by PAZ as well (Fig. 2B), with the lowest proportions (0.2-0.35) found in the four PAZs with the greatest mean C (Z = 1.93, P′ = 0.027). Proportions of total C in dead wood (snags and downed woody material (DWM)) were highest in the PICO and ABLA zones (0.31 and 0.24, respectively). In much of the region, Pinus contorta Dougl. ex Loud. (lodgepole pine) and Abies lasiocarpa (Hook.) Nutt. (subalpine fir) have experienced elevated mortality in recent years from mountain pine beetle (Dendroctonus ponderosae (Hopkins)) and western spruce budworm (Choristoneura occidentalis Freeman), respectively. Proportions of total stand C in dead wood ranged between 0.1 and 0.2 in the other PAZs.
The C stocks (Mg C/ha) of dead wood differed by stand age and management history, but the patterns were similar among PAZs, so we combined them for analysis (see Appendix Fig. A5 for chronosequences by PAZ and C pool). For plots with no record of tree cutting within 40 yr prior to measurement, C in standing dead trees (snags) was greater than in DWM for stands 0-20 yr old, but the reverse was true for stands 20-80 yr old as indicated by confidence intervals (Fig. 3). In contrast, in stands with any level of cutting within 40 yr, C in DWM was greater than in snags for all stand ages. C in DWM appeared similar between cut and uncut stands 60-300 yr of age, but C in snags was lower in cut than uncut stands for all age classes. Dead wood C (i.e., snags + DWM) was greater in stands >200 yr old than in younger stands (Z = 14.9, P < 0.0001). We expected that young stands <20 yr old that did not originate from harvest would have greater amounts of dead wood (snags and downed wood), and that they would be comparable to the live tree amounts in older stands from which they might have originated. When we removed sparse and non-stocked stands with no record of natural disturbance from the young-stand set, snag and DWM C stocks were higher than with the undisturbed young stands included (25 ± 3.3 and 9.5 ± 1.6 Mg C/ha for downed wood and snags, respectively), but the combined total (35 Mg C/ha) was still smaller than the mean live tree stocks for most PAZs.

Controls on carbon accumulation rates
On average, forests in all PAZs accumulated C on a net basis (Net∆C w > 0), except for JUOC and PICO (Z = 0.13, P = 0.45 and Z = 1.50, P = 0.067, respectively; Fig. 4). Net∆C w tended to be greater on PAZs with greater current C (Pearson r = 0.86), but there were notable exceptions. For example, Net∆C w was lower in the ABAM zone than in the TSHEPISI zone (Z = 2.58, P′ = 0.020), and was lower in the TSMEpark zone than in the ABCOGR zone (Z = 2.39, P′ = 0.026). Stands in the ABLA zone experienced a net loss in live tree C on average (Z = 3.20, P′ = 0.005). Some of the patterns are likely due to small-scale disturbance, with most C pools showing gains on forests that were not cut or burned ("undisturbed") between measurements (Fig. 4). The net change in live tree C was greater in the undisturbed forests than for all forests, and dominated the gains in all PAZs except for the PICO, ABLA, and ABCOGR zones, where large proportions of the gain in C (>0.4) was due to increases in downed wood. The net loss of snag C and gain of DWM C in the undisturbed PICO zone forests are suggestive of a prior disturbance pulse in a large number of stands (possibly due to mountain pine beetle) where snags fell and added to the downed wood pool.
C accumulation rates varied significantly with stand age. Rates of change with age for individual Fig. 3. Mean carbon in dead wood by stand age for standing dead trees ("snags") and down woody material (DWM), shown for forested plots with and without records of tree cutting ("harvest") within 40 yr prior to measurement. Bars are standard errors around the mean. PAZ by MAI class groups indicated differences in the age at which ∆Live and Net∆C w peaked and whether or how much they declined with age (Appendix Fig. A5). MAI class was the dominant effect on the patterns of Net∆C w with stand age, so we aggregated plots by MAI class and stand age for further analysis. NPP w increased to a plateau at the 80-100 yr age class on low MAI sites, rose more quickly to plateau in the 20-40 yr class on medium MAI sites, and peaked in the 20-60 yr ages and fell by ~33% in older stands on high-MAI sites (Fig. 5). Means with confidence intervals that did not overlap zero were also significantly different from zero in the Holm-Bonferroni tests except in the one case noted below. Mortality rates for the medium and high-MAI classes increased slowly but steadily as stand age increased, eventually matching the rates for NPP w . Consequently, ∆Live was not significantly different from zero for stands over 250 yr old for these two groups. For the low MAI class group, although the effect of the rate of C change due to mortality was more variable (plateauing in the 100 to 300 yr old stands, increasing sharply then dropping again), ∆Live was not significantly different from zero for most stand age classes over 175 yr old.
The rates of change in the dead wood C pools varied less with stand age than live trees, resulting in less effect on the pattern of Net∆C w with stand age. ∆Snag was negative in stands <20 yr old in all MAI classes but was variable and commonly not significantly different from zero in older stands. ∆DWM was significantly positive in most of the 75-200 yr ages for low and medium MAI site classes. On high-MAI site classes, ∆DWM was significantly negative for stands 0-60 yr of age; many stands likely originated from clearcuts, with little wood input from snags and small live trees to balance decomposition.
For all woody C pools combined, Net∆C w was positive in stands 20-175 yr old on all sites, but was mostly not different from zero in stands <20 and >200 yr of age as indicated by the confidence intervals (except for the high-MAI value for 300-400 yr old stands, where P′ = 0.074), though the mean values tended to be >0 and not less. For the combined set of undisturbed stands across all MAI classes, Net∆C w was significantly >0 for all age classes except for the youngest (<20 yr) and oldest (≥400 yr) stands. The Net∆C w rate for the ≥400 yr age class was 0.15 ± 0.64 Mg C·ha −1 ·yr −1 (95% confidence interval; Z = 0.47, P′ = 0.64). In age classes between 200-400 yr old, C accumulation was modest, but significant (0.34-0.70 Mg C·ha −1 ·yr −1 ; P′ < 0.001 for 200-250 yr, P′ = 0.044 for 250-300 yr, and P′ = 0.041 for 300-400 yr age classes).

Tree size effects
Most of the accumulation of C in undisturbed stands across the study region was in the smallest tree sizes, reflecting in part their greater abundance compared to larger trees. Trees <50 cm DBH at time 1 accounted for 69% of the NPP w and 87% of the ∆Live during the growth periods covered by the study (Fig. 6). NPP w of the largest-diameter trees was offset by mortality, with ∆Live significantly <0 for trees 100-150 cm DBH at time 1 (Z = 3.86, P′ = 0.002), and not different from zero for trees Fig. 4. Annual net change in C by pool and combined (Net∆C w ) by each plant association zone, for all forests and for undisturbed (not burned or cut during the measurement interval) forests only. Bars are 95% confidence intervals for Net∆C w . Positive numbers are net gain, negative are net loss. >150 cm DBH. Nevertheless, C stocks in large trees (and all size classes >25 cm DBH) increased in undisturbed stands due to recruitment from smaller size classes (P′ < 0.033 for each). The increase in C in large trees also held true when disturbed stands were included (P′ < 0.027 for each), although increases were proportionately lower, particularly in the smaller tree sizes.
Determining which size classes of trees within undisturbed stands account for most of the annu- al C accumulation depends on how occupancy, or the space used by each tree, is defined. Each relative diameter class made up 20% of the basal area (m 2 /ha) in each plot; summaries across groups of plots therefore describe the relative contributions of different within-stand sizes of trees. Across the study area, NPP w in undisturbed stands was evenly distributed among the first four relative diameter classes (RDCs) but was lowest for the largest size class (Z = 13.3, P′ < 0.001; Table 3). Patterns in ∆Live with RDC were similar, but ∆Live was approximately half (49-54%) of NPP w for the first four classes, but only a third (32%) of NPP w for the largest tree size class, indicating a greater effect of mortality on live C stocks in the largest class. To examine trends of NPP w with RDC using additional metrics of abundance (besides equal groupings of basal area), we calculated ratios of means using sums of growth divided by sums of tree density, C mass, crown area, and basal area per RDC (Fig. 7). As expected, NPP w increased significantly with RDC on an individual tree basis. On an estimated crown-area basis, NPP w also increased with tree size, but the differences among tree sizes were less than on a per tree basis and the largest tree sizes were more similar to each other (though still significantly different; Z = 2.07, P′ = 0.019) than were the smallest tree sizes with each other. As discussed above but expressed in different units to aid comparison, on a per unit basal area basis NPP w was lower for the largest trees than for the first 4 RDCs. Finally, on a mass basis (analogous to relative growth rate), smaller trees in stands had greater rates of NPP w per unit mass than large trees.  The patterns of NPP w by RDC we found for all stands (described above; Table 3) were similar to those among PAZs (not shown) and MAI classes on undisturbed plots (Table 4). ∆Live rates, however, were highest for the three intermediate RDCs on high-productivity MAI class sites (Z = 3.06, P′ = 0.001), but declined with increasing tree size on the low-productivity MAI class sites (RDCs 2-4; Z = 2.78, P′ = 0.003), indicating disproportionately higher mortality in large trees on low-productivity sites. Differences in NPP w by RDC varied among age classes as well. Except for stands <20 yr old, some of which were a combination of large old trees and younger regeneration, the majority of the NPP w per unit basal area shifted from large tree to small tree basal area groups as stand age increased (Fig. 8).

Chronosequences of carbon stocks
C stocks and accumulation rates varied significantly among PAZs and site productivity Fig. 7. Patterns of NPP w for different sizes of trees in within-stand relative diameter classes, using different variables as a basis for density, including individual trees, mass, basal area, and crown area, and showing ∆Live for the basal area metric. Only undisturbed stands were used. Bars are 95% confidence intervals around the mean. (MAI) classes. Predictably, the highest C stocks were in zones with abundant precipitation and moderate temperatures, but regional disturbance history had a role as well, with higher C stocks in the ABAM zone than the TSHEPISI zone, where site productivity was greater but stands were younger. Comparing mean current C stocks in the study area (157 Mg C/ha, or 1428 Tg total; not including mineral soil) with the estimated mean maximum from the oldest stands in each PAZ × MAI class (252 Mg C/ha, 2282 Tg total) suggests these forests currently store 63% of their potential maximum C. Of course even without harvest the maximum would never be reached across a region because of natural disturbances (e.g., severe wildfires, wind storms, and landslides; Smithwick et al. 2002), but under current management and disturbance regimes continued accumulation is likely. We do not consider C produced by the forest but stored off-site in buildings and landfills, or any energy or construction substitutions from use of wood products that might apply in some accounting approaches; comparisons with in-forest C storage are sensitive to assumptions about baselines, leakage, and disturbance regimes (e.g., McKinley et al. 2011). The maximum estimated C stocks for vegetation types found west of the Cascade crest from our study (<400 Mg C/ha for LIDE3 (Lithocarpus densiflorus), TSHEPISI, ABAM, and TSMEpark PAZs) are substantially lower than the >600 Mg C/ha (excluding mineral soil) from westside research plots in Smithwick et al. (2002), but similar to estimates for live trees by Raymond and McKenzie (2013), which was also based on forest inventory data. Across all forest PAZ × MAI classes, stands attained 75% of their estimated maximum C by age 127 yr. This and changes in ∆Live with stand age indicate that the speed and amounts of potential future annual C accumulation are greatest for forests with a large proportion of young stands, and that the potential for additional accumulation declines substantially as forests mature and stocks approach their maximum carrying capacity.
The contribution of dead wood to total C stocks varied among PAZs and did not match the expected chronosequence pattern well. The lowest proportions of total C in dead wood pools were found in the most productive PAZs. Differences in species' decomposition rate, susceptibility to insects and disease, or recent disturbance regimes may be responsible for higher dead wood proportions in less productive zones. Our chronosequence results of dead wood C did not follow a well-defined U-shaped curve of dead wood with stand age following stand-replacement fire in mature forests (Harmon et al. 1986). Dead wood C in young naturally disturbed stands is primarily determined by the live trees in the predisturbance stands, yet dead wood C in our study was lower in the youngest stands (<20 yr) than the live tree C in mature stands. Low C stocks in dead wood may be due to prior salvage-logging or the tendency for fires to burn at highest intensity in stands that are sparse to begin with (Azuma et al. 2004). Indeed, most of the unmanaged and burned stands in this study were in the low-productivity ABLA, PICO, PIPO, and PSME (Pseudotsuga menziesii) zones. Our analysis is probably not an appropriate chronosequence for dead wood because stands that began with different amounts of legacy wood were mixed together. The chronosequence did suggest a slow accumulation of dead wood in stands continuing past age 200 yr, which was also seen in the net change analysis. Stands with recent harvest had fewer snags than unharvested stands for all age classes, which could reflect less mortality due to reduced live tree density, or the felling of potentially hazardous snags in cutting units.
The value of a study based on a probabilistic forest inventory sample is that the scope of inference is well defined and results apply to all forests in the population of interest, albeit for a defined and relatively short period of time. Few studies have analyzed the components of change in live and dead C pools using repeat field measurements and design-based inference across the diverse set of ecosystems and stand conditions found here. One of the challenges of working with these data is that the large diversity of vegetation types and conditions does not always conform to simple conceptions of vegetation development. While the conceptual model of secondary succession is useful, many forests do not develop from a stand-replacing event without subsequent disturbance. Fire, cutting, or insect outbreaks can thin the original, or oldest tree cohort, allowing new trees to establish. These cohorts can rarely be distinguished by v www.esajournals.org GRAY ET AL.
diameter alone, and the concept of "stand age" becomes tenuous. In these situations, the ages of the dominant trees reflect neither the time of stand origin or the last major disturbance, but a combination of events ("uneven aged" in forestry terminology). While the stand age estimate reflects the importance of the oldest trees, it often does not mean development of an "even aged" stand experiencing few mortality events since stand origin (J. T. Stevens et al. in press). Indeed, the notion of undisturbed old-growth stands is somewhat artificial, given that both understory reinitiation and the old-growth stage are characterized by gap dynamics; i.e., periodic mortality events affecting the dominant trees in a stand. Stand age and even site index may also not be as meaningful in harsh environments where snow or rock may limit stand development more than time since disturbance. For example, stands in the environmentally harsh ABLA and JUOC zones in this study did not show much pattern of C stocks or accumulation rates with stand age. In addition, site index and MAI estimates are based on a few commercial species in "fully stocked" stands, and are therefore rough approximations for non-commercial forest types (Hanson et al. 2002), and sites where substrate or drought conditions limit stands from ever attaining "full stocking" (Cochran et al. 1994). We believe that stand age and MAI are still useful concepts, but additional information on disturbance and management history and site edaphic conditions could improve our understanding of forest C accumulation potential.

Carbon accumulation rates
Overall, we found stability and perhaps slight increases in measured aboveground C pools in undisturbed old-growth forests. Live tree C accumulation rates measured as NPP w were not lower in old stands than in young stands on low and moderate productivity sites, but were approximately one-third lower on the most productive sites. Speculation on why NPP w apparently declined with age on high-productivity sites and not lowerproductivity sites could include some of the mechanisms mentioned in the introduction, as well as the possibility that many of the productive stands <40 yr old were planted and tended and are more densely stocked than older stands. Our results suggest that declines in forest Net∆C w with stand age had little to do with declines in NPP w , and were primarily due to increases in heterotrophic respiration as stands became fully occupied and large trees began to die and transition to dead wood (Coomes et al. 2012, Xu et al. 2012, Foster et al. 2014).
We did not see the ~2 Mg C·ha −1 ·yr −1 of total C accumulation (Net∆C w ) in old-growth stands indicated by Luyssaert et al. (2008). Although dead wood pools seemed to increase in older stands, overall net increases in measured pools in our oldest (>400 yr old) stands were 0.15 Mg C·ha −1 ·yr −1 and not significantly different from zero. Our findings of a rough balance between change in live and dead C also do not support Luyssaert et al.'s (2008) suggestion that positive accumulation rates in old-growth were due to slow decomposition of dead wood and fast recovery of living vegetation from gap events. It is possible that some old-growth forests that are intact or recovering from prior gap disturbances are accumulating C at high rates, while others are experiencing gap events and losing C. Indeed, Wharton et al. (2012) found high variation in eddy-flux covariance measurements of C accumulation rates in an old-growth Pseudotsuga menziesii-Tsuga heterophylla forest around a mean of 0.5 Mg C·ha −1 ·yr −1 . This was in the range we found for 200-400 yr old stands (0.34-0.70 Mg C·ha −1 ·yr −1 ), suggesting a potentially large and sustained C sink as forests continue to age under current management and disturbance regimes.
There are a number of challenges to estimating C stocks and change accurately. It is possible that we would have more definitive estimates with longer remeasurement time spans and/ or measurements of changes in the forest floor and mineral soil C pools. Soil C research in the region has suggested that C in the mineral soil does not change greatly with stand development, but C does accumulate in the litter and humus layers (Homann et al. 2005, Giesen et al. 2008). Much of the uncertainty and discrepancy among different published estimates of total ecosystem C accumulation rates may be due to assumptions with respect to soil processes (McKinley et al. 2011). Although we do not have age-and PAZ-appropriate estimates of leaf and fine root turnover for a complete assessment of NPP, there v www.esajournals.org GRAY ET AL.
is little reason to think that we would find substantially different patterns between NPP w and overall NPP (e.g., Kira andShidei 1967, Ryan et al. 2004). The scope and applicability of available volume and biomass equations applied to regional analyses is unknown, though initial work is being done to address this (e.g., Temesgen et al. 2015). In particular, a single value is usually used to estimate the proportion of total tree biomass in coarse roots for a species, regardless of tree size or growing conditions, due to a lack of empirical information. Since proportional allocation of plants to below-ground structures usually increases as soil resources become more limiting (e.g., Santantonio and Hermann 1985), it is likely that the differences in above-ground NPP w among vegetation types we report are greater than if we had better information on below-ground allocation.

Role of tree size
Most of the C accumulation in the region was due to the NPP w of the smallest tree sizes (<50 cm DBH). Although large trees (>100 cm DBH) that remained alive did accumulate C (Stephenson et al. 2014), the change in live tree C due to mortality was greater for large trees than for smaller trees, with the net effect being a loss of live C or change not significantly different from zero. This should be expected since most living organisms grow large or old and eventually die. However, over this time period, the region was gaining C in large trees as trees from smaller sizes grew into larger sizes, more than replacing the live C lost from mortality. So despite apparent increases in mortality over time in the western USA (van Mantgem et al. 2009), and shifts in tree species' ranges on the west coast (Monleon and Lintz 2015), on balance these National Forest lands are gaining large trees. The change in large trees holds whether assessing all stands or only undisturbed stands, as mortality from fires and logging during this period disproportionately killed smaller trees. This is also not too surprising given the reduction in logging and shift in emphasis to restoration treatments on PNW national forests, and the low rates of intense wildfire on forestland (Moeur et al. 2011, Gray and Whittier 2014, T. R. Whittier and A. N. Gray 2016. Growth rates of large trees may be under-estimated, however, given the lack of data from large trees used in existing biomass equations (Temesgen et al. 2015), the use of diameter and height to predict biomass change, and the finding that some old trees accumulate more biomass on upper stems than at DBH (i.e., change taper; Sillett et al. 2010).
We found that larger individual trees do accumulate more C per tree than smaller trees, as has been found in other studies (e.g., Stephenson et al. 2014). Determining which sizes of trees within a stand accumulate more C per unit land area occupied turned out to depend greatly on the definition of occupancy. On a basal area basis (the most common density metric in forestry), NPP w per unit basal area in the largest tree class was 76% of the mean of the smaller size classes; the differences were even more pronounced on a mass basis. On an estimated crown-area basis, however, NPP w of the largest size class was similar to medium sizes and greater than the smallest sizes. The proportion of overall NPP w accumulated by smaller trees increased as stand age increased. This is counter to the expected decline in growth efficiency as more trees in a stand become non-dominant (Binkley 2004). The difference might be that many natural stands have multiple species with multiple resource tolerances and that over time stands become more spatially heterogeneous and multilayered and use space more efficiently, with large amounts of foliage in upper and lower layers of stands (Van Pelt and Nadkarni 2004). Indeed, in mixed-species stands it can be the understory trees that respond the most with growth after disturbance events (Gray et al. 2012).

Conclusions
Analyses of probabilistic inventories of large areas of forest land are not only important for estimating C stocks and flux (e.g., Coulston et al. 2012), but can also provide insights into the ecological drivers behind the vegetation patterns seen across a region. These systems of permanent plots capture most of the range in site conditions, stochastic events, and stand structure at a finer resolution than is possible through other means. While inventory plots are not designed for intensive process-level studies, any improvements in v www.esajournals.org GRAY ET AL.
site-specific information on disturbance and management history, microclimate, and edaphic conditions will improve our ability to interpret the vegetation changes that are measured. Revisiting our first hypothesis, we found that while forests in our region are accumulating woody C, the variation in accumulation rates is high and readily explained by differences in site productivity and stand age. Similarly, the potential C stocks in oldgrowth forests vary significantly by vegetation type and site productivity. Maintaining or increasing C stocks on the land is one of many management goals, but specific information on stocks and fluxes related to vegetation characteristics used in management could be helpful. With regard to the second hypothesis, our results suggest that the decline in C accumulation (Net∆C w ) with stand age was not due to a decline in tree growth (NPP), although some reduction in NPP on productive sites may have played a role; declines in Net∆C w were primarily related to increased decomposition-related losses of dead wood. We found that net C accumulation of wood in old (>200 yr) stands was low and indistinguishable from zero in the oldest (>400 yr) stands, but the dead wood pool may be where any increased accumulation is ending up in this age class. Finally, our analyses suggested that small trees play an important role in overall C accumulation, but that the C stored in large trees continues to increase in our region.