Effects of contemporary land-use and land-cover change on the carbon balance of terrestrial ecosystems in the United States

Changes in land use and land cover (LULC) can have profound effects on terrestrial carbon dynamics, yet their effects on the global carbon budget remain uncertain. While land change impacts on ecosystem carbon dynamics have been the focus of numerous studies, few efforts have been based on observational data incorporating multiple ecosystem types spanning large geographic areas over long time horizons. In this study we use a variety of synoptic-scale remote sensing data to estimate the effect of LULC changes associated with urbanization, agricultural expansion and contraction, forest harvest, and wildfire on the carbon balance of terrestrial ecosystems (forest, grasslands, shrublands, and agriculture) in the conterminous United States (i.e. excluding Alaska and Hawaii) between 1973 and 2010. We estimate large net declines in the area of agriculture and forest, along with relatively small increases in grasslands and shrublands. The largest net change in any class was an estimated gain of 114 865 km2 of developed lands, an average rate of 3282 km2 yr−1. On average, US ecosystems sequestered carbon at an annual rate of 254 Tg C yr−1. In forest lands, the net sink declined by 35% over the study period, largely a result of land-use legacy, increasing disturbances, and reductions in forest area due to land use conversion. Uncertainty in LULC change data contributed to a ~16% margin of error in the annual carbon sink estimate prior to 1985 (approximately ±40 Tg C yr−1). Improvements in LULC and disturbance mapping starting in the mid-1980s reduced this uncertainty by ~50% after 1985. We conclude that changes in LULC are a critical component to understanding ecosystem carbon dynamics, and continued improvements in detection, quantification, and attribution of change have the potential to significantly reduce current uncertainties.


Introduction
Changes in land use and land cover (LULC) and ecosystem disturbances play an important role in the global carbon budget. Since 1850, changes in LULC are estimated to have resulted in the loss of 145 Pg C from terrestrial ecosystems worldwide [1], mostly in response to land-clearing for agriculture and pasture, wood harvest, and fire. Over this period, emissions from LULC change represented approximately 33% of global CO 2 emissions [2]. More recently, global emissions from LULC change and disturbance were estimated at 1.11 ± 0.35 Pg C yr −1 [1]. Notably, carbon losses from LULC change have been offset by factors such as afforestation and expansion of woody plant communities [3], improvements in land management [4,5], climate warming [6,7], and CO 2 fertilization [8], resulting in an estimated net global land sink of 2.6 Pg C yr −1 [9]. Despite improvements in our understanding of the processes that affect carbon storage and flux, uncertainties in the magnitude of LULC impacts remain large [10], stemming from a lack of consistent measurements of land change processes at regional scales spanning sufficiently long time periods.
The United States (US) is one of the world's largest contributors of atmospheric greenhouse gases, releasing 6.3-7.4 Pg CO 2 eq. yr −1 between 1990-2015 [11]. Over the same period, terrestrial ecosystems of the US have been estimated as a persistent carbon sink, offsetting approximately 12% of these emissions annually [11,12]. However, large uncertainties exist in the estimation of the land sink, with estimates ranging from 35-630 Tg C yr −1 [13]. These uncertainties arise from variation among the types of land change processes and disturbances considered, the quality of LULC change data used, the spatial and temporal extent of the study, the types of ecosystems considered, and analytical method applied. A recent comparison of inventory methods, inversion models, and terrestrial biosphere models estimated net ecosystem exchange for the United States ranging from 302-685 Tg C yr −1 , with inversion models generally producing the largest sink estimates and inventory methods the lowest [14]. A similar comparison yielded net carbon sink estimates for North America ranging from 280-890 Tg C yr −1 [15]. In forest ecosystems, a review of disturbance impacts on forest carbon balance indicated that carbon losses due to harvest generally had the highest agreement, yet still ranged from approximately 90-150 Tg C yr −1 [16]. Emissions from fire in the conterminous US had higher variabilities, ranging from 10-133 Tg C yr −1 . While many studies have examined the impact of LULC change on ecosystem carbon balance, few include the effect of urbanization and most do not explicitly incorporate uncertainties in land change estimates, such as classification accuracy, into their simulations.
Multi-temporal remote sensing observations provide an opportunity for reducing uncertainties associated with LULC change and disturbance and their impacts on terrestrial carbon dynamics. The Landsat series of satellites provide a rich source of synopticscale data from which a wide range of land change processes can be identified on a consistent basis dating back to the early 1970s [17]. For example, recent studies have used Landsat data to produce time-series maps of forest disturbance [18][19][20][21][22], which have then been used to estimate corresponding ecosystem carbon fluxes with attribution to various land change processes [23,24]. To date these efforts have been limited to only certain land types (e.g. forests) spanning only a portion of the full Landsat archive. Moreover, most of these studies have focused on only a subset of land change processes (e.g. harvest and fire), and have not incorporated climatic variability or CO 2 fertilization effects.
In this study, we developed a fully coupled model of land change and ecosystem carbon dynamics to estimate the contributions and uncertainties of LULC change and disturbances on carbon stocks and fluxes across major terrestrial ecosystems for the conterminous United States (CONUS; i.e. not including Alaska and Hawaii) for the period 1973-2010. We used improved LULC histories derived from multi-temporal remote sensing, a gain-loss approach with full carbon accounting, and incorporated key uncertainties arising from LULC change data.

Modeling approach
We used the Land Use and Carbon Scenario Simulator (LUCAS) to estimate changes in LULC and ecosystem carbon dynamics for CONUS [25,26]. LUCAS uses a state-and-transition simulation model (STSM) approach to project changes in LULC, whereby a landscape is divided into a set of raster simulation cells, and the state and age of each cell is simulated forward stochastically in time in response to a suite of probabilistic transitions [27]. Changes in carbon (C) storage and flux were integrated into the STSM by adding a stock-flow (SF) model of carbon dynamics [28]. The result was an integrated model of LULC change and carbon gain-loss.
LUCAS was parameterized to run as a spatiallyreferenced model. The land area of CONUS (7 803 313 km 2 ) was stratified spatially into 84 ecoregions [29] and parameterized separately for each strata. However, no spatial interactions between cells within each stratum were considered. The model was simulated at a resolution of 5 × 5 km using an annual time-step for the period 1973-2010, and repeated for 100 Monte Carlo realizations. All of the simulations were generated using the version 3.1.2 of ST-Sim, running under SyncroSim software version 2.0.3 [30], with model inputs and outputs prepared using the R software version 3.3.1 [31], including the rsyncrosim package [30].
For the STSM portion of the model we tracked state type and age for each cell and timestep. A total of 12 possible state types were defined for each of the 84 ecoregions ( figure 1(a)), consisting of combinations of six LULC classes (Agriculture, Grassland, Shrubland, Forest, Developed, and Wetland), crossed with two levels of protection (Protected, Unprotected). Areas classified as Barren, Water, and Snow/Ice in 1973 were considered static for the duration of the simulation. Transitions were defined to represent possible changes in the state type of each cell over time, including Urbanization, Agricultural Expansion, Agricultural Contraction, Forest Harvest, Wildfire, and Protection (i.e. the establishment of new protected areas). Additional details on the LUCAS model, including validation of results, are reported in Sleeter et al [32,33]. The integrated SF model tracked eight carbon stocks as continuous state variables (figure 1(b)): Atmosphere, Living Biomass, Standing Deadwood, Down Deadwood, Litter, Soil, Aquatic, and three harvested products pools (Grain, Straw, and Harvested Wood). The Atmosphere, Aquatic, and three harvest product pools were included to maintain a carbon mass balance in the model simulations. We defined the following flows to transfer carbon between stock types: Growth, Mortality, Litterfall, Deadfall, Decay, Decomposition, Emission, Leaching, and Harvest. For this study, we define 'growth' as the amount of carbon fixed by photosynthesis added to a cell after plant (autotrophic) respiration but prior to any reductions resulting from mortality and/or litterfall. We classified flows as either base or transition-triggered. Base flows were carbon fluxes between stocks due only to the passage time (irrespective of changes in LULC), while transition-triggered flows were carbon fluxes between stocks that occurred only in response to LULC transitions.

Land use and land cover change
We used four remote sensing derived data sources to estimate changes in LULC classes over the Landsat record (table 1). All data sources draw primarily from archived Landsat imagery data spanning the last 45 years. A national database of protected areas was used to allocate lands into protected and unprotected categories and to identify establishment of new protected areas over the model simulation period. A spatially explicit map of forest age was used to initialize the forest classes (table 1).
For the period 1973-2000, Land Cover Trends (LCT) data were used to estimate all land-use change transitions (i.e. urbanization, agricultural expansion and contraction; figure 2). For each simulated timestep and Monte Carlo realization, a target area for each transition type was generated by sampling from a normal distribution based on the annualized reported mean rate of change and standard error for the corresponding period and ecoregion. Similarly, data from the National Land Cover Database (NLCD) was used to generate annual transition area targets from 2001-2010 (figure 2), again by sampling from a normal distribution based on the annualized reported mean rate of change (for each of two periods), and the standard error estimated from accuracies reported for each LULC change theme [40]. LCT change rates for the 1992-2000 period were extended to 2001 to fill the gap between LCT and NLCD datasets (figure 2).
Beginning in 1986, data from the North American Forest Dynamics (NAFD) dataset were used to estimate area targets for the 'forest to developed', 'forest to agriculture', and 'forest harvest' transitions. NAFD data were first spatially intersected with burn perimeters from the Monitoring Trends in Burn Severity (MTBS) data for each historical year. Disturbed forest cells in NAFD were removed if a wildfire was mapped in the same location ± 2 years of the forest disturbance. Next, the remaining NAFD disturbance cells were intersected with the 2011 NLCD land cover map. Any disturbed forest cells from NAFD classified as developed in the 2011 NLCD were assumed to be disturbances due to urbanization. These cells were summarized as 'forest to developed' transitions in the year the disturbance was first mapped. The same process was repeated for 'forest to agriculture' transitions, where any NAFD disturbed cells mapped as agriculture in the 2011 NLCD were assumed to be the result of 'forest to agriculture' transitions. All remaining forest disturbances were assumed to be caused from harvest activities. Based on an accuracy assessment of NAFD, we assumed all NAFD derived transitions were normally distributed ± 30% [43]. For the period 1973-1985, the rate of forest harvest in each timestep and Monte Carlo realization was sampled, with replacement, from the full distribution of ecoregion-based NAFD disturbance rates. Spatially explicit 30 meter maps of land cover for the conterminous United States [37][38][39] based on a hierarchical classification system. Accuracy assessments have been conducted for land cover and land cover change themes [40].

North American Forest Dynamics (NAFD)
Annual; 1985-2010 NAFD data consist of 30 meter spatially explicit maps of forest disturbances collected on an annual basis spanning the 1985-2010 period for the conterminous United States [18]. NAFD data only identify the location of a forest disturbance and not the causal agent, thereby requiring additional processing to identify the types of change occurring. Validation of the NAFD approach and general accuracy estimates are described in Thomas Figure 2. Sources of data used to parameterize the STSM transition rates. Land Cover Trends (LCT) data are from Sleeter et al [36], North American Forest Dynamics (NAFD) data from Goward et al [18], National Land Cover Database (NLCD) from Homer et al [39,44] and Fry et al [45], Monitoring Trends in Burn Severity (MTBS) from Eidenshink et al [22], Protected Areas Database (PADUS) from the US Geological Survey [41].
Monitoring Trends in Burn Severity (MTBS) burn perimeter and severity maps were used to estimate the total area of forest, grassland, and shrubland burned in each timestep, for each ecoregion, for each severity class (high, medium, and low), and for protected and unprotected land classifications. For the period 1973-1984 we sampled, with replacement, a year from the 1984-2000 MTBS time series for each simulated timestep and realization. Lastly, data from Protected Areas Database of the US (PADUS) was used to estimate the establishment of protected areas in each timestep. Additional details regarding the STSM for land change are found in section S1 available at stacks.iop.org/ERL/13/045006/mmedia.

Carbon stocks and flows
Similar to the approach used by Daniel et al [28], we used the Integrated Biosphere Simulator (IBIS), a dynamic global vegetation model [46,47], to estimate the base flow rates for Forest, Grassland, Shrubland, and Agriculture in each ecoregion [25,26]. The IBIS model was initiated with minimum vegetation (i.e. bare-ground), and then simulated forward for 250 years using the average climate for the period 1970-2000 and a fixed CO 2 level of 330 ppm (∼1970 level). IBIS was run for 250 years to allow vegetation to grow and generate baseline carbon stocks and fluxes by age, where age represents the time since the start of the bare-ground simulation. With the exception of the Growth flow, carbon flux rates for each LULC class and ecoregion were estimated by dividing the IBIS estimated fluxes by the size of the 'from stock' for each age. The Growth flow was estimated as the IBIS-generated carbon flux from Atmosphere to Living Biomass in each timestep (kg C m −2 ). For Grassland, Shrubland, and Agriculture, growth was estimated for each state type and ecoregion as the net primary production (NPP) averaged over the 250 year IBIS simulation; Growth for the Forest state type was estimated for each combination of ecoregion and age. For all other state types, transitiontriggered flow rates were set to zero and no fluxes of carbon were simulated. In addition to the baseline flow rates, transition-triggered flow rates [28] were used to model the transfer carbon between stocks in response to Urbanization, Agricultural Expansion, Forest Harvest, and all Wildfire transitions and were derived from literature values where available. Additional details on the SF model are in section S2.

Model initialization
The initial state type of each cell (i.e. in 1973) was based on 1973 LCT data (for the LULC class) and the 1973 protected areas map (for Protected vs Unprotected) using NLCD (2011) data to split the single Grassland/Shrubland class (see section S1). All realizations of the model were initialized with the same distribution of state type areas. Initial forest age was derived from a forest age map for 2006 [48], which was then modified to reflect age in 1973. The forest age in 1973 was estimated by subtracting 33 years from the 2006 age estimates on a cell by cell basis. Cells with a resulting forest age less than or equal to zero were imputed based on the mean forest age for the corresponding ecoregion (see section S3.2 for additional details). The forest age map was overlaid with the 1973 protected areas map to generate forest age distributions for protected and unprotected forest lands in each ecoregion.
Initial ecosystem carbon pools were based on a look-up table approach where each cell was assigned initial carbon stock value (live biomass, standing and down deadwood, litter, and soil) based on the combination of the strata of each cell (e.g. ecoregion), the state type of each cell (e.g. LULC class) and the age of each cell. All other stocks (e.g. Aquatic, Atmosphere, Harvest Products) and LULC classes (i.e. Developed, Barren, Snow/Ice, Water) were initialized to zero. Carbon stock look-up values were based on the estimates of carbon stock and age from the IBIS cold-start simulation.

Land change dynamics in the conterminous United States
Expansion and contraction of agriculture accounted for 317 830 km 2 of gross land area change over the study period, with losses of 5981 km 2 yr −1 and gains of 3799 km 2 yr −1 . Nationally, agriculture increased in area from 1973-1980 followed by large declines in the 1980s and 1990s, resulting in an estimated net decline of 80 729 km 2 . Grasslands experienced an estimated net increase of 28 453 km 2 over the study period. The total amount of land converted to development was estimated at 114 865 km 2 , with conversion from agriculture accounting for 38% of the total. Forest area, the largest proportion of land area in 1973, declined by an estimated 2.8% (−64 580 km 2 ) by 2010. Urbanization was the largest driver of net forest area change, with an estimated mean annual conversion rate of 1410 km 2 yr −1 . This was partially offset by the net balance of forest and agricultural land conversion; deforestation (forest converting to agriculture) occurred at a rate ∼40% higher than afforestation (agriculture converting to forest).
Forest harvest occurred on an estimated 22 612 km 2 yr −1 across the conterminous US, with the highest rates occurring in the Southeast and Pacific Northwest regions. Wildfire was estimated at 11 738 km 2 yr −1 , with a trend of increasing area over time ( figure 3). Wildfire was nearly ten times as likely on protected forests, which had an average disturbance rate of ∼1.1% yr −1 . When harvest was included, the annual disturbance rate of unprotected forests was ∼1.4% yr −1 . Additionally, new protected areas made up 2.6% of the CONUS land area. The majority of new establishments were shrublands, followed by forests and grasslands.

Carbon stocks and total ecosystem carbon storage
For the purposes of this study, total ecosystem carbon (TEC) is the combined carbon stored in terrestrial ecosystems (forest, grasslands, shrublands, and agriculture) in living biomass, standing and down deadwood, litter, and soil (2 meters in depth). In 1973, TEC is estimated to have been 100.0 Pg C, of which 73% was stored in soils and 17% in living biomass (figure 4). By 2010, TEC increased to 109.4 Pg C, with the largest increases associated with living biomass (+2.6 Pg C), deadwood (+3.0 Pg C), and soils (+1.7 Pg C) in forest ecosystems (figure 4). We also estimated an additional 4.3 Pg C was transferred to the harvested wood products (HWP) pool and  0.9 Pg C was transferred to the aquatic pools over the 37 year study period, which was roughly equivalent to the combined net increase in living biomass and soils (figure 4).  declined from 142 g C m 2 yr −1 to 98 g C m 2 yr −1 , resulting in a net decline in the strength of the forest carbon sink of 114 Tg C yr −1 . At the current sequestration rate (i.e. 2010), the estimated forest area lost due to land-use change represents a reduction in the forest carbon sink of 6 Tg C yr −1 , or approximately 5% of the total decline in the forest carbon sink. The remaining reduction is largely attributed to declines due to ageing, which can also be affected by changes in stand-replacing disturbance events, including a reduction in the rate of harvest (figure 3) and an increase in burn area. The net agriculture carbon sink increased over time, while grasslands and shrublands were generally carbon neutral, fluctuating between source and sink based on annual variability in productivity (figure 5). Two periods (1987)(1988)(1989)(1999)(2000)(2001)) experienced noticeably large declines in the size of the net land sink, with one year (1988) estimated as a net source of carbon to the atmosphere (figure 5). Both periods corresponded to exceptional drought conditions across large areas of the CONUS. LULC and disturbance impacts were also important during these periods, with forest harvest removing 172 Tg C yr −1 in 1987 and 159 Tg C yr −1 in 1988, the two highest years in our study. Wildfire emissions during this same two-year period were estimated to be two times higher than the average rate, exemplified by the 1988 fires in Yellowstone National Park. Despite drought conditions and amplified losses from LULC and disturbance, forests were still estimated as a small net sink. However, agriculture was estimated as a large source because of extreme drought conditions across much of the region.

Carbon balance of the conterminous United States
Harvest activities accounted for the largest annual carbon removals, totaling an estimated 523 Tg C yr −1 ( figure 6). Harvest of agricultural products accounted for 78% of these removals and were assumed to be returned to the atmosphere each year (table 2). The remainder of removals was the result of forest harvest activities, with a small amount attributed to clearing of land for development and agriculture ( figure 6). Estimates of total carbon emissions associated with LULC change and wood harvest averaged 73 Tg C yr −1 . Of these, urbanization had the largest effect, removing an average of 25 Tg C yr −1 . Wildfire emissions were of similar magnitude as those from forest harvest activities ( figure 6). However, it is important to note that emissions associated with the production and decay of HWP were not considered in this study, and that the total contribution of forest harvest activities to atmospheric carbon is likely much higher. Our approach also considered the lateral flux of carbon from terrestrial to aquatic systems (i.e. leaching, runoff), which was of the same magnitude as urbanization-related carbon losses ( figure 6).

Discussion
Since 1970, urbanization, the expansion and contraction of agriculture, forest harvest, and wildfire have resulted in significant changes to LULC composition in the CONUS. We estimated large net declines in area of both agriculture and forest along with relatively small increases in grasslands and shrublands. Table 2. Carbon balance for terrestrial ecosystems. Stock is the stock size in 1973. NPP is net primary production, Rh is heterotrophic respiration, NEP is net ecosystem productivity, HARV is carbon losses due to harvest activities, URB is losses due to urbanization transitions, FIRE is losses due to wildfire, AGR is losses due to agricultural expansion, DIST is losses due to background disturbance, LAT is the lateral flux to the aquatic pool, NBP is net biome productivity, TRF is for transfers of carbon between classes due to changes in state class types, and NECB is the net ecosystem carbon balance calculated as NBP minus TRF. All units are in Tg C yr −1 .  The largest net change in any class was an estimated 114 865 km 2 increase in developed lands. Taken together, we estimate terrestrial ecosystems of the CONUS sequestered an average of 254 Tg C yr −1 between 1973 and 2010. There was no noticeable trend in the land carbon sink over the 35 year study period, which averaged 282 Tg C yr −1 in the 1970s, 224 Tg C yr −1 in the 1980s, 284 Tg C yr −1 in the 1990s, and 241 Tg C yr −1 in the 2000s. However, there was considerable inter-annual variability in the land carbon sink driven by short-term climatic variations, including at least one year (1988) where terrestrial ecosystems of the CONUS were an overall net carbon source to the atmosphere. The most notable long-term trend was a decline in total forest carbon sink strength due to ageing and deforestation. Prior to 1985, uncertainty in LULC change data translated to a wide margin of CONUS NBP estimates (±40 Tg C yr −1 ). Uncertainty in NBP estimates due to LULC change was reduced by nearly 50% after 1985, largely due to increased spectral and spatial resolution, improved mapping capabilities, and increased data availability.

Comparison to other studies
Our estimates of the forest carbon sink are generally similar to a growing body of literature centered on a net annual rate of ∼200 Tg C yr −1 . Over the 35 year study period, we estimated an average sink in forests of 220 Tg C yr −1 , similar to a range of studies incorporating a variety of methods and spanning a range of temporal periods [11,[50][51][52][53][54]. A synthesis conducted for the First State of the Carbon Cycle Report (SOCCR; [3]) estimated a net forest sink of 220 Tg C yr −1 during the 1980s, which is similar to our estimate of 213 Tg C yr −1 over the same period. A separate analysis of inventory-based methods, atmospheric inversion modeling, and terrestrial biosphere models [14] estimated the net forest carbon sink was between 157-282 Tg C yr −1 for the period 2000-2006, which encapsulates our estimate of 171 Tg C yr −1 . Zheng et al [51] estimated a net decline in forest area and an average annual net carbon uptake of 180 Tg C yr −1 by CONUS forests between 1992-2001 using a combination of inventory and land-cover change data, which compared well with our estimate of 186 Tg C yr −1 . Process model estimates of the mean annual CONUS forest sink from 1901-2010 were slightly higher at 206 Tg C yr −1 [55], with NBP peaking in the 1970s in response to regrowth from afforestation in prior decades, and a subsequent decline in the 1980s and 1990s due to forest ageing. Our results identify this pattern as well; we estimated a net forest sink of 300 Tg C yr −1 during the 1970s with a steady decline to 191 Tg C yr −1 in the 2000s. A continuation of this decline in future decades would reduce the amount of anthropogenic CO 2 emissions offset by terrestrial ecosystems, which has important implications for climate change mitigation strategies.
Fewer studies have estimated the integrated carbon balance for forest, grasslands, shrublands and agriculture in the CONUS. A recent assessment conducted by the US Geological Survey (USGS) [56][57][58] estimated that terrestrial ecosystems sequestered 431 Tg C yr −1 from 2001-2005, nearly double our estimate for the same time period (237 Tg C yr −1 ). A large portion of the difference was due to the USGS estimate of 292 Tg C yr −1 being sequestered in forests, >100 Tg C yr −1 larger than our estimate for the same period. In the western US, from 2001-2005, Liu et al [56] estimated that grasslands and shrublands combined ranged from a net sink of 58.7 Tg C yr −1 to a net source of 7 Tg C yr −1 with the variation largely attributable to different estimation methods. Our results indicate the shruband grassland sink is towards the lower end of these estimates, with grassland estimated at 5 Tg C yr −1 and shrublands as a net source of −2.2 Tg C yr −1 . Hurtt et al [54] estimated that CONUS terrestrial ecosystems sequestered 330 Tg C yr −1 during the 1980s, which was higher than our estimate of 224 Tg C yr −1 . While our estimates of the forest sink are similar, Hurtt estimated an additional net sink of 130 Tg C yr −1 on non-forest lands and pastures, whereas we estimated these lands were a net source of −13 Tg C yr −1 over the same period. Given the large land area of non-forest ecosystems, reducing uncertainties in these sectors should be an important objective of national-scale carbon cycle research.

Component fluxes and attribution of carbon change
Forest harvest activities transfer carbon both within ecosystems (i.e. from live to dead wood), and outside of ecosystems to harvest product pools and to the atmosphere through direct CO 2 emissions. Our approach, using a time-series of forest disturbance maps and a causal attribution process, estimated that on average, 174 Tg C yr −1 was moved between various carbon pools due to harvest activities. Of this total, 54 Tg C yr −1 was transferred from live to dead carbon pools within ecosystems, 107 Tg C yr −1 was allocated to offsite products (HWP), and 13 Tg C yr −1 was emitted directly to the atmosphere. Fluxes of carbon due to harvest activities have varied over time, with the combined transfers due to emission, harvest, and mortality being reduced by nearly 50% since the mid-1980s. Our estimates of harvest compare favorably to a range of other studies, largely based on forest inventories. For example, the US Environmental Protection Agency estimated harvest removals at 143 Tg C in 1990, 134 Tg C in 2000, and 95 Tg C in 2009 [59,60], whereas this study estimated the combined harvest+emission removal at 135, 133, and 84 Tg C yr −1 for the same years. More generally, harvest removals have been estimated in a range of studies spanning different temporal periods and utilizing different modeling and accounting approaches, and generally range from 90-150 Tg C yr −1 [16], a range which supports the findings presented in this study.
We estimated wildfire-related transfers of live carbon to deadwood pools and direct atmospheric emissions due to combustion for three wildfire severity classes based on annual burn area maps from the MTBS project. Our study likely underestimated the effect of fire since MTBS data only accounted for events greater than 202 hectares in the eastern US and 404 hectares in the western US A new generation of spatially explicit annual burn area maps has recently been published which suggests the total burn area is more than double the estimate from MTBS, primarily because MTBS omits a large proportion of burned area in the Great Plains and the eastern US [61]. Furthermore, we used a generalized set of carbon flux parameters based on default combustion factors from the IPCC. We estimated that on average, wildfire resulted in the direct loss of 13 Tg C yr −1 with an additional 4 Tg C transferred from live to dead pools through mortality fluxes. Ghimire et al [62] used a process model and MTBS data to estimate carbon fluxes from fire for the western United States and estimated direct emissions of 4 Tg C yr −1 and mortality of 10.5 Tg C yr −1 for the period 1984-2008. Parameterization of carbon flux parameters was far more detailed than our approach, and while the total live biomass flux was similar, results suggest our study may overestimate combustion and underestimate mortality. However, other studies suggest carbon losses due to fire may be higher [51,63]. For example, Wiedinmyer and Neff used MODIS active fire maps and coarse resolution land cover data and estimated the average fire emission rate of 213 Tg CO 2 yr −1 (57 Tg C yr −1 ) for the period 2002-2006 [63]. However, it is important to note that wildfire burned >50% more area during this period as compared to the longterm average of our study period, suggesting that the period of analysis is an important consideration in understanding the broader contributions of wildfire to ecosystem carbon balance. Furthermore, reconciling estimates of burn area would likely lead to a significant reduction in uncertainty of ecosystem carbon balance estimates.

Consideration of uncertainty
While considerable uncertainty remains in our understanding of the US carbon budget, one potential source of uncertainty in our estimates stems from the absence of insect-induced mortality, such as bark beetles in the western United States. Berner et al [64] estimated bark beetles in the western US were responsible for the loss of 14.6 Tg C yr −1 (above ground live carbon only) between 2003 and 2012. However, annual emissions were estimated at only 2 Tg C yr −1 , indicating the true impact of these mortality events would be realized in future decades due to the time-lag associated with the decomposition of coarse woody debris. Storms pose another risk to ecosystem carbon storage. Zeng et al [65] estimated that tropical cyclones alone were responsible for the loss of 39 Tg C yr −1 of live biomass during the 20 th century while Dahal et al [66] estimated a loss of 18.2 Tg C yr −1 over a similar period. Widespread forest mortality events driven by extreme and persistent drought and increasing temperature have become more common and pose yet another risk to the terrestrial carbon sink [67]. While we did not explicitly model the direct effects of these disturbances (insects, storms, missing fires), we did include a background biomass disturbance rate parameterized for each ecoregion (average of 0.2% yr −1 ) which resulted in the combined annual loss of 58 Tg C yr −1 from terrestrial ecosystems. A focus of future research should be to refine this parameter and, where possible, attribute it towards explicit disturbances.
Remote sensing can play an important role in detecting and quantifying the magnitude and frequency of a wide range of disturbance events. Aside from programs designed specifically to map wildfire, no comprehensive multi-temporal national-scale data exist which identify a wide range of disturbance events, including harvest activities, insects, storms, disease, and a wide range of human management actions. Improving the ability of remote sensing to accurately attribute various disturbance types would greatly improve the ability to estimate changes in ecosystem carbon balance and reduce uncertainties associated with these transitions. Furthermore, improved characterization of land change processes through repeated remote sensing observations will help reconcile differences in the magnitude and trajectories of forest change estimates in the United States.
Another important uncertainty not explored in our study is the lateral transfer and interaction of carbon between terrestrial and aquatic systems. Research suggests that exclusion of lateral transfers of carbon may result in a significant overestimation of ecosystem NEP [68,69]. While we did include an estimate of carbon leached from soils, it did not include the effect of sediment burial or emissions from lakes and streams. Butman et al [68] estimated a total net lateral flux of 41.5 Tg C yr −1 , compared with our estimate of 26 Tg C yr −1 . Furthermore, Butman et al estimated an additional 85.3 Tg C yr −1 of emissions from water bodies which was offset by 20.6 Tg C yr −1 of burial in lakes and reservoirs. Overall, freshwater aquatic ecosystems were estimated as a net source of carbon at a rate of 106 Tg C yr −1 , an amount roughly equivalent to the annual carbon removal due to forest harvest activities. However, a thorough comparison of results is difficult due to inconsistencies in boundary conditions, definitions, and consideration of processes. Future research is needed to develop fully coupled terrestrial-aquatic ecosystem carbon linkages to reduce uncertainties and close the gap in the global carbon budget.
This study has explored important sensitivities and uncertainties associated with changes in LULC and disturbances on terrestrial ecosystem carbon balance by incorporating detailed land change histories derived from multi-temporal remote sensing observations and their associated estimates of uncertainty.
Uncertainties reported in this paper are the result of a specific set of key model parameters, namely the rates, types, and geographic and temporal patterns of LULC change. Uncertainties were derived by incorporating map accuracy information and statistical estimation data into the model and using the derived statistical distributions within a stochastic Monte Carlo sampling framework. While uncertainties in some LULC parameters were incorporated, other important LULC parameters were not, namely uncertainties associated with wildfire area and severity, and establishment of protected areas. In addition, due to a lack of observational data, some processes were not represented, such as insect disturbances and storm events. Furthermore, there are a number of other key model parameters for which uncertainty was not incorporated, including the distribution of initial conditions (e.g. LULC area, forest age structure), carbon flux rates, including both base rates and fluxes resulting from LULC change, the effects of climate variability and CO 2 fertilization, and carbon pool size. As such, the uncertainties attributed to historical changes in LULC reported here, represent a relatively small subset of the overall uncertainty attributable to model parameters on the terrestrial carbon balance of the United States.
The LUCAS model used for this study is built on a generalized modeling framework, utilizing a fullycoupled state-and-transition simulation model with stocks and flows [28]. One benefit of the generic modeling framework described here is the ability to estimate uncertainties resulting from model structure [70] which have been found to be important factors in estimating ecosystem carbon balance [71,72]. Unlike many models of LULC and ecosystem carbon dynamics, the structure of the LUCAS model can be modified to include virtually any configuration of LULC classes and transitions and carbon stocks and carbon fluxes. This flexibility allows the representation of a wide range of controlling processes. For example, we integrated parameters derived from a process-based dynamic global vegetation model (IBIS), an empirical model of NPP, and empirical data on the effect of CO 2 fertilization to track carbon through a generalized set of ecosystem pools. Future research should explore how sensitive our estimates are to key structural uncertainties by utilizing different models to derive base carbon parameters, modifying the number of carbon pools and fluxes, and/or change the number of LULC classes and transitions. Additional processes can also be incorporated, such as including a broader set of carbon stocks to track the fate of carbon as it moves through various production pools and end-use stages. Likewise, the LUCAS model can be modified to include a range of aquatic carbon pools, such as rivers, lakes and streams, and wetlands, which would allow an integrated estimate of the combined impact of terrestrial and aquatic systems on total ecosystem carbon balance.

Conclusion
In this study we used a range of remote sensing data to estimate the effects of LULC change and disturbances on terrestrial ecosystem carbon balance in the CONUS over 35 years while controlling for climatic variability and change. We found that CONUS ecosystems were, on average, a net carbon sink (254 Tg C yr −1 ), but that climatic variability had a strong impact on the magnitude of the annual carbon sink, ranging from a 496 Tg C yr −1 sink in 1983 to a net source of −152 Tg C yr −1 in 1988. Although forests represented the largest proportion of the CONUS land sink, we estimated the net forest sink declined by 35% over the study period (from 300 to 197 Tg C yr −1 ) due to ageing and an overall reduction in forest area from land use change. By incorporating uncertainty in LULC change data into our estimates, we found that the period prior to 1985 had roughly two-times the uncertainty as more recent years. Reduced uncertainties were largely attributable to the availability of annual maps of wildfire and forest disturbance which have not been developed using older Landsat sensors (i.e. multi-spectral scanner, MSS). We suggest incorporating uncertainty into estimates of ecosystem C balance can only improve the confidence in estimates going forward. The LUCAS model can be used to explore a wide range of additional uncertainties, including the relative impact of uncertainties in carbon flux rates on carbon sink strength estimates. Perhaps more importantly, the LUCAS model can be used to explore a wide range of future projections, such as establishing empirical baselines to evaluate future climate mitigation strategies.