Modeling climate-smart forest management and wood use for climate mitigation potential in Maryland and Pennsylvania

State and local governments are increasingly interested in understanding the role forests and harvested wood products play in regional carbon sinks and storage, their potential contributions to state-level greenhouse gas (GHG) reductions, and the interactions between GHG reduction goals and potential economic opportunities. We used empirically driven process-based forest carbon dynamics and harvested wood product models in a systems-based approach to project the carbon impacts of various forest management and wood utilization activities in Maryland and Pennsylvania from 2007–2100. To quantify state-wide forest carbon dynamics, we integrated forest inventory data, harvest and management activity data, and remotely-sensed metrics of landuse change and natural forest disturbances within a participatory modeling approach. We accounted for net GHG emissions across (1) forest ecosystems (2) harvested wood products, (3) substitution benefits from wood product utilization, and (4) leakage associated with reduced instate harvesting activities. Based on state agency partner input, a total of 15 management scenarios were modeled for Maryland and 13 for Pennsylvania, along with two climate change scenarios and two bioenergy scenarios for each state. Our findings show that both strategic forest management and wood utilization can provide substantial climate mitigation potential relative to business-as-usual practices, increasing the forest C sink by 29% in Maryland and 38% in Pennsylvania by 2030 without disrupting timber supplies. Key climate-smart forest management activities include maintaining and increasing forest extent, fostering forest resiliency and natural regeneration, encouraging sustainable harvest practices, balancing timber supply and wood utilization with tree growth, and preparing for future climate impacts. This study adds to a growing body of work that quantifies the relationships between forest growth, forest disturbance, and harvested wood products. 

Methods
Primary data inputs to CBM-CFS3 included a detailed forest inventory, growth-yield relations to estimate forest productivity, and estimates of harvest yields and intensity, land-use change, and natural disturbances. Inventory data are categorized by a series of forest classifiers defining relevant characteristics such as spatially referenced boundaries, ownership, forest type, site productivity, or reserve status. Allometric equations are used to predict tree volume-to-biomass relationships (Boudewyn et al., 2007). For this study, forest inventory, growth-yield curves, and harvest data were estimated from the USDA Forest Service Forest Inventory and Analysis (FIA) program, which we accessed through the FIA DataMart (USDA Forest Service, 2019) using the rFIA package (Stanke et al., 2020) in the R programming environment (R Core Team, 2020), which enables data exploration and user-defined spatio-temporal queries and estimation of the FIA database (FIADB). Methodologies derived from Bechtold & Patterson, 2005 and Pugh et al., 2018 were used to estimate each state’s forest inventory by a predetermined list of classifiers. Natural disturbance history was estimated from both the FIADB and LANDFIRE (USGS, 2016) datasets to better constrain initial belowground and soil carbon parameters during what the modeling framework refers to as the model spin-up period (Kurz et al., 2009). Estimates of merchantable volume and corresponding biomass from FIADB were used to calibrate the model’s allometric volume-to-biomass assumptions to match forest type groups and growth conditions in Maryland and Pennsylvania.

Harvest removals were estimated as average annual removal of merchantable timber in cubic feet between 2007 and 2019, converted to metric tons of carbon using methodologies and specific gravities reported by Smith et al., 2006. To assign a harvest type and intensity to each record of volumetric removal, stand age at the time of removal was calculated by taking the mid-point average between time t1 and t2 (Bechtold & Patterson, 2005) where t1 is the year the unharvested stand was measured and t2 is the repeat interval year measurement post-harvest. In collaboration with state partners, harvest type and intensity were determined heuristically for each forest type based upon state-level management documentation, peer-reviewed literature, and expert input. A complete list of harvest types and intensities prescribed to each forest type group as a product of stand age can be found in Table S3 of the manuscript.

Longer-term averages from 2007-2019 were used to estimate annual area targets for all land-use change (LUC) and natural disturbance events including wind, fire, disease, and insects. Annual LUC average rates by ownership and forest type group were derived by overlaying a geospatial forestland ownership dataset (Sass et al., 2020), the Protected Areas Database of the U.S. (PAD-US), a national geodatabase of protected areas (USGS, 2018), and the National Land Cover Database (NLCD), a remotely-sensed data product used to characterize land cover and land cover change (Wickham et al., 2021). Wind disturbance events were calculated using the LANDFIRE Historic Disturbance dataset (USGS, 2016), a remotely-sensed data product provided by the USGS that estimates annual disturbance events. Annual averages for wildfire disturbances were derived from the LANDFIRE Historic Disturbance dataset (USGS, 2016) and validated through annual reports from the National Interagency Fire Center (NIFC). Annual prescribed fire acres were estimated from reports provided by the Maryland DNR Forest Service and Pennsylvania DCNR Bureau of Forestry and scaled to represent treatments on forestlands only. Annual acreages of insect and disease disturbance were derived from the National Insect & Disease Detection Survey (USDA Forest Service, 2020), a spatial data product produced by USDA that collects and reports data on forest insects, diseases, and other disturbances. For more information on all input and activity data see Supplementary Materials 1.2. A complete list of BAU parameters can be found in Table S2.

State-specific trade and commodity data from Resource Planning Act (RPA) assessments (USDA Forest Service, 2021), US Commodity Flow Surveys (US Department of Transportation et al., 2020), US International Trade Commission export data (US International Trade Commission, 2021), and published peer-reviewed data (Howard & Liang, 2019) when available, or US averages from the same sources, were used to adapt and parameterize both HWP models. FAOSTAT data (FAO, 2021) were utilized to determine the commodity distributions of exported roundwood. Softwood products were parameterized and modeled separately from hardwood products, as the two wood types differ in exports and commodities produced as well as their associated product half-lives and displacement (Dymond, 2012; Howard et al., 2017). Published data were used to calculate softwood- and hardwood-specific half-lives for Maryland and Pennsylvania sawn wood and veneer products, while we relied on literature estimates for other products (Skog, 2008; J. E. Smith et al., 2006b). To calculate substitution benefits, we coupled region-specific data (USDA Forest Service, 2021), US consumption rates (Howard et al., 2017), product weights (C. Smyth et al., 2017), and LCA data (Bala et al., 2010; Dylewski & Adamczyk, 2013; Hubbard et al., 2020; Meil & Bushi, 2013; Puettmann, 2020a, 2020b; Puettmann et al., 2016; Puettmann & Salazar, 2018, 2019), following methods developed by Smyth et al. (2017). Landfill CO­2 and CH4 emissions rely on PICC defaults for methane generation (k) and landfill half-lives for wet, temperate climates (Pingoud et al., 2006). See Supplementary Materials 1.3 for more details on substitution and leakage calculation methods.