ABSTRACT Title of Dissertation: FINE RESOLUTION ASSESSMENT OF THE CARBON FLUXES FROM CONTEMPORARY FOREST DYNAMICS Weishu Gong, Doctor of Philosophy, 2022 Dissertation directed by: Professor Chengquan Huang, Department of Geographical Sciences Current estimation of the Earth?s carbon budget contains large uncertainties, with the largest ones in its terrestrial components. With an overarching goal to improve the understanding of carbon budget at regional to global scales, this study aimed to: 1. Develop a grid-based carbon accounting (GCA) model for estimating carbon fluxes from forest disturbance, tested over North Carolina; 2. Develop a consistent timber product output (TPO) record for a globally important timber production region, including seven states in the southeast U.S.; and 3. Further improve the GCA model based on results from objectives 1 and 2, and use it to derive carbon source/sink estimates for all forest land in North Carolina. The results show that several inputs/parameters such as pre-disturbance carbon density, disturbance intensity, allocation of removed carbon among slash and different wood product pools, and forest growth rates could have large impact on carbon estimates. The total emission between 1986 and 2010 from logging over North Carolina was reduced by one third and two thirds, respectively, when remote sensing-based disturbance intensity and biomass data were used to replace parameter values inherited from the original bookkeeping carbon accounting (BCA) model, and was reduced by over 70% when both were used. Use of the TPO data derived in Chapter 3 to partition the removed carbon among slash and different wood product pools resulted in noticeably higher emission estimates than those derived using the partitioning ratios provided by the original BCA model. In addition, without considering legacy effect from wood products harvested before 1986, the emission value derived using the prompt release method was 50% higher than that derived using the delayed release method. This study addresses multiple sources of uncertainties related to the terrestrial carbon budget. The TPO study demonstrated an approach for deriving consistent TPO records for large timber production regions. The GCA model produced state level carbon estimates comparable to those reported by the U.S. Forest Service while providing critical spatial details needed to support carbon management and advance forest-driven climate change mitigation initiatives. FINE RESOLUTION ASSESSMENT OF THE CARBON FLUXES FROM CONTEMPORARY FOREST DYNAMICS by Weishu Gong Dissertation submitted to the Faculty of the Graduate School of the University of Maryland, College Park, in partial fulfillment of the requirements for the degree of Doctor of Philosophy 2022 Advisory Committee: Professor Chengquan Huang, Chair Professor Shunlin Liang, Co-Chair Professor Matthew Hansen Doctor Richard A. Houghton Professor Joseph Sullivan ? Copyright by Weishu Gong 2022 Dedication To my parents, Jianya Gong and Jiansi Yang, who were always there for me when needed during the entire length of my PhD study. I wouldn?t be able to make it without your support, love, and understanding. To my beloved cat, Christoph, who, despite not knowing a thing about science or research, is still an important presence in my life just by being his usual adorable self. Thanks for being a clingy little fur ball that pesters me all the time for food and affection. AND Above all, to all my friends and family, who have been cheering me on for years. Thank you for your love and support. ii Acknowledgements This dissertation was written under the supervision of Dr. Chengquan Huang, who had provided me immense support during the entire process. Thank you for your guidance. This dissertation would not have been possible without you. I would also like to thank Dr. Houghton, whose bookkeeping carbon accounting model had served as the very foundation for this research, and who had engaged in so many discussions with me and my advisor to help improve my research. Thank you for all your feedback. The co-chair of my doctoral advisory committee, Dr. Shunlin Liang, has been mentoring me since I was an undergraduate student. Thank you for your years of advice on career planning, research, and life in general. The other members of my committee, Dr. Matthew Hansen and Dr. Joseph Sullivan, had both provided me highly valuable suggestions and feedback during the different stages of my dissertation research. Thank you both, Drs. Hansen and Sullivan. My previous and current coworkers, Dr. Pui-Yu Ling, Dr. Xin Tao, Dr. Feng ?Aaron? Zhao, and Ms. Jiaming Lu had provided me essential datasets that were used during the research. Dr. Ling assisted me with the TPO survey data and walked me through the basics of TPO modeling, Dr. Tao and Ms. Lu provided me the disturbance intensity dataset, while Dr. Zhao provided me the VCT disturbance dataset. Dr. Jack Ma, the MS program director, provided critical support in times of need. Thank you all for the help. The US Forest Service had provided the FIA inventory data that were crucial to many aspects of this research, including the calculation of forest growth rates and the derivation of the pre-disturbance carbon density map. My doctoral study was supported by the Department of Geographical Sciences at the University of Maryland, the Environmental Model and Data Optima (EMDO) Company, and the PIESAT Information Technology (Australia) Company. Additional assistantships were provided by projects funded by NASA, NOAA, the United States Geological Survey (USGS), and the United States Forest Service (USFS). iii Table of Contents Dedication ..................................................................................................................... ii Acknowledgements ...................................................................................................... iii Table of Contents ......................................................................................................... iv List of Tables ............................................................................................................... vi List of Figures ............................................................................................................. vii Chapter 1 Introduction .................................................................................................. 1 1.1 Forest Change and the Global Carbon Budget ................................................... 1 1.2 Carbon Estimation Approaches .......................................................................... 3 1.3 Remote Sensing of Forest Carbon Dynamics ..................................................... 6 1.4 Research Objectives and Dissertation Structure ............................................... 10 Chapter 2 Carbon fluxes from contemporary forest disturbances evaluated using a grid-based carbon accounting model .......................................................................... 14 2.1 Background ....................................................................................................... 14 2.2 Methods............................................................................................................. 17 2.2.1 Study Area ................................................................................................. 17 2.2.2 The original Bookkeeping Carbon Accounting (BCA) model .................. 17 2.2.3 Development of A Grid-Based Carbon Accounting (GCA) model ........... 20 2.2.4 Model inputs .............................................................................................. 22 2.2.5 Carbon modeling scenarios ........................................................................ 28 2.3 Results ............................................................................................................... 29 2.3.1 Model verification and sensitivity analysis................................................ 29 2.3.2 Carbon fluxes from disturbance and post-disturbance growth .................. 33 2.4 Discussion ......................................................................................................... 38 2.5 Conclusions ....................................................................................................... 40 Chapter 3 Estimation of Annual Harvested Wood Products based on Remote Sensing and TPO Survey Data ................................................................................................. 42 3.1 Introduction ....................................................................................................... 42 3.2 Study Area and Data Collection ....................................................................... 44 3.2.1 Study Site ................................................................................................... 44 3.2.2 Data collection ........................................................................................... 46 3.3 Methods............................................................................................................. 52 3.3.1 Preparation of predictor variables .............................................................. 52 3.3.2 Modeling approaches ................................................................................. 54 3.4 Results ............................................................................................................... 55 3.4.1 Limitations of the FELR method ............................................................... 55 3.4.2 Performance of the OLS and RF for TPO modeling ................................. 58 3.4.3 Spatial-Temporal Variability of TPO in the Southeast .............................. 60 3.5 Discussion and Conclusions ............................................................................. 65 Chapter 4 GCA Model Improvement and Implications for Carbon Management and Flux Estimation over North Carolina .......................................................................... 69 4.1 Introduction ....................................................................................................... 69 4.2 GCA Model Improvement ................................................................................ 70 4.2.1 Partitioning of Removed C ........................................................................ 70 iv 4.2.2 C accumulation from forest growth ........................................................... 72 4.2.3 Model verification and adjustment ............................................................ 73 4.3 North Carolina?s Forest Carbon Dynamics....................................................... 78 4.3.1 Source and Sink Analysis .......................................................................... 78 4.3.2 Impact of management approaches and emission calculation methods ..... 81 4.4 Summary and Conclusions ............................................................................... 86 Chapter 5 Conclusions ................................................................................................ 88 5.1 Summary of Major Findings ............................................................................. 88 5.2 Significance of Original Contributions ............................................................. 91 5.3 Limitations and Future Research Directions ..................................................... 93 Bibliography ............................................................................................................... 96 v List of Tables Table 2-1 Pools tracked by the BCA model for LULCC types considered in this study (indicated by ?X?). ...................................................................................................... 18 Table 2-2 Ecozone-specific parameters used by the BCA model ............................... 18 Table 2-3 Growth rate statistics from selected FIA plots. The last two columns show the values provided with the original BCA model...................................................... 27 Table 2-4 Scenarios for assessing model sensitivity to spatial carbon density and disturbance intensity data. ........................................................................................... 29 Table 3-1 Years in which survey data are available for each state, "x" indicates the data is available for that year in that state. .................................................................. 48 Table 3-3-2 Average annual timber production volume (in thousand cubic meters) in each state based on the survey data listed in Table 3-1. ............................................. 49 Table 3-3 IFZ Disturbance magnitude (IFZ-DM) levels at which county level disturbance areas were calculated. .............................................................................. 53 Table 3-4 R2 values of the relationships between TPO values predicted by OLS and survey data derived through 10-fold cross-validation. ............................................... 59 Table 3-5 R2 values of the relationships between TPO values predicted by RF and survey data derived through 10-fold cross-validation. ............................................... 59 Table 3-6 Estimated 30-year (1986-2015) total TPO (thousand m3, or kcume) for the 7 states in the southeast United States. ....................................................................... 61 Table 4-1 Scenarios for assessing the impact of logging intensity and emission expression methods on forest C outcome estimates. .................................................. 82 vi List of Figures Figure 1-1 Conceptual diagram of global carbon dioxide budget (Le Qu?r? et al., 2018). ............................................................................................................................ 2 Figure 1-2 The overarching goal and three specific research objectives of this dissertation study. ....................................................................................................... 11 Figure 1-3 Study regions of the three research objectives of this dissertation study. . 13 Figure 2-1 The process of determining disturbance attribution. ................................. 24 Figure 2-2 Comparison of carbon density values between plot measurements and derived maps. After filtering the original map with 95th or 99th percentile in a 6 km by 6 km window, the range of values are closer to those of undisturbed plots. ......... 26 Figure 2-3 Distribution of forest growth rates calculated from FIA plot data. ........... 27 Figure 2-4 Comparison of carbon fluxes by disturbance type between different scenarios. Using both RS-based carbon density map and disturbance intensity data has turned the 25-year net total C flux from source to sink. ....................................... 31 Figure 2-5 Annual total carbon absorbed and released by disturbance types during the study period under the four different scenario settings. The trendlines show large differences among the sources under the four scenarios. ........................................... 33 Figure 2-6 A browse image of the 30-m net flux map derived by summing up all source and sink terms over the study period for each pixel. The full resolution zoon-in examples show the flux patterns over areas with (a) extensive forest-to-urban conversion, (b) a large fire occurred in 2008, (c) active wood harvesting followed by strong growth throughout the study period, (d) minimal disturbances within national park/national forests, and (e) storm damage/salvage logging followed by recovery. 34 Figure 2-7 Total carbon absorption and emission by disturbance types at county level during the 25-year study period. Wood harvesting activities are prevalent in the eastern part of the state, resulting in both large sources and sinks. Only a handful of counties are affected by fire. Emission from conversion to urban land is prominent in counties with major urban centers. Overall, the emission is less than or almost equal to the absorption in all counties. ................................................................................. 37 Figure 3-1 Study area located in the US southeastern States, including North Carolina, South Carolina, Alabama, Florida, Georgia, Mississippi, Tennessee. ........ 46 Figure 3-2 A forest disturbance event typically results in an abrupt increase in IFZ, which decreases gradually due to post-disturbance recovery (Huang et al., 2015). Disturbance year is determined by the timing of the abrupt change and the change in the IFZ value provides a spectral measure of the disturbance intensity. .................... 51 Figure 3-3 Temporal mismatch between survey year and potential disturbance event (Huang et al. 2015). .................................................................................................... 54 Figure 3-4 Comparison of model estimates (y-axis) derived using the fixed effects linear regression (FELR) method over South Carolina (top) and Tennessee (bottom) and TPO survey data. Each point represents a county. The orange lines are 1:1 lines. ..................................................................................................................................... 57 Figure 3-5 Abnormally high or low values were predicted in many years by the FELR method in South Carolina (left) and Tennessee (right), which were obviously wrong. Such obvious errors were also found in Florida, Georgia, and Mississippi. .............. 58 vii Figure 3-6 Coefficient of determination (R2) values derived by comparing survey data and values predicted by OLS (y-axis) and Random Forest (x-axis) through 10-fold cross validation. .......................................................................................................... 60 Figure 3-7 County-level total TPO values predicted by the Random Forest for the 30- year study period (1986-2015). ................................................................................... 62 Figure 3-8 Comparison of the average annual TPO values calculated from the 30-year record derived through this study and available TPO survey data. Each point represents one wood product type in one state. The two right most points were from Georgia. ....................................................................................................................... 63 Figure 3-9 Annual timber production volume in softwood, hardwood, and total for seven states: (a) North Carolina; (b) South Carolina; (c) Alabama; (d) Florida; (e) Georgia; (f) Mississippi; (g) Tennessee. ..................................................................... 64 Figure 3-10 Total C influx to the P10 (C released in 10 years) and P100 (C released in 100 years) pools calculated using the derived TPO record over the 30-year study period. The blue ellipsoids highlight the counties that had more harvested C allocated to the P10 pool than the P100 pool, or vice versa. ...................................................... 67 Figure 4-1 Comparison of annual C changes in the aboveground (including harvested wood products) pool derived using the GCA model and from the USFS inventory data over forest land remaining forest land. The slope of the increasing trend in the annual flux estimates by the GCA model was reduced after adjusting the estimates to account for legacy emission from wood products harvested during 100 years before each model year (see the text for details on the adjustment method). ........................ 74 Figure 4-2 Because the GCA model did not include legacy emissions from wood products harvested before 1986 and started to accumulate the wood product pool in 1986, model estimates of legacy emission from wood products increased as the number of years since 1986 increased, roughly following a logarithmic curve. ........ 75 Figure 4-3 Comparison of annual C changes in the aboveground (including harvested wood products) pool derived using the GCA model with growth rates adjusted based on FIA plot data and from the USFS inventory data over forest land remaining forest land. After adjusting the estimates to account for legacy emission from wood products harvested during 100 years before each model year, the annual sink estimated by the GCA model became slightly smaller than the inventory-based estimates. ..................................................................................................................... 77 Figure 4-4 Compared to the initial GCA model developed in Chapter 2, improvements made in this chapter resulted in higher source and lower sink estimates for forest land disturbed by logging. The same was true for conversion to urban land (no sink) and forest land subject to fire disturbances. ................................................ 78 Figure 4-5 Contributions of different source and sink processes to the total source and sink over North Carolina from 1986 and 2010. .......................................................... 80 Figure 4-6 Annual carbon emissions derived using the delayed (top) and prompt (bottom) emission methods at the four disturbance intensity levels shown in Table 4- 1................................................................................................................................... 84 Figure 4-7 Total fluxes over the 25-year study period derived using the delayed (top) and prompt (bottom) emission methods at the four disturbance intensity levels shown in Table 4-1. ................................................................................................................ 85 viii Chapter 1 Introduction 1.1 Forest Change and the Global Carbon Budget Numerous studies have pointed out that recent anthropogenic carbon emissions are the highest in human history and are extremely likely the main cause of the rise in atmospheric CO2 concentration that leads to the observed global warming since the mid-20th century (IPCC, 2014). Monitoring and assessing the global carbon cycle is therefore of utmost importance for national and international decision makers when considering climate mitigation strategies (Griscom et al., 2017; Marland et al., 2003; Noormets et al., 2015). However, there are many uncertainties in the current estimation of the Earth?s carbon budget, with the largest ones in its terrestrial components (Houghton, 2013; Antonarakis, 2014). According to a recent study (Figure 1-1, Le Qu?r? et al., 2018), the global carbon budget for the decade of 2007-2016 had a budget imbalance (or missing sink) of 0.6 GtC yr?1. However, uncertainties in estimating the land sink and emissions from land use change were 0.8 GtC yr?1 and 0.7 GtC yr?1, respectively. The combined uncertainties from these two flux estimates were 250% of the missing sink. Reducing these uncertainties is crucial for improved characterization of the missing sink and the global carbon budget. Understanding the mechanisms responsible for this missing sink is crucial for accurate projections of future concentrations of CO2 in the atmosphere (Houghton & Davidson, 1998) and the change in the Earth?s climate. 1 Figure 1-1 Conceptual diagram of global carbon dioxide budget (Le Qu?r? et al., 2018). Covering over 4.1 billion hectares, forests contain over 80% of aboveground and 40% of belowground terrestrial carbon, which is more than any other terrestrial ecosystem (Bradford et al., 2008; Gray & Whittier, 2014; Liu et al., 2008). Relatively minor alterations to carbon storage or cycling in forest ecosystems may have substantial impact on atmospheric carbon dioxide concentrations. Forest carbon dynamics are governed in large part by disturbance and subsequent regrowth processes. While forest growth provides a mechanism for transferring atmospheric carbon to the forest ecosystem, through disturbance forest carbon is released to the atmosphere or transferred to other pools (e.g., slash, wood products) where carbon is released 2 gradually. Along with other forest characteristics, disturbances regulate the amount of carbon stored in forests and influence C fluctuations overtime (Dugan et al., 2015; Goetz et al., 2012; Pan et al., 2011; Powell et al., 2010). A major goal of forest carbon management is to increase carbon sequestration and storage by forests and related carbon pools (e.g., wood products) while reducing and/or delaying carbon emission from those pools (Canadell and Raupach, 2008; Denise et al., 2011), which has been proposed as an important climate change mitigation strategy in many climate initiatives, including the Kyoto Protocol and the United Nations Framework Convention on Climate Change REDD+ (UNFCCC-REDD+, with REDD standing for ?reducing emissions from deforestation and forest degradation?, and the + standing for ?the role of conservation, sustainable management of forests and enhancement of forest carbon stocks?) program (Agrawal et al., 2011; Schulze et al., 2002). Implementing this strategy requires robust carbon accounting systems to provide reliable estimates of carbon credits for carbon trade and to support the measurement, reporting, and verification (MRV) of carbon pools and fluxes (Birdsey et al., 2006; Fahey et al., 2010; Lamb et al., 2021). 1.2 Carbon Estimation Approaches Many methods have been developed for deriving carbon estimates at regional to global scales, which can be grouped into four categories. The first is the inverse models used with variations in atmospheric concentrations of CO2 to infer sources and sinks (Fan et al., 1998; Pacala et al., 2001; Crevoisier et al., 2007). These studies have generally estimated larger net sinks than other approaches, but the uncertainties are 3 high. The second is to use process-based ecosystem models or Terrestrial Biosphere Models (TBMs) to simulate changes in carbon storage on land (Tian et al., 1999; Hurtt et al., 2002). The estimates of annual Net Ecosystem Productivity (NEP) from TBMs are also highly variable, with major components of NEP (GPP and Rs) having even larger among-models variability (Huntzinger et al., 2012). The third approach relies on data from forest inventories to calculate carbon budgets (Birdsey & Heath, 1995; Turner et al., 1995a; Smith et al., 2007; Pan et al., 2011a). Ground-based inventory of the nation?s forests has been provided by the US Forest Service Forest Inventory and Analysis (FIA) program. FIA forest monitoring includes three phases of measurements conducted approximately every five years (Smith, 2002). Measurements at each forested plot include forest type, tree species, tree size (diameter, height), and stand age, among others. Each plot represents an area of about 1 ha. The biomass of each individual tree is calculated from tree size measurements using allometric equations (e.g. Jenkins et al., 2003). Tree level biomass estimates are then summarized to calculate plot level biomass density. Since 2000, FIA plots in the eastern US have been inventoried once every 5 years and those in western US once every 10 years. Estimates from inventory data are less variable than estimates from inverse studies and ecosystem models, and generally the estimated net carbon sink in forests is also smaller. However, inventory data are based on sampling strategies designed to estimate stocks (growth and mortality) at the national and state levels, not to monitor changes at sub-state levels (Bradford et al., 2010). Use of inventory data alone, therefore, is difficult to support forest carbon management and the measurement, reporting, and verification of carbon pools and fluxes at scales where local management 4 decisions are made typically. The fourth group utilizes land use change data, including changes in agricultural and urban lands, forest harvest, and wildfires to estimate the net carbon fluxes attributable to these changes (Houghton et al., 1999; Zheng et al., 2011). The bookkeeping model (Houghton, 1999) is a carbon estimation approach that takes land use / land cover change (LULCC) and biomass into consideration. The model keeps track of the carbon in four major pools: living aboveground and belowground biomass; dead biomass, including coarse woody debris; harvested wood products; and soil organic carbon (Houghton and Nassikas, 2017). Four types of land use / land cover changes are considered: forests disturbance by fire, industrial wood harvest, conversion from forest to cropland, and conversion from forest to urban land. While this model directly addresses the land sink and fluxes driven by land use change, it can only produce estimates at the ecozone level, which lack critical spatial details needed to support carbon management decision making by local authorities and individual land owners. A few synthesis studies have reviewed a combination of different methods and data sets (Schimel et al., 2000; Pacala et al., 2001; CCSP, 2007; Hayes et al., 2012). In addition, several C budget models have been developed to model ecosystem responses to climate drivers and other disturbances, and these models represent an established method for estimating C fluxes on a national to regional scale. For example, the Canadian national forest carbon monitoring accounting and reporting system (NFCMARS) uses the carbon budget model of the Canadian forest sector (cbm-cfs3) and is used as a decision support tool for forest managers to quantify forest carbon 5 dynamics on a landscape scale. The Australian National Academy of Sciences has developed the National Carbon Accounting System (NCAS), which calculates annual Carbon fluxes (Richards, 2020). 1.3 Remote Sensing of Forest Carbon Dynamics Understanding of the land sink and carbon fluxes driven by land use change requires reliable information on key land surface characteristics, including vegetation carbon and land use change history. Lack of such information can lead to large uncertainties in modeling disturbance rates and biomass densities, which are among the major factors contributing to the widely varying results derived from different studies (Huntzinger et al., 2012). With the ability to cover large areas in relatively short time, satellite remote sensing has been used to observe the dynamics of the vegetation cover on Earth. In the past few decades, remotely sensed data have not only contributed to increasing the speed, cost efficiency, precision, and timeliness associated with inventories, but they have facilitated construction of maps of forest attributes with high spatial resolutions (McRoberts & Tomppo, 2007). Among the various earth observing systems, the Landsat satellite series is arguably the most valuable for monitoring land surface characteristics. With its first satellite launched in 1972, Landsat has produced a fine resolution imagery record of the Earth?s surface for half a century. Landsat data has been used to map land cover and various surface characteristics in numerous studies (Loveland and Dwyer, 2012; Wulder et al., 2019). The multi-decadal time series observations provided by Landsat are especially valuable for understanding the carbon dynamics related to land use change over the last 6 few decades. The opening of the Landsat archive in 2008 for public access ushered new opportunities for mapping land change history using time series Landsat observations. LandTrendr (Kennedy et al., 2007; Kennedy et al., 2010), vegetation change tracker (VCT) (Huang et al., 2010), and continuous change detection and classification (CCDC) are among the various algorithms designed to leverage time series Landsat observations for land change studies. Since 2003, Goward and colleagues (Huang, Masek, Cohen, Moisen and others) have led efforts to characterize US forest disturbances using time series of Landsat observations (Goward et al., 2002; Goward et al., 2005; Goward et al., 2007; Goward et al., 2010). These studies became known as the North American Forest Dynamics (NAFD) study and was identified as a core project of the North American Carbon Program (Goward et al., 2008). Major NAFD products include forest disturbances mapped at an annual time step. These products were derived using the VCT algorithm and a 30-year surface reflectance record (Zhao et al., 2018). VCT detects anomalous events in the per-pixel spectral time series caused by forest disturbances, including harvest/logging, fire, storm damages and insect outbreaks. With biennial Landsat observations, VCT detected stand clearing and partial disturbances with an accuracy of 92% and 60% respectively (Thomas et al., 2011). The Landsat time series also provided opportunities for deriving other important biophysics related to forest carbon dynamics. Building on the NAFD disturbance products, Schleeweis et al. (Schleeweis et al., 2020) mapped forest disturbance types/causal agents across the conterminous U.S. (CONUS) from 1986 to 2010. By using field measurements collected by the Forest Inventory and Analysis (FIA) 7 program of US Forest Service as calibration data, Tao et al. (Tao et al., 2019) estimated the percentage of basal area removal (PBAR) as a measure of disturbance intensity for disturbance events detected by the VCT algorithm over North Carolina. Building on that study, a CONUS-wide disturbance intensity dataset was developed recently (Lu et al. in review). Further, the VCT disturbance products made it possible to derive an annual, multi-decadal (1986-2015) record on the timber product output (TPO) from forest harvest (Huang et al., 2015), and the influx of harvested carbon to wood products (Ling et al., 2016). Compared to the TPO survey data that were only available for a few selected years (Huang et al. 2015), the Landsat-based TPO record was much longer and provided estimates for every year. In addition to the NAFD products, several other geospatial datasets have been developed to provide information on specific disturbance types at national or sub- national scales. In particular, the Monitoring Trends in Burn Severity (MTBS) project (Eidenshink et al., 2007), a collaboration between the US Forest Service and USGS, has mapped the extent and severity of large fires across the United States using field records and Landsat images acquired from 1984 to present. The US Forest Service have been producing Aerial Detection Survey (ADS) sketch maps recording the location (polygons) of insect outbreaks, which can be used to produce consolidated data products on insect related mortality (Williams & Birdsey, 2003; Meddens et al., 2012). Hurricane and tornado tracks have been recorded by NOAA as early as 1851. Combining these tracks with wind model allows assessment of wind damages from tropical storms (Zeng et al., 2009). Biomass is one of the most valuable products for carbon studies, because it can be 8 directly related to landscape carbon, and can provide a constraint for both growth models and calculations of emissions from disturbance (Houghton, 2005; Houghton et al., 2009). Of the three remote sensing instrument types ? optical, radar, and LiDAR, LiDAR can provide metrics that are directly related to forest structure and height (Dubayah and Drake, 2000; Lefsky et al., 1999), and hence has the best potential for biomass estimation (Lefsky et al., 2002; Nelson et al., 2017). However, current spaceborne LiDAR systems, including ICESAT-2 and the Global Ecosystem Dynamics Investigation (GEDI), can only sample along their tracks. They cannot provide wall- to-wall observations needed to create spatially contiguous map products. In general, radar is more sensitive to vegetation structure than optical systems. While radar offers promise for predicting forest biomass and for mapping general forest types and tree species in floristically simple landscapes (Neeff et al., 2005; Pulliainen et al., 2003; Rauste, 2005a; Rignot et al., 1994), radar signal saturates at mid- to high- biomass levels, with the location of the saturation point being wavelength dependent (Austin et al., 2003; Balzter et al., 2003a, b; Castel et al., 2002; Fransson et al., 2000; Patenaude et al., 2005; Rauste, 2005a, b; Santos et al., 2003a; Santos et al., 2003b; Schroder et al., 2005). Although optical remote sensing data may not be as sensitive to forest structure and biomass as lidar and radar data (Franklin et al., 2003; Donoghue & Watt, 2006), they proved to be useful for mapping forest biomass over large areas (Bradford et al., 2008; (Powell et al., 2010). Calibrated using field measurements collected by the FIA program or LiDAR measurements from the ICESAT-1 mission, Landsat and MODIS observations were one of the primary inputs in developing biomass maps for the US 9 (Blackard et al., 2008; Wilson et al., 2013) and the tropics (Saatchi et al., 2011). Given the large quantities of LiDAR samples being collected by ICESAT-2 and GEDI, it can be expected that more and better biomass map products will be derived through fusion of these samples with field measurements and optical and/or radar observations soon. 1.4 Research Objectives and Dissertation Structure Given the large uncertainties in current estimation of the land sink and carbon fluxes driven by land use change, this study is designed with an overarching goal to improve the understanding of carbon budget at regional to global scales. The limitations of existing carbon modeling approaches call for new methods and data products that can leverage the rapid advances in satellite remote sensing of forest dynamics and provide more reliable carbon estimates with better spatial details. Such estimates are needed to answer key carbon budget questions at global, national, and local scales. In particular, estimates with local details are needed to determine whether a specific region is a carbon source or sink during a specific time period, and whether and to what degree different forest management practices, including those prescribed by various national and international climate change initiatives, may alter the source and sink estimates from forest land. The specific research objectives of this dissertation include (Figure 1-2): 10 Figure 1-2 The overarching goal and three specific research objectives of this dissertation study. - Given the advantages of the bookkeeping carbon accounting (BCA) model over other modeling approaches for tracking carbon fluxes driven by land use changes, the first goal is to develop a grid-based carbon accounting (GCA) model that can take advantage of many newly available remote sensing products and use the BCA model to produce carbon flux estimates arising from forest disturbance and post-disturbance recovery at sub-ha spatial resolutions. Such fine resolution estimates will provide critical management relevant details required by local or individual land managers. - Because large portions of harvested carbon can be stored in wood products for decades or longer before being completely released to the atmosphere, and only limited survey-based timber product output (TPO) estimates are available for a 11 few selected years, the second goal is to derive an annual TPO record spanning over three decades that are consistent among states and over time. This record will shed light on the spatial and temporal variable of the influx to and emission from wood products over the entire study region. - The third objective is to further improve the GCA model such that it can ingest the annual TPO data record developed through Objective 2 and produce carbon estimates for all forest areas. With this capability, the model is then used to determine whether a region is a carbon source or sink, and whether and to what degree the source and sink estimates are affected by different forest management practices. The study region is located in the southeast United States where large portions of the land base are dedicated to forestry and forest harvest rates are high (Figure 1-3). After the clear cutting and agriculture expansion that occurred in the 18th and 19th centuries, the region had extensive commercial wood harvest during the 20th century (Birdsey et al., 2006). Previous studies show that the regrowth of the forests has likely turned the region into a net carbon sink since the mid-20th century (e.g. Chen et al., 2006). The GCA model and TPO data developed through this study will be used to determine whether the region has been a source or sink over the last few decades. Due to constraints by data availability and computing resources, North Carolina is selected as the study area for the first and third objectives. For Objective 2, the study region is expanded to include 7 states in the Southeast, including Alabama, Georgia, North Carolina, South Carolina, Florida, Tennessee, and Mississippi. 12 Figure 1-3 Study regions of the three research objectives of this dissertation study. Following this introduction chapter, studies for the three specific research objectives will be presented in Chapters 2, 3, and 4. A summary will be provided in Chapter 5 along with a discussion of the significance and broad impact of this study. 13 Chapter 2 Carbon fluxes from contemporary forest disturbances evaluated using a grid-based carbon accounting model 2.1 Background The carbon exchange between the atmosphere and terrestrial ecosystems, especially forests, is an important component in the global carbon cycle. Forests take up carbon dioxide (CO2) through photosynthesis, and release back a large portion to the atmosphere via respiratory processes while store some carbon in plant biomass for decades or even centuries; meanwhile, a disturbance event, such as fire, disease, insect outbreaks, drought, and harvesting, can trigger accelerated release of stored carbon back to the atmosphere (Goward et al. 2008). However, large uncertainties exist in current estimates of carbon fluxes between the biosphere and the atmosphere (Houghton, 2013; Antonarakis, 2014). According to Le Qu?r? et al. (2016) in their recent Global Carbon Budget report, current emission from land-use change is 1.3?0.7 GtC/yr: the uncertain is more than 50%. Reducing such uncertainties is of great interest in the scientific community, especially in quantifying carbon fluxes at more local and regional scales (Turner et al. 2016). There are many ways to estimate forest carbon pools and fluxes. The first one is the inverse models used with variations in atmospheric concentrations of CO2 to infer sources and sinks (Fan et al., 1998; Pacala et al., 2001; Crevoisier et al., 2007). These studies have generally estimated larger net sinks than other approaches, but the uncertainties are high. The second approach is to use process-based ecosystem models or terrestrial biosphere models (TBMs) to simulate changes in carbon storage on land (e.g. Tian et al., 1999; Hurtt et al., 2002). The estimates of annual NEP from TBMs are 14 also highly variable; for example, the estimates of annual NEP for temperate North America to vary from a sink of 1600 to a source of 500 TgC/yr for the period 2000- 2005 (Huntzinger et al., 2012). The major components of NEP (GPP and Rs) varied even more than NEP among the models. The third approach relies on data from forest inventories to calculate carbon budgets (Birdsey & Heath, 1995; Turner et al., 1995a; Smith et al., 2007; Pan et al., 2011a). Estimates from inventory data are less variable than estimates from inverse studies and ecosystem models, and generally the estimated net carbon sink in forests is also smaller. But inventory data are based on sampling strategies designed to estimate stocks (growth and mortality), not to monitor changes (Bradford et al., 2010). The fourth approach is to utilize carbon accounting methods to track carbon fluxes arising from land use conversion, forest harvest, wildfire, and other land change types (Houghton et al., 1999; Zheng et al., 2011; Brack and Richards 2002). Driven by tabular statistics on land use change and related carbon pools, this approach typically produces results with spatial characteristics no better than those of the tabular input data. Such results might be useful at national or regional levels, but lack the spatial details needed to support carbon management decision makings by local agencies or individual landowners. Deriving carbon estimates with needed spatial details was difficult in the past, partly because spatial products on the required model inputs did not exist. With rapid advances in remote sensing technology, however, it has become increasingly more feasible to map many variables important for tracking track carbon fluxes with increasingly better quality. Landsat images have been used to map land cover and land 15 cover change for decades (Powell et al., 2010; Schroeder et al., 2014). Following the opening of the entire Landsat archive for no-cost access, the multi-decade Landsat record has been used to reconstruct forest disturbance history (Huang et al., 2010), map disturbance agent (Schroeder et al., 2014), and quantify disturbance intensity (Tao et al., 2019). Integration of remote sensing observations with field plot data and/or other measurements allowed derivation of biomass products at local, national, to continental scales (Hall et al., 2006; Santi et al., 2017). With increasingly more optical, radar, and lidar observations provided by existing and forthcoming satellite missions including S2, S1, IceSAT-2, and GEDI, the ability to produce high quality data products on land change, biomass density, and other biophysical variables needed to calculate terrestrial carbon fluxes will continue to improve (Xiao et al., 2019). Effective use of the rich, remote sensing-based datasets to advance carbon management decision makings requires a framework to integrate these datasets with models to produce carbon estimates with required spatial-temporal details. A major goal of this study was to develop a grid-based framework where disturbance, forest carbon, and other remote sensing products can be used as input to a well-established carbon accounting model (Houghton et al. 1999) to produce spatially detailed map products of carbon pools and fluxes. We have tested this framework over North Carolina where detailed forest carbon and disturbance products are available and used it to produce 30-m map products of carbon fluxes arising from forest disturbances occurred between 1986 and 2010. These fine resolution map products are valuable for understanding the spatial-temporal patterns of carbon sources from forest disturbances and sinks from post-disturbance growth. More importantly, they can provide much 16 needed details that are mostly unavailable thus far for carbon accounting, management, and related decision support at individual property owner, municipal, county, or even state levels. Note that any latency effect from before 1986 was not calculated in this study, nor were the absorption from undisturbed pixels included in the model. The total fraction of disturbed forest pixels is about 44.7% in North Carolina. 2.2 Methods 2.2.1 Study Area North Carolina is located in the southeastern United States. Its 100 counties are distributed from the Atlantic coast in the east to the Great Smoky Mountains in the west. About 60% (75199.3 km2) of the state?s 139,390 km2 land base is forest land (Brown et al., 2014), most of which is classified as timberland (Bardon et al., 2010). Major forest type groups include Oak-Hickory, Loblolly-shortleaf pine, Oak-pine, and Oak-Gum-Cypress. More than half of the state?s forests were disturbed at least one between 1985 and 2010 (Huang et al., 2015). While timber harvest is the dominant disturbance type, damages from hurricane, insect outbreak, snow/ice, fire, and other natural disturbances are also common. The state?s total area subject to stand clearing disturbances was relatively stable, but the impact of partial disturbances had large inter- annual variability (Tao et al., 2019). 2.2.2 The original Bookkeeping Carbon Accounting (BCA) model The Bookkeeping Carbon Accounting (BCA) model (Houghton et al. 1999) is a well-established model for tracking carbon fluxes arising from land use and land cover change. This model divides the globe into ecozones following the Food and Agriculture 17 Organization of the United Nations (FAO) Global Ecological Zones (GEZ, second edition). For each ecozone, the BCA model keeps track of the carbon in four major pools: living aboveground and belowground biomass; dead biomass, including coarse woody debris; harvested wood products; and soil organic carbon (Houghton and Nassikas, 2017). These major pools are divided into smaller categories for calculation purposes: soil release and uptake, slash (woody debris left on site), carbon burned on site, regrowth after disturbance, and decays in various industrial wood products. Carbon pools and the fluxes for different disturbance types are tracked at the ecozone level. Table 2-1 lists the forest change processes and carbon pools considered in this study, and the key parameters used to calculate these pools and fluxes are listed in Table 2-2. Table 2-1 Pools tracked by the BCA model for LULCC types considered in this study (indicated by ?X?). Fire Wood Harvesting Conversion to Urban Land Soil Release Soil Uptake Slash X X Burned X X Regrowth X X Wood X Products Table 2-2 Ecozone-specific parameters used by the BCA model Name Unit Description 18 CPrim g/Ha Carbon density in undisturbed, mature/primary forest CMin g/Ha Minimum carbon density after disturbance CSec g/Ha Carbon density in recovered, secondary forest FSlash Fraction of carbon ended up in slash pool DRSlash Decay rate coefficient for slash pool FP1 Fraction of carbon ended up in 1-year decay pool FP10 Fraction of carbon ended up in 10-year decay pool FP100 Fraction of carbon ended up in 100-year decay pool Tms Year Time for forest to grow into secondary forest from stand-clearing disturbance Tsp Year Time for forest to grow into mature forest from secondary forest For a disturbance event that resulted in a disturbance area of LC, the release from the carbon ended up in the slash in any given year X is Slashr in X = (Slash Pool + LC * C * F ) * (1 ? e -DRslash X-1 Prim Slash ) (2-1) And the carbon pool after each year?s release is Slash Pool in X = (Slash PoolX-1 + LC * CPrim * FSlash) * e -DRslash (2-2) The carbon burned on site is released immediately in the year of disturbance, and can be calculated as P1r = LC * FP1 * (CPrim ? CMin) (2-3) For regrowth after disturbance, it is assumed that the growth rate is fast when the disturbed land is recovering from minimum to secondary forest, then the rate becomes slow until the forest reaches mature stage, at which point the growth stops. During each regrowth stage, it is assumed that the growth rate is constant. The annual carbon uptake or sink by regrowth from minimum to secondary forest and from secondary forest to mature forest are calculated using equations (2-4) and (2-5) Regrowthmin2sec = LC * (CSec ? CMin) / Tms (2-4) Regrowthsec2mat = LC * (CMat ? CSec) / Tms (2-5) 19 The total annual uptake by regrowth is the sum of annual uptake by regrowth from all previous disturbances. For the wood products, based on different product types, the pool is divided into a fast decay pool (10-year decay pool, or P10, e.g., paper products) and a slow decay pool (100-year decay pool, or P100, e.g., furniture, wood used in buildings). The calculations for these pools and their annual releases are similar to Eq. (2-3) and Eq. (2-4), except that the exponent is replaced by a constant decay rate. 2.2.3 Development of A Grid-Based Carbon Accounting (GCA) model A major limitation of the current BCA model is its inability to provide results at a resolution required for project level or landscape level carbon management. To overcome this problem, we reimplemented the BCA model within a gridded framework to leverage the increasingly more available remote sensing products that could be used to derive the parameters required by the model. In this framework, a study area is divided into even-sized grids (e.g., 30-m pixels). Carbon pools and fluxes are calculated for individual grid cells instead of ecozones. For each cell, the parameters that can be derived from available remote sensing products will have cell-specific values derived from those products. The remaining parameters will inherit the ecozone-based values from the original BCA model according to which ecozone that cell belongs to. As will be discussed in section 2.3, several remote sensing products were available over the study area, including 1) a carbon content map, and 2) a suite of disturbance products providing details on the timing/year, intensity, as well as type/attribution of each disturbance event. With these products, the GCA model was implemented such 20 that it could account for carbon fluxes arising from multiple disturbance events that occurred within the same grid cell in different years. Further, we used these products to improve the calculation of carbon fluxes arising from disturbance and post- disturbance regrowth through the following steps: First, the carbon density map was used to replace CPrim in determining the initial carbon density of a grid cell right before the first disturbance event detected at that location (pre-disturbance carbon density). Second, the amount of carbon removed from the living biomass within that grid cell due to the first disturbance event was calculated as the product of pre-disturbance carbon density and the percentage of carbon removed (PCR) by that event, which was assumed to be 100% (stand clearing) without using the disturbance intensity data. As will be discussed in section 2.4.1, the disturbance intensity products provided estimates of percent basal area removal (PBAR), which was the percentage of the total basal area of live trees removed by a disturbance event. Based on the allometric equations of Jenkins et al. (2003), which were developed to convert tree diameter measurements to biomass, the percentage of carbon removed (PCR) by a disturbance event can be calculated from PBAR using the following equation: P 6/5 CR = PBAR Third, whether to use equation (2-4) or (2-5) to calculate the initial growth immediately following that disturbance event was determined according to the remaining carbon in live biomass after that event, which was the difference between pre-disturbance density and the amount of carbon removed by that event as calculated 21 above. If the remaining carbon is below CSec, equation (2-4) is used. Otherwise, use equation (2-5). Finally, if more disturbances were detected after the first disturbance over a pixel location, the pre-disturbance carbon density for each subsequent disturbance event was calculated as the sum of the remaining carbon after the previous disturbance and the carbon accumulated through the growth calculate in the third step. The fluxes arising from disturbance and post-disturbance growth were then calculated following the second and third steps. 2.2.4 Model inputs Since the GCA model was in essence a reimplementation of the BCA model within a gridded framework, most of its parameters would have values equivalent to those used in the original BCA model. That is, for each of those parameters, all grid cells in the GCA model that were located within the same ecozone of the BCA model had the same ecozone-based values used by the original BCA model. For the parameters that could be derived from available remote sensing products, their values over each grid cell will be derived from remote sensing products and can vary within the ecozones. The remote sensing products used in this study included a suite of forest disturbance maps and a pre-disturbance carbon density map. 2.2.4.1 Forest disturbance data Several products were used to produce a disturbance dataset with information on the timing (year), type (attribution), and intensity of disturbances that occurred in North Carolina from 1985 to 2010. The first were annual forest disturbance maps showing 22 which pixels had disturbances in what years (Huang et al., 2015). These maps were derived using Landsat time series stacks (LTSS) (Huang et al., 2009a) and the vegetation change tracker (VCT) algorithm (Huang et al., 2010). For each disturbance event detected by the LTSS-VCT approach, if it was also mapped as a burned pixel by the annual fire maps derived through the MTBS project (Eidenshink et al., 2007), then it was a fire disturbance. Otherwise it was a wood harvesting event unless it was classified as a conversion to urban as follows ? if at least one disturbance was detected at a pixel location and that pixel was mapped as an urban pixel in 2011 by the National Land Cover Dataset (NLCD 2011) (Homer et al., 2015), then the last disturbance detected over that pixel was a conversion to urban. It should be noted that if two or more disturbance events were detected over the same pixel location, each event could be of a different type. However, if one of the types was a conversion to urban, that type must be the last one, because wood harvesting or fire could not happen over an area already converted to urban. The decision process used to attribute the detected disturbances is shown in Figure 2-1. 23 Figure 2-1 The process of determining disturbance attribution. For each disturbance event determined above, its intensity was derived based on the PBAR dataset produced by (Tao et al., 2019). As mentioned earlier, PBAR was a measure of the percentage of total live tree basal area that was removed by a disturbance event. The PBAR value for a disturbance event was derived based on the spectral change values associated with that disturbance event and a model calibrated using reference PBAR data calculated from pre- and post-disturbance field measurements collected through the USFS Forest Inventory and Analysis (FIA) program. While in the Tao et al. (2019) study, PBAR was calculated without considering the disturbance type, for this study, the PBAR values were set to 100% for all conversion-to-urban disturbances. 24 2.2.4.2 Pre-disturbance carbon density data As mentioned in section 2.2, pre-disturbance carbon density ? the carbon density at a pixel location right before a disturbance event, is required for calculating the carbon change arising from that event. Since disturbances could and did occur throughout the entire observing period, this would require annual carbon density maps for the entire study period, which unfortunately do not exist. Although there were a couple of CONUS-wide carbon density maps (Kellndorfer et al., 2013; Wilson et al. 2013), the values from those maps could not be used as the pre-disturbance carbon density for disturbances occurred before the derivation of those maps. Given the fact that harvesting is more likely to occur over older forests than younger ones, we applied a spatial filtering method to an existing map to derive a pre-disturbance carbon density map that might provide more realistic values than the values in the original maps or the ecozone level values used by the original BCA model. Specifically, for any pixel location, we calculated the 95 or 99 percentile value with a relatively large area surrounding that pixel, say a 6 km by 6 km window, and used that value as the pre- disturbance carbon density value for that location. To evaluate how realistic the pre-disturbance carbon density values derived this way were, we selected the FIA plots that had field measurements right before they were disturbed and plotted the FIA-based pre-disturbance carbon values against those from the original carbon density map as well as the pre-disturbance carbon density maps derived above. Figure 2-2 shows that the median values from the pre-disturbance maps derived using both 95 and 99 percentiles were much closer to FIA measurements than 25 that derived from the original carbon map, demonstrating that the spatial filtering method did provide more realistic pre-disturbance carbon density values. Figure 2-2 Comparison of carbon density values between plot measurements and derived maps. After filtering the original map with 95th or 99th percentile in a 6 km by 6 km window, the range of values are closer to those of undisturbed plots. 2.2.4.3 Other inputs We used the pre-set parameter values provided with the original bookkeeping model for all other inputs. Attempts were made to find more realistic forest growth rates based on FIA plot data, however they varied greatly among different plots, even in the same ecozone (Figure 2-3). The median values of the plot-based growth rates were close to the fast growth rates provided in the original BCA for forest prior to reaching the secondary forest stage (Table 2-3). Therefore, no change was made to the growth rates of the original model setting. 26 Figure 2-3 Distribution of forest growth rates calculated from FIA plot data. Table 2-3 Growth rate statistics from selected FIA plots. The last two columns show the values provided with the original BCA model. Ecozone No. Median Average Min Max Model Model code FIA (fast) (slow) plots 21 3326 2.09 2.99 0.001 35.31 2.97 0.99 35 919 1.77 2.51 0.002 64.14 1.58 0.53 27 2.2.5 Carbon modeling scenarios The ability of the GCA model to use spatially explicit remote sensing products made it possible to assess the impact of those inputs on carbon estimates. In this study, four scenarios were designed to evaluate the impact of the products described in section 2.2.4 (Table 2-4). In the first scenario, neither the disturbance intensity dataset nor the pre-disturbance carbon density dataset was used. Instead, pre-defined mature forest carbon density values were used as pre-disturbance carbon density and all disturbances were considered stand clearing according to the original model setting (OMS) of the BCA model. In scenario 2, only the spatial disturbance intensity (SDI) dataset was used. Pre-defined mature forest carbon density values were used as pre-disturbance carbon density values. In scenario 3, only the spatial pre-disturbance carbon density (SPC) dataset was used. All disturbances were assumed to be stand clearing. In the last scenario, both spatial datasets (SDI+SPC) were used. Assuming the remote sensing products would produce more realistic carbon estimates, the final carbon fluxes arising from the mapped disturbances were calculated using scenario 4 (SDI+SPC). But the differences among the four scenarios could be used to examine the impact of the remote sensing products on the estimation of carbon fluxes arising from forest disturbances and recover in recent decades. 28 Table 2-4 Scenarios for assessing model sensitivity to spatial carbon density and disturbance intensity data. Scenario Pre-disturbance C value used Disturbance intensity data used? Original model Ecozone-specific constant No setting (OMS) Spatial Ecozone-specific constant Yes disturbance intensity only (SDI) Spatial Pre- Spatial pre-disturbance C density No disturbance C map density only (SPC) SPC + SDI Spatial pre-disturbance C density Yes map 2.3 Results 2.3.1 Model verification and sensitivity analysis The model was first coded following the logic and assumptions of the original BCA model, with every disturbance assumed to be stand-clearing, the initial carbon density right before the disturbance always at the mature forest level, and conversion to agricultural or urban land can happen multiple times on a pixel. With this intermediate model, not only the total net fluxes over the entire observing period but also the annual source and sink estimates for all major carbon pools derived using the two models were the same, demonstrating that we have successfully deployed the BCA model in a gridded framework for producing spatially detailed carbon estimates. Then, disturbance attribution process was improved to account for multiple disturbance events on the same pixel (see Figure 2-1), disturbance intensity and starting carbon 29 density were replaced by gridded data derived from satellite remote sensing products, and the carbon density throughout the disturbance-regrowth process was tracked in the finalized GCA model. Model results derived with and without using the spatially explicit pre-disturbance carbon density map and the disturbance intensity dataset as inputs to the GCA model show that these inputs had large impact on both the total carbon estimates over the 25- year study period (Figure 2-4) and annual flux estimates (Figure 2-5). Compared to benchmark results derived using original model settings (the OMS scenario), the net carbon flux arising from mapped disturbances over North Carolina during the study period was reduced from 242.9 MT (net source) to 136.4 MT (net source) when only the spatial disturbance intensity (the SDI scenario) was used. The flux value became negative (net sink) when the pre-disturbance carbon density map was used but all disturbances were assumed stand clearing (the SPC scenario). The net flux was further reduced to -47.2 MT when both spatial maps were used (the SPC+SDI scenario). 30 Figure 2-4 Comparison of carbon fluxes by disturbance type between different scenarios. Using both RS-based carbon density map and disturbance intensity data has turned the 25-year net total C flux from source to sink. Because wood harvesting was the dominant disturbance type in the study region, the two spatial datasets had the largest impacts on harvest source estimation. Emission from wood harvesting was estimated at 366.9 MT for the study region using the OMS scenario. It was reduced to 247.9 MT and 133.4 MT under scenarios SDI and SPC respectively, and further down to 88.5 MT when both datasets were used (the SPC_SDI scenario). The impacts of the two datasets on fire source estimates were similar but at much smaller scales because the total area affected by fire was only a small fraction of that subject to harvesting. Driven by post-disturbance growth rates, sink calculation for both post-harvest and post-fire growth should not be directly affected by the two spatial datasets. However, because the model assumes a fast growth rate at relatively low carbon density 31 and changes the growth rate to a lower value when a pixel?s carbon density exceeds a threshold value, use of the two spatial datasets can affect a pixel?s carbon value at any given time and hence can have indirect impact on the sink estimates for these post- disturbance growth processes. While the impact of the pre-disturbance C dataset on source calculation for forest- to-urban conversion is similar to that for wood harvesting and fire events, the impact of the disturbance intensity dataset is more complicated. In theory, source calculation for forest-to-urban should not be affected by the disturbance intensity data, because the intensity should always be 100% for this change process. However, a pixel converted to urban could have one or more disturbances (mostly harvesting) prior to the final conversion to urban. For such a pixel, if a mature forest carbon value was assigned to it prior to its earliest disturbance event, a low disturbance intensity would result in a small carbon source and trigger the model to start calculating carbon increase from post-disturbance regrowth and the pixel?s carbon density could exceed the assumed value for mature forest before the pixel was converted to urban. As a result, the source estimation for that forest-to-urban conversion event under the SDI scenario could be higher than that calculated under the OMS scenario. This could only happen to pixels that had other disturbances prior to the forest-to-urban conversion. This is why the annual forest-to-urban source values calculated under the OMS and SDI scenarios were similar in the first couple of years but differed in later years (Figure 2-5(C)). 32 Figure 2-5 Annual total carbon absorbed and released by disturbance types during the study period under the four different scenario settings. The trendlines show large differences among the sources under the four scenarios. 2.3.2 Carbon fluxes from disturbance and post-disturbance growth Carbon fluxes arising from disturbances and post-disturbance growth over North Carolina were modeled using the SPC+SDI scenario where both the disturbance intensity and pre-disturbance carbon density datasets were used as inputs to the GCA model. Sources arising from harvesting, fire, and forest-to-urban conversion, as well as sinks from post-harvest and post-fire regrowth were calculated on an annual basis for every 30-m pixel that had at least one disturbance mapped. Figure 2-6 shows the net flux images derived by summing up all source and sink terms over the study period for 33 each pixel, along with a few full resolution zoom-in examples showing the fluxes driven by different change processes. Figure 2-6 A browse image of the 30-m net flux map derived by summing up all source and sink terms over the study period for each pixel. The full resolution zoon-in examples show the flux patterns over areas with (a) extensive forest-to-urban conversion, (b) a large fire occurred in 2008, (c) active wood harvesting followed by strong growth throughout the study period, (d) minimal disturbances within national park/national forests, and (e) storm damage/salvage logging followed by recovery. 34 These 30-m maps can be aggregated for any geographic or administrative regions (e.g., areas affected by individual disturbance events, properties of individual landowners, districts, counties, and state) to derive estimates needed for addressing specific carbon management and/or decision support needs. At the state level, forest harvesting and fire from 1986 to 2010 released 88.5 MT and 1.6 MT carbon respectively. During the same period, regrowing trees over the logged area absorbed 142.7 MT carbon while those over burned area added 1.6 MT more. The net flux from harvesting, fire, and post-disturbance growth was -52.5 MT. Conversion of forest to urban resulted in a net source of 5.3 MT. Overall, the areas subject to the three types of disturbances and post-disturbance growth was a net sink of 47.2 MT carbon over the entire study period. The source and sink estimates differed substantially among the counties within the state (Figure 2-7). Most counties with high carbon releases due to harvesting are located in the eastern side of the state, with Bertie and Beaufort having the highest releases. Except for Rutherford County and Wilkes County, harvesting emissions in most counties on the west side were only small fractions of many counties on the east side. Emissions due to the other two disturbance types were small in general. However, each of the two disturbance types had some hotspots. As expected, counties surrounding population centers had high emissions from urbanization, including the Mecklenburg County (the Charlotte metropolitan area), Wake County and Durham County (the Raleigh-Durham-Chapel Hill area, Figure 2-6(a)), and Guilford County (suburban Greensboro). The few other counties that had sizable emissions from urbanization also had mid-sized cities. 35 Large fire is not common in NC. In 2008, however, a fire ignited by lightening burned more than 41,500 acres, mostly within the Pocosin Lakes National Wildlife Refuge, resulting in high emissions in Tyrrell and Hyde counties. Several other counties, including Pender, Craven, and to lesser degrees, Burke and Washington, also had sizable emissions due to fire. The carbon sink related to wood harvesting and fire is driven by post-disturbance growth. Given the fast forest growth rates in the study region (especially in ecozone 21), it does not take many decades for the carbon accumulation through post- disturbance growth to surpass the amount released by a disturbance event, especially when the pre-disturbance carbon is low or the intensity of that disturbance event is low, or both. As a result, most counties with high carbon releases due to wood harvesting had even higher sink values over the entire study period. Except for Pender County where the carbon sink from post-fire growth by 2010 exceeded the release from a fire occurred in 1986, the few counties that had fire disturbances had minimal sinks from post-fire growth. Overall, the total C sink arising from post-disturbance recovery exceeded the total release from all three disturbance types in most counties in east and central North Carolina, resulting in negative net fluxes in those counties during the study period. The relatively large emissions from forest-to-urban conversion in Mecklenburg County and the 2008 Pocosin Lakes fire that spilled into Tyrrell County were offset by sinks from post-disturbance growth. As a result, the net fluxes from disturbance and post- disturbance growth were near 0 in both counties. Three other counties in the east along with more than a dozen counties in the west also had near 0 net fluxes. Most of these 36 counties had very little emissions from disturbances and hence not much post- disturbance growth. The slightly positive net fluxes in the Swain, Jackson, and Macon Counties were likely due to slightly higher disturbances towards the end of study period. Figure 2-7 Total carbon absorption and emission by disturbance types at county level during the 25- year study period. Wood harvesting activities are prevalent in the eastern part of the state, resulting in both large sources and sinks. Only a handful of counties are affected by fire. Emission from conversion to urban land is prominent in counties with major urban centers. Overall, the emission is less than or almost equal to the absorption in all counties. 37 2.4 Discussion The BCA model provides a useful tool for estimating carbon pools and fluxes arising from land use change. It has been used to calculate carbon estimates at national and global scales. However, because the smallest modeling unit of the model is an ecozone, it is difficult to use it to derive estimates with sub-ecozone details. To address this limitation, we reimplemented the BCA model within a gridded framework. The resultant grid-based carbon accounting (GCA) model provides a flexible framework for integrating remote sensing products into the carbon accounting process to produce spatially disaggregated carbon estimates, which are increasingly needed to support fine scale carbon management decision makings and related applications. As with the original BCA model, the estimates produced by the GCA model are attributed to different pools and change processes at an annual time step. While the BCA model has been widely used to estimate carbon fluxes arising from historical land use changes (Houghton et al. 1999; Houghton 2003b; Houghton et al. 2012; Houghton & Nassikas 2017), the GCA model is intended for modeling the carbon dynamics of forest changes that occurred in recent decades during which key datasets such as disturbance history and carbon density could be mapped using RS data. The sensitivity analysis conducted in this study revealed large differences between the carbon estimates derived based on RS products and those derived using parameter values of the original BCA model, which were fine tuned for modeling historical land use changes. Compared to estimates derived by assuming the forests had a mature forest carbon density value and all disturbances were stand clearing events, the carbon release from wood harvesting was reduced by more than 30% and 60%, respectively, when RS 38 based disturbance intensity and carbon density data were used in deriving those estimates, and by more than 75% when both were used. As a result, the net flux arising from the observed disturbances and post-disturbance growth in North Carolina changed from 242.9 MT (source) to -47.2 MT (sink) when the disturbance intensity and pre- disturbance carbon density values used by the original BCA model were replaced by values derived based RS products. Like any other science data, RS products, including those used in this study, typically have varying levels of uncertainties. However, as RS technology advances rapidly and increasingly more and better calibration data are becoming available, more RS products with better quality will be developed. The gridded framework of the GCA model allows rapid integration of these new products to improve carbon estimation. This framework can be adapted for use with remote sensing products of any other spatial resolutions. While methods have been developed for producing forest disturbance products on an annual basis, forest carbon maps over large areas are available only for very few years (Wilson et al. 2013). It should be noted that the carbon value from a map developed for a specific year cannot be used directly as the pre-disturbance carbon for calculating carbon emissions from a disturbance event that occurred in a different year. For areas where wood harvesting is the dominant disturbance type, however, it is reasonable to assume that older forests with higher carbon values are more likely to be logged than younger ones. Based on this assumption, we developed a method to create a pre-disturbance carbon density map based on a map produced for a specific year. The 39 values from the derived map were much closer to pre-disturbance carbon values measured by FIA field crew. It should be noted that GCA model implemented in this study only tracks the carbon accumulated through growth after a disturbance event. It does not calculate carbon changes from growth for the years before the first disturbance was detected over a pixel location, nor does it do so over forest areas that have no detected disturbances. While we intend to add modules to track these carbon changes, the carbon estimates derived through this study do not provide a complete picture of forest carbon dynamics over North Carolina, and hence should not be used as evidence as to whether harvesting reduces or enhances carbon sequestration. More comprehensive experiments are needed to address this important question. 2.5 Conclusions The Bookkeeping Carbon Account (BCA) model was implemented within a gridded framework to allow it to produce carbon estimates with spatial details beyond the ecozone-based modeling unit of the original BCA model. The resultant grid-based carbon accounting (GCA) model has been calibrated and parameterized such that it produces the same state level flux estimates as the original BCA model over North Carolina when the inputs to both models are equivalent. As with the original BCA model, the GCA model calculates carbon fluxes between major carbon pools arising from several disturbance types and post-disturbance growth over forestlands, including timber harvesting, fire, and conversion to urban area. An important feature of the GCA model is that it can integrate remote sensing products into the carbon accounting process to produce spatially disaggregated carbon 40 estimates, which are increasingly needed to support fine scale carbon management decision makings and related applications. In North Carolina, the model was used to estimate carbon fluxes arising from forest disturbances mapped using historical Landsat observations over 1986-2010. The net flux from those disturbances and post- disturbance growth over the study period was 242.9 MT (source) when derived based on the original parameter values of the BCA model, which were tuned for modeling fluxes from historical land use change over several centuries. The flux was reduced to -47.2 MT (sink) when remote sensing based disturbance attribution, disturbance intensity, as well as pre-disturbance carbon density products were used as model inputs, demonstrating that use of remote sensing products could have large impact on modeling the carbon fluxes arising from forest disturbance and post-disturbance growth. Therefore, improved estimates of those fluxes will rely heavily on improving RS products using better observations and more accurate and more representative calibration data. Currently, the GCA model only tracks carbon accumulated through growth after a disturbance event. It does not calculate carbon absorption by forests not disturbed during the study period. The net carbon sink estimation derived through this study does not suggest that harvesting enhances or reduces carbon sequestration. To address this question, it?s necessary to estimate the carbon fluxes over the disturbed areas by assuming those areas were not disturbed during the study period, which will be done in Chapter 4. 41 Chapter 3 Estimation of Annual Harvested Wood Products based on Remote Sensing and TPO Survey Data 3.1 Introduction Timber harvest is a dominant anthropogenic method for removing forest carbon from the standing biomass pool. In the US, it is the number one cause of anthropogenic disturbances to forests, especially in major timber production regions (Zheng et al. 2011; Harris et al. 2016). As a measure of the primary output of this activity, Timber Product Output (TPO) can be used to estimate carbon pools and fluxes in many carbon budget models. For example, Harris et al. (2016) used TPO survey data from the US to calculate the net carbon change from timber harvest, which showed that timber harvesting activities accounted for 86% ~ 92% of total C loss. On the other hand, the carbon in many timber products is not completely released immediately after timber harvest. While paper and some other wood products can store carbon for only a few years, furniture and building materials typically last for decades to centuries before they are completely decomposed (Ruddell et al. 2007; Zeng 2008; Lippke et al. 2011). Detailed information on different wood product type is therefore important for accurate estimation of the carbon pools and emissions from wood products. In the United States, timber product output (TPO) reports are produced by the Forest Service Forest Inventory and Analysis (FIA) program based on survey data collected from wood processing mills (Woodbury et al., 2007). While these reports have been produced for all 48 states in the conterminous U.S., they are not available for every year for any given state. Further, the number of available reports, as well as the reporting frequencies and reporting years are not consistent among states (Huang 42 et al. 2015), making it difficult to derive consistent estimates of carbon pools in wood products and related fluxes among states (Birdsey, 2004; Zhu et al., 2010). With the ability to image the Earth?s land surface repeatedly, satellite remote sensing can provide more consistent observations of forest disturbances than allowed by using ground-based methods, which are often labor intensive and time consuming, and hence difficult to deploy with adequate frequency at the national scale. While direct estimation of timber volume and other forest attributes using optical data may be challenging (e.g. Makela and Pekkarinen, 2001; Trotter et al., 1997), the multi-decadal Landsat record has enabled annual mapping of forest disturbances at regional (Huang et al., 2009b; Kennedy et al., 2012), national (White et al., 2017; Zhao et al., 2018), and global scales (Hansen et al., 2013). In addition to detecting disturbance events, time series Landsat data have also been used to attribute disturbance agents (Kennedy et al., 2015; Schleeweis et al., 2020) and quantify disturbance intensity (Tao et al., 2019). Because Landsat-based disturbance products are derived from a globally consistent observational record that has been calibrated over multiple decades (Chander et al., 2009; Mishra et al., 2016), they may provide a basis for deriving annual TPO estimates that are spatially more consistent and temporally more complete than available survey data. For example, Huang et al. (2015) demonstrated that in North Carolina, the total TPO can be estimated annually from 1986 to 2010 using the Ordinary Least Square (OLS) regression method and Landsat-based annual forest disturbance products, whereas the survey data was available only for 10 of the years from 1992 to 2009. Building on that study, Ling et al. (2016) developed a suite of linear models for estimating the carbon influxes to different wood product types for 43 each year between 1986 and 2010. By accounting for the specific effects in both spatial (county) and temporal (year) domains, those models provided annual estimates of 6 harvested wood products for each county in North Carolina, including saw log, pulpwood, and fuelwood from softwood and hardwood, respectively. The purpose of this chapter is to further explore the feasibility to derive multi- decadal TPO records for a much larger timber production region, namely, the southeast US, which includes North Carolina, South Carolina, Alabama, Florida, Georgia, Mississippi, and Tennessee (Figure 3-1). Given the demonstrated robustness of many advanced machine learning algorithms in a wide range of remote sensing applications (Huang and Jensen, 1997; Shao and Lunetta, 2012; Verrelst et al. 2012), this chapter seeks to evaluate whether one of the most commonly used algorithms ? Random Forest (RF), can provide more robust TPO estimates across this large region. RF is an integrated learning method that can solve both classification and regression problems (Ho, 1998; Breiman, 2001). It has been used in numerous studies to generate classification maps or derive quantitative estimates of vegetation biophysical variables (Du et al., 2015; Pal, 2005; Rodriguez-Galiano et al., 2012; Pearse et al., 2017; Ploton et al., 2017; Liu et al 2019). 3.2 Study Area and Data Collection 3.2.1 Study Site Located in the southeastern of the United States, the seven states of the study area cover a total area of 354,053 square miles that include both lowlands and highlands. While the region is dominated by rolling hills, plateaus and rich river valleys in the 44 north, the south is characterized by beaches, swamps, and wetlands. The region has a subtropical monsoon climate in the south transitioning to a temperate climate to the north. The summer is long, hot, and humid, while the winter is short and warm. The average annual precipitation in this area ranges from about 43 to 54 inches. The forest coverage of the seven states is between 50.68% and 70.57% respectively, of which Alabama, South Carolina, Georgia, and Mississippi have the largest forest coverage, all above 65%. Southern forests have the highest planted timberland rates in the US (71% of all planted timberland), of which Alabama has 33% planted timberland rate, 32% for Georgia, 32% for Mississippi, 31% for Florida, and 31% for Louisiana, all have the highest proportions of planted to total timberland, nationally (Oswalt et al. 2017). In the southern planted timberland, the main forest types are Loblolly pine (Pinus taeda) and shortleaf pine (Pinus echinate), which account for 71% of all planted forests in the South. 45 Figure 3-1 Study area located in the US southeastern States, including North Carolina, South Carolina, Alabama, Florida, Georgia, Mississippi, Tennessee. 3.2.2 Data collection 3.2.2.1 TPO Survey Data The TPO survey data (USDA Forest Service, 2019) were acquired from Forest Inventory and Analysis (FIA) national program of the United States Forest Service (USFS). FIA carried out state level TPO studies to investigate round wood uses for industrial and non-industrial purposes in a given state using ground-based survey methods. In order to determine the geographic origin, harvest date, volume, species, and use of harvested roundwood products, FIA canvasses all primary wood-using mills, harvest sites, residential users, and commercial producers that harvest and sell wood products for selected study years (Johnson, 1996, Woodbury, et al. 2006). For each study year the collected data are compiled to produce a TPO report that provides county 46 level wood product estimates for that year. The survey data includes county level estimation of industrial roundwood, fuelwood, chips, post, poles, and so on. The industrial roundwood estimation mainly consists of the estimates of industrial hardwood (sawlog, pulpwood) and industrial softwood (sawlog, pulpwood). In this study, industrial roundwood was considered only, as it is the main component of total roundwood production in the study region. Table 3-1 lists the TPO reports available from 1995 to 2009 over the study region. All of the 7 states had data for 1995, 1999, 2005, 2007, and 2009. Except for Mississippi, the other 6 states also had data for 1997 and 2003. In addition, TPO reports were available for Georgia, North Carolina, South Carolina, and Tennessee in 2001. While Mississippi had TPO data for only 5 years, it?s the only state that had data for 2002. These reports provided county level estimates for six major wood product types, including softwood saw log (SSL), softwood pulpwood (SPW), softwood fuelwood (SFW), hardwood saw log (HSL), hardwood pulpwood (HPW), and hardwood fuelwood (HFW). 47 Table 3-1 Years in which survey data are available for each state, "x" indicates the data is available for that year in that state. State 1995 1997 1999 2001 2002 2003 2005 2007 2009 Alabama x x x x x x x Florida x x x x x x x Georgia x x x x x x x x Mississippi x x x x x x No. Carolina x x x x x x x x So. Carolina x x x x x x x x Tennessee x x x x x x x x As the only data source that can provide survey-based TPO estimates for individual counties during a specific year, the TPO reports are highly valuable for calibrating and validating models designed to produce county level TPO estimates on an annual basis (Huang et al., 2015). The average annual timber production volumes for hardwood, softwood, and both calculated from the available TPO reports are listed in Table 3-2. Except for Tennessee, all other states had substantially higher softwood production than hardwood production. With a much lower total wood production than that of each of the other six states, Tennessee?s hardwood production was over 3 times of its softwood production. 48 Table 3-3-2 Average annual timber production volume (in thousand cubic meters) in each state based on the survey data listed in Table 3-1. State Hardwood Softwood Total Alabama 9733.46 23986.03 33719.49 Florida 1574.83 13006.01 14580.83 Georgia 7831.17 27986.87 35818.03 Mississippi 7675.47 19330.68 27006.15 North Carolina 9217.77 14760.85 23978.62 South Carolina 4299.35 13804.73 18104.08 Tennessee 7409.33 2328.98 9738.31 3.2.2.2 Forest disturbance and land cover data Forest disturbance data used in this study included 30m disturbance maps derived using Landsat time series stacks (LTSS) and the vegetation change track (VCT) algorithm (Huang et al., 2011; Zhao et al., 2018). An LTSS is a stack of Landsat images assembled for a World Reference System (WRS) path/row tile to provide clear view observations at a regular time step (Huang et al., 2009a). An annual LTSS consists of an individual Landsat acquisition or a composited image for each year to provide cloud free or near cloud free (< 5% cloud cover) time series observations for the summer leaf- on season of a multi-year study period. VCT is an automated algorithm designed to analyze an LTSS to produce disturbance products (Huang et al., 2010). This is achieved through two major steps. 49 First, each image in the LTSS is analyzed to mask water, cloud, and cloud shadow pixels and to select forest samples from densely forested areas. Based on these forest samples, an integrated forest z-score (IFZ) index is calculated for each pixel using the following formula: ? ? b ? b ? IFZ = ? ? i i ? ? ?? SD ? i ? band?,?,? (3-1) where bi denotes the pixel?s value in band i; and bi and SDi are the mean reflectance value and its standard deviation of the forest samples selected above. IFZ index is a non-negative and inverse indicator of the likelihood of a pixel being a forest pixel (Figure 3-2). The higher IFZ value, the more likely the pixel being a non-forest pixel. On the contrary, the closer IFZ is to zero, the more likely the pixel is to be a forest pixel. In step 2, the IFZ time series of each pixel is analyzed to determine whether that pixel had forest cover during at least part of the study period. If yes, the pixel is further evaluated to determine whether disturbances occurred during the study period. If yes, the disturbance year is determined for each detected disturbance event and its disturbance magnitude is calculated as the change in IFZ as shown in Figure 3-2. 50 Figure 3-2 A forest disturbance event typically results in an abrupt increase in IFZ, which decreases gradually due to post-disturbance recovery (Huang et al., 2015). Disturbance year is determined by the timing of the abrupt change and the change in the IFZ value provides a spectral measure of the disturbance intensity. The VCT algorithm has been used to produce 30m disturbance products for the conterminous U.S. from 1986 to 2010 (Zhao et al., 2018). The disturbance products used in this study was generated as part of that effort, which were extended to 2015. Because the TPO data provide details on softwood and hardwood products, a land cover map derived as part of the National Land Cover Database (NLCD) was used in this study to provide information on the distribution of deciduous and evergreen forests over the study region. NLCD is an ongoing US-wide land cover mapping effort that has resulted in NLCD products for many epochs since 1990s. The 1992 map should provide a better characterization of pre-disturbance forest distributions than later NLCD products, and hence was used in this study. 51 3.3 Methods 3.3.1 Preparation of predictor variables A key hypothesis in modeling TPO from disturbance data is that timber harvest volume should be correlated with Landsat based disturbance estimates (Huang et al., 2015). Given the disturbance products, harvested timber volume (HTV) could be calculated from stock density, harvest area and harvest intensity as follows: HTV = stock density x (logging area x logging intensity) (3-2) And TPO would be the difference between HTV and slash: TPO = HTV ? slash (3-3) However, logging is a dominant but not the only disturbance agency in the southeast. As a result, the disturbance data can only provide proxy information for logging. Also, no good data is available on the slash amount from logging and the stock density before each logging event. Therefore, direct calculation of TPO using the above equations was not possible. Statistical models calibrated using TPO survey data were used to establish relationships between TPO and the disturbance data. Since the TPO survey data were only available at the county level, the 30m disturbance data were aggregated to produce county level estimates. Based on lessons learned from the Huang et al. (2015) study, the disturbance areas within each county in any given year of the study period were calculated at four IFZ disturbance magnitude (IFZ-DM, Table 3-3). Further, these areas were calculated for disturbances over deciduous, evergreen, and mixed forest pixels as mapped by the 1992 NLCD land cover map. Therefore, a total of 12 disturbance area values were calculated for each county in each disturbance year. 52 Table 3-3 IFZ Disturbance magnitude (IFZ-DM) levels at which county level disturbance areas were calculated. Distubance Magnitude Threshold Values Level 1 IFZ-DM ? 3 2 3 < IFZ-DM ? 6 3 6 < IFZ-DM ? 9 4 IFZ-DM ? 9 The effective date range of the TPO survey data for any given year does not match that of the disturbances mapped by the VCT for that year (Figure 3-3). The survey data is collected from the beginning to the end of a year (hereby called survey year), while most disturbance events mapped for a year by the VCT algorithm could occur any time between the acquisition dates of the images used for that year and the previous year (Figure 3-3). Further, due to residual clouds and other data quality issues, the actual disturbance year could also be one year before or after the mapped disturbance year (Huang et al., 2010; Thomas et al., 2011). To reduce the influence of such temporal mismatches, for each TPO modeling/prediction year it was necessary to consider two more years ? one immediately before and one immediately after that year (Huang et al. 2015, Ling et al. 2016). Therefore, a total of 36 predictor variables (12 per year x 3 years) were used to model and predict TPO for any given year. 53 Figure 3-3 Temporal mismatch between survey year and potential disturbance event (Huang et al. 2015). 3.3.2 Modeling approaches Three modeling algorithms were examined in this study, including the ordinary least square (OLS) linear regression and the fixed effects linear regression (FELR) modeling approach used in previous studies (Huang et al. 2015; Ling et al. 2016), as well as the Random Forest algorithm. The OLS linear regression model can be written as: ??? = ? + ???? + ??? (3-4) where ???, ???, and ??? is are the TPO value for a specific wood product type, a vector of the 36 predictor variables as described in section 3.3.1, and model residual for county i in year t, and ? and ? are the regression intercept and the slope vector for the predictor variables, respectively. Both ? and ? are independent of the county index i and year index t. By considering both county specific and year specific effects, the fixed effects regression model can be written as: ??? = ?? + ???? + ?? + ??? (3-5) where ?? and ?? denote county specific and year specific effects, respectively. Random Forest is an ensemble learning method based on classification and regression trees (Breiman, 2001). The algorithm ensembles a large number of trees 54 where each tree is developed using a random subset of the predictor variables and is trained using a random subset of the training data (Breiman, 1996; Breiman et al, 2017). With demonstrated robustness for addressing both classification and regression problems (Du et al., 2015; Pal, 2005; Rodriguez-Galiano et al., 2012), it has been used in many forest remote sensing applications (Pearse et al., 2017; Ploton et al., 2017; Liu et al., 2019). The three approaches were used to model each of the 6 individual product types, including softwood saw log (SSL), softwood pulpwood (SPW), softwood fuelwood (SFW), hardwood saw log (HSL), hardwood pulpwood (HPW), and hardwood fuelwood (HFW), as well as the total of the 6 types. The robustness of OLS and Random Forest was evaluated through a ten-fold cross-validation method. In this method, the reference samples were divided into 10 random subsets. Each subset was used to evaluate a model calibrated using the remaining 9 subsets. This was repeated 10 times such that each time a different subset was used as test data. Results from all 10 experiments were combined to produce representative accuracy estimates for the final model to be derived using all reference samples. As will be discussed in section 3.4.1, the FELR method produced abnormally large and erroneous TPO values. Therefore, there was no need to further assessed that algorithm using cross validation. 3.4 Results 3.4.1 Limitations of the FELR method Similar to the results reported by the Ling et al. (2016), the FELR method appeared to be able to model the individual wood product types for all states in the Southeast. A 55 comparison of model estimates and the survey data over South Carolina and Tennessee is shown in Figure 3-4. However, the derived models predicted abnormally high or low values that appeared to be obvious errors (Figure 3-5). A detailed examination of the model inputs revealed that in many of years where the model predictions were abnormally high or low, there were some counties that had TPO values but their disturbance areas were 0, or had disturbances mapped but their TPO values were 0. Obviously, these were errors in either the disturbance data or the survey data, or both. It?s like that the FELR method was extremely sensitive to such errors. Abnormally high or low predictions by the FELR method were also found in Florida, Georgia, and Mississippi. Given these large errors in the predictions, the FELR method was not considered any further in this study. 56 Figure 3-4 Comparison of model estimates (y-axis) derived using the fixed effects linear regression (FELR) method over South Carolina (top) and Tennessee (bottom) and TPO survey data. Each point represents a county. The orange lines are 1:1 lines. 57 Figure 3-5 Abnormally high or low values were predicted in many years by the FELR method in South Carolina (left) and Tennessee (right), which were obviously wrong. Such obvious errors were also found in Florida, Georgia, and Mississippi. 3.4.2 Performance of the OLS and RF for TPO modeling Several wood product types were modeled reasonably well using both OLS and RF. Comparison of the predicted values against the survey data resulted in coefficient of determination (R2) values of 0.5 or higher for multiple wood product types in several states (Tables 3-4 and 3-5). Except for Tennessee, both OLS and RF resulted in better R2 values on the other 6 states for the softwood types than for the hardwood types, indicating that relationships between timber output for hardwood types and the disturbance data were not as good as those for softwood types. Among the 7 states, South Carolina?s TPO appeared to be more difficult to model. Except for the hardwood sawlog (HS) in Florida, the R2 values derived for this state were lower than those for the other 6 states for all 6 wood types. 58 Table 3-4 R2 values of the relationships between TPO values predicted by OLS and survey data derived through 10-fold cross-validation. AL FL GA MS TN NC SC SP 0.61 0.41 0.56 0.39 0.53 0.55 0.26 SS 0.61 0.47 0.40 0.53 0.09 0.57 0.27 SF 0.64 0.62 0.51 0.54 0.34 0.60 0.36 HP 0.46 0.06 0.11 0.23 0.52 0.09 0.00 HS 0.22 0.00 0.16 0.16 0.28 0.35 0.05 HF 0.56 0.14 0.20 0.55 0.22 0.17 0.02 Table 3-5 R2 values of the relationships between TPO values predicted by RF and survey data derived through 10-fold cross-validation. AL FL GA MS TN NC SC All 7 States SP 0.63 0.57 0.64 0.60 0.62 0.56 0.37 0.62 SS 0.63 0.66 0.48 0.66 0.18 0.60 0.41 0.59 SF 0.59 0.71 0.56 0.57 0.44 0.62 0.43 0.37 HP 0.55 0.23 0.23 0.28 0.54 0.18 0.04 0.41 HS 0.34 0.05 0.20 0.21 0.46 0.41 0.18 0.35 HF 0.55 0.24 0.29 0.44 0.34 0.30 0.10 0.29 Overall, RF performed better in modeling TPO than the OLS. Its R2 values were higher than those derived using OLS for most of the states and wood types (Figure 3- 6). While the differences were small in some cases, on average use of RF instead of the OLS resulted in an improvement of 0.08 in the R2 value, and the improvements exceeded 0.1 for at least one of the 6 wood types in each of the 7 states. OLS had substantially better performance than RF for only two cases: softwood fuelwood (SF) in Alabama and hardwood fuelwood (HF) in Mississippi, suggesting that in these two cases the RF models might have overfitting issues. 59 0.80 0.70 0.60 SPW 0.50 SSL 0.40 SFW 0.30 HPW 0.20 HSL 0.10 HFW 0.00 1:1 line 0.00 0.20 0.40 0.60 0.80 R2 from Random Forest Figure 3-6 Coefficient of determination (R2) values derived by comparing survey data and values predicted by OLS (y-axis) and Random Forest (x-axis) through 10-fold cross validation. 3.4.3 Spatial-Temporal Variability of TPO in the Southeast Given the overall better performance of the Random Forest as compared to OLS for TPO modeling, it was used to produce annual TPO estimates for all 7 states in the southeast. From 1986 to 2015, the region produced more than 5 billion m3 wood products, including 3.7 billion m3 softwood products and 1.6 billion m3 hardwood products (Table 3-6). Of the 7 states, Georgia and Alabama had the highest TPO values, while Tennessee produced the least amount of wood products, although the state was larger and had more forest land than South Carolina. Tennessee was also the only state that produced more hardwood products than softwood products. All other 6 states had higher total TPO values for the softwood types than for the hardwood types. 60 R2 from OLS Softwoods are mostly used by the construction industry or to produce paper pulp and card products (Ryan, 2019). They account for about 80% of the world's timber production, with traditional production centers being located in the Baltic region, North America, and China (Europe, U.N.E.C,2007). Table 3-6 Estimated 30-year (1986-2015) total TPO (thousand m3, or kcume) for the 7 states in the southeast United States. State SP SS SF HP HS HF All Types AL 330487 404207 7269 175158 87575 43301 1065097 FL 235876 215251 4439 45939 16294 18238 559074 GA 439448 496884 11946 148388 104004 62604 1309772 MS 238282 362725 6186 146179 76195 35335 878470 TN 45456 58719 3152 72831 105075 36378 324335 NC 159241 239199 6613 111021 94753 50262 675047 SC 206721 215908 4521 61986 36193 24301 567153 All 1655511 1992892 44126 761500 520089 270420 5378947 States The derived TPO estimates had large variations among the counties within each state. Counties that had high TPO values were mostly in Alabama, South Carolina, the central and southern part of Mississippi, northern Florida, and the east part of North Carolina (Figure 3-7). Most counties in Tennessee, southern Florida, western Mississippi, north central Georgia, and the northern most part of Alabama had low TPO values. 61 Figure 3-7 County-level total TPO values predicted by the Random Forest for the 30-year study period (1986-2015). Overall, the average state-level annual TPO estimates calculated from the 30-year record were comparable to those derived from the TPO survey data, with the 30-year averages being higher than the survey data for the softwood types in some states (Figure 3-8). In particular, the 30-year average TPO values for softwood sawlog (SS) and softwood pulpwood (SP) in Mississippi were 14.6 and 16.6 million m3/year, whereas the averages from the survey data were 12.7 and 12.9 million m3/year, respectively. Alabama?s 12.5 million m3/year 30-year average value for softwood sawlog was ~2 million m3/year higher than the 11.6 million m3/year average value calculated from survey data. A detailed examination of the annual TPO estimates over the 30-year 62 study period revealed that the TPO values had large variations over time (Figure 3-9). The highest TPO values in Georgia were found before 1990 and after 2010, but those years were not represented by the TPO survey data. This may explain why the 30-year averages were substantially higher than the averages from the survey data. The slightly higher 30-year average values than averages from the survey data in Alabama, Florida, Mississippi, and South Carolina were likely due to similar reasons. 18000 16000 14000 12000 SP SS 10000 SF 8000 HP 6000 HS 4000 HF 1:1 line 2000 0 0 5000 10000 15000 30-year Average TPO (x 1000 m3) Figure 3-8 Comparison of the average annual TPO values calculated from the 30-year record derived through this study and available TPO survey data. Each point represents one wood product type in one state. The two right most points were from Georgia. 63 Annual Average TPO from Survey Data (x 1000 m3) Figure 3-9 Annual timber production volume in softwood, hardwood, and total for seven states: (a) North Carolina; (b) South Carolina; (c) Alabama; (d) Florida; (e) Georgia; (f) Mississippi; (g) Tennessee. The temporal dynamics of the TPO the 30-year study period had similarities and differences among the 7 states (Figure 3-9). The 2007-2008 economic downturn appeared to have large impact on wood product output across the region. All 7 states had the lowest TPO values around that time, which were followed by sharp increases in the next few years. Several states also had high TPO values before or around early 64 1990s that were followed relatively lower TPO values throughout much of the 1990s and 2000s, including Florida, Georgia, and South Carolina, whereas Tennessee?s TPO values were relatively high during that period. Both Alabama and Mississippi had high TPO values over multiple years that peaked in early to mid-1990s, which were followed by a smaller peak centered around 2005. 3.5 Discussion and Conclusions This chapter explored three approaches for deriving annual TPO estimates, including a fixed effects linear regression (FELR) model, the ordinary least square (OLS) linear regression, and the Random Forest (RF) regression tree algorithm. The FELR method produced abnormally large values that were obvious errors in 5 of the 7 states, indicating that this method was not appropriate for TPO modeling, likely due to its sensitive to noises in the predictor variables and/or calibration data. Both OLS and RF produced more stable TPO estimates. Cross-validation results revealed that on average, the coefficient of determination (R2) between predicted and survey data was improved by 0.08 when TPO was modeled using RF instead of the OLS. Therefore, the RF was used to derive annual TPO estimates from remote sensing-based disturbance products over three decades (1986 to 2015). Predictions by the RF algorithm revealed that from 1986 to 2015, the 7 states in the Southeast produced more than 5 billion m3 wood products, of which over 40% were contributed by Georgia and Alabama. The average annual TPO estimates calculated from the 30-year record were comparable to those derived from the TPO survey data, with the 30-year averages being higher than the survey data in some states. This was 65 mainly due to large temporal variations of the TPO estimates over the 30-year study period, with many of the years in the first and/or last 5-10 years of the study period having high TPO values that were not represented by the survey data. The TPO dataset derived through this study represents multiple improvements over the available survey data. First, the time span of the derived data record (1986- 2015) is twice as long as that covered by the TPO survey data used in this study (1995- 2009). Second, the annual time step of the derived data can better capture the temporal variability of the TPO than the survey data available for each state, which typically skipped one or more years (Table 3-1). Third, the years during which survey data are available in one state are not always the same as those in other states, making it difficult to assemble a complete TPO dataset that covers all 7 states for some years. As an annual dataset derived consistently across the region and over time, the TPO record derived through this study can provide a complete TPO dataset for the entire region in any given year between 1986 and 2015. Modeled based on forest disturbances mapped using satellite imagery, the TPO data derived through this study can provide an observational basis for calculating the amount of C transferred from the standing biomass to the wood product through logging and for partitioning of harvested C among different wood product pools. As an example, Figure 3-10 shows the total amount of C calculated from the derived TPO record that were transferred to the P10 (C released in 10 years) and P100 (C released in 100 years) pools over the 30-year study period. While the amount transferred to the two product pools were comparable in most counties, many counties had higher shares of the harvested C in the P10 pool than in the P100 pool, or vice versa. The TPO record 66 derived through this study will be used to improve the estimates of C fluxed related to wood products in Chapter 4. Figure 3-10 Total C influx to the P10 (C released in 10 years) and P100 (C released in 100 years) pools calculated using the derived TPO record over the 30-year study period. The blue ellipsoids highlight the counties that had more harvested C allocated to the P10 pool than the P100 pool, or vice versa. 67 With demonstrated robustness over 7 states in the Southeast, the RF-based modeling framework developed in this study might be applicable in other regions of the United States. Annual forest disturbances have already been mapped for the conterminous U.S. (CONUS) using time series Landsat observations from 1986 to 2010 (Zhao et al., 2018). As of the Huang et al. (2015) study, every state in CONUS had TPO survey data for at least two years. Therefore, it is possible to expand the geographic domain of this study to the entire CONUS. It should be noted that several data products that became available recently might be better than the data used in this study and hence could be used to improve TPO modeling. One is a disturbance intensity dataset developed by Tao et al. (2019), which quantifies the percentage of the basal area removed by disturbance events mapped by VCT. Calibrated using basal area removal data calculated from field observations, this dataset should provide a more reliable measure of logging intensity that in turn should result in better estimation of the timber output from logging. Another dataset is the disturbance attribution data developed by Schleeweis et al. (2020), which explicitly mapped forest harvest and several other disturbance agents across CONUS. Excluding disturbances that are unlikely to result in timber output likely will improve TPO estimation. Another desirable dataset is annual aboveground biomass. Such a dataset could be used to calculate pre-harvest C stock, which is needed together with a measure of disturbance intensity to calculate harvested C. 68 Chapter 4 GCA Model Improvement and Implications for Carbon Management and Flux Estimation over North Carolina 4.1 Introduction The Grid-based Carbon Accounting (GCA) model developed in Chapter 2 provided a framework for estimating C fluxes over areas affected by forest disturbances at the 30m resolution, which is needed to calculate the C fluxes arising from harvest or other disturbances over individual forest patches. Such patch level C estimates are crucial for supporting various carbon management, carbon trade, and other climate change mitigation initiatives at individual project and/or individual property owner levels (Birdsey et al. 2006; Corbera and Schroeder 2011). In order to derive C estimates for all forest lands, however, the fluxes over forests not affected by contemporary disturbances should also be included. Further, the annual TPO records derived in Chapter 3 provided an opportunity to improve C flux estimates related to different wood product pools. Following the original bookkeeping carbon accounting (BCA) model, the GCA model uses fixed fractions when dividing harvested wood products into different release pools. In reality, the release pool fractions of carbon sources in different grids could be different. Derived using satellite-based annual forest disturbance data and TPO survey data, the TPO records from Chapter 3 provide year and county specific estimates of fuelwood, pulpwood, and sawlog volumes. These records could be used to derive more realistic partitioning of the total biomass removed by disturbance events into the P1, 69 P10 and P100 release pools (i.e., pools with C assumed to release in 1, 10, and 100 years, respectively). This chapter had two major goals. The first was to further improve the GCA model such that it could produce flux estimates for pre-disturbance growth over disturbed areas as well as the growth of undisturbed forests, and to use the TPO records developed in Chapter 2 to derive more realistic partitioning of the total biomass removed by disturbance events into the P1, P10 and P100 release pools. Following a comparison of model outputs to available inventory-based flux estimates, the growth rates used by the model were also adjusted according to rates calculated using FIA field plot data. The second goal was to use the improved model to evaluate how forest C flux estimates were affected by forest management practices as well as emission calculation methods. The spatial and temporal domains of this chapter were exactly the same as in Chapter 2: the study area included the entire state of North Carolina and the years ranged from 1986 to 2010. Except for those described in the following sections, all other model inputs were exactly the same as those used in Chapter 2. 4.2 GCA Model Improvement 4.2.1 Partitioning of Removed C The aboveground carbon removed by disturbance events can be divided into four pools from which the C is released through different processes at different rates ? slash, P1, P10, and P100 (Table 2-2). The original BCA model provided the partitioning ratios among the four pools for each ecozone, which were assumed to remain the same across counties and over time within each ecozone. In reality, however, these ratios 70 could differ from county to county and vary from year to year. The TPO data provided an opportunity to derive county and year specific partitioning among the P1, P10, and P100 pools. Specifically, softwood and hard fuelwoods should go to the P1 pool because they likely will be burned in one year. The pulpwood might store C for up to a decade and should be in the P10 pool. Sawlogs could store C in furniture and buildings for over a century, and hence should be in the P100 pools. No estimate of slash was provided by the TPO survey data. In theory, however, slash should be the difference between removed C and the amount of C stored in wood products: Cslash = Cvct-Ctpos (4-1) Where, Cvct, Cslash, and Ctpos are the total C removed by disturbance events mapped by the VCT algorithm and the amounts stored in slash and wood products, respectively. In any given year, Cvct for a county was calculated as the sum of the carbon removed at individual pixel locations (Cpj) for all pixels within that county (4-2): Cvct=? Cpj (4-2) For each TPO survey year, Ctpos was provided by the TPO survey data. Therefore, a slash ratio (Rslash) could be calculated for each TPO survey year for every county: Rslash = Cslash/Cvct (4-3) The county specific Rslash values calculated this way, however, varied greatly among counties and TPO survey years. Assuming the Rslash value was mostly determined by harvest methods, which were unlikely to differ that much among counties in the same state, the large among-county and among-TPO survey year differences were more likely caused by uncertainties in both the disturbance products 71 and the pre-disturbance C data. To minimize the impact of such uncertainties, an average Rslash was calculated based on the county and TPO survey year specific values. The resultant Rslash value, 0.4672, was then used to calculate the influx of C to the slash pool from the total C removed by each year?s disturbances for each county (CVslash): CVslash = Cvct x Rslash (4-4) After subtracting this slash pool from the total removed C, the P1, P10, and P100 pool fractions were calculated annually for each county by ratios of C in different wood products based on the annual TPO record derived in Chapter 2. Finally, the C removed at each pixel location within a county was then partitioned according to the above derived county and year specific fraction values. Release of C stored in each product pool was calculated the same way as described in Chapter 2. 4.2.2 C accumulation from forest growth Forests transfer C from the atmosphere to the biomass pool through tree growth. Over disturbed areas, the GCA model developed in Chapter 2 used pre-defined growth rates to calculate the C accumulated by regrowing trees reestablished after each logging and/or fire disturbance event. The growth before the first disturbance event, however, was not tracked, nor was the growth over undisturbed forest areas. In this chapter, both of these two cases were included in calculating the C fluxes arising from forest growth. 72 4.2.3 Model verification and adjustment While the GCA model could produce C estimates at the 30m spatial resolution, validation of fine resolution flux estimates remains a challenge. Flux estimates with spatial details close to the 30m resolution are difficult to find in the published literature. Such estimates are often reported at global, national, or regional scales (Houghton and Nassikas, 2017; Woodall et al., 2015). However, state level details have been reported in some studies (Han et al., 2007; Houghton and Nassikas, 2017; Woodall et al., 2015; Zheng et al., 2011), which could be useful for evaluating the results derived in this study. In particular, the U.S. Forest Service published an inventory of greenhouse gas emissions and removals from forest land, woodlands, and urban trees (Domke et al. 2020). This report provided annual flux estimates with state level details from 1990 to 2018. Estimates were derived for both forest land remaining forest land and forest land conversions, and included fluxes among aboveground biomass (including harvested wood products), belowground biomass, and other major carbon pools. Those related to aboveground biomass over forest land remaining forest land and forest land conversion to urban land between 1990 and 2010 were used to evaluate the results derived through this study. Figure 4-1 shows a comparison between the GCA model and inventory data on carbon changes in the aboveground biomass pool over forest land remaining forest land. Overall, the aboveground biomass pool had large negative fluxes (sink) during the study period. But the GCA estimates were much larger. The annual fluxes were twice as much as the inventory data in early 1990s and remained about 50% larger by 73 2010. As a result, the 21-year total sink estimate from the GCA model was 80% larger than the inventory-based estimate. 0.00 -2.00 -4.00 -6.00 -8.00 -10.00 -12.00 -14.00 -16.00 -18.00 Year GCA USFS Inventory GCA-Adjusted for Legacy Emission Figure 4-1 Comparison of annual C changes in the aboveground (including harvested wood products) pool derived using the GCA model and from the USFS inventory data over forest land remaining forest land. The slope of the increasing trend in the annual flux estimates by the GCA model was reduced after adjusting the estimates to account for legacy emission from wood products harvested during 100 years before each model year (see the text for details on the adjustment method). Compared to the relatively stable inventory data, the annual sink derived from the GCA model had an obvious decreasing trend over the 21-year period. This is because the GCA model did not include legacy emissions from wood products harvested before 1986. As the model started to accumulate the wood product pool in 1986, the pool increased over time, and so did the model estimate of the legacy emission from that pool, roughly following a logarithmic curve (Figure 4-2). 74 C Flux (Tg) 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Figure 4-2 Because the GCA model did not include legacy emissions from wood products harvested before 1986 and started to accumulate the wood product pool in 1986, model estimates of legacy emission from wood products increased as the number of years since 1986 increased, roughly following a logarithmic curve. Assuming all factors controlling the influx of C to the wood product pools and emission from the pools between 1986 and 2010 would not change that much in future years, the fitted curve could be used to estimate the legacy emission in any year after 1986 from wood products harvested between 1986 and that year. For example, by the 100th year after 1986, the total annual legacy emission from wood products harvested in the 100 years after 1986 would be 6.8Tg. Assuming the same method could also be used to estimate the legacy emission in any year between 1986 and 2010 from wood products harvested 100 years before that year, then it would be possible to adjust the flux value estimated by the GCA model for each year such that the value would include legacy emission from wood products harvested during the 100 years before that year. Figure 4-1 shows that after this adjustment, the annual C sink values from the GCA model were still 2.5Tg ? 4.5Tg larger than the inventory data, indicating that other 75 factors, e.g., inflated growth rates, could have also contributed to the larger sink estimates by the GCA model. As shown in Table 2-3, the growth rate used by the model for young forests (fast growth) in ecozone 21 was close to the mean growth rate calculated from FIA plot data. However, that mean growth rate likely was overestimated, because a good number of FIA plots appeared to have abnormally high (5-10 times higher than the mean) growth rates (Figure 2-3). Such high values had disproportionately large impact on the mean value, but they were errors because no trees could grow that fast. Because medians are typically less affected by extremely large or small values, the median growth rates from Table 2-3 were used to replace the original growth rates used by the GCA model for young forests (fast). After this adjustment, annual sink estimates from the GCA model were still mostly larger than inventory data, but the differences were much smaller than those derived using the original growth rates (Figure 4-3). After compensating for legacy emissions from wood products harvested during the 100 years before each model year, the annual sink estimates became slightly smaller than the inventory data, with differences ranging from 1Tg and 2Tg (Figure 4-3). These differences were less than half of the 2.5Tg-4.5Tg differences derived above using the original growth rates in the GCA model, suggesting that the growth rates of forests across North Carolina were better represented by the FIA-based median values than the original growth rates used by the model. Therefore, the FIA-based median growth rates were used in the GCA model in the remaining analysis in this chapter. 76 0.00 -2.00 -4.00 -6.00 -8.00 -10.00 -12.00 Year GCA Inventory GCA-Adjusted for Legacy Emission Figure 4-3 Comparison of annual C changes in the aboveground (including harvested wood products) pool derived using the GCA model with growth rates adjusted based on FIA plot data and from the USFS inventory data over forest land remaining forest land. After adjusting the estimates to account for legacy emission from wood products harvested during 100 years before each model year, the annual sink estimated by the GCA model became slightly smaller than the inventory-based estimates. Compared to the initial GCA model developed in Chapter 2, the improved model developed in this chapter not only included fluxes from growth of pre-disturbance forests and undisturbed forests, but also adjusted forest growth rates and the method for calculating the influx of harvest carbon to slash and the wood product pools. Therefore, both sources and sinks derived from the improved model were different from those derived using the initial GCA model developed in Chapter 2, even over the disturbed areas (Figure 4-4). 77 C Flux (Tg) 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 Logging Source 10.0 8.0 6.0 4.0 2.0 0.0 Year Initial model (Chapter 2) Improved model (Chapter 4) Logging Sink 12.0 10.0 8.0 6.0 4.0 2.0 0.0 Year Initial model (Chapter 2) Improved model (Chapter 4) Figure 4-4 Compared to the initial GCA model developed in Chapter 2, improvements made in this chapter resulted in higher source and lower sink estimates for forest land disturbed by logging. The same was true for conversion to urban land (no sink) and forest land subject to fire disturbances. 4.3 North Carolina?s Forest Carbon Dynamics 4.3.1 Source and Sink Analysis The improved GCA model was used to produce C flux estimates at the 30m resolution for North Carolina. The results showed that without considering soil carbon 78 Carbon (Tg) Carbon (Tg) 1986 1986 1987 1987 1988 1988 1989 1989 1990 1990 1991 1991 1992 1992 1993 1993 1994 1994 1995 1995 1996 1996 1997 1997 1998 1998 1999 1999 2000 2000 2001 2001 2002 2002 2003 2003 2004 2004 2005 2005 2006 2006 2007 2007 2008 2008 2009 2009 2010 2010 as well as the legacy emissions from wood products harvested before 1986 and the emissions from wood products harvested from 1986 to 2010 that will be released after 2010, North Carolina?s forest land was a net C sink of 218.1Tg over the 25-year (1986- 2010) study period. Major sources included logging (136.7Tg), conversion of forest to urban land (6.9Tg), and fire (2.1Tg). About 55% of the forests in North Carolina remained undisturbed during the study period. The growth of these forests resulted in the largest sink (179.1Tg). While logging resulted in the largest source, post-logging regrowth and the C removed by logging that was transferred to the harvested wood product pool and was not released by 2010 created a sink of 103.2Tg, and the sink from post-fire growth was 0.7Tg. Further, pre-disturbance forest growth over the disturbed areas resulted in a sink of 80.8Tg. Figure 4-5 shows the contributions of each source and sink process to the total source and sink estimates over the 25-year study period. 79 Source Converstion to urban land Fire 5% 1% Logging 94% Logging Fire Converstion to urban land Pre-disturbance growth Sink 22% Post-logging growth 29% Post-fire growth Undisturbed forest 0% 49% Post-logging growth Post-fire growth Undisturbed forest Pre-disturbance growth Figure 4-5 Contributions of different source and sink processes to the total source and sink over North Carolina from 1986 and 2010. On average, the forest land absorbed 14.6Tg per year between 1986 and 2010, including 7.2Tg by undisturbed forests, 4.1Tg from post-disturbance growth and C stored in harvested wood products not released by 2010, and 3.2Tg from pre- disturbance growth. The average annual release from logging was 5.5Tg, while fire and conversion to urban land created annual sources of 0.08Tg and 0.3Tg, respectively. 80 Overall, the forest land in North Carolina had an average net C flux of -8.7Tg/year (sink) between 1986 and 2010. This sink estimate seems to be comparable to those reported in a few studies that provided flux estimates for North Carolina. For example, Zheng et al. (2011) reported an average net C change of -9.7Tg/year for forest areas in North Carolina that remained forested between 1992 and 2001. Without specifying a year range, Han et al. (2007) estimated that the ?current annual? C change in forest biomass was -8.5Tg/year. While it is not clear how legacy emissions from wood products were analyzed in those studies, the net sink estimate from this study should be smaller if legacy emissions from all wood products harvested over 100 years were to be fully accounted for (see section 4.2.3). 4.3.2 Impact of management approaches and emission calculation methods While the GCA model provides a way to derive high resolution source and sink estimates for forest land, those estimates could be affected by many factors. In chapter 2, it was demonstrated that the use of remote sensing-based biomass and disturbance intensity data had large impact on the estimation of fluxes driven by disturbance and post-disturbance recovery. The use of newly available TPO records developed in Chapter 3 and field-based growth rates also resulted in estimates that were different from the earlier model (Figure 4-4) and were closer to inventory data (Figures 4-1 and 4-3). While forest management in general can improve carbon storage (Canadell and Raupach, 2008), different harvest and management practices have different C outcomes (Naudts et al., 2016; Nunery and Keeton, 2010; Van Deusen, 2010). Other important factors to consider in derive flux estimates include land-cover dynamics and methods for expressing emissions from forest disturbances (Ramankutty et al., 2007). 81 The GCA model provided an opportunity to examine the C outcome of different logging intensity over North Carolina, and how the estimation of the outcome could be affected by the use of different emission expression methods. Four harvest intensity methods were considered, including clear cut (disturbance intensity = 100%), observed disturbance intensity represented by the disturbance intensity data developed by Tao et al. (2019), which had already been used in Chapter 2, half of the observed disturbance intensity, and zero disturbance (i.e., no logging). The emission expression methods considered included ?prompt emission? (C released in the year of clearing) and ?delayed emission? (C released following decay processes in multiple years) (Fearnside, 1997). The combination of logging intensity and emission expression methods considered in this study resulted in 8 scenarios (Table 4-1): Table 4-1 Scenarios for assessing the impact of logging intensity and emission expression methods on forest C outcome estimates. Disturbance Intensity Delayed Release Prompt Release Clear Cut (100%) CD CP Mapped Intensity (Tao et al. MD MP 2019) Half of Mapped Intensity HD HP Zero Intensity ZD ZP As expected, clear cut resulted in the highest fluxes from logging, followed by logging at the mapped intensity, half of the mapped intensity, and zero intensity (no 82 logging). Except for the zero intensity case, flux estimates derived using the prompt emission method were always higher than those derived using the delayed emission method (Figure 4-6) for each of the other three logging intensities. The annual fluxes derived using the delayed emission method for each of these three logging intensity levels had an increasing trend over time. As discussed earlier, this was mainly because the GCA model started to accumulate the wood product pools in 1986, and did not include legacy emissions from wood products harvested before 1986. The increasing trends were mostly the result of a growing wood product pool, which then resulted in higher emissions from that pool as the number of years since 1986 increased. Fluxes derived using the prompt emission method had larger temporal variability than those derived using the delayed emission method but did not have obvious temporal trends. 25.00 20.00 15.00 10.00 5.00 0.00 Year MD CD HD ZD 83 C Flux (Tg) 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 25.00 20.00 15.00 10.00 5.00 0.00 Year MP CP HP ZP Figure 4-6 Annual carbon emissions derived using the delayed (top) and prompt (bottom) emission methods at the four disturbance intensity levels shown in Table 4-1. When fluxes were calculated using the delayed emission method, the logging source estimate over the logging areas for the 25-year study period would be ~60Tg more if those areas had clear cut instead of having the logging intensities mapped by Tao et al. (2019). On the other hand, that estimate would be reduced by ~60Tg if the observed logging intensities were reduced by half, and further down to 0 if there logging intensity were also reduced to 0 (Figure 4-7). Logging intensity had little impact on post-logging growth. The flux from logging at the observed intensity (scenario MD) and post-logging growth over the logging areas was a net source of 33.5Tg over 25 years. That source would become 95.8Tg if the logging practice were changed to clear cut. However, if the observed logging intensity were reduced by half, the sink from post-logging growth would exceed emissions from logging, resulting in a net sink of 28.6Tg. The sink from the disturbed areas would reach 94.5Tg should logging be eliminated during the study period. 84 C Flux (Tg) 1986 1987 1988 1989 1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008 2009 2010 350.00 250.00 150.00 50.00 -50.00 -150.00 -250.00 -350.00 Logging Source Sink from post-logging Net Flux over Logging Net Flux over All Forest growth Areas Land CD MD HD ZD 350.00 250.00 150.00 50.00 -50.00 -150.00 -250.00 -350.00 Logging Source Sink from post-logging Net Flux over Logging Net Flux over All Forest growth Areas Land CP MP HP ZP Figure 4-7 Total fluxes over the 25-year study period derived using the delayed (top) and prompt (bottom) emission methods at the four disturbance intensity levels shown in Table 4-1. The emission expression methods had no impact on the calculation of sink from post-logging growth. However, because the logging source estimates derived using the prompt emission method were substantially higher than those derived using the delayed emission method, logging source estimates derived using the prompt emission method 85 C Flux (Tg) C Flux (Tg) -344.05 -343.90 -279.41 -237.27 -218.06 -141.23 -156.82 -48.48 -94.54 -94.54 -28.63 13.36 33.49 110.18 95.78 203.99 94.54 94.54 101.85 101.85 103.21 103.21 103.98 103.98 0.00 0.00 73.22 115.22 136.70 213.39 199.76 307.96 for the three non-zero intensity cases exceeded estimates of the sink from post-logging growth, resulting in positive net flux estimates for all three logging intensity levels. Due to the large sink estimates from pre-disturbance growth and the growth of undisturbed forests, North Carolina?s forest land had a net sink estimate over the 25- year study period regardless of the logging intensity levels or the emission expression method used. The estimated size of the sink, however, decreased substantially as logging intensity increased from no logging to clear cut, and the decreasing pace was faster when flux estimates were derived using the prompt emission method than using the delayed emission method. 4.4 Summary and Conclusions Building on results derived in Chapters 2 and 3, a number of improvements were made to the GCA model such that it tracks the fluxes among all major aboveground C pools, including pre-disturbance growth and the growth of undisturbed forests, which were not included in the initial model developed in Chapter 2. Based on available TPO survey data and the annual TPO records derived in Chapter 3, efforts were also made to improve the partitioning of removed C among slash and different wood product pools. Specifically, an average slash ratio was calculated for North Carolina based on the difference between removed C and the influx of C to wood products, which was calculated from the TPO survey data. The annual TPO records derived in Chapter 3 were then used to calculate county and year specific ratios for partitioning the rest of the removed C among the P1, P10, and P100 pools, which were then applied to each 30m pixel based on which county and year that pixel belongs to. 86 While the GCA model was designed to use remote sensing-based products to estimate fluxes arising from contemporary disturbances and hence could not track legacy emissions from wood products harvested before the study period, the model did track emissions from wood products harvested during the study period, which provided a basis for estimating the amount of legacy emissions from wood products that were not accounted for by the GCA model. After compensating for these missed legacy emissions and adjusting forest growth rates based on FIA field plot data, the model produced estimates that were close to the greenhouse gas emission and sink inventory data provided by the USFS for North Carolina and comparable to estimates reported in several other studies. Without considering soil carbon as well as the legacy emissions from wood products harvested before 1986 and the emissions from wood products harvested from 1986 to 2010 that will be released after 2010, North Carolina?s forest land was a net C sink of 218.1Tg over the 25-year (1986-2010) study period. That sink would be smaller if all logging areas had clear cut, and could be larger if the logging intensities could be reduced by half or down to 0. In general, the prompt emission method resulted in larger source estimates from disturbances than the delayed emission method, and hence reduced the sink estimates at all disturbance intensity levels. 87 Chapter 5 Conclusions Designed with an overarching goal to improve the understanding of the terrestrial carbon budget, this study sought to improve the estimation of forest carbon dynamics in two fronts. One was to develop an improved carbon modeling capability ? a grid- based carbon accounting (GCA) model was developed based on the bookkeeping carbon accounting (BCA) model developed by Houghton et al. (1999). The initial GCA model only included modules for estimating fluxes arising from forest disturbance and post-disturbance recover (Chapter 2). Additional improvements were made in Chapter 4, which allowed the model to produce carbon estimates for all forest land. The other front was to derive new remote sensing-base data products, including an improved forest carbon density map that better represented pre-disturbance forest carbon than the original map, a series of annual forest attribution maps, and an annual timber product output (TPO) record that was much longer and more consistent than available TPO survey data. Together with a suite of disturbance datasets, these products were incorporated into the GCA model to produce carbon estimates over a quarter century (1986-2010) across North Carolina. 5.1 Summary of Major Findings Previous inter-model comparison studies showed that flux estimates derived using different carbon models are highly variable (Huntzinger et al., 2012). This study demonstrated that given the same model, different inputs can also result in vastly differences flux estimates. Among the major inputs/parameters that could have large impact on carbon estimates included pre-disturbance carbon density, disturbance 88 intensity, allocation of removed carbon among slash and different wood product pools, and forest growth rates. The first three inputs/parameters primarily affect emission estimation from disturbances. Compared to estimates derived using an initial GCA model that inherited the parameters provided by the original BCA model, the total emission between 1986 and 2010 from logging over North Carolina was reduced by almost one third and two thirds when those parameters were replaced by remote sensing-based disturbance intensity and biomass data, respectively, and when both products were used, the emission was reduced by more than three quarters. Use of the TPO data derived in Chapter 3 to partition the removed carbon among slash and different wood product pools resulted in noticeably high emission estimates than those derived using the partitioning ratios provided by the original BCA model. The emission expression method used also had considerable impact on the source estimation for logging. Without considering legacy effect from wood products harvested before 1986, the emission value derived using the prompt release method was 50% higher than that derived using the delayed release method. Of course, much of these differences would be gone if the legacy emissions from wood products harvested over a sufficiently long period (e.g. 100 years) were considered. Forest growth rates mainly affect the sink estimates from the growth of undisturbed forests as well as from pre- and post-disturbance growth. Compared with flux estimates reported by the greenhouse gas inventory program of the U.S. Forest Service, the growth rates prescribed by the original BCA appeared to be too high. After adjusting those rates based on FIA inventory data and accounting for legacy emissions from wood products harvested before 1986, the GCA model produced flux estimates 89 comparable to those reported by the USFS greenhouse gas inventory program. Without considering soil carbon, the forest land in North Carolina had an average annual sink of 5.5Tg C/Year from 1990 to 2010, compared to the 6.7Tg C/Year sink reported by the USFS. Alternative logging practices could have significant impact on the size of this sink. Without considering the legacy emissions from wood products harvested before 1986 and the emissions from wood products harvested from 1986 to 2010 that will be released after 2010, the sink would be reduced by ~30% if all logging areas had clear cut, but would be increased by ~30% or ~60%, if the logging intensities could be reduced by half or down to 0, respectively. Of the three approaches explored in this study for TPO modeling, the fixed effects linear regression (FELR) model produced abnormally large values that were obvious errors, and hence was deemed inappropriate for TPO modeling. Both the ordinary least square (OLS) linear regression and the Random Forest (RF) regression tree algorithm produced more stable TPO estimates, with the coefficient of determination (R2) values derived using RF being 0.08 higher than those derived using the OLS on average. Predictions by the RF algorithm revealed that from 1986 to 2015, the 7 states in the Southeast produced more than 5 billion m3 wood products, over 90% of which were pulpwood and sawlogs that were assumed to release carbon in 10 years (P10) and 100 years (P100), respectively. While the influxes to the P10 and P100 were roughly the same for the entire Southeast over the 30-year study period, the ratio of the influxes to the two pools had substantial variations from year to year and among counties and states. 90 5.2 Significance of Original Contributions A major contribution of this research is the development of the GCA model. By integrating the BCA model with spatially detailed remote sensing products within a grid-based framework, this model made it possible to utilize a well-established carbon accounting method designed for use at regional to global scales to produce fine resolution carbon estimates. As with the original BCA model, the GCA model tracks carbon fluxes arising from natural disturbances (e.g., fire), conversion, and forest management (e.g. logging). Given that about half of the forest area in the U.S. is disturbed each decade (Pan et al., 2011), being able to account for the fluxes from different disturbance types should greatly reduce the uncertainties with the flux estimates over forest land (Turner et al. 2016). A key input for calculating the carbon flux from a given disturbance event is the carbon density before that disturbance event. Because disturbances can happen in any year within a relatively large study region (e.g., a county or state), ideally annual carbon maps for an entire study period should be required by the GCA model. Unfortunately, no such annual carbon maps exist for most areas. Use of a carbon map developed for a specific year in other years can result in substantial uncertainties, because the carbon density in most areas changes from one year to another due to growth or disturbance. Assuming that larger trees are more likely to be harvested than smaller ones, a local window-based method was developed in this study to derive more realistic pre- disturbance carbon values for any year using a carbon map available for a specific year. The pre-disturbance carbon values derived using this method were much closer to those derived from FIA plot data than values from the carbon map available in a specific year 91 or the values prescribed in the original BCA model. This method should be useful for other carbon studies where pre-disturbance or relatively mature forest carbon values in different years are needed but carbon maps are only available for one or a few selected years. While both the ordinary least square (OLS) linear regression and fixed effects linear regression (FELR) have been used to model TPO from satellite-based disturbance data in previous studies, it was demonstrated in this study the Random Forest algorithm was more robust for TPO modeling over 7 states in the southeast United States. The TPO record derived through this study represents multiple improvements over the survey data used in this study. With a time span (1986-2015) twice as long as that covered by the survey data (1995-2009), it provides a more complete picture of the influx to and emissions from wood product pools over three decades across the entire Southeast. The annual time step of the derived data revealed large temporal variations that were not captured by the survey data and makes it straightforward to assemble a complete TPO dataset for any given year of the study period that covers all 7 states, which was not possible by using the survey data alone for most years. Derived based on ground survey data and time series satellite observations, this TPO record made it possible to develop observation-based partitioning of harvested carbon among different release pools for each county in each year, which should be more realistic than the fixed partitioning ratios prescribed in the original BCA model. Therefore, it should be highly valuable for future efforts aimed to improve carbon estimates in the southeast region. Given that Landsat-based forest change maps can be and have been produced across the globe or for any land area of 92 the earth, the modeling approach should be applicable to any other regions where at least some survey data on wood products are available. Given that carbon management decisions are often made at local, community, or individual property levels, the 30m carbon maps produced by the GCA model should provide adequate spatial details to support decision makings at these levels. The grid size of the GCA model can be adjusted for use with remote sensing products of any spatial resolutions. As remote sensing technology continues to evolve rapidly, new products with better quality and finer resolutions will continue to emerge. The GCA model provides a framework for integrating new products as they become available to support various carbon management applications, including calculating carbon credits for carbon trade, and supporting various climate change initiatives regarding the measurement, reporting, and verification of carbon pools and fluxes (Birdsey et al., 2006; Fahey et al., 2010; Lamb et al., 2021). 5.3 Limitations and Future Research Directions Despite the many advantages of the GCA model as compared to the original BCA model, there is plenty of room for further improvement. For example, the new model still uses a single set of growth rates for each ecozone. In reality, however, forest growth rates likely are more variable across space and over time. An integrated analysis of the FIA plot data with time series satellite observations and climate data may provide growth rate estimates that are more variable and likely more realistic. While a method was developed to determine the pre-disturbance carbon value for any disturbance year based on an available forest carbon map developed for a specific year, 93 annual forest carbon maps, if derived properly, likely will provide more realistic representation of pre-disturbance carbon values, and hence are more desirable. With the launch of Sentinel-1 back in 2014, global fine resolution radar observations have become routinely available, which in general should be more sensitive to forest biomass than optical data. Since 2018, global samples of LiDAR measurements have been collected by GEDI and ICESAT-2. As LiDAR measurements in general are better correlated with biomass than optical and radar data, these samples could be used together with ground measurements to calibrate optical and/or radar images to produce annual biomass maps. While the Random Forest algorithm in general performed better than the OLS algorithm in TPO modeling, relationships between the predicted values and survey data were not always that good. In particular, the hardwood types appeared to be more difficult to model for most states. Improvements might be achievable using two products. One is the disturbance intensity dataset developed by Tao et al. (2019). The other is a pre-disturbance carbon density map, like the one developed in Chapter 2. Because logging is the dominant disturbance agent to the forests in the Southeast, the logging emission in any given year depends to a large degree on the amount of carbon harvested in that year transferred to wood products and the amount of legacy emission from wood products harvested in previous years. While the GCA model cannot account for all legacy emissions by itself, it does track the portion from wood products harvested during a study period. For a sufficiently long study period (e.g., 25 years), results from the model could be used to estimate the total legacy emission from wood products harvested over a much longer period (e.g., 100 years), assuming the 94 wood product pools were relatively stable over time. Knowing this assumption likely is not true, historical TPO or TPO-like data, if available, should be used to calculate the legacy emission from wood products harvested before the beginning of the study period. Despite these limitations and some desirable new improvements, the GCA model is ready for use in other regions, at least other states in the Southeast, to produce high resolution carbon estimates. Such estimates should provide more definitive (likely with less uncertainty) answer to questions on the carbon dynamics of the forests in this region. More importantly, the fine spatial resolution of the derived products should provide sufficient details to support carbon management applications at a broad range of geographical scales, including the measurement, reporting, and verification of carbon pools and fluxes at local, community, or even individual property owner levels. 95 Bibliography AGO (Australian Greenhouse Office) 2000a. National Carbon Accounting System: Phase 1 Implementation Plan for the 1990 Baseline. NCAS Technical Report No. 10. Australian Greenhouse Office, Canberra,56pp. AGO (Australian Greenhouse Office) 2000b. International Review of the Implementation Plan for the 1990 Baseline. NCAS Technical Report No. 11. Australian Greenhouse Office, Canberra, 16pp. AGO (Australian Greenhouse Office) 2000c. System Design. NCAS Technical Report No. 21. Australian Greenhouse Office, Canberra, 25pp. AGO (Australian Greenhouse Office) 2000d, THE FULLCAM CARBON ACCOUNTING MODEL: DEVELOPMENT, CALIBRATION AND IMPLEMENTATION FOR THE NATIONAL CARBON ACCOUNTING SYSTEM, NCAS Technical Report. Agrawal, A., Nepstad, D., & Chhatre, A. (2011). Reducing Emissions from Deforestation and Forest Degradation. Annual Review of Environment and Resources, 36, 373-396. Ahmed, O. S., Franklin, S. E., Wulder, M. A., & White, J. C. (2015). Characterizing stand- level forest canopy cover and height using Landsat time series, samples of airborne LiDAR, and the Random Forest algorithm. ISPRS Journal of Photogrammetry and Remote Sensing, 101, 89?101. https://doi.org/10.1016/j.isprsjprs.2014.11.007 Anand, A. (2016). Effect of Vegetation Structure on Undercanopy Solar Radiation Using LiDAR Remote Sensing.Anderson, J., Martin, M.E., Smith, M.L., Dubayah, R.O., Hofton, M.A., Hyde, P., Peterson, B.E., Blair, J.B. & Knox, R.G. (2006). The use of 96 waveform lidar to measure northern temperate mixed conifer and deciduous forest structure in New Hampshire. Remote Sensing of Environment, 105, 248-261. Antonarakis, A. S. ?Uncertainty in Initial Forest Structure and Composition When Predicting Carbon Dynamics in a Temperate Forest.? Ecological Modelling 291 (2014): 134?41. doi:10.1016/j.ecolmodel.2014.07.030. Austin, J.M., Mackey, B.G., & Van Niel, K.P. (2003). Estimating forest biomass using satellite radar: an exploratory study in a temperate Australian Eucalyptus forest. Forest Ecology and Management, 176, 575-583. Baccini, A., Goetz, S. J., Walker, W. S., Laporte, N. T., Sun, M., Sulla-Menashe, D., ? Houghton, R. A. (2012). Estimated carbon dioxide emissions from tropical deforestation improved by carbon-density maps. Nature Climate Change, 2(3), 182? 185. https://doi.org/10.1038/nclimate1354 Balzter, H. (2001). Forest mapping and monitoring with interferometric synthetic aperture radar (InSAR). Progress in Physical Geography, 25, 159-177. Balzter, H., Rowland, C.S. & Saich, P. (2007). Forest canopy height and carbon estimation at Monks Wood National Nature Reserve, UK, using dual-wavelength SAR interferometry. Remote Sensing of Environment, 108, 224-239. Balzter, H., Skinner, L., Luckman, A., & Brooke, R. (2003). Estimation of tree growth in a conifer plantation over 19 years from multi-satellite L-band SAR. Remote Sensing of Environment, 84, 184-191. Bardon, R.E., Megalos, M.A., New, B., & Brogan, S. (2010). North Carolina?s Forest Resources Assessment. In (p. 489). Raleigh, North Carolina: NC Division of Forest Resources. 97 Barford, C.C., Wofsy, S.C., Goulden, M.L., Munger, J.W., Pyle, E.H., Urbanski, S.P., Hutyra, L., Saleska, S.R., Fitzjarrald, D. & Moore, K. (2001). Factors controlling long- and short-term sequestration of atmospheric CO2 in a mid-latitude forest. Science, 294, 1688-1691. Bari, M.A., Smettem, K.R.J., (2006). A conceptual model of daily water balance following partial clearing from forest to pasture. Hydrol. Earth Syst. Sci. 10, 321?337. Bengtsson, J., Nilsson, S.G., Franc, A., Menozzi, P., (2000). Biodiversity, disturbances, ecosystem function and management of European forests. For. Ecol. Manag. 132, 39? 50. Bigsby, H. (2009). Carbon banking: Creating flexibility for forest owners. Forest Ecology and Management, 257(1), 378?383. https://doi.org/10.1016/j.foreco.2008.09.018 Birdsey, R. (2004). Data Gaps for Monitoring Forest Carbon in the United States: An Inventory Perspective. Environmental Management, 33, S1-S8. Birdsey, R.A. & Heath, L.S. (1995). Carbon changes in U.S. forests. In: L.A. Joyce (Ed.), Productivity of America's forests and climate change (pp. 56-70), Fort Collins, Colorado USDA Forest Service. Birdsey, R.A., Lewis, G.M., 2003. Current and historical trends in use, management, and disturbance of U.S. forestlands. In: Kimble, J.M., Heath, L.S., Birdsey, R.A., Lal, R. (Eds.), The Potential of U.S. Forest Soils to Sequester Carbon and Mitigate the Greenhouse Effect. CRC Press, New York, pp. 15?33. Birdsey, R., Pregitzer, K., & Lucier, A. (2006). Forest Carbon Management in the United States. Journal of Environmental Quality, 35(4), 1461?1469. https://doi.org/10.2134/jeq2005.0162 98 Biswas, S. (2016). Contemporary Forest Cover Dynamics in Myanmar.Blackard, J.A., Finco, M.V., Helmer, E.H., Holden, G.R., Hoppus, M.L., Jacobs, D.M., Lister, A.J., Moisen, G.G., Nelson, M.D., Riemann, R., Ruefenacht, B., Salajanu, D., Weyermann, D.L., Winterberger, K.C., Brandeis, T.J., Czaplewski, R.L., McRoberts, R.E., Patterson, P.L. & Tymcio, R.P. (2008). Mapping US forest biomass using nationwide forest inventory data and moderate resolution information. Remote Sensing of Environment, 112, 1658-1677. Boudreau, J., Nelson, R.F., Margolis, H.A., Beaudoin, A., Guindon, L. & Kimes, D.S. (2008). Regional aboveground forest biomass using airborne and spaceborne LiDAR in Qu?bec. Remote Sensing of Environment, 112, 3876-3890. Bradford, J. B., R. A. Birdsey, L. A. Joyce, and M. G. Ryan. ?Tree Age, Disturbance History, and Carbon Stocks and Fluxes in Subalpine Rocky Mountain Forests.? Global Change Biology 14, no. 12 (2008): 2882?97. doi:10.1111/j.1365-2486.2008.01686.x. Bradford, J.B., Weishampel, P., Smith, M.-L., Kolka, R., Birdsey, R.A., Ollinger, S.V. & Ryan, M.G. (2010). Carbon pools and fluxes in small temperate forest landscapes: Variability and implications for sampling design. Forest Ecology and Management, 259, 1245-1254. Breiman, L. Bagging predictors. Machine Learn. 1996, 24, 123?140. Breiman, L., Random forest. Mach. Learn. 2001,45. Breiman, L.; Cutler, A. Random Forests. Available online: http://www.stat.berkeley.edu/~breiman/Random Forests/cc_home.htm#inter January 2017. 99 Brown, M.J., New, B.D., Johnson, T.G., & Chamberlain, J.L. (2014). North Carolina's forests, 2007. Asheville, NC: USDA-Forest Service, Southern Research Station (p. 112): Resour. Bull. SRS-RB-199. Cairns, M.A., Brown, S., Helmer, E.H. & Baumgardner, G.A. (1997). Root biomass allocation in the world's upland forests. Oecologia, 111, 1-11. Cameron, A. Colin; Trivedi, Pravin K. (2005). Microeconometrics: Methods and Applications. Cambridge University Press. pp. 717?19. Canadell, J.G., & Raupach, M.R. (2008). Managing forests for climate change mitigation. Science, 320, 1456-1457. Castel, T., Guerra, F., Caraglio, Y., & Houllier, F. (2002). Retrieval biomass of a large Venezuelan pine plantation using JERS-1 SAR data. Analysis of forest structure impact on radar signature. Remote Sensing of Environment, 79, 30-41. CCSP (2007). The First State of the Carbon Cycle Report (SOCCR): The North American Carbon Budget and Implications for the Global Carbon Cycle, Asheville, NC: National Oceanic and Atmospheric Administration, National Climatic Data Center. Chander, G., Markham, B.L., & Helder, D.L. (2009). Summary of current radiometric calibration coefficients for Landsat MSS, TM, ETM+, and EO-1 ALI sensors. Remote Sensing of Environment, 113, 893-903. Cheng, D.-L. & Niklas, K.J. (2007). Above- and Below-ground Biomass Relationships Across 1534 Forested Communities. Annals of Botany, 99, 95-102. Chen, Jing M., (2017), Mapping leaf chlorophyll content: a new way to improve global carbon cycle estimation, International Conference of GeoInformatics? 2017, Buffalo, US. 100 Chen, Jing M. and Cihlar, Josef (1996.), Retrieving leaf area index of boreal conifer forests using Landsat TM images. Remote Sensing of Environment, 55(2):153-162 Chen, Y., Huang, C., Ticehurst, C., Merrin, L., & Thew, P. (2013). An evaluation of MODIS daily and 8-day composite products for floodplain and wetland inundation mapping. Wetlands, 33(5), 823?835. https://doi.org/10.1007/s13157-013-0439-4 Cohen, W.B., Spies, T.A. & Fiorella, M. (1995). Estimating the age and structure of forests in a multi-ownership landscape of western Oregon, U.S.A. International Journal of Remote Sensing, 16, 721-746. Coops, N.C. (2002). Eucalypt forest structure and synthetic aperture radar backscatter: a theoretical analysis. Trees, 16, 28-46. Crevoisier, C., Shevliakova, E., Gloor, M., Wirth, C. & Pacala, S. (2007). Drivers of fire in the boreal forests: Data constrained design of a prognostic model of burned area for use in dynamic global vegetation models. Journal of Geophysical Research-Atmospheres, 112. Corbera, E., & Schroeder, H. (2011). Governing and implementing REDD+. Environmental Science and Policy, 14(2), 89?99. https://doi.org/10.1016/j.envsci.2010.11.002 Crevoisier, C., Shevliakova, E., Gloor, M., Wirth, C., & Pacala, S. (2007). Drivers of fire in the boreal forests: Data constrained design of a prognostic model of burned area for use in dynamic global vegetation models. Journal of Geophysical Research Atmospheres, 112(24). https://doi.org/10.1029/2006JD008372 Crevoisier, C., Sweeney, C., Gloor, M., Sarmiento, J.L. & Tans, P.P. (2010). Regional US carbon sinks from three-dimensional atmospheric CO2 sampling. Proceedings of the National Academy of Sciences, 107, 18348-18353. 101 Dang, An Thi Ngoc;Nandy, Subrata; Srinet, Ritika;Luong, Nguyen Viet;Ghosh,Surajit; Kuma, A.Senthil .Forest aboveground biomass estimation using machine learning regression algorithm in Yok Don National Park.Vietnam, Ecological Informatics. 2019, Volume 50, Pages 24-32. Denise, G., M.A., S., & Stephen, C. (2011). Forest carbon management and carbon trading: A review of Canadian forest options for climate change mitigation. The Forestry Chronicle, 87, 625-635. Diggle, Peter J.; Heagerty, Patrick; Liang, Kung-Yee; Zeger, Scott L. (2002). Analysis of Longitudinal Data (2nd ed.). Oxford University Press. pp. 169?171. ISBN 0-19- 852484-6. Dolan, K. A. (2015). Spatial and Temporal Dynamics of Disturbance Within and Between Forest Regions of the US.Dolan, K., Masek, J.G., Huang, C. & Sun, G. (2009). Regional forest growth rates measured by combining ICESAT GLAS and Landsat data. Journal of Geophysical Research: Biogeosciences, 114, G00E05, doi:10.1029/2008JG000893. Domke, Grant M.; Walters, Brian F.; Nowak, David J.; Smith, James, E.; Ogle, Stephen M.; Coulston, J.W.; Wirth, T.C. 2020. Greenhouse gas emissions and removals from forest land, woodlands, and urban trees in the United States, 1990-2018. Resource Update FS-227. Madison, WI: U.S. Department of Agriculture, Forest Service, Northern Research Station. 5 p. https://doi.org/10.2737/FS-RU-227. Donoghue, D.N.M. & Watt, P.J. (2006). Using LiDAR to compare forest height estimates from IKONOS and Landsat ETM+ data in Sitka spruce plantation forests. International Journal of Remote Sensing, 27, 2161-2175. 102 Du, P.; Samat, A.; Waske, B.; Liu, S.;Li, Z. Random forest and rotation forest for.fully polarized SAR image classification using polarimetric and spatial features.ISPRS J. Photogramm. Remote Sens. 2015, 105, 38?53. Dubayah, R.O., & Drake, J.B. (2000). Lidar remote sensing for forestry. Journal of Forestry, 98, 44-46. Dubayah, R.O., Sheldon, S.L., Clark, D.B., Hofton, M.A., Blair, J.B., Hurtt, G.C. & Chazdon, R.L. (2010). Estimation of tropical forest height and biomass dynamics using lidar remote sensing at La Selva, Costa Rica. Journal of Geophysical Research- Biogeosciences, 115. Dugan, Alexa J, Richard A Birdsey, Sean P Healey, Christopher Woodall, Fangmin Zhang, M Jing, Alexander Hernandez, and James B Mccarter.(2015),?UTILIZING FOREST INVENTORY AND ANALYSIS DATA, REMOTE SENSING, AND ECOSYSTEM MODELS FOR NATIONAL FOREST SYSTEM CARBON ASSESSMENTS.? New Directions in Inventory Techniques & Applications Forest Inventory & Analysis (FIA) Symposium, 2015, 147?52. Eidenshink, J., Schwind, B., Brewer, K., Zhu, Z.-L., Quayle, B. & Howard, S. (2007). A Project For Monitoring Trends In Burn Severity. Fire Ecology, 3, 3-21. Europe, U.N.E.C.f. Forest Products Annual Market Review 2006-2007; United Nations Publications: 2007. Fahey, T.J., Woodbury, P.B., Battles, J.J., Goodale, C.L., Hamburg, S.P., Ollinger, S.V., & Woodall, C.W. (2010). Forest carbon storage: ecology, management, and policy. Frontiers in Ecology and the Environment, 8, 245-252. 103 Fan, S., Gloor, M., Mahlman, J., Pacala, S., Sarmiento, J., Takahashi, T. & Tans, P. (1998). A large terrestrial carbon sink in North America implied by atmospheric and oceanic carbon dioxide data and models. Science, 282, 442-446. Fearnside, P.M. (1997). GREENHOUSE GASES FROM DEFORESTATION IN BRAZILIAN AMAZONIA: NET COMMITTED EMISSIONS. Climatic Change, 35, 321-360. Fitzmaurice, Garrett M.; Laird, Nan M.; Ware, James H. (2004). Applied Longitudinal Analysis. Hoboken: John Wiley & Sons. pp. 326?328. ISBN 0-471-21487-6. Flanagan, S. A. (2016). Plant Migrations Impact on Potential Vegetation and Carbon Redistribution in Northern North America from Climate Change. https://doi.org/10.16953/deusbed.74839 Foody, G.M., Palubinskas, G., Lucas, R.M., Curran, P.J. & Honzak, M. (1996). Identifying terrestrial carbon sinks: Classification of successional stages in regenerating tropical forest from Landsat TM data. Remote Sensing of Environment, 55, 205-216. Franklin, S.E., Hall, R.J., Smith, L. & Gerylo, G.R. (2003). Discrimination of conifer height, age and crown closure classes using Landsat-5 TM imagery in the Canadian Northwest Territories. International Journal of Remote Sensing, 24, 1823-1834. Fransson, F., WaLter F., Ulander, L., (2000). Estimation of forest parameters using CARABAS-II VHF SAR data, IEEE Transactions on Geoscience and Remote Sensing, Volume: 38, Issue: 2. Friedlingstein, P., Jones, M. W., O?Sullivan, M., Andrew, R. M., Hauck, J., Peters, G. P., Peters, W., Pongratz, J., Sitch, S., Le Qu?r?, C., DBakker, O. C. E., Canadell1, J. G., 104 Ciais1, P., Jackson, R. B., Anthoni1, P., Barbero, L., Bastos, A., Bastrikov, V., Becker, M., ? Zaehle, S. (2019). Global carbon budget 2019. Earth System Science Data, 11(4), 1783?1838. https://doi.org/10.5194/essd-11-1783-2019 Gardiner, Joseph C.; Luo, Zhehui; Roman, Lee Anne (2009). "Fixed effects, random effects and GEE: What are the differences?". Statistics in Medicine. 28: 221?239. doi:10.1002/sim.3478. Garestier, F. & Dubois-Fernandez, P.C. (2008). Forest Height Inversion Using High- Resolution P-Band Pol-InSAR Data. IEEE Transactions on Geoscience and Remote Sensing, 46, 3544-3559. Gillespie, A.J.R., (1999.) Rationale for a National Annual Forest Inventory program. J. For. 97, 16?20. Globeland 30((http://www.globeland30.cn/GLC30Download/index.aspx). Goetz, S., Baccini, A., Laporte, N., Johns, T., Walker, W., Kellndorfer, J., Houghton, R. & Sun, M. (2009). Mapping and monitoring carbon stocks with satellite observations: a comparison of methods. Carbon Balance and Management, 4, 2. Goetz, S. J., B. Bond-Lamberty, B. E. Law, J. A. Hicke, C. Huang, R. A. Houghton, S. McNulty, et al. ?Observations and Assessment of Forest Carbon Dynamics Following Disturbance in North America.? Journal of Geophysical Research: Biogeosciences 117, no. 2 (2012): 1?17. doi:10.1029/2011JG001733. Goward, S.N., Huang, C., Masek, J.M., Cohen, W. & Moisen, G. (2007). Role of North America Forest Disturbance and Regrowth In NACP: Integrated Analyses Of Landsat and U.S. Forest Service FIA Data - Phase 2 National Aeronautics and Space Adminstration, pp. 50. 105 Goward, S.N., Huang, C., Masek, J.M., Cohen, W., Moisen, G. & Nemani, R. (2010). US forest disturbance history from Landsat. National Aeronautics and Space Administration, pp. 44 Goward, S.N., Masek, J. & Cohen, W. (2002). Assessing Eastern North America Forest Disturbance and Regrowth: Potential from Passive Optical Remote Sensing Evaluated in the US Mid-Atlantic Region National Aeronautics and Space Adminstration, pp. 32. Goward, S.N., Masek, J., Cohen, W. & Moisen, G. (2005). North American Forest Disturbance and Regrowth since 1972: Empirical Assessment with Field Measurements and Satellite Remotely Sensed Observations. National Aeronautics and Space Adminstration, pp. 78. Goward, S.N., Masek, J.G., Cohen, W., Moisen, G., Collatz, G.J., Healey, S., Houghton, R., Huang, C., Kennedy, R., Law, B., Powell, S., Turner, D. & Wulder, M.A. (2008). Forest disturbance and North American carbon flux. Eos Transactions, 89, 105-116. Gray, Andrew N., and Thomas R. Whittier. ?Carbon Stocks and Changes on Pacific Northwest National Forests and the Role of Disturbance, Management, and Growth.? Forest Ecology and Management 328 (2014): 167?78. Gray, A. N., & Whittier, T. R. (2014). Carbon stocks and changes on Pacific Northwest national forests and the role of disturbance, management, and growth. Forest Ecology and Management, 328, 167?178. https://doi.org/10.1016/j.foreco.2014.05.015 Greene, W.H., 2011. Econometric Analysis, 7th ed., Prentice Hall Griscom, B.W., Adams, J., Ellis, P.W., Houghton, R.A., Lomax, G., Miteva, D.A., Schlesinger, W.H., Shoch, D., Siikam?ki, J.V., Smith, P., Woodbury, P., Zganjar, C., Blackman, A., Campari, J., Conant, R.T., Delgado, C., Elias, P., Gopalakrishna, T., 106 Hamsik, M.R., Herrero, M., Kiesecker, J., Landis, E., Laestadius, L., Leavitt, S.M., Minnemeyer, S., Polasky, S., Potapov, P., Putz, F.E., Sanderman, J., Silvius, M., Wollenberg, E., & Fargione, J. (2017). Natural climate solutions. Proceedings of the National Academy of Sciences, 114, 11645-11650. Grossman, Robert; Seni, Giovanni; Elder, John; Agarwal, Nitin; Liu, Huan. Ensemble Methods in Data Mining: Improving Accuracy Through Combining Predictions. Morgan & Claypool. 2010. Guo, L.B. & Gifford, R.M. (2002). Soil carbon stocks and land use change: a meta analysis. Global Change Biology, 8, 345-360. Han, F.X., Plodinec, M.J., Su, Y., Monts, D.L., & Li, Z. (2007). Terrestrial carbon pools in southeast and south-central United States. Climatic Change, 84, 191-202. Hansen, M.C., Egorov, A., Roy, D.P., Potapov, P., Ju, J.C., Turubanova, S., Kommareddy, I. & Loveland, T.R. (2011). Continuous fields of land cover for the conterminous United States using Landsat data: first results from the Web-Enabled Landsat Data (WELD) project. Remote Sensing Letters, 2, 279-288. Hansen, M.C., Potapov, P.V., Moore, R., Hancher, M., Turubanova, S.A., Tyukavina, A., Thau, D., Stehman, S.V., Goetz, S.J., Loveland, T.R., Kommareddy, A., Egorov, A., Chini, L., Justice, C.O., & Townshend, J.R.G. (2013). High-Resolution Global Maps of 21st-Century Forest Cover Change. Science, 342, 850-853. Harmon, M.E., Ferrell, W.K. & Franklin, J.F. (1990). Effects on Carbon Storage of Conversion of Old-Growth Forests to Young Forests. Science, 247, 699-702. Harris, N. L., Brown, S., Hagen, S. C., Saatchi, S. S., Petrova, S., Salas, W., ? Lotsch, A. (2012). Baseline map of carbon emissions from deforestation in tropical regions. 107 Science, 336(6088), 1573?1576. https://doi.org/10.1126/science.1217962Harris, N. L.; Hagen, S. C.; Saatchi, S. S.; Pearson, T. R.H.; Woodall,C. W. ;Domke, G. M.; Braswell,B. H. et al. Attribution of Net Carbon Change by Disturbance Type across Forest Lands of the Conterminous United States. Carbon Balance and Management.2016, 11, no. 1. doi:10.1186/s13021-016-0066-5. Hayes, D.J., Turner, D.P., Stinson, G., McGuire, A.D., Wei, Y.X., West, T.O., Heath, L.S., Dejong, B., McConkey, B.G., Birdsey, R.A., Kurz, W.A., Jacobson, A.R., Huntzinger, D.N., Pan, Y.D., Mac Post, W. & Cook, R.B. (2012). Reconciling estimates of the contemporary North American carbon balance among terrestrial biosphere models, atmospheric inversions, and a new approach for estimating net ecosystem exchange from inventory-based data. Global Change Biology, 18, 1282-1299. Healey, S.P., Yang, Z., Cohen, W.B., Pierce, D.J., (2006). Application of two regression based methods to estimate the effects of partial harvest on forest structure using Landsat data. Remote Sens. Environ. 101, 115?126. Ho, Tin Kam. Random Decision Forests. Proceedings of the 3rd International Conference on Document Analysis and Recognition.Montreal, QC, 14?16. August 1995. pp. 278? 282. Ho, TK. The Random Subspace Method for Constructing Decision Forests. IEEE Transactions on Pattern Analysis and Machine Intelligence. 1998, 20(8): 832?844. doi:10.1109/34.709601. Homer, C., Dewitz, J., Yang, L., Jin, S., Danielson, P., Xian, G., Coulston, J., Herold, N., Wickham, J., & Megown, K. (2015). Completion of the 2011 National Land Cover 108 Database for the Conterminous United States ? Representing a Decade of Land Cover Change Information. Photogrammetric Engineering & Remote Sensing, 81, 345-354. Houghton, R.A. (2003a). Why are estimates of the terrestrial carbon balance so different? Global Change Biology, 9, 500-509. Houghton, R.A. (2003b). Revised estimates of the annual net flux of carbon to the atmosphere from changes in land use and land management 1850-2000. Tellus Series B-Chemical and Physical Meteorology, 55, 378-390. Houghton, R.A. (2005a). Aboveground Forest Biomass and the Global Carbon Balance. Global Change Biology, 11, 945-958. Houghton, R. A.(2005b), Tropical deforestation as a source of greenhouse gas emissions, Tropical Deforestation and Climate Change, Edited by Paulo Moutinho & Stephan Schwartzman. Houghton, R.A. (2007). Balancing the Global Carbon Budget. Annual Review of Earth and Planetary Sciences, 35, 313-347. Houghton, Richard A. ?Keeping Management Effects Separate from Environmental Effects in Terrestrial Carbon Accounting.? Global Change Biology 19, no. 9 (2013): 2609?12. doi:10.1111/gcb.12233. Houghton, R A, and E A Davidson. Missing Sinks, Feedbacks, and Understanding the Role of Terrestrial Ecosystems in the Global Carbon Balance. Global Biogeochemical ? 12, no. 1 (1998): 25?34. doi:10.1029/97GB02729. Houghton, R.A., Hackler, J.L. & Lawrence, K.T. (1999). The U.S. carbon budget: Contributions from land-use change. Science, 285, 574-578. 109 Houghton, R.A. & Hackler, J.L. (2000a). Changes in terrestrial carbon storage in the United States. II: The role of fire and fire management. Global Ecology & Biogeography, 9, 145-170. Houghton, R.A. & Hackler, J.L. (2000b). Changes in terrestrial carbon storage in the United States. 1: The roles of agriculture and forestry. Global Ecology and Biogeography, 9, 125-144. Houghton, R.A., Hall, F. & Goetz, S.J. (2009). Importance of biomass in the global carbon cycle. Journal of Geophysical Research-Biogeosciences, 114, G00E03. Houghton, R.A., House, J.I., Pongratz, J., van der Werf, G.R., DeFries, R.S., Hansen, M.C., Le Quere, C. & Ramankutty, N. (2012). Carbon emissions from land use and land-cover change. Biogeosciences, 9, 5125-5142. Houghton, R A, J L Hackler, K T Lawrence, R A Houghton, J L Hackler, and K T Lawrence. ?The U.S. Carbon Budget: Contributions from Land-Use Change? 285, no. 5427 (2017): 574?78. Houghton,R A, Nassikas, Alexander A (2017). Global and regional fluxes of carbon from land use and land cover change 1850?2015. AGU:Global Biogeochemical Cycles,456- 472 Houghton, R.A., Skole, D.L., Nobre, C.A., Hackler, J.L., Lawrence, K.T. & Chomentowski, W.H. (2000). Annual fluxes or carbon from deforestation and regrowth in the Brazilian Amazon. Nature, 403, 301-304. Huang, C., Goward, S.N., Masek, J.G., Gao, F., Vermote, E.F., Thomas, N., Schleeweis, K., Kennedy, R.E., Zhu, Z., Eidenshink, J.C. & Townshend, J.R.G. (2009a). 110 Development of time series stacks of Landsat images for reconstructing forest disturbance history. International Journal of Digital Earth, 2, 195-218. Huang, C., Goward, S.N., Schleeweis, K., Thomas, N., Masek, J.G. & Zhu, Z. (2009b). Dynamics of national forests assessed using the Landsat record: case studies in eastern U.S. Remote Sensing of Environment, 113, 1430-1442. Huang, C., Goward, S.N., Masek, J.G., Thomas, N., Zhu, Z. & Vogelmann, J.E. (2010). An automated approach for reconstructing recent forest disturbance history using dense Landsat time series stacks. Remote Sensing of Environment, 114, 183-198. Huang, Chengquan, Ling, Pui-Yu and Zhu, Zhiliang (2015), North Carolina?s forest disturbance and timber production assessed using time series Landsat observations, International Journal of Digital Earth. 1-23. http://dx.doi.org/10.1080/17538947.2015.1034200 Huang, C., Schleeweis, K., Thomas, N., & Goward, S.N. (2011). Forest dynamics within and around the Olympic National Park assessed using time series Landsat observations. In Y. Wang (Ed.), Remote Sensing of Protected Lands (pp. 71-93). London: Taylor & Francis. Huang, C. & Townshend, J.R.G. (2003). A stepwise regression tree for nonlinear approximation: applications to estimating subpixel land cover. International Journal of Remote Sensing, 24, 75-90. Huang, W. (2015). Monitoring and Assessing Forest Biomass from Disturbance and Recovery with LiDAR and SAR.Huang, X., & Jensen, J. R. (1997). A Machine- Learning Approach to Automated Knowledge-Base Building for Remote Sensing 111 Image Analysis with GIS Data. Photogrammetric Engineering & Remote Sensing, 63(10), 1185?1194. Huntzinger, D.N., Post, W.M., Wei, Y., Michalak, A.M., West, T.O., Jacobson, A.R., Baker, I.T., Chen, J.M., Davis, K.J., Hayes, D.J., Hoffman, F.M., Jain, A.K., Liu, S., McGuire, A.D., Neilson, R.P., Potter, C., Poulter, B., Price, D., Raczka, B.M., Tian, H.Q., Thornton, P., Tomelleri, E., Viovy, N., Xiao, J., Yuan, W., Zeng, N., Zhao, M. & Cook, R. (2012). North American Carbon Program (NACP) regional interim synthesis: Terrestrial biospheric model intercomparison. Ecological Modelling, 232, 144-157. Hurtt, G.C., Pacala, S.W., Moorcroft, P.R., Caspersen, J., Shevliakova, E., Houghton, R.A. & Moore, B. (2002). Projecting the future of the U.S. carbon sink. PNAS, 99, 1389- 1394. Idol, T.W., Figler, R.A., Pope, P.E. & Ponder Jr, F. (2001). Characterization of coarse woody debris across a 100 year chronosequence of upland oak???hickory forests. Forest Ecology and Management, 149, 153-161. IPCC, (2014), AR5 Synthesis Report: Climate Change 2014 ? IPCC, https://www.ipcc.ch/report/ar5/syr/ Jenkins, J.C., Chomnacky, D.C., Heath, L.S. & Birdsey, R.A. (2003). National scale biomass estimators for United States tree species. Forest Science, 49, 12-35. Jensen, J.R., Qiu, F. & Ji, M. (1999). Predictive modelling of coniferous forest age using statistical and artificial neural network approaches applied to remote sensor data. International Journal of Remote Sensing, 20, 2805-2822. 112 Johnson, E.W. & Wittwer, D. (2008). Aerial detection surveys in the United States. Australian Forestry, 71, 212-215. Johnson, T.G. (1996). United States timber industry-an assessment of timber product output and use, 1996. Gen. Tech. Rep. SRS-45. Asheville, NC: US Department of Agriculture, Forest Service, Southern Research Station. 145 p. 2001, 45. Kasischke, E.S., Amiro, B.D., Barger, N.N., French, N.H.F., Goetz, S.J., Grosse, G., Harmon, M.E., Hicke, J.A., Liu, S. & Masek, J.G. (2013). Impacts of disturbance on the terrestrial carbon budget of North America. Journal of Geophysical Research: Biogeosciences, 118, 303-316. Kauffman, J. S. (2016). Spatiotemporal Informatics for Sustainable Forest Production Utilizing Forest Inventory and Remotely Sensed Data. Retrieved from https://vtechworks.lib.vt.edu/handle/10919/74974 Kellndorfer, J.M., Walker, W.S., LaPoint, E., Kirsch, K., Bishop, J. & Fiske, G. (2010). Statistical fusion of lidar, InSAR, and optical remote sensing data for forest stand height characterization: A regional-scale method based on LVIS, SRTM, Landsat ETM plus, and ancillary data sets. Journal of Geophysical Research-Biogeosciences, 115, doi:10.1029/2009JG000997. Kellndorfer, J., Walker, W., LaPoint, E., Cormier, T., Bishop, J., Fiske, G. & Kirsch, K. (2013). Vegetation height, biomass, and carbon stock for the conterminous United States: A high-resolution dataset from Landsat ETM+, SRTM-InSAR, National Land Cover Database, and Forest Inventory and Analysis data fusion. 113 Kennedy, R.E., Cohen, W.B. & Schroeder, T.A. (2007). Trajectory-based change detection for automated characterization of forest disturbance dynamics. Remote Sensing of Environment, 110, 370-386. Kennedy, R.E., Yang, Z., Braaten, J., Copass, C., Antonova, N., Jordan, C., & Nelson, P. (2015). Attribution of disturbance change agent from Landsat time-series in support of habitat monitoring in the Puget Sound region, USA. Remote Sensing of Environment, 166, 271-285. Kennedy, R.E., Yang, Z.G., & Cohen, W.B. (2010). Detecting trends in forest disturbance and recovery using yearly Landsat time series: 1. LandTrendr - Temporal segmentation algorithms. Remote Sensing of Environment, 114, 2897-2910. Kennedy, R.E., Yang, Z., Cohen, W.B., Pfaff, E., Braaten, J., & Nelson, P. (2012). Spatial and temporal patterns of forest disturbance and regrowth within the area of the Northwest Forest Plan. Remote Sensing of Environment, 122, 117-133. Kimes, D.S., Nelson, R.F., Skole, D.L. & Salas, W.A. (1998). Accuracies in mapping secondary tropical forest age from sequential satellite imagery. Remote Sensing of Environment, 65, 112-120. Laird, Nan M.; Ware, James H. (1982). "Random-Effects Models for Longitudinal Data". Biometrics. 38 (4): 963?974. JSTOR 2529876. Lamb, R.L., Hurtt, G.C., Boudreau, T.J., Campbell, E., Sep?lveda Carlo, E.A., Chu, H.- H., de Mooy, J., Dubayah, R.O., Gonsalves, D., Guy, M., Hultman, N.E., Lehman, S., Leon, B., Lister, A.J., Lynch, C., Ma, L., Martin, C., Robbins, N., Rudee, A., Silva, C.E., Skoglund, C., & Tang, H. (2021). Context and future directions for integrating 114 forest carbon into sub-national climate mitigation planning in the RGGI region of the U.S. Environmental Research Letters, 16, 063001. Law, B.E., Van Tuyl, S., Thornton, P.E., Irvine, J. & Anthoni, P.M. (2001). Carbon storage and fluxes in ponderosa pine forests at different developmental stages. Global Change Biology, 7, 755-777. Lefsky, M.A., Acker, S.A., Gower, S.T., Cohen, W.B., Harding, D.J., & Parker, G.G. (2002). Lidar remote sensing of above-ground biomass in three biomes. Global Ecology and Biogeography, 11, 393-399. Lefsky, M.A., Cohen, W.B., Acker, S.A., Parker, G.G., Spies, T.A., & Harding, D. (1999). Lidar remote sensing of the canopy structure and biophysical properties of douglas-fir western hemlock forests. Remote Sensing of Environment, 70, 339-361. Lefsky, M.A., Carabajal, C.C., Del Bom Espirito-Santo, F., Hunter, M.O., de Oliveira Jr, R., Harding, D.J., Keller, M. & Cohen, W.B. (2005). Estimates of forest canopy height and aboveground biomass using ICESat. Geophysical Research Letters, 32, 1-4 Article Number L22S02. Le Qu?r?, C., Raupach, M.R., Canadell, J.G., Marland, G., Bopp, L., Ciais, P., Conway, T.J., Doney, S.C., Feely, R.A., Foster, P., Friedlingstein, P., Gurney, K., Houghton, R.A., House, J.I., Huntingford, C., Levy, P.E., Lomas, M.R., Majkut, J., Metzl, N., Ometto, J.P., Peters, G.P., Prentice, I.C., Randerson, J.T., Running, S.W., Sarmiento, J.L., Schuster, U., Sitch, S., Takahashi, T., Viovy, N., van der Werf, G.R. & Woodward, F.I. (2009). Trends in the sources and sinks of carbon dioxide. Nature Geoscience, 2, 831-836. 115 Le Qu?r?, C., Andrew, R. M., Canadell, J. G., Sitch, S., Ivar Korsbakken, J., Peters, G. P., Manning, A. C., Boden, T. A., Tans, P. P., Houghton, R. A., Keeling, R. F., Alin, S., Andrews, O. D., Anthoni, P., Barbero, L., Bopp, L., Chevallier, F., Chini, L. P., Ciais, P., ? Zaehle, S. (2016). Global Carbon Budget 2016. Earth System Science Data, 8(2), 605?649. https://doi.org/10.5194/essd-8-605-2016 Le Qu?r?, C., Andrew, R. M., Friedlingstein, P., Sitch, S., Pongratz, J., Manning, A. C., Ivar Korsbakken, J., Peters, G. P., Canadell, J. G., Jackson, R. B., Boden, T. A., Tans, P. P., Andrews, O. D., Arora, V. K., Bakker, D. C. E., Barbero, L., Becker, M., Betts, R. A., Bopp, L., ? Zhu, D. (2018). Global Carbon Budget 2017. Earth System Science Data, 10(1), 405?448. https://doi.org/10.5194/essd-10-405-2018 Lesiv, Myroslava; Shvidenko, Anatoly; Schepaschenko, Dmitry; See,Linda;Fritz, Steffen, A spatial assessment of the forest carbon budget for Ukraine, Mitigation and Adaptation Strategies for Global Change. 2019, Volume 24, Issue 6, pp 985?1006 Li, A., Huang, C., Sun, G., Shi, H., Toney, C., Zhu, Z., Rollins, M.G., Goward, S.N. & Masek, J.G. (2011). Modeling the growth of young forests regenerating from recent disturbances in Mississippi using Landsat time series observations and ICESat/GLAS lidar data. Remote Sensing of Environment, 115, 1837?1849. Lippke, B.; Oneil, E.; Harrison, R.; Skog, K.; Gustavsson, L.; Sathre, R. (2011). Life cycle impacts of forest management and wood utilization on carbon mitigation: knowns and unknowns. Carbon Management, 2, 303-333. Ling, P.-Y. (2016). Integrating Forest Ecological Processes, Management, and Carbon Payments To Assess Ecosystem Service Trade-Offs. 116 Ling, P.-Y., Baiocchi, G., & Huang, C. (2016). Estimating annual influx of carbon to harvested wood products linked to forest management activities using remote sensing. Climatic Change, 134, 45-58. Liu, J., Liu, S., Loveland, T. R., & Tieszen, L. L. (2008). Integrating remotely sensed land cover observations and a biogeochemical model for estimating forest ecosystem carbon dynamics. Ecological Modelling, 219(3?4), 361?372. https://doi.org/10.1016/j.ecolmodel.2008.04.019 Liu, S., Loveland, T. R., & Kurtz, R. M. (2004). Contemporary carbon dynamics in terrestrial ecosystems in the Southeastern Plains of the United States. Environmental Management, 33(SUPPL. 1). https://doi.org/10.1007/s00267-003-9152-z Liu, Shanshan;Wei, Xinliang; Li ,Dengqiu and Lu,Dengsheng.Examining Forest Disturbance and Recovery in the Subtropical Forest Region of Zhejiang Province Using Landsat Time-Series Data. Remote Sensing. 2017, 9, 479; doi:10.3390/rs9050479 Liu, W., Song, C., Schroeder, T. A., & Cohen, W. B. (2008). Predicting forest successional stages using multitemporal Landsat imagery with forest inventory and analysis data. International Journal of Remote Sensing, 29(13), 3855?3872. https://doi.org/10.1080/01431160701840166 Liu, Yanan; Gong,Weishu;Hu. Xiangyun et al., Forest Type Identification with Random Forest Using Sentinel-1A, Sentinel-2A, Multi-Temporal Landsat-8 and DEM Data. Remote Sens. 2018, 10, 946; doi:10.3390/rs10060946 Liu, Y.; Gong, W.; et al., 2019, Estimation of the forest stand mean height and aboveground biomass in Northeast China using SAR Sentinel-1B, multispectral 117 Sentinel-2A, and DEM imagery, ISPRS Journal of Photogrammetry and Remote Sensing, Volume 151, May 2019, Pages 277-289 Liu. Yanan; Hu, Xiangyun;Wu, Hao; Zhang, Anqi et al., Spatiotemporal Analysis of Carbon Emissions and Carbon Storage Using National Geography Census Data in Wuhan, China. ISPRS Int. J. Geo-Inf. 2019, 8, 7; doi:10.3390/ijgi8010007. Loveland, T.R., & Dwyer, J.L. (2012). Landsat: Building a strong future. Remote Sensing of Environment, 122, 22-29. Makela, H., & Pekkarinen, A. (2001). Estimation of timber volume at the sample plot level by means of image segmentation and Landsat TM imagery. Remote Sensing of Environment, 77, 66-75. Marland, G., Pielke, R.A., Apps, M., Avissar, R., Betts, R.A., Davis, K.J., Frumhoff, P.C., Jackson, S.T., Joyce, L.A., Kauppi, P., Katzenberger, J., MacDicken, K.G., Neilson, R.P., Niles, J.O., Niyogi, D.d.S., Norby, R.J., Pena, N., Sampson, N., & Xue, Y. (2003). The climatic impacts of land surface change and carbon management, and the implications for climate-change mitigation policy. Climate Policy, 3, 149-157. Masek, J., Goward, S., Kennedy, R., Cohen, W., Moisen, G., Schleeweis, K. & Huang, C. (2013). United States Forest Disturbance Trends Observed Using Landsat Time Series. Ecosystems, DOI: 10.1007/s10021-013-9669-9, 1-18. Masek, J, Vermote ,E., et al.,2006, A Landsat Surface Reflectance Dataset for North America, 1990?2000, IEEE Geoscience and Remote Sensing Letters, Volume: 3 , Issue: 1 . McKinley, D.C., Ryan, M.G., Birdsey, R.A., Giardina, C.P., Harmon, M.E., Heath, L.S., Houghton, R.A., Jackson, R.B., Morrison, J.F., Murray, B.C., Pataki, D.E. & Skog, 118 K.E. (2011). A synthesis of current knowledge on forests and carbon storage in the United States. Ecological Applications, 21, 1902-1924. McRoberts, Ronald E., and Erkki O. Tomppo. ?Remote Sensing Support for National Forest Inventories.? Remote Sensing of Environment 110, no. 4 (2007): 412?19. doi: 10.1016/j.rse.2006.09.034. Meddens, A.J.H., Hicke, J.A. & Ferguson, C.A. (2012). Spatiotemporal patterns of observed bark beetle-caused tree mortality in British Columbia and the western United States. Ecological Applications, 22, 1876-1891. Michalak, A.M., Jackson, R.B., Marland, G., Sabine, C.L., Anderson, R.F., Bronk, D., Davis, K.J., DeFries, R.S., Denning, A.S., Dilling, L., Jacobson, A., Lohrenz, S., McGuire, A.D., McKinley, G.A., Miller, C., III, B.M., Ojima, D.S., O?Neill, B., Randerson, J.T., Running, S.W., Sohngen, B., Tans, P.P., Thornton, P.E., Wofsy, S.C. & Zeng, N. (2011). A U.S. Carbon Cycle Science Plan: University Corporation for Atmospheric Research. Mishra, N., Helder, D., Barsi, J., & Markham, B. (2016). Continuous calibration improvement in solar reflective bands: Landsat 5 through Landsat 8. Remote Sensing of Environment, 185, 7-15. MRLC, 2019, https://www.mrlc.gov/ Mokany, K., Raison, R.J. & Prokushkin, A.S. (2006). Critical analysis of root : shoot ratios in terrestrial biomes. Global Change Biology, 12, 84-96. Montesano, P. M. (2015). The Uncertainty of Spaceborne Observation of Vegetation Structure in the Taiga-Tundra Ecotone: A Case Study in Northern Siberia. ProQuest Dissertations and Theses, 147. Retrieved from 119 https://search.proquest.com/docview/1733293169?accountid=27233%0Ahttp://techni on-primo.hosted.exlibrisgroup.com/openurl/972TEC/972TEC_SP? Mougin, E., Proisy, C., Marty, G., Fromard, F., Puig, H., Betoulle, J.L. & Rudant, J.P. (1999). Multifrequency and Multipolarization Radar Backscattering from Mangrove Forests. IEEE Transactions on Geoscience and Remote Sensing, 37, 94-102. Murty, D., Kirschbaum, M.U.F., McMurtrie, R.E. & McGilvray, H. (2002). Does conversion of forest to agricultural land change soil carbon and nitrogen? a review of the literature. Global Change Biology, 8, 105-123. Naudts, K., Chen, Y., McGrath, M.J., Ryder, J., Valade, A., Otto, J., & Luyssaert, S. (2016). Europe's forest management did not mitigate climate warming. Science, 351, 597-600. NBCD,2000, NACP Aboveground Biomass and Carbon Baseline Data, V.2 (NBCD 2000), U.S.A., 2000, https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1161 Neeff, T. ?Dutra,L. ?Santos, J. ? Freitas?C., Araujo, L., 2003, Tropical forest biomass measurement by backscatter and DEM information as derived from airborne SAR, IEEE International Geoscience & Remote Sensing Symposium, Volume: 4, DOI:10.1109/IGARSS.2003.1294512 Neeff, T., de Alencastro Graca, P.M., Dutra, L.V., & da Costa Freitas, C. (2005). Carbon budget estimation in Central Amazonia: Successional forest modeling from remote sensing data. Remote Sensing of Environment, 94, 508-522. Nelson, R., Boudreau, J., Gregoire, T.G., Margolis, H., Naesset, E., Gobakken, T. & Stahl, G. (2009). Estimating Quebec provincial forest resources using ICESat/GLAS. 120 Canadian Journal of Forest Research-Revue Canadienne De Recherche Forestiere, 39, 862-881. Nelson, R., Margolis, H., Montesano, P., Sun, G., Cook, B., Corp, L., Andersen, H.-E., deJong, B., Pellat, F.P., Fickel, T., Kauffman, J., & Prisley, S. (2017). Lidar-based estimates of aboveground biomass in the continental US and Mexico using ground, airborne, and satellite observations. Remote Sensing of Environment, 188, 127-140. Nilson, T. & Peterson, U. (1994). Age dependence of forest reflectance: Analysis of main driving factors. Remote Sensing of Environment, 48, 319-331. Noormets, A., Epron, D., Domec, J.C., McNulty, S.G., Fox, T., Sun, G., & King, J.S. (2015). Effects of forest management on productivity and carbon sequestration: A review and hypothesis. Forest Ecology and Management, 355, 124-140. Novick, K. A., Oishi, A. C., Ward, E. J., Siqueira, M. B. S., Juang, J. Y., & Stoy, P. C. (2015). On the difference in the net ecosystem exchange of CO2 between deciduous and evergreen forests in the southeastern United States. Global Change Biology, 21(2), 827?842. https://doi.org/10.1111/gcb.12723 Nunery, J.S., & Keeton, W.S. (2010). Forest carbon storage in the northeastern United States: Net effects of harvesting frequency, post-harvest retention, and wood products. Forest Ecology and Management, 259, 1363-1375. O?Connell, Barbara M.; Conkling, Barbara L.; Wilson, Andrea M.; Burrill, Elizabeth A.; Turner, Jeffery A.; Pugh, Scott A.; Christiansen, Glenn; Ridley, Ted; Menlove, James 2016. The Forest Inventory and Analysis Database: Database description and user guide version 6.1.1 for Phase 2. U.S. Department of Agriculture, Forest Service. 870 121 p. [Online]. Available at web address: http://www.fia.fs.fed.us/library/database? documentation/. Oswalt, Sonja N.; Smith, W. Brad; Miles, Patrick D.; Pugh, Scott A. coords.(2017). Forest Resources of the United States, 2017: a technical document supporting the Forest Service 2020 RPA Assessment. Gen. Tech. Rep. WO-97. Washington, DC: U.S. Department of Agriculture, Forest Service, Washington Office.2019, 223 p. https://doi.org/10.2737/WO-GTR-97. Pacala, S.W., Hurtt, G.C., Baker, D., Peylin, P., Houghton, R.A., Birdsey, R.A., Heath, L., Sundquist, E.T., Stallard, R.F., Ciais, P., Moorcroft, P., Caspersen, J.P., Shevliakova, E., Moore, B., Kohlmaier, G., Holland, E., Gloor, M., Harmon, M.E., Fan, S.M., Sarmiento, J.L., Goodale, C.L., Schimel, D. & Field, C.B. (2001). Consistent Land- and Atmosphere-Based U.S. Carbon Sink Estimates. Science, 292, 2316-2320. Pal, M., Random forest classifier for remote sensing classification. Int. J. Remote Sens. 2005, 26, 217?222. Pan, Y., Birdsey, R.A., Fang, J., Houghton, R., Kauppi, P.E., Kurz, W.A., Phillips, O.L., Shvidenko, A., Lewis, S.L., Canadell, J.G., Ciais, P., Jackson, R.B., Pacala, S., McGuire, A.D., Piao, S., Rautiainen, A., Sitch, S. & Hayes, D. (2011a). A Large and Persistent Carbon Sink in the World's Forests. Science, 333, 988-993. Pan, Y., Chen, J.M., Birdsey, R., McCullough, K., He, L. & Deng, F. (2011b). Age structure and disturbance legacy of North American forests. Biogeosciences, 8, 715- 732. 122 Patenaude G.L?Briggs B.D.J?R Milne?Rowland C.S?Dawson T.P?Pryor S.N, (2003). The carbon pool in a British Semi-Natural Woodland, Forestry 76(1) Patenaude, G., Milne, R., & Dawson, T.P. (2005). Synthesis of remote sensing approaches for forest carbon estimation: reporting to the Kyoto Protocol. Environmental Science & Policy, 8, 161-178. Pearse, Grant; Morgenrotha,Justin; Watt, Michael; Dash, Jonathan;2017, Optimising prediction of forest leaf area index from discrete airborne lidar, Remote Sensing of Environment, Volume 200, October 2017, Pages 220-239. Petranka, J.W., Eldridge, M.E., Haley, K.E., (1993). Effects of timber harvesting on southern Appalachian salamanders. Conserv. Biol. 7, 363?377. Pflugmacher, D., Cohen, W.B. & E. Kennedy, R. (2012). Using Landsat-derived disturbance history (1972-2010) to predict current forest structure. Remote Sensing of Environment, 122, 146-165. Pflugmacher, Dirk, Warren B. Cohen, Robert E. Kennedy, and Zhiqiang Yang. ?Using Landsat-Derived Disturbance and Recovery History and Lidar to Map Forest Biomass Dynamics.? Remote Sensing of Environment 151 (2014): 124?37. doi: 10.1016/j.rse.2013.05.033. Ploton, P.; Barbiera, N.; et al, 2017, Toward a general tropical forest biomass prediction model from very high resolution optical satellite images, Remote Sensing of Environment, Volume 200, October 2017, Pages 140-153 Post, W.M. & Kwon, K.C. (2000). Soil carbon sequestration and land-use change: processes and potential. Global Change Biology, 6, 317-327. 123 Powell, S.L., Cohen, W.B., Healey, S.P., Kennedy, R.E., Moisen, G.G., Pierce, K.B. & Ohmann, J.L. (2010). Quantification of live aboveground forest biomass dynamics with Landsat time-series and field inventory data: A comparison of empirical modeling approaches. Remote Sensing of Environment, 114, 1053-1068. Pretzsch, H. (2001). Models for Pure and Mixed Forests. In: J. Evans (Ed.), The Forest Handbook (pp. 210-228): Vol 1, London: Blackwell Science Ltd. Pulliainen, J. ?Engdahl, M. ?Hallikainen, M. (2003). Feasibility of multi-temporal interferometric SAR data for stand-level estimation of boreal forest stem volume, Remote Sensing of Environment 85(4):397-409. Ramankutty, N., Gibbs, H.K., Achard, F., Defries, R., Foley, J.A. & Houghton, R.A. (2007). Challenges to estimating carbon emissions from tropical deforestation. Global Change Biology, 13, 51-66. Ramsey, F., Schafer, D., 2002. The Statistical Sleuth: A Course in Methods of Data Analysis, 2nd ed. Duxbury Press Rauste, Yrjo,( 2005). Multi-temporal JERS SAR data in boreal forest biomass mapping, Remote Sensing of Environment, 97, 263 ? 275. Reybold, W.U. & TeSelle, G.W. (1989). Soil geographic data bases. Journal of Soil and Water Conservation, 44, 28-29. Richards, Gary,2020, https://www.globalcarbonproject.org/ global/ presentations/ 2_Terrestrial/Richards.pdf Richardson, Heather J.; Hill, David J.; Denesiuk, Dan R. & Fraser, Lauchlan H. A comparison of geographic datasets and field measurements to model soil carbon using 124 random forests and stepwise regressions (British Columbia, Canada) ?Journal of GIScience & Remote Sensing.2017, Volume 54, Issue 4, Pages 573-591 Rignot, E., Way, J.B., McDonald, K., Viereck, L., Williams, C., Adams, P., Payne, C., Wood, W., & Shi, J. (1994). Monitoring of environmental conditions in Taiga forests using ERS-1 SAR. Remote Sensing of Environment, 49, 145-154. Rodriguez-Galiano, V.F., Ghimire, B., Rogan, J., Chica-Olmo, M., Rigol-Sanchez,J.P. (2012). An assessment of the effectiveness of a random forest classifier for land cover classification. ISPRS J. Photogramm. Remote Sens., 67, 93?104. Rollins, M.G. (2009). LANDFIRE: a nationally consistent vegetation, wildland fire, and fuel assessment. International Journal of Wildland Fire, 18, 235-249. Ruddell, S., Sampson, R., Smith, M., Giffen, R., Cathcart, J., Hagan, J., Sosland, D., Godbee, J., Heissenbuttel, J., Lovett, S., Helms, J., Price, W., & Simpson, R. (2007). The role for sustainably managed forests in climate change mitigation. Journal of Forestry, 105(6), 314?319. https://doi.org/10.1093/jof/105.6.314 Ruefenacht, B., Finco, M.V., Nelson, M.D., Czaplewski, R., Helmer, E.H., J.A. Blackard, Holden, G.R., Lister, A.J., Salajanu, D., Weyermann, D. & Winterberger., K. (2008). Conterminous U.S. and Alaska Forest Type Mapping Using Forest Inventory and Analysis Data. Photogrammetric Engineering & Remote Sensing, 74, 1379-1388. Running, S.W. (2008). Climate change - Ecosystem disturbance, carbon, and climate. Science, 321, 652-653. Ryan, V. REVISION CARDS - SOFTWOODS. Availabe online: http://www.technologystudent.com/joints/softwoods1.html (accessed on Accessed 2019.12). 125 Saatchi, S.S., Harris, N.L., Brown, S., Lefsky, M., Mitchard, E.T.A., Salas, W., Zutta, B.R., Buermann, W., Lewis, S.L., Hagen, S., Petrova, S., White, L., Silman, M. & Morel, A. (2011). Benchmark map of forest carbon stocks in tropical regions across three continents. Proceedings of the National Academy of Sciences of the United States of America, 108, 9899-9904. Santos, J.R., Freitas, C.C., Araujo, L.S., Dutra, L.V., Mura, J.C., Gama, F.F., Soler, L.S., & Sant'Anna, S.J.S. (2003). Airborne P-band SAR applied to the aboveground biomass studies in the Brazilian tropical rainforest. Remote Sensing of Environment, 87, 482- 493. Sarmiento, J.L., C.Wofsy, S., Shea, E., Denning, A.S., Easterling, W., Field, C., Fung, I., Keeling, R., McCarthy, J., Pacala, S., Post, W.M., Schimel, D., Sundquist, E., Tans, P., Weiss, R. & Yoder, J. (1999). A U.S. Carbon Cycle Science Plan, Washington, DC: U.S.Global Change Research Program. Shah, Syed Haleem;Angel,Yoseline;Houborg, Rasmus;Ali, Shawkat and McCabe, Matthew F., A Random Forest Machine Learning Approach for the Retrieval of Leaf Chlorophyll Content in Wheat, Remote Sensing. 2019, 11(8), 920; https://doi.org/10.3390/rs11080920 Shao, Y., & Lunetta, R. S. (2012). Comparison of support vector machine, neural network, and CART algorithms for the land-cover classification using limited training data points. ISPRS Journal of Photogrammetry and Remote Sensing, 70, 78?87. https://doi.org/10.1016/j.isprsjprs.2012.04.001 Schimel, D., Melillo, J., Tian, H.Q., McGuire, A.D., Kicklighter, D., Kittel, T., Rosenbloom, N., Running, S., Thornton, P., Ojima, D., Parton, W., Kelly, R., Sykes, 126 M., Neilson, R. & Rizzo, B. (2000). Contribution of increasing CO2 and climate to carbon storage by ecosystems in the United States. Science, 287, 2004-2006. Schleeweis, K., Goward, S.N., Huang, C., Masek, J.G., Moisen, G., Kennedy, R.E. & Thomas, N.E. (2013). Regional dynamics of forest canopy change and underlying causal processes in the contiguous U.S. Journal of Geophysical Research: Biogeosciences, 118, 1-19. Schleeweis, K.G., Moisen, G.G., Schroeder, T.A., Toney, C., Freeman, E.A., Goward, S.N., Huang, C., & Dungan, J.L. (2020). US National Maps Attributing Forest Change: 1986?2010. Forests, 11, 653. Schofield, N.J., Ruprecht, J.K., (1989). Regional-analysis of stream salinization in southwestWestern Australia. J. Hydrol. 112, 19. Schroder, Jonathan; Carpena, Rafael Mu?oz; Dukes,,Michael; Li , Yuncong, 2005. The Development and Testing of an Automated Drip Fertigation System for Sustainable Agriculture, [American Society of Agricultural and Biological Engineers 2005 Tampa, FL July 17-20, 2005] Schroeder, T. A., Healey, S. P., Moisen, G. G., Frescino, T. S., Cohen, W. B., Huang, C., Kennedy, R. E., & Yang, Z. (2014). Improving estimates of forest disturbance by combining observations from landsat time series with U.S. Forest Service Forest Inventory and Analysis data. Remote Sensing of Environment, 154(1), 61?73. https://doi.org/10.1016/j.rse.2014.08.005 Schroder, R., Puls, J., Hajnsek, I., Jochim, F., Neff, T., Kono, J., Renato Paradella, W., Marcos Quintino da Silva, M., de Morisson Valeriano, D., & Pereira Farias Costa, M. 127 (2005). MAPSAR: a small L-band SAR mission for land observation. Acta Astronautica, 56, 35-43. Schroeder, T.A., Wulder, M.A., Healey, S.P., Moisen, G.G., (2011). Mapping wildfire and clear cut harvest disturbances in boreal forests with Landsat time series data. Remote Sens. Environ. 115, 1421?1433. Schulze, E.D., Valentini, R., & Sanz, M.J. (2002). The long way from Kyoto to Marrakesh: Implications of the Kyoto Protocol negotiations for global ecology. Global Change Biology, 8, 505-518. Sexton, J.O., Song, X.-P., Feng, M., Noojipady, P., Anand, A., Huang, C., Kim, D.-H., Collins, K.M., Channan, S., DiMiceli, C. & Townshend, J.R. (2013). Global, 30-m resolution continuous fields of tree cover: Landsat-based rescaling of MODIS continuous fields and lidar-based estimates of error. International Journal of Digital Earth, DOI:10.1080/17538947.2013.786146. Smith JE, Heath LS, Skog KE, Birdsey RA (2006). Methods for calculating forest ecosystem and harvested carbon with standard estimates for forest types of the United States. Gen. Tech. Rep. NE-343. Newtown Square, PA, 216 Smith, B.W., Miles, P.D., Perry, C.H. & Pugh, S.A. (2007). Forest resources of the United States: a technical document supporting the Forest Service 2010 RPA assessment. General Techincal Report, Washington, D.C.: U.S. Department of Agriculture, Forest Service, Washington Office. Smith, W.B. (2002). Forest inventory and analysis: a national inventory and monitoring program. Environmental Pollution, 116, S233-S242. 128 Song, X. (2015). Improved Quantification of Forest Cover Change and Implications for the Carbon Cycle. Steininger, M. (2004). Net Carbon Fluxes from Forest Clearance and Regrowth in the Amazon. Ecological Applications, 14(4). Retrieved from https://www.jstor.org/stable/4493648 Stockmann, K.D.; Anderson, N.M.; Skog, K.E.; Healey, S.P.; Loeffler, D.R.; Jones, G.; Morrison, J.F. Estimates of carbon stored in harvested wood products from the United States forest service northern region, 1906-2010. Carbon balance and management 2012, 7, 1 Sun, G., Ranson, K.J., Khairuk, V.I. & Kovacs, K. (2003). Validation of surface height from shuttle radar topography mission using shuttle laser altimeter. Remote Sensing of Environment, 88, 401-411. Tang, H. (2015). LiDAR Remote Sensing of Vertical Foliage Profile and Leaf Area Index. Tao, X., Huang, C., Zhao, F., Schleeweis, K., Masek, J., & Liang, S. (2019). Mapping forest disturbance intensity in North and South Carolina using annual Landsat observations and field inventory data. Remote Sensing of Environment, 221, 351-362. Thomas, N., Huang, C., Goward, S.N., Powell, S., Rishmawi, K., Schleeweis, K. & Hinds, A. (2011). Validation of North American Forest Disturbance dynamics derived from Landsat time series stacks. Remote Sensing of Environment, 115, 19-32. Thornton, P.E., Running, S.W. & White, M.A. (1997). Generating surfaces of daily meteorological variables over large regions of complex terrain. Journal of Hydrology, 190, 214-251. 129 Tian, H., Melillo, J.M., Kicklighter, D.W., McGuire, A.D. & Helfrich, J. (1999). The sensitivity of terrestrial carbon storage to historical climate variability and atmospheric CO2 in the United States. Tellus Series B-Chemical and Physical Meteorology, 51, 414-452. Tian, S.; Zhang, X.; Tian, J.; Sun, Q. Random forest classification of wetland land covers from multi-sensor data in the arid region of Xinjiang, China. Remote Sens. 2016, 8, 954. Tong, X., Brandt, M., Yue, Y., Ciais, P., Rudbeck Jepsen, M., Penuelas, J., ? Fensholt, R. (2020). Forest management in southern China generates short term extensive carbon sequestration. Nature Communications, 11(1), 1?10. https://doi.org/10.1038/s41467- 019-13798-8 Trotter, C.M., Dymond, J.R., & Goulding, C.J. (1997). Estimation of timber volume in a coniferous plantation forest using Landsat TM. International Journal of Remote Sensing, 18, 2209-2223. Turner, D.P., Koerper, G.J., Harmon, M.E. & Lee, J.J. (1995a). A Carbon Budget for Forests of the Conterminous United-States. Ecological Applications, 5, 421-436. Turner, D.P., Koerper, G.J., Harmon, M.E. & Lee, J.J. (1995b). Carbon Sequestration by Forests of the United-States - Current Status and Projections to the Year 2040. Tellus Series B-Chemical and Physical Meteorology, 47, 232-239. Turner, D. P., Ritts, W. D., Kennedy, R. E., Gray, A. N., & Yang, Z. (2016). Regional carbon cycle responses to 25 years of variation in climate and disturbance in the US Pacific Northwest. Regional Environmental Change, 16(8), 2345?2355. https://doi.org/10.1007/s10113-016-0956-9 130 Van Deusen, P. (2010). Carbon sequestration potential of forest land: Management for products and bioenergy versus preservation. Biomass and Bioenergy, 34, 1687-1694. Verrelst, J., Mu?oz, J., Alonso, L., Delegido, J., Rivera, J. P., Camps-Valls, G., & Moreno, J. (2012). Machine learning regression algorithms for biophysical parameter retrieval: Opportunities for Sentinel-2 and -3. Remote Sensing of Environment, 118, 127?139. https://doi.org/10.1016/j.rse.2011.11.002 von Gadow, K. & Hui, G. (1999). Modelling Forest Development (pp. 213). Forestry Science, 57, Dordrecht: Kluwer Academic Publishers. Wang?Li?ai; Zhou?Xudong; Zhu? Xinkai ; Dong?Zhaodi ; Guo? Wenshan ; Estimation of biomass in wheat using random forest regression algorithm and remote sensing data. The Crop Journal. 2016, Volume 4, Issue 3, Pages 212-219 Walker, W.S., Kellndorfer, J.M., LaPoint, E., Hoppus, M. & Westfall, J. (2007). An empirical InSAR-optical fusion approach to mapping vegetation canopy height. Remote Sensing of Environment, 109, 482-499. White, J.C., Wulder, M.A., Hermosilla, T., Coops, N.C., & Hobart, G.W. (2017). A nationwide annual characterization of 25years of forest disturbance and recovery for Canada using Landsat time series. Remote Sensing of Environment, 194, 303-321. Williams, C.A., Collatz, G.J., Masek, J.G. & Goward, S.N. (2012). Carbon consequences of forest disturbance and recovery across the conterminous United States. Global Biogeochemical Cycles, 26. Williams, D.W. & Birdsey, R.A. (2003). Historical patterns of spruce budworm defoliation and bark beetle outbreaks in North American conifer forests: an atlas and 131 description of digital maps. Gen. Tech. Rep. NE-308. In: F.S. U.S. Department of Agriculture, Northeastern Research Station (Editor). Newtown Square, PA. Wilson, B.T., Woodall, C. & Griffith, D. (2013). Imputing forest carbon stock estimates from inventory plots to a nationally continuous coverage. Carbon Balance and Management, 8, 1. Wofsy, S.C. & Harris, R.C. (2002). The North American Carbon Program Plan (NACP): A Report of the Committee of the U.S. Carbon Cycle Science Steering Group, Washington, DC: US Global Change Research Program. Woodall, C.W., Coulston, J.W., Domke, G.M., Walters, B.F., Wear, D.N., Smith, J.E., Andersen, H.-E., Clough, B.J., Cohen, W.B., Griffith, D.M., Hagen, S.C., Hanou, I.S., Nichols, M.C., Perry, C.H., Russell, M.B., Westfall, J.A., & Wilson, B.T. (2015). The U.S. forest carbon accounting framework: stocks and stock change, 1990-2016. In (p. 49). Newtown Square, PA: U.S. Department of Agriculture, Forest Service, Northern Research Station. Woodall, C. W., Walters, B. F., Coulston, J. W., D?Amato, A. W., Domke, G. M., Russell, M. B., & Sowers, P. A. (2015). Monitoring Network Confirms Land Use Change is a Substantial Component of the Forest Carbon Sink in the eastern United States. Scientific Reports, 5(October), 1?9. https://doi.org/10.1038/srep17028 Woodbury, P.B., Smith, J.E., & Heath, L.S. (2007). Carbon sequestration in the U.S. forest sector from 1990 to 2010. Forest Ecology and Management, 241, 14-27. Wulder, M.A., Loveland, T.R., Roy, D.P., Crawford, C.J., Masek, J.G., Woodcock, C.E., Allen, R.G., Anderson, M.C., Belward, A.S., Cohen, W.B., Dwyer, J., Erb, A., Gao, F., Griffiths, P., Helder, D., Hermosilla, T., Hipple, J.D., Hostert, P., Hughes, M.J., 132 Huntington, J., Johnson, D.M., Kennedy, R., Kilic, A., Li, Z., Lymburner, L., McCorkel, J., Pahlevan, N., Scambos, T.A., Schaaf, C., Schott, J.R., Sheng, Y., Storey, J., Vermote, E., Vogelmann, J., White, J.C., Wynne, R.H., & Zhu, Z. (2019). Current status of Landsat program, science, and applications. Remote Sensing of Environment, 225, 127-147. USDA Forest Service. ?Region 6 - About the Region,? n.d. https://www.fs.usda.gov/main/r6/about-region. USDA Forest Service. ?Pacific Northwest Region Almanac 2017,? 2017. USDA Forest Service, 2019, https://www.nrs.fs.fed.us/fia/topics/tpo/ Zhang?Huan; Wu?Pengbao,; Yin?Aijing; Yang? Xiaohui ;Zhang?Ming ; Gao? Chao. Prediction of soil organic carbon in an intensively managed reclamation zone of eastern China: A comparison of multiple linear regressions and the random forest model. Science of The Total Environment. 2017, Volume 592, 15 , Pages 704-713 Zhao, F. (2015). Impact of Recent Forest Management and Disturbances on Carbon Dynamics in the Greater Yellowstone Ecosystems, (1). Zhao, F., Huang, C., Goward, S.N., Schleeweis, K., Rishmawi, K., Lindsey, M.A., Denning, E., Keddell, L., Cohen, W.B., Yang, Z., Dungan, J.L., & Michaelis, A. (2018). Development of Landsat-based annual US forest disturbance history maps (1986?2010) in support of the North American Carbon Program (NACP). Remote Sensing of Environment, 209, 312-326. Zhao, Feng R.; Meng, Ran; Huang, Chengquan; Zhao, Maosheng ; Zhao, Feng ; Gong, Peng; Yu, Le and Zhu, Zhiliang. Long-Term Post-Disturbance Forest Recovery in the 133 Greater Yellowstone Ecosystem Analyzed Using Landsat Time Series Stack. Remote Sens. 2016, 8(11), 898; https://doi.org/10.3390/rs8110898. Zhao, S., Liu, S., Sohl, T., Young, C., & Werner, J. (2013). Land use and carbon dynamics in the southeastern United States from 1992 to 2050. Environmental Research Letters, 8(4). https://doi.org/10.1088/1748-9326/8/4/044022 Zheng, D., Heath, L. S., Ducey, M. J., & Smith, J. E. (2011a). Corrigendum: Carbon changes in conterminous US forests associated with growth and major disturbances: Environmental Research Letters, 6(1), 1992?2001. https://doi.org/10.1088/1748- 9326/6/1/019502 Zheng, D., Heath, L. S., Ducey, M.J. & Smith, J.E. (2011b). Carbon changes in conterminous US forests associated with growth and major disturbances: 1992-2001. Environmental Research Letters, 6, 014012. Zheng, D., Heath, L. S., & Ducey, M. J. (2013). Carbon benefits from protected areas in the conterminous United States. Carbon Balance and Management, 8(1), 1?14. https://doi.org/10.1186/1750-0680-8-4 Zheng, D., Heath, L. S., Ducey, M. J., & Butler, B. (2010). Relationships between major ownerships, forest aboveground biomass distributions, and landscape dynamics in the New England region of USA. Environmental Management, 45(2), 377?386. https://doi.org/10.1007/s00267-009-9408-3 Zheng, D., Rademacher, J., Chen, J., Crow, T., Bresee, M., Le Moine, J., & Ryu, S. R. (2004). Estimating aboveground biomass using Landsat 7 ETM+ data across a managed landscape in northern Wisconsin, USA. Remote Sensing of Environment, 93(3), 402?411. https://doi.org/10.1016/j.rse.2004.08.008 134 Zheng, G., Chen, J. M., Tian, Q. J., Ju, W. M., & Xia, X. Q. (2007). Combining remote sensing imagery and forest age inventory for biomass mapping. Journal of Environmental Management, 85(3), 616?623. https://doi.org/10.1016/j.jenvman.2006.07.015 Zhu, Z., Bergamaschi, B., Bernknopf, R., Clow, D., Dye, D., Faulkner, S., Forney, W., Gleason, R., Hawbaker, T., Liu, J., Liu, S., Prisley, S., Reed, B., Reeves, M., Rollins, M., Sleeter, B., Sohl, T., Stackpoole, S., Stehman, S., Striegl, R., & Wein, A. (2010). A method for assessing carbon stocks, carbon sequestration, and greenhouse-gas fluxes in ecosystems of the United States under present conditions and future scenarios. In (p. 188). Reston, VA: US Geological Survey. Zhu, K., Zhang, J., Niu, S., Chu, C., & Luo, Y. (2018). Limits to growth of forest biomass carbon sink under climate change. Nature Communications, 9(1). https://doi.org/10.1038/s41467-018-05132-5 Zeng, H.C., Chambers, J.Q., Negron-Juarez, R.I., Hurtt, G.C., Baker, D.B. & Powell, M.D. (2009). Impacts of tropical cyclones on US forest tree mortality and carbon flux from 1851 to 2000. Proceedings of the National Academy of Sciences of the United States of America, 106, 7888-7892. Zeng, N. (2008), Carbon sequestration via wood burial. Carbon Balance and Management, 3, 1-12. 135