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-m2 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, Pregitzer et al. 2019b) 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 carbon pools (Supplemental Table 1). 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, 2020).
Live Tree C Stocks
We performed calculations separately on the overstory (trees DBH ≥ 10 cm) and mid-story (trees between 2 and 9.99-cm 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):
Equation 1: ABMt = Exp(β0 + β1 ln DBHt), where:
ABMt = aboveground biomass of tree t (kg dry weight)
Exp = exponential function
β0 and β1 = species group coefficients
ln = natural log
DBHt* = 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:
Equation 2: Rct = Exp(β0 + β1/DBHt), where:
Rct = ratio of component c to aboveground biomass of tree t
Exp = exponential function
β0 and β1 = component coefficients (e.g. coarse root, foliage, etc.)
DBHt = 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:
Equation 3: TBMt = ABMt + Rfrt*ABMt + Rcrt*ABMt, where:
TBMt = total biomass of tree t (kg dry weight)
ABMt = aboveground biomass of tree t (kg dry weight) – from Eq. 1
Rfrt = fine root component ratio of tree t – from Eq. 2
Rcrt = 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% (Woodall et al. 2011).
Equation 4: Cx = (0.50*∑(TBMx/1000))/0.0314ha, where:
Cx = carbon stored in trees for plot x (Mg/ha)
TBMx = 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):
Equation 8: NWBx = (b1 ∗ PCx)/(b2 + PCx), where:
NWBx = aboveground nonwoody biomass in plot x (Mg/ha)
PCx = percent cover in plot x (ratio)
b1 and b2 = forest type coefficients
Equation 9: SBx = b1 + b2 * SDx, where:
SBx = aboveground seedling biomass in plot x (Mg/ha)
SDx = seedling density (stems/ha) in plot x – from Eq. 10
b1 and b2 = 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:
Equation 10: SDx = (STx/SPx) * 10000, where:
SDx = seedling density in plot x (stems/ha)
STx = number of stems in plot x (total of all subplots)
SPx = number of 1 m2 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).
Equation 11: THCx = (NWBx + SBx)*1.11*0.50, where:
THCx = total herbaceous carbon in plot x (Mg/ha)
NWBx = aboveground nonwoody biomass in plot x (Mg/ha)
SBx = 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 (Equations 1 and 2). To account for volume loss, we reduced component biomass using structural loss adjustment (SLA) factors (Domke 2011) 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 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:
Equation 12: ACBMct = ((ABMt*Rct)/BD)*SLAc*BDs, where:
ACBMct = adjusted biomass of component c of SDT t (kg dry weight)
ABMt = total aboveground biomass of SDT t (kg dry weight) – from Eq. 1
Rct = ratio of component c to total biomass of SDT t – from Eq. 2
SLAc = structural loss adjustment factor for component c
BDs = bulk density of species s (ratio) – from Harmon 2008, see CWM calculations
To account for density loss, we reduced total biomass further using the following equation:
Equation 13: TBMt = DRF*∑ACBMct, where:
TBMt = total biomass of SDT t adjusted for structural and decay loss (kg dry weight)
DRF = decay reduction factor (ratio) – from Harmon 2011
∑ = sum all components of SDT t
ACBMct = adjusted biomass of component c of SDT t (kg dry weight) – from Eq. 12
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 stocks per unit area, we used the same methods as the live tree pool (see Live Tree Stocks).
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 plot-level estimator - see Eq. 18 below - as opposed to a plotless approach to estimate final CWM stocks). This method is modified from Woodall et al. (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):
Equation 14: VDWl = (Ll/12)*(5Abl + 5Aul + 2√ Abl Aul ), where:
VDWl = volume of downed log l (cm3)
Ll = length of downed log l (cm)
Abl = cross-sectional area at the base of downed log l (cm2)
Aul = cross-sectional area at the upper end of downed log l (cm2)
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. 16) 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 et al. (2008):
Equation 15: Vp = PπHpWpLp/4, where:
Vp = volume of pile p (cm3)
P = packing ratio = 0.15 (Hardy, 1996)
Hp = height of pile p (cm)
Wp = width of pile p (cm)
Lp = length of pile p (cm)
We calculated volume of the one fence in the dataset using Huber’s formula (Fraver et al., 2007):
Equation 16: Vf = Lf*Amf, where:
Vf = volume of fence f (cm3)
Lf = length of fence f (cm)
Amf = cross-sectional area at the longitudinal midpoint of fence f (cm2) – 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 et al., 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).
Equation 17: Cix = Vix*Dix*CPd, where:
Cix = carbon in piece i in plot x (g)
Vix = volume of piece i in plot x (cm3) – from Eq. 15 or 16
Dix = density of piece i in plot x (g/cm3) – absolute for downed logs and fences (Harmon et al., 2008), FWM bulk multiplied by 0.8 for piles (Woodall et al., 2008)
CPd = percent carbon for decay class d – Harmon et al., 2008
We calculated C stock per unit area for each plot using the plot-level line-intersect sampling (LIS) estimator (adapted from Woodall et al. 2008):
Equation 18: Cx = 100(π/2TL) ∑Cix/Lix, where:
100 = convert from g/cm2 to Mg/ha
Cx = CWM carbon in plot x (Mg/ha)
TL = transect length = 2000 (cm)
∑ = sum all CWM pieces in plot x
Cix = carbon of CWM piece i in plot x (g) – from Eq. 17
Lix = length of CWM piece i in plot x (cm)
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 et al. (2008):
Equation 19: Vsx = ∑ (π2/8)*(S*nsx*QMDIs2)/TLs, where:
Vsx = volume of pieces of size class s in plot x (m3/ha)
∑ = sum for all size classes in plot x
S = slope correction factor (1.13 default)
nsx = number of pieces of size class s in plot x
QMDIs = quadratic mean diameter for size class s (cm)
TLs = length of transect sampled for size class s (m) – see Table 9
We obtained quadratic mean diameter (QMDI) from Woodall et al. (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 et al. (2008). As with QMDI, we used an average if the plot did not have an FIA forest type.
Equation 20: Cx = Vx*BDf*DRF*0.50/1E6, where:
Cx = FWM carbon in plot x (Mg/ha)
Vx = volume of all pieces in plot x (m3/ha) – from Eq. 15
BDf = FWM bulk density for forest type f (g/m3) – from Woodall et al. (2008)
DRF = decay reduction factor; average = 0.8 – from Harmon et al. (2008)
0.50 = conversion to carbon (50% of dry-weight biomass is carbon)
1E6 = conversion from g to Mg
Litter & Duff Stocks
We calculated average depth for both layers in each plot using the LIS estimator for litter/duff (Woodall et al. 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:
Equation 21: Cx = (Dlx*BDlf + Ddx*BDdf)*100*0.50, where:
Cx = litter & duff carbon in plot x (Mg/ha)
Dlx = average litter layer depth for plot x (cm)
BDlf = litter bulk density for FIA forest type f (g/cm3) – from Woodall 2008
Ddx = average duff layer depth for plot x (cm)
BDdf = duff bulk density for FIA forest type f (g/cm3) – from Woodall 2008
100 = conversion from g/cm2 to Mg/ha
0.50 = conversion to carbon (50% of dry-weight biomass is carbon)
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:
Equation 22: Vsx = (Vh-Vxr)*(1-CFx), where:
Vsx = volume of soil in plot x to 10 cm depth (cm3/ha)
Vh = total volume of a hectare to 10 cm depth = 1E9 cm3
Vxr = volume of coarse roots in plot x to 10 cm depth (cm3) = 65% of belowground coarse root biomass (g) ÷ bulk density (g/cm3)
CFx = 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/explore/soilgrids) 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 used an equation from Al-Shammary et al. (2018):
Equation 23: BDx = 1.177 + 0.00263*SAx − 0.0439*log(SIx) + 0.00208*SIx, where:
BDx = bulk density of soil in plot x (g/cm3)
SAx = percent sand for plot x (ratio)
log = base 10 log function
SIx = percent silt for plot x (ratio)
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).
Equation 24: Cx = (Vsx*BDx*CPx)/0.74/1E6, where:
Cx = SOC in plot x to 30 cm (Mg/ha)
Vsx = volume of soil in plot x to 10 cm (cm3/ha) – from Eq. 22
CPx = percent organic carbon (or 0.58 LOI) for plot x (ratio)
BDx = bulk density for plot x (g/cm3) – from Soil Grid or Eq. 23
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