Estimating carbon storage in urban forests of New York City

Forests play an important role in mitigating many of the negative effects of climate change. One of the ways trees mitigate impacts of climate change is by absorbing carbon dioxide and storing carbon in their wood, leaves, roots, and soil. Field assessments are used to quantify the carbon storage across different forested landscapes. The number of trees, their size, and total area inform estimates of how much carbon they store. Urban forested natural areas often have greater tree density compared to trees planted in designed cityscapes suggesting that natural area forests could be an important carbon stock for cities. We report a carbon budget for urban forested natural area using field-collected data across an entire city and model carbon stock and annual stock change in multiple forest pools. We find that natural area forests in New York City store a mean of 263.04 (95% CI 256.61, 270.40) Mg C ha-1 and we estimate that 1.86 Tg C (95% CI 1.60, 2.13 Tg C) is stored in the city’s forested natural areas. We provide an upper estimate that these forests sequester carbon at a mean rate of 7.42 (95% CI 7.13, 7.71) Mg C ha-1 y-1 totaling 0.044 Tg (95% CI 0.028, 0.055) of carbon annually, with the majority being stored in trees and soil. Urban forested natural areas store carbon at similar and in some cases higher rates compared to rural forests. Native oak-dominated forests with large mature trees store the most carbon. When compared to previous estimates of urban-canopy carbon storage, we find that trees in natural area forests in New York City account for the majority of carbon stored despite being a minority of the tree canopy. Our results show that urban forested natural areas play an important role in localized, natural climate solutions and should be at the center of urban greening policies looking to mitigate the climate footprint of cities.


Introduction
Cities are net sources of carbon emissions (Rosenzweig et al. 2010) and human development and urbanization has been linked to loss in tree cover globally (Crowther et al. 2015). As a part of addressing these negative impacts, cities have been enacting policies and programs to reduce net emissions (Kabisch et al. 2016). Natural climate solutions are one of the approaches cities have promoted and can include increasing urban tree canopy cover and planting trees within municipal boundaries (Fargione et al. 2018). Urban tree cover is found to be associated with multiple positive outcomes that can provide benefits to city life, including reducing urban heat islands (Melaas et al. 2016), stormwater capture (Berland et al. 2017) and in the US urban trees are estimated to store 708 Tg of carbon (C) and sequester 25 Tg C y -1 (Nowak and Crane 2002). However, these estimates depend on accurate characterization of tree species, size, density, and growing conditions associated with urban tree canopies. Tree density, basal area and biomass vary markedly among tree canopy types in cities, especially between designed and natural area ecosystems (Pregitzer et al. 2019b).
Forested natural areas are ecosystems embedded within the boundaries of a city. They differ from street trees, yard trees, and landscaped parks in that they usually have higher tree density and more native canopy (Pregitzer et al. 2019b). They also require different management strategies that include invasive species management and promoting natural regeneration at the stand-scale, rather than individual tree care (Steenberg et al. 2019). Urban forested natural areas occur in patchy spatial arrangements in areas of preserved parkland or private property, rather than along linear streetscapes, plazas or backyards (Zipperer et al. 1997). Due to this patchy spatial arrangement and low relative land area as a proportion of city-scale tree-canopy cover, forested natural areas often get missed or underestimated in assessments of urban forests. For example, in New York City forested natural areas make up one quarter of the total tree canopy but occupy only 5.5% of the city land area, yet are estimated to contain approximately 70% of the total number of trees in the entire city (Pregitzer et al. 2019b). This comparison suggests that urban natural areas have potential to contribute to city-level natural climate change solutions. However, most estimates of the carbon mitigation potential for urban forests and trees do not account for spatial variation in tree densities and composition within urban canopies (Nowak et al. 2008(Nowak et al. , 2013. When estimates do account for differences in canopy type, urban natural areas are often lacking data and thus assumed to be more similar to rural forests than other urban canopy types, and hence estimates of carbon storage are based on rural values (Fargione et al. 2018).
To test the veracity of such assumptions requires empirical measurement of carbon stocks and fluxes in urban natural area forests, which would then account for any differences in forest structure, composition and processes due to the urban context (e.g. increased invasive species, older trees) that might impact estimated carbon budgets.
To better understand and characterize the role that urban natural area forests could play in mitigating carbon emissions through city-level natural climate solutions, we applied field-collected data and remote sensing to ask the following questions in New York City, NY, USA: 1) How much carbon is stored in urban natural area forests? 2) How much carbon is stored in different carbon pools in urban natural area forests? 3) What is the range in carbon storage across different vegetation types, including native and invaded forests? and 4) What is the potential for carbon sequestration in natural area forests? We then compare and discuss these findings in relation to existing estimates of carbon storage in cities and in rural forests.

Study area and design
New York City (NYC,40.7128°N,74.0060°W) is located on the edge of the Atlantic Ocean between the New England and Mid-Atlantic US regions. NYC is the most populous city in the US, but despite high population density, 40% of NYC's land cover is greenspace of which 8,006 ha are upland natural areas. These natural areas have been under the care and management of the NYC Department of Parks and Recreation (Parks Department) since 1984 to conserve native ecosystems, and the benefits they provide, through management and restoration activities such as invasive species removal and tree planting. We sampled a subset of upland natural areas (2,947 ha) that are publicly owned by the Parks Department, that have been designated as "Forever Wild". These Forever Wild natural areas are primarily forest but include some open grassland, shrubland and bare soil. Using ArcGIS, a 2-ha grid was clipped to the Forever Wild Parkland boundary, and then within each grid cell one random point was generated and designated as plot center.

Field measurements
Each point was visited in the field in the 2013 or 2014 growing season. If >50% of the plot area was impermeable surface, landscaped, wetland or was unsafe to access we did not sample it (n=200). A total of 1,124 plots (10-m radius) were measured. Field measurements for vegetation are reported in Pregitzer et al. (2019a). Briefly, to estimate the species composition, tree density, and basal area in each plot, the diameter of each overstory tree was measured at 1.37 m from the ground (diameter at breast height; DBH). Overstory trees were defined as woody species >10 cm DBH and included both live and standing-dead individuals. Midstory trees were classified as those between 2-10 cm DBH and tallied. Four 1-m 2 subplots were established 5 m from plot center in each cardinal direction. In these subplots, tree seedlings (<2 cm DBH) were counted and all vegetation <1 m high including herbaceous and woody plants were recorded by species; areal cover was estimated for each species to nearest 1%. In addition to the vegetation measures, for the C budget we also measured coarse woody material (CWM >10-cm dia.), fine woody material (FWM) (0.1-9.9 cm dia.) and litter and duff depth. These pools were sampled at different points with a line-intercept sampling method using a 20-m line going in a random azimuth through the center point of the plot. In addition, four soil samples were collected 5 m from plot center in cardinal directions to a depth of 10 cm and pooled. Samples were analyzed in the lab for several properties, including percent organic carbon and loss on ignition (LOI).

Assigning vegetation types and extrapolation
After sampling, all plots were assigned a vegetation association by ecologists in the New York State Natural Heritage Program (NYSNHP), resulting in 62 different vegetation types. The 62 types were classified into ten vegetation groups for broader comparison for this analysis. We then used the proportion of these plots to estimate the hectares for each group across NYC based on a spatially-explicit, remotely-sensed map (O'Neil-Dunne et al. 2014) that classified hierarchical vegetation classes across the entire extent of NYC. Because the plot classifications and map classifications are independent data sets and were not aligned at the NYSNHP vegetation type, we assigned each of the 10 vegetation groups to the most similar vegetation mapping class to determine the total hectares for each of the 10 vegetation groups. Across NYC, 8,006 ha of natural areas fit within these classifications. We provide estimates of total carbon for each vegetation group in the supplemental materials based on their estimated hectares per group. These vegetation groups have been previously reported at a coarser scale (Pregitzer et al. 2019a, b) but not in relation to their estimated C stock and accrual rates, which we determined as described below.

Overview and assumptions for the C accounting
To determine the C stock in NYC natural area forests, we calculated the C stock in different ecosystem pools (i.e. live trees, downed woody material) and vegetation groups within the forest and uplands using a combination of field-collected data and published estimates of C that could be used to model C from our field measurements. Aside from soils, estimates for C stock and stock change for each pool are based on models and published equations that are applied to field-collected forest plot data. Table 1 lists each forest C pool and provides a brief description of how the pool was defined for this analysis, published benchmark estimates of pool sizes (both urban and rural), net atmosphere exchange rates per hectare from the scientific literature, as well as factors to consider during interpretation of these estimates due to the urban context. We then took our field-based measurements (i.e. volume of downed wood) and applied the most appropriate published estimates for C stock (i.e. 50% of wood biomass).
To model a rate of annual C stock change for each pool we projected a year of growth for the trees and based other stock change estimates on atmosphere exchange alone (except soil), which requires the assumption of no transfers between C pools. As such, our net C accrual rates do not account for tree mortality, which could result in an overestimate of live-tree C accrual but also likely underestimates forest C uptake in other pools for which we report net annual losses, such as downed woody debris and soil, which both receive new inputs over time. However, given a paucity of data for litter and woody debris inputs to urban forests, we did not feel confident estimating these inputs and so our C budget should be seen as conservative as we only account for exchange with the atmosphere and assume no transfer between pools. All equations for estimating annual C stock change can be found in the supplemental materials.

Estimating uncertainty and confidence
All carbon calculation methods used in this budget require data, equations, and assumptions obtained from field assessments, scientific literature, and online databases. Sources of uncertainty vary between pools but in general include: 1) Measurement error from data collected in the field from our field staff, and from field staff in previous published studies (e.g. species identification, DBH, sample processing, etc.); 2) Model assumptions and misspecifications in key equations such as the allometric equations or the percent carbon in biomass (see Martin et al. 2018 who has disputed the common 50% rule) and the application to our study area, which sometimes have no rural forest equivalents from which most allometric estimates are derived; 3) Our study design to capture an appropriate sample size across pools and vegetation groups; 4) Assumptions about biotic and abiotic influences on forest carbon dynamics and translations to the urban environment. In this budget we did our best to minimize uncertainty but given reasons 1) to 4), quantifying individual sources of uncertainty is challenging. Modeling error propagation and overall uncertainty is even more complex as variables are not always independent and errors may be compounded, as is the case with calculations that entail multiple assumptions based on tree species. Our analysis therefore raises questions about confidence in process-based knowledge and parameter estimates which are required for accurate estimates of urban natural-area forest carbon budgets. As such, we view our budget as an initial effort that uses the best-available data and assumptions, and throughout this paper highlight key knowledge and parameter gaps that when addressed will build confidence in carbon budgets for urban natural-area forests.
Here, we estimated uncertainty for each pool using margin of error calculations based on a t-distribution, which we selected because all our data is right-skewed. We recognize there is error and uncertainty in coefficients and values applied from previously published work, as such our uncertainty only reflects variability across plots and measurement error from our study, not error associated with previous published work. We also did not calculate uncertainty around assumptions about biotic and abiotic factors related to our estimates, but acknowledge that many aspects of the urban environment, including heat island, invasive species pressure and preserved and hence older trees, raise important questions about the validity of using conversions developed for rural forests. We believe the most logical estimation approach given our empirical data and these uncertainties is to calculate carbon by pool for each vegetation type. Pools were summed per plot to get a total for stock and annual stock change. Below we detail our calculation steps for each pool. We report mean, median, confidence interval, standard error, and standard deviation across all vegetation type and Table 1 A selection of the published estimates of different carbon pools sampled and urban considerations to contextualize the results from New York City natural area C stocks. The benchmark estimates are intended to provide a point of reference to help contextualize our calculations for carbon pools in NYC's forests. Forest carbon is highly variable and dependent on microclimatic conditions such as moisture, microbial communities, and nutrient availability, all of which can be impacted by human activity in urban or altered environments. Given that most carbon pool estimates are based on rural forests and thus do not account for differences in urban areas, such as the urban heat island effect, we expect to see some differences in these estimates. In addition, some differences in sampling techniques, estimation methods, and/or carbon pool attributes vary between studies that could contribute to some differences between our estimates and others. Our NYC estimates represent the mean across all vegetation types. Standard errors and 95% confidence intervals can be found in the paper text and in supplemental materials Pool considered in NYC Natural Area

Published Estimates Carbon Stock (Mg C ha -1 ) NYC Estimated C Stock (Mg C ha-1) Urban Considerations
Live Trees: All trees (>2 cm DBH) including above and below ground 87.1 -Northeastern US (Smith et al. 2013) 73.3-NYC assuming 100% cover (Nowak et al. 2013) 135.4 Lower ozone levels, higher CO 2 , warmer temperatures, and higher nutrient deposition could lead to increase growth rates and annual carbon sequestration. However, pollutants in soil (e.g. heavy metals), increased pests, and greenhouse gasses in the atmosphere (e.g. NoX and So2) could decrease annual tree growth and C sequestration (Gregg et al. 2003) Groundcover: All vegetation growing <1m height Decomposition increases with temperature (Hanson et al. 2003); decreased ozone levels facilitate litter decay (Carreiro et al. 2009)
Soil could be compacted. carbon pools (Supplemental Table S1). To arrive at the total stock and stock change across all natural-area forests we extrapolated the estimates per pool to the appropriate hectares based on vegetation group. We report all extrapolated estimates based on the mean and report the 95% confidence intervals (lower, upper). We used R statistical software version 4.0.3 for all statistical analysis (R Core Team 2021).

Live tree c stocks
We performed calculations separately on the overstory (trees DBH ≥ 10 cm) and mid-story (trees between 2 and 9.99cm DBH) datasets and combined them. All overstory trees were visually assessed for their vigor class based on missing branches and fine twig die back. Vigor class 1 is over 90% healthy to vigor class 5, which is dead. Dead trees (vigor class of 5) were removed from both datasets and analyzed separately (see standing dead trees pool). Overstory trees with missing DBH (n=5) were also excluded. Mid-story tree DBH was not measured individually but tallied as a size class between 2-10 cm DBH, so we assumed a median DBH of 6 cm for each observation. Given the minimal contribution small trees have on the overall C budget we feel confident this assumed median DBH has minimal impact on the overall results of the C accounting for live trees. We calculated dry-weight aboveground biomass (i.e. bole, branches, twigs, and foliage) using the allometric equation from Jenkins et al. (2003): where: ABM t = aboveground biomass of tree t (kg dry weight) Exp = exponential function β 0 and β 1 = species group coefficients ln = natural log DBH t * = diameter at breast height of tree t (cm) *Note that this equation is intended for use with trees ≥ 2.5 cm DBH We used species group coefficients from Chojnacky et al. (2013) and Woodall et al. (2011), and adhered to the following process to match each species in the dataset to the appropriate species group. First, we matched species to the appropriate combination of taxa (genus or family) and median specific gravity in Chojnacky et al. (2013). For species with both hardwood and woodland coefficients, hardwood was selected. Then if the taxa were not included in Chojnacky et al. (2013), we obtained coefficients from the "REF_SPECIES" spreadsheet in the Woodall et al. (2011) supplemental documents which uses coefficients from Jenkins et al. (2003). Then, if the species was not listed in "REF_SPECIES" table, we substituted coefficients for mixed hardwoods from Jenkins et al. (2003). For trees with missing species information (n=10), we substituted a weighted average based on the proportion of stems in each of the existing species' coefficients. We calculated belowground biomass using the Jenkins et al. (2003) ratio equation: where: R ct = ratio of component c to aboveground biomass of tree t Exp = exponential function β 0 and β 1 = component coefficients (e.g. coarse root, foliage, etc.) DBH t = diameter at breast height of tree t (cm) To get the biomass of each component, we multiplied the ratio of each belowground biomass component (coarse and fine roots) by the tree's aboveground biomass. Total biomass of the tree was obtained by summing all components: where: TBM t = total biomass of tree t (kg dry weight) ABM t = aboveground biomass of tree t (kg dry weight) -from Eq. (1) R frt = fine root component ratio of tree t -from Eq. (2) R crt = coarse root component ratio of tree t -from Eq. (2) To determine the biomass per unit area in megagrams per hectare, total biomass of the plot (the sum of all trees) was converted from kilograms to megagrams and divided by the plot area (0.0314 ha). Carbon content was assumed to be 50% .
where: C x = carbon stored in trees for plot x (Mg/ha) TBM x = total biomass of all trees in plot x (kg) -from Eq. (3) 0.50 = conversion to carbon (50% of dry-weight biomass is carbon) 1000 = conversion from kg to Mg 0.0314 = hectares in plot (circle with 10 m radius)

Herbaceous stocks
We calculated nonwoody plant and seedling biomass using the following formulas from Johnson et al. (2017): where: NWB x = aboveground nonwoody biomass in plot x (Mg/ha) PC x = percent cover in plot x (ratio) b 1 and b 2 = forest type coefficients where: SB x = aboveground seedling biomass in plot x (Mg/ha) SD x = seedling density (stems/ha) in plot x -from Eq. (7) b 1 and b 2 = forest type coefficients Johnson et al. (2017) only provides one set of nonwoody biomass coefficients for all forest types. For seedling biomass coefficients, there are only two forest types available for the Northeast: maple/beech/birch and spruce/fir. Maple/ beech/birch coefficients were used for all plots. We used the following formula to calculate seedling density: where: SD x = seedling density in plot x (stems/ha) ST x = number of stems in plot x (total of all subplots) SP x = number of 1 m 2 subplots in plot x (4 or 10) We calculated aboveground C stocks for each plot by adding together the nonwoody and seedling biomass and multiplying the total biomass by 0.50. Belowground C stocks were assumed to be 11% of aboveground stocks (Smith et al. 2013).
where: THC x = total herbaceous carbon in plot x (Mg/ha) NWB x = aboveground nonwoody biomass in plot x (Mg/ha) SB x = aboveground seedling biomass in plot x (Mg/ha) 1.11 = adds belowground biomass 0.50 = conversion to carbon (50% of dry-weight biomass is carbon)

Standing dead tree stocks
We calculated standing dead tree (SDT) biomass for each tree component using the Jenkins et al. (2003) method (Eqs. 1 and 2). To account for volume loss, we reduced component biomass using structural loss adjustment (SLA) factors ) with the following modifications. (1) We assumed limbs and branches all present and all bark remaining in decay class 1 with 75% reduction from live tree.
(2) To account for height loss, we further reduced stem wood and bark structural loss adjustment factor (the Domke et al. 2011 SLAs are for biomass calculated with the Component Ratio Method, which includes height). For decay classes 2 through 4, we reduced stem wood and bark structural loss adjustment factors by 50% of the previous decay class structural loss adjustment factors. Using this method, total structural loss adjustment factor weighted by decay class is 54%, which is close to the ratio of SDT height to live tree height (0.53) that we calculated using tree observations with height data. For foliage and fine root components, we applied the structural loss adjustment factor directly to the component biomass. For all other components (stem wood, stem bark, above and belowground coarse roots, and branches), we applied the structural loss adjustment factor to the component volume using the following formula: where: ACBM ct = adjusted biomass of component c of SDT t (kg dry weight) ABM t = total aboveground biomass of SDT t (kg dry weight) -from Eq. (1) R ct = ratio of component c to total biomass of SDT t -from Eq.
(2) SLA c = structural loss adjustment factor for component c BD s = bulk density of species s (ratio) -from Harmon et al. (2008), see CWM calculations To account for density loss, we reduced total biomass further using the following equation: where: TBM t = total biomass of SDT t adjusted for structural and decay loss (kg dry weight) DRF = decay reduction factor (ratio) -from We obtained decay reduction factors (DRFs) from Harmon et al. (2011) by matching each observation to the appropriate combination of decay class and species classification (hardwood or softwood). Since 99% of live trees in the dataset were hardwoods, we assumed unknown species to be hardwood. We also assumed SDT carbon to be 50% of biomass (Harmon et al. 2008). To calculate C

Coarse woody material stocks
We field-measured coarse woody material (CWM) pieces >10-cm dia. in 894 of the 1,124 forest plots using a line intercept sampling method (20 m) that crossed the center of the circular 10-m radius plot (and note that we used a plotlevel estimator -see Eq. (15) below -as opposed to a plotless approach to estimate final CWM stocks). This method is modified from Woodall and Monleon (2008). A subset of plots (n=230) did not have consistent sampling so were listed as NA for those plots. For each piece that intercepted with the line we measured the diameter at both ends and the length of the piece as well as rated the structural decay (1-5, with 1 being the most intact and 5 being the most decayed). We then calculated the volume of all downed coarse woody material (CWM) pieces with different formulas depending on the piece type (downed log, pile, or fence). Since we calculated SDT C using the tree dataset, we removed SDTs from the coarse woody material (CWM) dataset. The tree dataset provided greater accuracy than the CWM dataset because all SDTs were recorded, as opposed to only SDTs that intercepted the transect line. For downed logs-the vast majority of pieces-we used the conic-paraboloid formula from Fraver et al. (2007): where: VDW l = volume of downed log l (cm 3 ) L l = length of downed log l (cm) A bl = cross-sectional area at the base of downed log l (cm 2 ) A ul = cross-sectional area at the upper end of downed log l (cm 2 ) To calculate cross-sectional areas, we used the two diameter measurements recorded for each piece using a measuring tape as well as the length and the area formula for a circle. For observations missing a diameter (n=2) measurement, we calculated volume with Huber's formula (Eq. 13) substituting the one recorded diameter measurement for the midpoint diameter. Most class 5 pieces (n=5 out of 6 across all plots) did not have diameter measurements and were thus excluded from the calculations. The CWM dataset noted whether a piece was hollow inside, and if so, provided the diameter of the hollow area. However, only a small number of CWM pieces were hollow (n=31 out of 1884 across all plots), and a study on CWM C stocks found that hollowness has a minor impact on uncertainty (Campbell et al. 2019). Therefore, we did not subtract the hollow area from the piece volume. For piles (n=3 across all plots), we calculated volume using the half-elliptical cylinder formula in Woodall and Monleon (2008): where: V p = volume of pile p (cm 3 ) P = packing ratio = 0.15 (Hardy 1996) H p = height of pile p (cm) W p = width of pile p (cm) L p = length of pile p (cm) We calculated volume of the one fence in the dataset using Huber's formula (Fraver et al. 2007): where: V f = volume of fence f (cm 3 ) L f = length of fence f (cm) A mf = cross-sectional area at the longitudinal midpoint of fence f (cm 2 )assumed to be DBH To account for structural loss in pieces with advanced decay, we multiplied the piece volume by a structural reduction factor (SRF). SRFs were obtained for classes 4 (0.800) and 5 (0.412) from Fraver et al. (2013). Downed logs were the only piece type with decay classes 4 and 5.
To determine the biomass of downed logs and the fence, we multiplied the adjusted piece volume by its absolute density. We obtained absolute density from Harmon et al. (2008) based on the species and decay class. For pieces with unknown decay class (n=5), we substituted decay class 3, the model decay class of CWM pieces. We followed this protocol to match species with appropriate density values. Note that (1) If the species value was unavailable, we used the genus value; (2) If the genus value was unavailable, we substituted the value from a species with a similar wood specific gravity (Jenkins et al. 2003); (3) If wood specific gravity was not available (e.g. invasive shrubs), we used the Ailanthus species value; or (4) If the species was unknown, we used the weighted average for the plot's NY community type (based on relative abundance of live tree species). For piles, we used the FWM bulk density for the plot's forest type instead of absolute density (Woodall and Monleon 2008). We multiplied the volume of the pile by the FWM bulk density and a decay reduction factor of 0.8 (see FWM calculations). To determine the C content of each piece, we multiplied the biomass by percent C (0.48 to 0.51) based off of decay class (Harmon et al. 2008). where:

Fine woody material stocks
We tallied FWM pieces by three size classes, small (0.02-0.6 cm), medium (0.61-2.5 cm) and large (2.51-9.9 cm) using calipers along a portion of the line transect. Small and medium pieces were tallied up to 5 m and large pieces were tallied up to 8 m. We calculated the volume of pieces in all three size classes at the plot level using the following formula from Woodall and Monleon (2008): where: V sx = volume of pieces of size class s in plot x (m 3 /ha) ∑ = sum for all size classes in plot x S = slope correction factor (1.13 default) n sx = number of pieces of size class s in plot x QMDI s = quadratic mean diameter for size class s (cm) TL s = length of transect sampled for size class s (m) We obtained quadratic mean diameter (QMDI) from Woodall and Monleon (2008) for the appropriate FWM size class based on FIA (Forest Inventory and Analysis) forest type comparisons to our plot data. We estimated plot-level C with the following formula. We obtained FWM bulk density for the plot's FIA forest type from Woodall and Monleon

Litter and duff stocks
We calculated average depth for both layers in each plot using the LIS estimator for litter/duff (Woodall and Monleon 2008). This formula is simply the sum of depth measurements divided by the number of depth measurements. We used the following formula to determine C per unit area for each plot: where:

Soil organic carbon stocks
At each plot, four soil samples were collected 5 m from plot center in cardinal directions to a depth of 10 cm, pooled, air dried for 30 days and analyzed in the laboratory (n=1,017, missing 107 soil samples). Soil texture was measured using the hydrometer method (Ashworth et al. 2001) and soil organic matter was calculated by mass loss of soils heated at 440° C for 18 h in a muffle furnace, loss on ignition (LOI) (Nelson and Sommers 1996). Total organic C was analyzed in the lab for a subset of the same plots (n=230) using a Thermo CHN elemental analyzer (Thermo Electron, Milan, Italy). For plots missing % C data, we estimated % C to be 0.58 LOI (Pribyl 2010). We tested other methods (De Vos et al. 2005) for modeling %C, including a linear regression model, on observations with both percent carbon and LOI data; the conversion factor of 0.58 produced the most accurate results when compared to the subset of measured estimates of %C. We calculated the volume of mineral soil for each plot using the following formula: where: V sx = volume of soil in plot x to 10 cm depth (cm 3 /ha) V h = total volume of a hectare to 10 cm depth = 1E9 cm 3 V xr = volume of coarse roots in plot x to 10 cm depth (cm 3 ) = 65% of belowground coarse root biomass (g) ÷ bulk density (g/cm 3 ) CF x = coarse fragment of plot x (ratio)estimated from NRCS soil series description Because we used a composite soil sample, we did not analyze soil samples for bulk density, so we obtained bulk density from SoilGrids (Retrieved from https:// www. isric. org/ explo re/ soilg rids) using the plot center coordinates (Soil Survey Staff 2019). SoilGrids provides bulk density for depths of 5 and 15 cm, so we averaged the two values to obtain bulk density to 10 cm. To estimate the coarse root volume in the soil, we used belowground biomass results from the live tree pool calculations. We assumed that 65% of coarse roots are in the top 10 cm of soil (estimated from figures in Yanai et al. 2006). To get volume, we divided 65% of the root biomass by the bulk density of the tree species. We summed the root volume for all trees in the plot and subtracted this from the total soil volume in the plot. The coarse fragment of report coarse fragments for a typical pedon with different estimates for each soil horizon. However, if the first 10 cm of soil included multiple horizons, we used an average of the horizons or if a range was provided instead of a single value, we used the median of the range. However, for soils designated as "rocky," "stony," or any variation of the two, we selected the higher end of the range. We assumed rock outcrop to be 100% coarse fragment. If a soil series description was not located, we used the average of other soil types. If we could not estimate bulk density from SoilGrids, we We used the following formula to calculate SOC per unit area for each plot to a depth of 10 cm and also to 30cm. Most carbon budgets report soil stocks to 30 cm or 1 m (i.e. IPCC GHC inventories measure to 30 so many people want to compare that number). Our assumptions to extrapolate down to 30 cm, include that 74% of soil C is in the first 10 cm. This assumption was based on results from a study of soils using national soil survey data, which found that NYC woodland soils have 104 Mg C/ha in the top 30 cm of soil (Cambou et al. 2018). Our SOC plot average to 10 cm is 77 Mg C/ha, which is 74% of the 30 cm estimate from Cambou et al. (2018). This percentage is similar to the Gaudinski et al. (2000) results, which show that, to a total depth of 30 cm, 71% of carbon is in the first 10 cm of soil of the Harvard Forest (percentage estimated from charts). Below is the equation we used to estimate soil C to 10 cm and 30 cm. All main findings present soil to 30 cm and the results to 10 cm can be found in the supplemental materials (but for reference are 0.74 of the results to 30 cm per plot).
where: C x = SOC in plot x to 30 cm (Mg/ha) V sx = volume of soil in plot x to 10 cm (cm 3 /ha) -from Eq. (19) CP x = percent organic carbon (or 0.58 LOI) for plot x (ratio) BD x = bulk density for plot x (g/cm 3 ) -from Soil Grid or Eq. (20) 0.74 = ratio of SOC in top 10 cm to SOC in top 30 cm **did not use to estimate to 10cm 1E6 = conversion from g to Mg

Carbon stocks in nyc natural area forests
We estimate that natural area forests in NYC store 1.86 Tg C (95% CI, 1.60 and 2.13) (Fig. 1). Across all plots, C stocks ranged from 2.67 to 1344.16 Mg C ha -1 with a mean value of 263.5 (95% CI, 256.61, 270.40) Mg C ha -1 (Fig. 2). While we have several putative outliers in terms of plot C stock ha -1 in the northeast (n=6, > 800 Mg C ha -1 ), each estimate was verified as having either especially high coarse woody material C or live tree C. Even so, we evaluated the impact of omitting these values from our estimates and the mean decreased marginally to 259.66 (95% CI 253.26, 266.05) Mg C ha -1 suggesting these play a small role in our overall estimates. We also analyzed the impact of the total stock (21) = ( * * )∕ . ∕ accounting using just the values of soil C to 10cm (rather than the 30cm typically used), this resulted in an overall decrease in the estimate of carbon to be 1.66 Tg C (95% CI, 1.43 and 1.91). Carbon stocks varied by vegetation type. Overall, native Oak-hickory forests stored the most C per hectare with a mean of 311.46 (95% CI 301.17, 321.44) Mg C ha -1 (Fig. 2).

Carbon stocks by pool in natural area forests
Overall, the majority of C was found in the tree and soil stocks. We estimate that 0.83 Tg C (95% CI 0.73, 0.92) is stored in live trees accounting for approximately 45% of the total. We estimate that approximately 41% or 0.76 Tg C (95% CI 0.66, 0.86) is stored in soil, with the rest being stored across the other pools (Fig. 1, Supplemental materials). Across all plots, the amount of C in live trees (above and belowground) ranged from 0.0 to 948.4 Mg C ha -1 with a mean of 135.4 (95% CI 130.13, 140.60) Mg C ha -1 . The amount of C in soil (to 30 cm) ranged from 0.42 to 646.8 Mg C ha -1 , with a mean of 105.11 (95% CI 101.58, 108.64) Mg C ha -1 (Fig. 3). Downed woody material and the litter layer each made up ~4% of the total budget. Downed woody material had stock values ranging from 0.0 to 1206.1 Mg C ha -1 with a mean of 15.25 (95% CI 12.73, 17.77) Mg C ha -1 , and the litter layer had stock values ranging from 0.0 to 83.62 Mg C ha -1 with a mean of 10.95 (95% CI 10.59, 11.32) Mg C ha -1 . The amount of C stored in standing dead trees ranged from 0.0 to 262.8 Mg C ha -1 , with a mean of 5.8 (95% CI 5.08, 6.66) Mg C ha -1 , approximately half that of downed woody material. An approximately similar sized stock of C was found for understory vegetation, whose stock values ranged from 0.0 to 16.0 Mg C ha -1 with a mean of 5.5 (95% CI 5.37, 5.73) Mg C ha -1 . Estimates for all pools by vegetation type can be found in Supplemental materials.

Projected annual c stock change
We project an estimated potential for annual C change in NYC's natural area forest to be positive with a mean value of 7.42 (95% CI 7.13, 7.71) Mg ha -1 y -1 (Supplemental materials), totaling an increase in C across natural area forests of 0.044 Tg y -1 (95% CI 0.028, 0.055) (Supplemental materials). However, individual forest plots varied, with a potential annual change in C ranging from -78.1 to 30.9 Mg C ha -1 y -1 (Fig. 4). Carbon stocks and rates of change varied by vegetation type. The maple lowland vegetation group had the highest mean C change, estimated at 11.42 (95% CI 10.59, 12.25) Mg ha -1 y -1 and the grassland shrubland group and marsh had the lowest net C change, with a mean estimate of 1.60 Mg (95% CI 0.86, 2.34) ha -1 y -1 and -0.57 (95% CI-3.73, 1.09) Mg ha y -1 , respectively (Fig. 2). Native-dominated vegetation types have higher annual rate of net C change than invasive-dominated vegetation, with a mean stock change of 7.9 (95% CI 7.64, 8.22) compared to 5.05 (95% CI 4.17, 5.94) Mg ha -1 y -1 , respectively (Fig. 4). In live trees, annual net C change ranged from 0.0 to 29.4 Mg C ha -1 y -1 , with a mean of 8.98 (95% CI 8.74, 9.21) Mg C ha -1 y -1 (Fig. 4).

Discussion
Our analysis provides one of the first comprehensive C accounting estimates for urban forest using field-collected data (similar studies include Jo 2002;Hutyra et al. 2011). Our C stock budget, which addresses our first objective to estimate the C stock of natural area forests in NYC, relies most on our field, empirical measures and so we have the greatest confidence in this estimate. Notably, our estimate of total C stock in natural area forests is higher than previous estimates for all urban canopy in NYC, even when we consider the uncertainty in our estimate and the low value of 1.58 Tg C stock. Specifically, it has been estimated that all trees in NYC (including forested natural areas) have a stock value of 1.2 Tg C (Nowak et al. 2018). Whereas our mean estimate of the stock for just urban natural area forest is 1.5-times greater, at 1.86 Tg C, and at the high end of our estimate is 2.13 Tg C. We note that prior estimates of C stocks for urban forests likely undercounted natural areas given their patchy distribution, meaning they were not enumerated reliably with previous, non-stratified approaches to sampling trees in urban areas (Pregitzer et al. 2019b). Notably, when just comparing tree stocks, our mean estimate is that live trees in natural area forests store 0.83 Tg C, which is a majority (~70%) of the previous total despite representing only 25% of all estimated tree canopy. As such, natural areas have disproportionally high carbon per hectare compared to other types of urban tree canopy (i.e. street trees). Prior estimates of citywide live-tree C stocks therefore are likely an underestimate of the C stock of live trees in cities, highlighting the need to account for natural area forests to provide robust urban canopy C budget numbers, which our data suggest have been underbudgeted for NYC.
Given their focus on live trees, prior estimates of urban forest C stocks presumably additionally underestimate the C stocks of urban forests given that they do not measure other C pools common to forests, such as downed dead wood (Nowak and Crane 2002). When other forest C pools are considered, our C stock mean estimate more than doubles. For example, the amount of C we found in soils was nearly equivalent to the mass in live trees, just as is also observed in managed temperate forests in the eastern US (Environmental Protection Agency 2017). In addressing our second objective to understand the C stock distribution among different pools, our estimates for live trees and soils therefore suggests a similar distribution between the live tree and soil pools to rural, managed forests. Notably, our estimates are similar to estimates for equivalent forest types in rural New York State, which is not surprising given the similarity in tree density and native species composition between the two (Pregitzer et al. 2019a). We also found that our estimates of downed woody material (mean 12.14 Mg C ha -1 for DWM), were within the range for rural forests, but on the high end of some estimates (Table 1, note our estimate is both FWM and CWM). The high-end estimate may be because tree damage and mortality in natural area forests is high, which will build the down woody material pool. In addition, NYC experienced tropical storms in 2011 and 2012, which may have augmented downed wood carbon stocks in our plots. Nevertheless, damage and mortality may be high in the absence of tropical storms because urban natural area forests are often preserved remnant forest patches, in which timber harvesting has not occurred for many decades, which may lead to some "over-mature" forests. Supporting this possibility, our estimate of C stored in live trees is 20% higher than an estimate in a benchmark study of a similar forest type (Northeastern hardwood, Smith et al. 2013), supporting our observations of multiple trees in many of our plots of age classes >150 years.
The high density of trees in natural area forests compared to other urban canopy types is undoubtedly one factor that contributes to the higher estimated C stocks in our study compared to previous work for NYC (Nowak et al. 2018). Yet in addressing our third objective to understand the range Fig. 4 Estimates of modeled annual carbon stock change (Mg C ha -1 y -1 ) in different vegetation groups in New York City natural area forests. Estimates for stock change that are negative represent net emission, values that are positive represent net accumulation. Each dot represents the value for one plot. Boxes display the first quartile, median and third quartile. Mean values are represented by red dot and are estimated for the vegetation group in C stocks across the different forest types in NYC, we identified an additional potential reason. Specifically, the fact that the oak-hardwood forest type is one of the most abundant in NYC natural areas likely contributes to the high citywide forest C stock that we estimate. This forest type is characterized by large and long-lived, native hardwood tree species, with relatively dense wood when compared with the non-native trees that comprise up to 50% of street tree canopy in NYC (Pregitzer et al. 2019b). The presence of these oak-hardwoods helps to create variation in C stocks across natural areas, and certainly the high wood density and long-life spans of the canopy trees, and perhaps the different land use of many of these stands (i.e. that many of them are likely relic forests as opposed reforested areas), presumably explains their high C stocks relative to other forest types. A threat to these carbon-rich, oak-hardwood forests in NYC is non-native invasive understory species which have become prevalent (Pregitzer et al. 2019a). The abundance of these species makes our estimate of herbaceous understory C stocks ~5-times those of rural forest estimates. Non-native herbaceous species are common in cities, and in NYC non-native herbaceous plants account for ~50% of total understory cover and many of these include species with high areal cover (Pregitzer et al. 2019a). However, despite the high understory abundance, the herbaceous stocks were a relatively low proportion of the total C in our budget. Yet they have the potential to reduce future C stocks because many of the non-natives can compete with native tree seedlings. As such, this dynamic of large overstory trees with invaded understory, which is common in urban areas, could result in a future forest with a different canopy-species composition and structure. Our C budget data suggest that such a shift to non-native forest types (Fig. 2), or a loss of large trees without replacement, would reduce C stocks and annual C increments in urban forests. As such, if C in trees and forests is a climate-mitigation goal for cities, management to maintain native forest types will be an important goal, otherwise C-rich forest types could be jeopardized and existing stocks could decline. Forest stand and invasive species management and natural area protection will therefore be an important priority to consider for C mitigation and climate action plans in urban areas.
Our fourth and last objective was to assess the potential for C sequestration in natural area forests. Our estimate for potential annual C accrual is the most uncertain of all our budget numbers, given that we derive estimates of net C accrual using literature-based equations and assumptions, and from a single field campaign of standing C stock. Many of these limitations are detailed in the Supplemental Methods and in Table 1 we highlight these results are a first estimate of C stock change and knowledge and data gaps that must be filled to develop more confident estimates of the potential for urban natural area forests to accrue C. Probably the most important action is to re-sample our plots to test how the C stocks have changed since their initial measurement in 2013 and 2014. In acknowledging the high number of assumptions that went into our estimate of C annual stock change, we note that others have used similar approaches (Woodbury et al. 2007;Smith et al. 2019;Nowak and Crane 2002) and that our estimate does, at the very least, reveal the potential for urban natural area forests to help contribute to city-level efforts to mitigate their C emission footprint, but further research is needed to verify this. To put this potential into the context of local climate solutions, we estimate that our estimates of annual C stock change in natural area forests could be equal to offsetting emissions from of 35,087 cars (assuming an average emission of 4.6 metric tons of CO 2 per year), offsetting emissions from the estimated 13,000 taxis active in NYC two and a half times over (Environmental Protection Agency 2021).
Despite the fact that urban natural areas may store and show potential to accrue C at similar (and sometimes higher) rates to non-urban forests, urban land accounts for just 3% of the total land cover in the US and a similar proportion globally (United Nations 2018). Urban forested natural areas therefore account for a mere fraction of the C mitigation potential of forests at national, international and global scales. Yet urban forests do offer a direct and local strategy for cities to meet their C emission reduction goals and, further, provide many additional benefits such as cooling, storm water capture, and recreational opportunities. Further, 80% of people globally live in cities, making urban natural areas of high value as educational venues to learn about the C cycle and its reliance on healthy forests. These healthy forests are the most accessible form of high-quality nature to many urban residents, making their protection and preservation a vital component of livable cities for humans and many other species. Nevertheless, they are under threat as increasing numbers of people move to cites. For example, between 2014 and 2019, 4% of natural area parkland, or ~38,000 ha, was converted to other land use types in the 100 most populous cities in the US (Pregitzer et al. 2021). As cities look towards new areas to put housing, and invest in installing trees within the built environment to reach tree canopy goals, natural areas may be overlooked as a vital part of urban infrastructure. Accurate estimates of the C stock potential of natural area forests in cities, along with a fuller appreciation of their value to urban populations and wildlife, will help to inform decisions about natural area preservation as cities seek to make land-use decisions that meet multiple goals, from limiting C emissions to land development to urban livability.
Forgione and Sophie Plitt, Tessa O'Connell and Hunter Armstrong. Thank you to the field staff who helped to collect the data including Pitor Bartoszewski, Hannah Buck, Becca Carden, Kevin Corrigan, Jean Epiphian, Rebecca Gorney, Emory Griffin-Noyes, Catherine Molanphy, Jesse Moy, Devon Nemire-Pepe, Beth Nicholls, Nathan Payne, Hayden Ripple, Aaron Rogers, Stephanie Smith, Brian Tarpinian, Kimberly Thompson, Lillis Urban, Alec Wong, and to Mina Kim and Rachel Charrow for GIS and database support. Thanks to NYC Parks Natural Resources Group, and the NYC Urban Field Station for assistance and support throughout the assessment. Thank you to Rich Hallett and Chris Woodall for reviews of an earlier version of this work. Thank you to Jen Shin for making the illustrations.
Funding This research was supported by the JM Kaplan Fund and a doctoral scholarship to C. Pregitzer from the Yale School of the Environment.
Availability of data and material Raw data tables will be made available on The Dryad Digital Repository.

Conflicts of interest/competing interests
The authors declare no competing interests.