Abstract
The inclusion of forest carbon in climate change mitigation planning requires the development of models able to project potential future carbon stocks—a step beyond traditional monitoring, reporting and verification frameworks. Here, we updated and expanded a high-resolution forest carbon modelling approach previously developed for the state of Maryland to 11 states in the Regional Greenhouse Gas Initiative (RGGI) domain, which includes Connecticut, Delaware, Maine, Maryland, Massachusetts, New Hampshire, New Jersey, New York, Pennsylvania, Rhode Island, and Vermont. In this study, we employ an updated version of the Ecosystem Demography (ED) model, an improved lidar initialization strategy, and an expanded calibration/validation approach. High resolution (90 m) wall-to-wall maps of present aboveground carbon, aboveground carbon sequestration potential, aboveground carbon sequestration potential gap (CSPG), and time to reach sequestration potential were produced over the RGGI domain where airborne lidar data were available, including 100% of eight states, 62% of Maine, 12% of New Jersey, and 0.65% of New York. For the eight states with complete data, an area of 228 552 km2, the contemporary forest aboveground carbon stock is estimated to be 1134 Tg C, and the forest aboveground CSPG is estimated to be larger at >1770 Tg C. Importantly, these estimates of the potential for added aboveground carbon sequestration in forests are spatially resolved, are further partitioned between continued growth of existing trees and new afforested/reforested areas, and include time estimates for realization. They are also assessed for sensitivity to potential changes in vegetation productivity and disturbance rate in response to climate change. The results from this study are intended as input into regional, state, and local planning efforts that consider future climate mitigation in forests along with other land-use considerations.
Export citation and abstract BibTeX RIS
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 license. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
1. Introduction
Over the past decade, terrestrial ecosystems have sequestered one-third of global fossil fuel emissions, with forests especially critical for mitigating climate change (Pan et al 2011, Friedlingstein et al 2019). Initiatives to increase the capacity of the terrestrial carbon sink through afforestation and reforestation are actively being implemented from local to global scales (Griscom et al 2017, Fargione et al 2018). Such efforts utilize a range of scientific approaches to quantify forests' carbon capture benefit, often relying on national or continental level estimates (Roxburgh et al 2006, Depro et al 2008, Rhemtulla et al 2009). Higher-resolution information on present and potential future forest carbon stocks could further support policy development and management activities and provide important insights for climate mitigation planning (Cook-Patton et al 2020, Goldstein et al 2020). In particular, there is a growing need for advanced forest carbon models to quantify potential future options for additional carbon storage in forests by integrating scientific advances in remote sensing and modelling (Achard et al 2007, Gibbs et al 2007, Hurtt et al 2014, Reinmann et al 2020).
To be most useful, forest models should contain multiple features. First, they should be highly accurate. While the accuracy of future projections is difficult to ascertain by itself, abundant data exists for the present across a wide range of heterogeneous conditions that can be used in calibration and validation. Moreover, mapping present and modelling future carbon within the same framework can help ensure data consistency and maintain logical consistency between present and future ecosystem condition estimates. This capability may be particularly relevant for carbon budgets, which seek to track emission reductions relative to established mitigation goals and baselines.
Second, models should be able to quantify forest carbon at high spatial resolution. Land-use decisions must ultimately be implemented at local scales, as climate mitigation activities for terrestrial carbon relate directly to land ownership. Furthermore, historical and ongoing deforestation and degradation has intensively fragmented forests and resulted in high spatial heterogeneity in forest carbon distribution (Turner et al 2003, Haddad et al 2015, Ordway and Asner 2020).
Third, models should operate across a range of policy relevant spatial domains including local, state, regional, and ultimately national and global scales. Modelling across large domains supports consistent carbon estimates over space and minimizes discrepancies between jurisdictional estimates, offering greater flexibility for data users. A system crossing large spatial domains could particularly benefit multi-actor collaboratives, which may require a common scientific approach in support of forest carbon trading and reporting schemes.
Advances in forest ecosystem modelling and remote sensing over the past two decades offer the potential to begin to meet these needs. A new generation of forest ecosystem models have been developed to efficiently track more detailed ecological processes over large spatial scales (Hurtt et al 1998, Moorcroft et al 2001, Fisher et al 2018). In remote sensing, optical data have increased in resolution, and in combination with Light Detection and Ranging (lidar) has enabled accurate and high-resolution measurements of forest vertical structure (Dubayah and Drake 2000, Drake et al 2002, O'Neil-Dunne et al 2014a). Together, these advances have been combined to provide accurate initial conditions for forest models (Hurtt et al 2004, 2010, 2016, Thomas et al 2008, Antonarakis et al 2011, 2014). In parallel, these same advances have enabled high-resolution empirical biomass mapping which can provide important baseline maps and be used to validate forest model estimates (Huang et al 2015, 2019, Tang et al 2020). Operationally, a prototype of this integration has been developed for and implemented over the state of Maryland (USA), where the Ecosystem Demography model (ED) was initialized with 1 m airborne lidar forest height and optical imagery to generate 90 m 'wall-to-wall' maps of contemporary forest aboveground carbon stocks, future aboveground carbon sequestration potential, and the time to achieve it in years (Hurtt et al 2019).
This study aims to improve and expand the approach to integrated forest modelling and lidar remote sensing another order of magnitude, to an important multi-state region. Here, we focus on forest aboveground carbon stocks and future forest aboveground carbon sequestration potential over an 11-state region consisting of ten members of the Regional Greenhouse Gas Initiative (RGGI), the nation's first mandated cap-and-trade program for CO2 emissions, plus Pennsylvania, which is expected to join RGGI by 2022. Although Virginia joined RGGI in 2021, it is not included in this study. We refer to these states collectively as the 'RGGI' region. The RGGI region encompasses a land area nearly ten times that of Maryland and includes a large gradient of topography, temperature and precipitation.
2. Data and methods
2.1. Study area
The 11-state RGGI domain of this study encompasses about 281 695 km2 of total land area, including 12 542 km2 in Connecticut, 5047 km2 in Delaware, 23 187 km2 in New Hampshire, 2362 km2 in New Jersey, 804 km2 in New York, 49 977 km2 in Maine, 25 142 km2 in Maryland, 20 202 km2 in Massachusetts, 115 883 km2 in Pennsylvania, 2678 km2 in Rhode Island, and 23 871 km2 in Vermont, respectively. Maine, New Jersey and New York land area is only partially included because of limited lidar coverage over these states.
There are five US EPA level II ecoregions including Mixed Wood Plains, Southeastern USA Plains, Ozark, Ouachita-Appalachian Forests, Mississippi Alluvial and Southeast USA Coastal Plain, and Atlantic Highlands. Across the domain, the annual average temperature ranges from −1.4 °C to 14.4 °C, total annual precipitation ranges from 795 mm yr−1 to 2178 mm yr−1. Most of this area belongs to the Appalachian Highlands with dense forest, and the Atlantic Coastal Plain, which is generally flat and low in elevation. Elevation ranges from sea level to 1917 m at Mount Washington in the state of New Hampshire. Dominant forest types are deciduous forests (i.e., maple-beech-birch, oak-hickory and aspen-birch).
Similar to most of the eastern US, forests in the RGGI region were intensively cleared for agricultural expansion by the mid-19th century, with partial reforestation and restoration after agricultural abandonment and westward relocation. In the mid to late 20th century, significant loss and fragmentation occurred due to urbanization. As a result, land-cover in the RGGI region is comprised of 59% forests and 16% cropland and pasture according to the 2011 National Land Cover Dataset (NLCD).
All states in the region have greenhouse gas (GHG) reduction goals based on state legislation or gubernatorial directives. While the states are actively pursuing efforts in support of climate mitigation, the degree to which forest carbon is currently included in climate planning is highly variable (Lamb et al 2021). Ten of the eleven states in the study domain are members of the US Climate Alliance (USCA) (excepting New Hampshire) and currently participate in the USCA's National and Working Lands Challenge, which commits states to integrating priority actions and pathways for land-based carbon into state GHG mitigation plans. Furthermore, as members of the NASA Carbon Monitoring System's Multi-State Working Group, all 11 states have signalled an interest in better utilizing advanced tools and technologies to better estimate forest carbon stocks and fluxes.
2.2. Model
This study uses the ED model which is an individual-based prognostic ecosystem model. By integrating submodules of growth, mortality, hydrology, carbon cycle, and soil biogeochemistry, ED can track plant dynamics including growth, mortality, and reproduction. Along with plant dynamics, ED can track the carbon cycle, including carbon uptake by leaf photosynthesis, carbon allocation to biomass growth in leaves, roots and stems, carbon redistribution from plants to soil based on plant tissue turnover from dead plants due to mortality and disturbance, carbon decomposition in various pools (metabolic litter pool, structural litter pool, soil slow pool, soil passive pool, wood product pool, harvested crop pool, etc) as well as carbon combustion from fire. Over the last two decades, ED has been continuously developed and combined with lidar and land-use change data to predict ecosystem dynamics and associated water and carbon fluxes across spatial scales (e.g., site, regional, and global) and temporal scales (short-term seasonal to long-term decadal and century) (Hurtt et al 2002, 2004, 2010, 2016, Fisk et al 2013, Flanagan et al 2019). ED distinguishes itself from most other ecosystem models by explicitly tracking vegetation structure and scaling fine-scale physiological processes to large-scale ecosystem dynamics (Hurtt et al 1998, Moorcroft et al 2001, Fisher et al 2018). In ED, vegetation structure (e.g., height and diameter at breast height) and physiological processes (e.g., leaf photosynthesis and phenology) are modelled at the individual scale, where individual plants compete mechanically for light, water, and nutrients. Explicitly modelling vegetation height facilitates a potential connection to lidar data. The most advanced version of ED was used in this study and it has been recently calibrated and evaluated globally by various benchmarking datasets such gross primary productivity, leaf area index, aboveground biomass (AGB), and net biome productivity (NBP) (Ma et al 2020).
2.3. Data
2.3.1. Model drivers
Meteorological variables used to drive ED come from NASA Daymet, available from 1980 to 2017 at daily temporal resolution and 1 km spatial resolution (Thornton et al 2016), and MERRA2, available from 1980 to 2017 at hourly temporal resolution and 0.5° spatial resolution (Gelaro et al 2017). As ED requires hourly meteorological variables to compute leaf carbon assimilation and transpiration, the climatology data from Daymet was interpolated to hourly metrics using the MERRA2 climatology. Specifically, hourly air temperature, humidity and shortwave radiation were calculated using following equations:
where 
, 
and 
are MERRA2 air temperature, air humidity and shortwave radiation at a given hour (hr), respectively. 
and 
are MERRA2 daily average air humidity and shortwave radiation, respectively. 
and 
are MERRA2 daily maximum and minimum air temperature, respectively. 
and 
are Daymet daily average air humidity and shortwave radiation, respectively, and 
and 
are Daymet daily maximum and minimum air temperature, respectively.
CO2 concentration was held constant at 360 ppm both spatially and temporally, a value near the middle of the CO2 range between 1981 and 2017 (Keeling 2008). Physical soil and hydraulic property inputs are from Probabilistic Remapping of SSURGO (POLARIS) (Chaney et al 2016) and CONUS-SOIL (Miller and White 1998). The POLARIS dataset remaps the SSURGO database and fills gaps using a machine learning model and high-resolution geospatial environmental data. The soil water module of ED calculates water content and percolation rate based on saturated hydraulic conductivity, saturated water content and Van Genuchten water retention curve parameters from POLARIS and depth to bedrock data from CONUS-SOIL.
The annual average air temperature, annual total precipitation from NASA Daymet, and soil depth from CONUS-SOIL are shown in figure S2 (available online at stacks.iop.org/ERL/16/045014/mmedia).
2.3.2. Canopy cover and height
ED was initialized with canopy height and tree canopy cover maps to generate AGB. The canopy height map was derived from a 1 m lidar Canopy Height Model (CHM) (O'Neil-Dunne et al 2014a, 2014b). Utilizing suggested height metrics from Hurtt et al (2019) and ED's native 10 diameter canopy scale, the 1 m CHM was first aggregated to 10 m by extracting the max canopy height and then averaging to 90 m (hereafter referred to as lidar canopy height). The tree canopy cover map was derived using an object-based approach that integrated the lidar data and multi-spectral optical images from the National Agricultural Imagery Program (NAIP). NAIP optical images with 1 m spatial resolution were first classified as tree cover in conjunction with lidar canopy height, and then aggregated to 90 m (O'Neil-Dunne et al 2014a, 2014b). The lidar canopy height and NAIP tree canopy cover at 90 m resolution are shown in figure 1, with lidar acquisition year shown in figure S24. Lidar and NAIP data were available for 0.65% of New York, 12% of New Jersey, 62% of Maine and 100% of Connecticut, Delaware, Maryland, Massachusetts, New Hampshire, Pennsylvania, Rhode Island and Vermont. Specifically, lidar and NAIP data cover the New York counties of Bronx, Kings, New York, Queens, and Richmond as well as the New Jersey counties of Morris, Sussex, and Warren.
Figure 1. Ninety meter lidar canopy height (a) and NAIP tree canopy cover (b) over the RGGI region. Using a sample region in New Hampshire, (c)–(h) illustrate the process of 90 m lidar canopy height and NAIP tree canopy cover generation, where (e) and (h) are lidar canopy height and NAIP aerial imagery at 1 m resolution; (d) utilizes (e) to identify the maximum lidar canopy height over 10 m land cells; (g) is NAIP tree canopy classification of (h) at 1 m; (c) and (f) are derived by averaging (d) and (g) respectively to 90 m resolution.
Download figure:
Standard image High-resolution image2.3.3. Land cover
Land cover of non-forested wetland, inland water, and impervious surface was excluded from the analysis. Specifically, land-cover types of open water and herbaceous wetland in NLCD 2011 were aggregated from 30 m to 90 m by counting the percentage of each land-cover type, and Percent Developed Imperviousness in NLCD 2011 was aggregated to 90 m by averaging the original 30 m fraction (Xian et al 2011, Jin et al 2013). Then the aggregated 90 m data was used to proportionally exclude water, wetland and impervious surface from each land cell.
2.4. Model initialization, projection and evaluation
This study generally follows the initialization and projection approach used in Hurtt et al (2019), but proposes a modification to the initialization method (here defined as weighting-based initialization method) to improve AGB estimates where ED-modelled canopy height saturates. Full details of this initialization and projection process can be found in supplement S1. We simulated ecosystem succession from bare ground to mature forest by running the ED model for 500 years with meteorology, CO2 and soil properties as inputs (described in section 2.3.1). This model run generated a gridded lookup table that stores a time-series trajectory of AGB and canopy height over the course of succession for each 90 m land cell. In post-processing, ED's stored AGB-height lookup table was then initialized and projected with maps of lidar canopy height, NAIP canopy tree cover and NLCD impervious fraction.
ED initialization combined the stored lookup table with lidar height and tree canopy cover to estimate contemporary AGB (hereafter referred as ED initialized AGB). The ED initialized AGB represents the present aboveground carbon stock of existing trees at the time of lidar acquisition, and it sets the successional baseline for future carbon sequestration. ED then projections AGB growth from ED initialized AGB to the maximum potential AGB based on current meteorology, CO2 and soil properties for every 90 m land cell. This projection does not reduce current impervious surface area, nor does it consider land-use change projections or local laws restricting areas for afforestation or reforestation. Future carbon sequestration was calculated from both the continued growth of existing trees, as identified by the NAIP canopy cover map (hereafter referred as continued growth), and newly reforested or afforested trees, as simulated on any proportion of the 90 m grid cell not otherwise covered by impervious surface, open water, or herbaceous wetland (hereafter referred as regrowth). Following Hurtt et al (2019), several projection-related metrics were defined: 95% of the future maximum AGB the land can reach was defined as the carbon sequestration potential (CSP), the difference between CSP and ED initialized AGB was defined as carbon sequestration potential gap (CSPG), the time in years required to reach CSP from present AGB was defined as carbon sequestration potential time gap (CSPTG). ED initialization and projection processes were only completed for areas where lidar canopy height and NAIP tree canopy cover data were available (section 2.3.2). Note that this study's projections of CSP, CSPG and CSPTG only involve aboveground carbon, because it is observable component of forest carbon from lidar remote sensing.
ED initialized AGB was directly validated using aboveground live biomass estimates from US Forest Service Forest Inventory and Analysis (FIA) plots (hereafter referred to as FIA plot AGB). For this study, 4540 FIA plots were extracted from the FIA database based on three criteria: (a) geographically located within RGGI regional domain; (b) estimates were reported within period of 2004–2015; and (c) sites are in forest condition and with AGB larger than zero. For each FIA plot, all overlapping 90 m land cells of ED initialized AGB were averaged to compare against FIA plot AGB. Spatial mismatches between FIA plots and ED initialized AGB may affect evaluation, especially in highly fragmented forest, because the FIA plot footprint (about 0.4 ha) is smaller than the land cell size of ED initialized AGB maps (about 0.8 ha).
ED initialized AGB was also compared to wall-to-wall AGB maps including AGB estimated from lidar-informed empirical models (Huang et al 2019, Tang et al 2020), hereafter referred to as Lidar empirical AGB, and to other existing AGB products (Blackard et al 2008, Saatchi et al 2012, Kellndorfer et al 2013, Wilson et al 2013, Santoro et al 2018). These wall-to-wall AGB maps vary widely in terms of input data, modelling method and spatial resolution. However, comparison to these maps allowed for further evaluation of ED initialized AGB over regions with a lack of FIA plots with trees measured, such as non-forested areas, and allowed us to examine the efficacy of using high-resolution (1 m) lidar and NAIP imagery data to capture fine-scale AGB heterogeneity.
3. Results
3.1. ED initialization and evaluation
ED initialized AGB is shown in figure 2(a). The AGB spatial pattern corresponds well to landcover, lidar canopy height and topography. Generally, high AGB occurs in mountainous areas along the Appalachians Mountains in western Maryland, northcentral Pennsylvania, and southern Vermont and New Hampshire. Low AGB occurs in populated and cultivated areas such as eastern Maryland, south-eastern Pennsylvania, most of the Atlantic Coastal Plain, eastern Massachusetts, and the Connecticut River Valley. Table 1 summarizes average AGB density and aboveground carbon stocks from ED by state. Southern states such as Maryland, Delaware and Pennsylvania have relatively lower AGB densities (<100 Mg ha−1) than the New England states; this pattern was also identified by inspection of the NAIP imagery. Aboveground carbon stock in the RGGI region is estimated as 1134 Tg C for 228 552 km2, excluding the states of Maine, New York and New Jersey because of their partial coverage of lidar and NAIP data.
Figure 2. RGGI region maps of (a) ED initialized AGB; (b) CSP; (c) CSPG, defined as the difference between CSP and initialized AGB; (d) CSPTG, defined as time in years to reach to carbon sequestration potential from initialized AGB.
Download figure:
Standard image High-resolution imageTable 1. Statewide average NAIP tree cover, average AGB density and carbon stocks of ED initialized AGB, as well as CSP, CSPG and average CSPTG for the states of Connecticut, Delaware, Maryland, Massachusetts, New Hampshire, Pennsylvania, Rhode Island and Vermont.
| NAIP tree canopy cover (%) | ED initialized AGB | CSP | CSPG | |||||
|---|---|---|---|---|---|---|---|---|
| Density (Mg ha−1) | Stocks (Tg C) | Density (Mg ha−1) | Stocks (Tg C) | Density (Mg ha−1) | Stocks (Tg C) | CSPTG (year) | ||
| Connecticut | 68.8 | 119.2 | 75.4 | 247.2 | 156.3 | 139.8 | 81.3 | 231 |
| Delaware | 35.8 | 43.7 | 11.1 | 251.8 | 64.0 | 212.6 | 53.0 | 292 |
| Maryland | 48.8 | 80.1 | 101.7 | 270.0 | 343.1 | 203.5 | 242.1 | 262 |
| Massachusetts | 69.7 | 111.2 | 114.1 | 221.3 | 227.0 | 121.6 | 113.6 | 229 |
| New Hampshire | 85.5 | 133.8 | 156.5 | 224.6 | 262.7 | 100.9 | 107.2 | 207 |
| Pennsylvania | 63.7 | 90.7 | 528.1 | 265.8 | 1547.5 | 183.8 | 1022.0 | 279 |
| Rhode Island | 69.5 | 111.8 | 15.3 | 237.0 | 32.4 | 134.0 | 17.2 | 242 |
| Vermont | 79.6 | 109.7 | 131.8 | 225.6 | 270.9 | 122.4 | 139.7 | 246 |
ED initialized AGB correlates with AGB from FIA plots and lidar wall-to-wall maps at grid-level comparison (figure 3). Using over 4000 plots, figures 3(a) and (b) suggest that ED initialized AGB explains moderate variations in FIA plot AGB (R2 of 0.35), and also shows close alignment with the 1:1 line with a bias of 7.22 Mg ha−1 and RMSE of 61.87 Mg ha−1. Relatively stronger agreement is shown in figures 3(c) and (d), with comparison of ED to lidar empirical AGB; there was a higher R2 of 0.85 and smaller RMSE of 29.54 Mg ha−1. The comparison of ED to wall-to-wall lidar empirical AGB includes an expansive spatial domain (about 34.9 million 90 m land cells) as well as non-forest area not otherwise sampled by FIA plots with measured trees. Despite the overall high correlation between these two maps, ED initialized AGB differs from lidar empirical AGB for larger ED-based AGB (>200 Mg ha−1), where there are fewer lidar-based AGB estimates above 200 Mg ha−1. Further comparison of lidar empirical AGB to FIA plot AGB in figure S3 indicates that lidar empirical AGB is likely to underestimate AGB where estimates exceed 250 Mg ha−1. ED initialized AGB using the mid-point method is also compared to FIA plot AGB in figure S4, which shows larger RMSE and bias than the weighting-based method in figures 3(a) and (b). Comparison between figures 3(a) and (b) and figures S4(a) and (b) highlights the improvements gained by using the weighting-based method, such as increased correlation between ED initialized AGB and FIA plot AGB and correction for the overestimation of AGB after height saturation.
Figure 3. Density scatter plots and histograms comparing ED initialized AGB to FIA plot AGB in (a) and (b), and to lidar empirical AGB in (c) and (d) for all 90 m land cells. For (a) and (b), the corresponding ED initialized AGB is obtained by averaging original 90 m ED initialized AGB over overlapping land cells within the bounded circle area of four FIA subplots (about 40 m in radius).
Download figure:
Standard image High-resolution imageED initialized AGB was also compared to the lidar empirical AGB and existing wall-to-wall AGB maps at both the state and county levels. Figure 4 illustrates that, at county level, ED initialized AGB has relatively more agreement (i.e., RMSE and bias) with lidar empirical AGB, NBCD (Kellndorfer et al 2013), Saatchi et al (2012) and GlobBiomass (Santoro et al 2018) than with the Blackard et al (2008) and Wilson et al (2013) maps. Stronger correlation among datasets can be found when comparing carbon stocks rather than average densities at the county level. Aboveground carbon stocks at the state level (table 2) also show closer agreement among ED initialized AGB, lidar empirical AGB, NBCD, Saatchi et al (2012), and GlobBiomass. Aboveground carbon stock from ED initialized AGB is estimated as 1134 Tg C in the RGGI region (excl. Maine, New Jersey and New York because of incomplete lidar and NAIP data coverage), which is within range of other AGB maps (939–1152 Tg C).
Figure 4. Comparison of ED initialized AGB to lidar empirical AGB and existing AGB products, including NBCD (Kellndorfer et al 2013), Saatchi et al (2012), GlobBiomass (Santoro et al 2018), Blackard et al (2008), and Wilson et al (2013), for county-wide (a) average AGB density and (b) carbon stocks.
Download figure:
Standard image High-resolution imageTable 2. Statewide aboveground carbon stocks (Tg C) estimated by ED initialized AGB, lidar empirical AGB, and existing AGB products, including NBCD, Saatchi et al (2012), GlobBiomass, Blackard et al (2008) and Wilson et al (2013).
| State | ED initialized AGB (Tg C) | Lidar empirical AGB (Tg C) | NBCD (FIA) (Kellndorfer et al 2013) (Tg C) | Saatchi et al (2012) (Tg C) | GlobBiomass (Santoro et al 2018) (Tg C) | Blackard et al (2008) (Tg C) | Wilson et al (2013) (Tg C) |
|---|---|---|---|---|---|---|---|
| Connecticut | 75.4 | 71.9* | 71.9* | 53.5*** | 66.4** | 50.6**** | 54.3*** |
| Delaware | 11.1 | 12.8++ | 10.8* | 11.2+ | 13.0++ | 8.1*** | 11.9+ |
| Maryland | 101.7 | 101.7+ | 81.3*** | 87.3** | 98.7* | 64.9**** | 79.7*** |
| Massachusetts | 114.1 | 110.3* | 111.7* | 82.0*** | 111.0* | 80.9*** | 87.8*** |
| New Hampshire | 156.5 | 149.3* | 150.0* | 120.8*** | 151.8* | 138.2** | 126.6** |
| Pennsylvania | 528.1 | 560.8+ | 548.8+ | 561.4+ | 550.4+ | 453.2** | 474.0** |
| Rhode Island | 15.3 | 15.1* | 12.0*** | 10.6**** | 13.2** | 9.0***** | 9.5**** |
| Vermont | 131.8 | 130.4* | 139.6+ | 116.6** | 138.0+ | 134.0+ | 119.3* |
Superscripts represent deviation degree between ED and other AGB products. The + indicates that the estimate is greater than ED initialized AGB, the * indicates that estimate is lower than ED initialized AGB. The number of +/* symbols next to each estimate represents the relative difference at the intervals of 0%–10%, 10%–20%, 20%–30%, 30%–40%, 40%–50%. For example, ***/+++ represents either a −20% to −30% or 20% to 30% difference from ED initialized AGB.
3.2. ED projection of carbon sequestration potential
The ED-projected CSP for the RGGI region is shown in figure 2(b). The spatial pattern of CSP more generally reflects the heterogeneity of environmental conditions rather than current landcover. Areas of higher sequestration potential tend to appear in regions with deep soil or warmer air temperatures (figure S2). For example, warmer air temperatures together with deeper soil in central Maryland and south-eastern Pennsylvania result in larger sequestration potential estimates than in other areas of the region. Relatively cooler air temperatures in the White Mountains lead to lower sequestration potentials in northeastern Vermont, northern New Hampshire and western Maine. The statewide total aboveground CSP is estimated at 156.3 Tg C in Connecticut, 64.0 Tg C in Delaware, 343.1 Tg C in Maryland, 227.0 Tg C in Massachusetts, 262.7 Tg C in New Hampshire, 1547.5 Tg C in Pennsylvania, 32.4 Tg C in Rhode Island and 270.9 Tg C in Vermont. Statewide CSP of Maine, New Jersey and New York is not reported here because of their incomplete lidar and NAIP data coverage.
The projected maps of CSPG and CSPTG are shown in figures 2(c) and (d), indicating strong spatial variation across the RGGI region. As expected, large sequestration gaps and longer sequestration time gaps generally appear where present AGB is low, such as in eastern and central Maryland, southeastern Pennsylvania, and southeastern Maine. Relatively smaller gaps are located in northcentral Pennsylvania and western Maryland, most of Rhode Island, Connecticut, Massachusetts, Vermont and New Hampshire. Statewide average CSPTG ranges from 207 to 292 years, and correlates to the relative fraction of CSPG to CSP. The longest and shortest CSPTG appear in Delaware and New Hampshire, respectively. Statewide, Delaware is currently at 1/5 of its aboveground carbon sequestration potential, Maryland and Pennsylvania are at about 1/3, Rhode Island, Connecticut, Massachusetts, Vermont and New Hampshire are at about 1/2. The relatively larger CSPG in Delaware is likely due to its low tree cover and young forest. As table 1 indicates, Delaware has the lowest tree cover (i.e. 35.8%) in the region, less than half that of highly forested states such as New Hampshire and Vermont. Despite substantial crop abandonment in 2008 and 2016 (Lark et al 2020), tree cover in DE is only 30% lower than the adjacent state of Maryland, but its AGB density is only half that of Maryland. This difference implies that Delaware forests are younger. The CSPG is also stratified by continued growth areas and regrowth area for each state (shown in figure 5) and each county (shown in figures S8–S15). The stratification in figure 5 indicates that large gaps are primarily located in regrowth area for Maryland, Delaware and Pennsylvania, but primarily in continued growth areas for the other six states.
Figure 5. CSPG over areas with continued growth (green) vs that over regrowth (red) in Maryland, Delaware, Pennsylvania and Rhode Island, Connecticut, Massachusetts, Vermont, New Hampshire and part of Maine.
Download figure:
Standard image High-resolution imageIn addition to aboveground carbon sequestration potential, we also project the potential future growth trajectory of aboveground carbon. Figure 6(a) shows the potential annual aboveground carbon stock from present to 300 years in the future for each state, and Figure 6(b) are corresponding maps of potential AGB in years 2050, 2100, 2200 and 2300, respectively. The relative contribution of AGB from continued growth and regrowth varies by state. For example, the contribution of regrowth to newly gained AGB by 2300 is 46.3% in Connecticut, 70.7% in Delaware, 65.9% in Maryland, 43.0% in Massachusetts, 26.1% in New Hampshire, 51.1% in Pennsylvania, 35.9% in Rhode Island and 34.9% in Vermont. Delaware has the highest contribution from regrowth among all states; the second and third highest contributions are in Maryland and Pennsylvania, respectively. In contrast, continued growth is the largest contributor to the aboveground carbon stocks in all other states. Annual estimates of aboveground carbon stocks from present to 300 years for each county in the RGGI region (except Maine, New Jersey and New York) can be found in figures S16–S23, and county-level aboveground carbon stocks in 2020, 2030, 2040 and 2050 are summarized in tables S1–S8.
Figure 6. ED potential AGB from present to 300 years in the future. Blue line in (a) represents statewide annual aboveground carbon stocks for Maryland, Delaware, Pennsylvania and Rhode Island, Connecticut, Massachusetts, Vermont and New Hampshire. The four numbers along each curve correspond to the stock value at years 2050, 2100, 2200 and 2300. Corresponding maps of AGB density are also mapped in (b). Green and yellow lines in (a) represent the relative contributions of continued growth and regrowth to the carbon stocks.
Download figure:
Standard image High-resolution imageProjections of CSP and CSPTG are further assessed for sensitivity to changes in NPP and disturbance rate. As an example, figure 7 examines how the average CSP and CSPTG across Maryland, Delaware and Pennsylvania, which together account for 50% the land area in the RGGI region, respond to different NPP and disturbance rates. As expected, averaged CSP and CSPTG of Maryland, Delaware and Pennsylvania are generally proportional to NPP and inversely proportional to disturbance. High CSP and CSPTG appear at high NPP but low disturbance, representing conditions where a forest may gain carbon quickly and lose less of it over time to disturbance. In contrast, low NPP and high disturbance result in low CSP and CSPTG because of a slowing carbon sequestration rate and high losses due to disturbance.
Figure 7. Sensitivity of average CSP and CSPTG over the states of Maryland, Delaware and Pennsylvania in response to percent changes in NPP and disturbance rate. NPP and disturbance rates are changed from 50% to 150% at an increment of 20%.
Download figure:
Standard image High-resolution image
4. Discussion
Forests play a crucial role in climate mitigation. Avoided deforestation, improved forest management and reforestation could provide two-thirds of the cost-effective nature-based climate mitigation needed to hold warming to below 2 °C (Griscom et al 2017), with the regrowth of natural forest the single largest natural climate solution both globally and within the United States (Fargione et al 2018, Cook-Patton et al 2020). In this context, accurate high-spatial resolution estimates of the potential for additional aboveground carbon sequestration in forests is needed. This work combined advances in ecosystem modelling and remote sensing to estimate present forest aboveground carbon stocks and project future forest aboveground carbon sequestration potential at 90 m resolution across the 11 states in RGGI region, including 34.9 million 90 m land cells over an area of 283 000 km2.
The RGGI region has large aboveground CSPG compared to its present aboveground carbon stocks. We found that present AGB stocks in Delaware are at approximately one-fifth of their potential, Maryland and Pennsylvania are at one-third, Connecticut, Massachusetts, New Hampshire, Rhode Island and Vermont are at about half. The significant gap between present AGB and aboveground carbon sequestration potential provides opportunities for regional climate mitigation. Maximum potential gains in AGB by 2050 would be 21.2 Tg C in Connecticut, 13.8 Tg C in Delaware, 28.0 Tg C in Massachusetts, 64.2 Tg C in Maryland, 26.6 Tg C in New Hampshire, 257.7 Tg C in Pennsylvania, 4.2 Tg C in Rhode Island, and 32.0 Tg C in Vermont. Together, these eight states have the potential to gain 209.1 Tg C through continued growth of existing trees and 238.6 Tg C through regrowth of new trees. The average annual sequestration rate (i.e. 11.5 Tg C yr−1) is equivalent to 9.6% of these eight states' average annual fossil fuel emissions between 2011 and 2017 (i.e. 119 Tg C yr−1) according to US EPA fossil fuel combustion inventories (U.S. Environmental Protection Agency 2019). The high spatial resolution underlying these region-wide estimates may help decision-makers prioritize areas for reforestation. For example, this work indicates that counties of Lancaster, Crawford, Bradford, Washington and York in Pennsylvania, Sussex in Delaware, and Coos in New Hampshire, where the CSPG is over 20 Tg C, have exceptionally high potential to gain additional forest aboveground carbon. These counties with the high CSPG vary in terms of dominate land cover type. The six Pennsylvania and Delaware counties are predominately pasture and cropland, whereas the New Hampshire county is mostly forest.
Historical and present land-use activities and natural disturbance have resulted in a diverse landscape of heavily fragmented forest, mixed with non-forest patches. While this land cover matrix poses a great challenge for biomass mapping, this study incorporates forest structure data as detailed as 1 m spatial resolution to address this challenge. To illustrate, three areas representing typical forest, residential and agricultural area were chosen as examples. NAIP aerial imagery in figures S5, S6 and S7, respectively show a forested area as a mixture of trees and gaps (black shadows), a residential area as a mixture of houses and trees along roadsides and backyards, and an agricultural area as a mixture of cropland and scattered trees. Across all three heterogeneous landscapes, treed versus non-tree areas are easily identified using high-resolution (1 m) lidar canopy height and tree cover classification. The resulting variations in AGB are consequently well captured in ED initialized AGB and lidar empirical AGB estimates. These fine-scale detailed maps also identify a substantial amount of land carbon being stored in trees that are scattered across non-forest areas, which emphasizes the importance of trees outside of forests to regional and global carbon budgets (Huang et al 2015, Zomer et al 2016, Chapman et al 2020, Spawn et al 2020). These high-resolution capabilities are incredibly valuable for carbon modelling and climate mitigation planning across regions with strong AGB spatial heterogeneity.
This study expands the previous carbon modelling system prototyped for the state of Maryland (Hurtt et al 2019) to a larger domain using an updated version of the ED model and an improved initialization method. We compared new results to prior results for MD, where both studies used the same lidar canopy height and NAIP tree cover data inputs for initialization. This study yields comparable estimates to Hurtt et al (2019). In particular, our results concur with Hurtt et al (2019)'s spatial pattern of sequestration potential and sequestration gap. For example, the Maryland counties of Frederick and Baltimore have the largest CSP (>25 Tg C) in both studies. Further, this study shows agreement with the location of Maryland's largest aboveground carbon stocks, first in Garrett County, followed by Charles, Baltimore, Allegany and Frederick counties, respectively. Although our study includes woody wetland, we have excluded it here to provide a direct comparison to Hurtt et al (2019). This study projects total CSP in Maryland as 304.0 Tg C, which is within 3.5% of CSP estimates in Hurtt et al (2019) (314.8 Tg C). There are, however, differences in contemporary aboveground carbon stock estimates between this study and Hurtt et al (2019); this study estimates Maryland's present carbon stocks at 84.5 Tg C, which is ∼26 Tg C lower than Hurtt et al (2019). The possible causes for this difference in present aboveground carbon stocks include different model driver inputs (e.g. soil depth, air temperature, etc) as well as the change in initialization method. The weighting-based initialization method in this study corrects overestimation at high AGB area and yields improved correlation with FIA plots. For example, switching the initialization method used in this study back to the prior method results in an estimate of aboveground carbon stocks of 94.7 Tg C, implying about 40% of the difference may be due to the initialization method.
We also compared our estimates of natural forest regrowth rate with a recent global estimate from Cook-Patton et al (2020), which combines field measurements across multi-biomes and environmental covariates to produce a global 1 km map of potential AGB growth rates for the first 30 years of natural forest regrowth. We averaged the potential aboveground carbon growth rate between 5 and 30 years to align with Cook-Patton et al. Average growth rates across the RGGI region are similar between this study and Cook-Patton et al, at 0.96 Mg C ha yr −1 and 1.05 Mg C ha yr−1, respectively. Spatial comparisons are shown in figure S25. We find that absolute differences for more than half of all land area (59.2%) are below 0.30 Mg C ha year−1, with 40.5% within 0.20 Mg C ha yr −1 and 20% of land area within 0.10 Mg C ha year−1. Figure S26 (a) shows the average potential growth rates stratified by soil depth, suggesting the strong dependence of AGB growth on soil depth in ED. Most area has a soil depth between 100 cm and 160 cm, and the two carbon estimates are closest at this depth range. The figure S26 (b) shows ED potential growth rates as a function of stand age, suggesting the growth rate varies with age during the first 30 years of natural forest regrowth. This nonlinear age dependence is included in ED estimates but not in Cook-Patton et al (2020).
Understanding uncertainty sources is essential to applications of our estimates of climate mitigation potential and also to provide insights into subsequent modelling development. Uncertainties regarding the estimation of present aboveground carbon stocks arise from ED drivers and initialization inputs. First, due to lack of available meteorology datasets with both high spatial and temporal resolution, the meteorological driver is interpolated by both fine spatial resolution (Daymet) and coarse spatial resolution (MERRA2) data products. Imperfect interpolation, inherent mismatches between two meteorology products or underlying uncertainties within them may propagate errors to ED estimates. Second, in this study, ED is run with average meteorology and CO2 level between 1981 and 2017 despite temporal variability and positive increasing trends. Some states may also have lower disturbance rates than the spatially constant value used in ED (1.2% yr−1). These factors may jointly lower AGB growth rates and aboveground carbon sequestration potential, especially if these average rates change over time. Third, the uncertainty associated with initialization inputs is inconsistent acquisition times of the lidar canopy height data. Our lidar data comes from acquisition between 2004 and 2018, and some archived lidar data may have less sampling density than others (Huang et al 2019, Tang et al 2020).
Uncertainties regarding the estimation of future forest carbon stocks arise from future climate change, CO2 levels, and the disturbance regime. Future climate change and continued increases in CO2 levels will likely have complicated impacts on vegetation growth. For example, field experiments found that increasing CO2 could enhance tree growth indicating more carbon will be sequestered (Norby and Zak 2011, Walker et al 2019). However, faster growth may also accelerate harvesting events, reducing tree longevity, and resulting in less carbon residence time in trees (Körner 2017). Warming could also lengthen growing seasons and in turn enhance annual vegetation growth (Menzel and Fabian 1999); however, it may also inversely increase evapotranspiration and autotrophic and heterotrophic respiration (White et al 1999, Piao et al 2008, Crowther et al 2016), decrease soil moisture and increase fire risk (Westerling et al 2006). Additionally, future disturbance is likely to increase as warming continues. Warmer and drier conditions facilitate disturbance related to fire, drought and insect outbreaks and decrease disturbance related to snow and ice, while warmer and wetter conditions increase disturbances from wind and pathogens (Seidl et al 2017). Given the unclear net impact of climate change and lack of climate change scenarios harmonized at high spatial resolutions, this study does not explicitly consider changes in future climate, CO2 and disturbance, but alternatively investigates how potential alterations in NPP and disturbance rate could propagate to estimates of future carbon sequestration.
Finally, this study has focused on forest aboveground carbon stocks from present to future. The work has capitalized on the qualitative advances in data and modelling of aboveground forest structure, and focused on producing new high-resolution estimates of above ground carbon that are most directly related to these measurements and also important for policy. This work does not yet include other forest carbon pools such as soil carbon or wood products, nor the dynamics related to forest management (Birdsey et al 2006, Lippke et al 2011). Our projections of CSP and CSPG are thus conservative underestimates of ecosystem's total carbon sequestration potential (Domke et al 2020a, 2020b), and future work could build off these advances above ground to address these other components. Looking ahead, the forest carbon modelling approach used in this study is currently being expanded to US CONUS and global domains and to include soil carbon and associated carbon fluxes. Key to expanding this work's spatial coverage is to switch from airborne remote sensing to orbital observations. Despite the demonstrated capability of 1 m airborne lidar data to capture fine-scale heterogeneity, leveraging this existing data at national and global scales encounters several challenges, including inconsistency of instrument quality and acquisition time as well as limited spatial coverage. New airborne lidar collection can also be expensive. Future work should utilize satellite-based forest structure measurements from the ongoing NASA missions of Global Ecosystem Dynamics Investigation (Dubayah et al 2020), Ice, Cloud, and land Elevation Satellite-2 (Markus et al 2017) and Landsat-based tree cover from Global Forest Change dataset (Hansen et al 2013). Key to accounting soil carbon will be sufficient data to constrain model estimates and account for the effects of disturbance and land use history. Relevant soil carbon inputs include the Harmonized World Soil Database (Wieder et al 2014) and SoilGrids250m dataset (Hengl et al 2017). Relevant products for disturbance and land use history include the Global Forest Change dataset (Hansen et al 2013), the North American Forest Dynamics dataset (Goward et al 2016), and the Land Use Harmonization 2 dataset (Hurtt et al 2020).
Acknowledgments
This work was funded by NASA Carbon Monitoring System (NASA-CMS) project (80NSSC17K0710), NASA New Investigator Program (NASA-NIP) (80NSSC18K0708), and also supported by the University of Maryland Ann GWylie Dissertation Fellowship. We also acknowledge in-kind support from the US Forest Service, Forest Inventory and Analysis Program.
Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.