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

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 spatially-referenced model. The land area of CONUS (7 803 313 km²) 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].

Zoom In Zoom Out Reset image size

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.

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).

Zoom In Zoom Out Reset image size

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.

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.

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₂ 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⁻²). 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, transition-triggered 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.

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.

Expansion and contraction of agriculture accounted for 317 830 km² of gross land area change over the study period, with losses of 5981 km² yr⁻¹ and gains of 3799 km² yr⁻¹. 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². Grasslands experienced an estimated net increase of 28 453 km² over the study period. The total amount of land converted to development was estimated at 114 865 km², 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²) by 2010. Urbanization was the largest driver of net forest area change, with an estimated mean annual conversion rate of 1410 km² yr⁻¹. 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² yr⁻¹ across the conterminous US, with the highest rates occurring in the Southeast and Pacific Northwest regions. Wildfire was estimated at 11 738 km² yr⁻¹, 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⁻¹. When harvest was included, the annual disturbance rate of unprotected forests was ~1.4% yr⁻¹. 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.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

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).

Average annual NPP from 1973–2010 was estimated at 2811 Tg C yr⁻¹. Forest and agricultural lands accounted for nearly 85% of total NPP while accounting for 56% of the total land area. Since 1973, NPP increased by ~200 Tg C despite overall reductions in photosynthetic area from land-use change. Total litter and soil emissions (i.e. heterotrophic respiration, Rh) were estimated at 1876 Tg C yr⁻¹. Mean annual net ecosystem productivity (NEP), calculated by subtracting Rh from NPP, was estimated at 935 Tg C yr⁻¹. Despite large inter-annual variation caused by climate variability, NEP remained relatively stable over the 37 year period. Net biome productivity (NBP) was calculated by subtracting removals and emissions of carbon due to LULC change, disturbance, harvest activities, and lateral transfers due to leaching, from NEP. We estimate that ecosystems in CONUS represented a net sink of 254 Tg C yr⁻¹ between 1973 and 2010. However, annual variability in the strength of the sink was large, with NBP estimates ranging from +496 Tg C yr⁻¹ in 1983 to a net source of −152 Tg C yr⁻¹ in 1988 (figure 5).

Zoom In Zoom Out Reset image size

Forests accounted for the largest proportion of the net sink, but forest sink strength declined over time (figure 5), primarily as a result of forest ageing. Forest NBP was, on average, 308 Tg C yr⁻¹ between 1973 and 1983 and declined to 191 Tg C yr⁻¹ from 2000–2010. The carbon sequestration rate of 'forest remaining forest' declined from 142 g C m² yr⁻¹ to 98 g C m² yr⁻¹, resulting in a net decline in the strength of the forest carbon sink of 114 Tg C yr⁻¹. 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⁻¹, 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–1989, 1999–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⁻¹ in 1987 and 159 Tg C yr⁻¹ 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.

Harvest activities accounted for the largest annual carbon removals, totaling an estimated 523 Tg C yr⁻¹ (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⁻¹. Of these, urbanization had the largest effect, removing an average of 25 Tg C yr⁻¹. 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).

Zoom In Zoom Out Reset image size

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⁻¹. Over the 35 year study period, we estimated an average sink in forests of 220 Tg C yr⁻¹, similar to a range of studies incorporating a variety of methods and spanning a range of temporal periods [11, 50–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⁻¹ during the 1980s, which is similar to our estimate of 213 Tg C yr⁻¹ 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⁻¹ for the period 2000–2006, which encapsulates our estimate of 171 Tg C yr⁻¹. Zheng et al [51] estimated a net decline in forest area and an average annual net carbon uptake of 180 Tg C yr⁻¹ 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⁻¹. Process model estimates of the mean annual CONUS forest sink from 1901–2010 were slightly higher at 206 Tg C yr⁻¹ [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⁻¹ during the 1970s with a steady decline to 191 Tg C yr⁻¹ in the 2000s. A continuation of this decline in future decades would reduce the amount of anthropogenic CO₂ 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–58] estimated that terrestrial ecosystems sequestered 431 Tg C yr⁻¹ from 2001–2005, nearly double our estimate for the same time period (237 Tg C yr⁻¹). A large portion of the difference was due to the USGS estimate of 292 Tg C yr⁻¹ being sequestered in forests, >100 Tg C yr⁻¹ 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⁻¹ to a net source of 7 Tg C yr⁻¹ with the variation largely attributable to different estimation methods. Our results indicate the shrub- and grassland sink is towards the lower end of these estimates, with grassland estimated at 5 Tg C yr⁻¹ and shrublands as a net source of −2.2 Tg C yr⁻¹. Hurtt et al [54] estimated that CONUS terrestrial ecosystems sequestered 330 Tg C yr⁻¹ during the 1980s, which was higher than our estimate of 224 Tg C yr⁻¹. While our estimates of the forest sink are similar, Hurtt estimated an additional net sink of 130 Tg C yr⁻¹ on non-forest lands and pastures, whereas we estimated these lands were a net source of −13 Tg C yr⁻¹ 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.

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₂ emissions. Our approach, using a time-series of forest disturbance maps and a causal attribution process, estimated that on average, 174 Tg C yr⁻¹ was moved between various carbon pools due to harvest activities. Of this total, 54 Tg C yr⁻¹ was transferred from live to dead carbon pools within ecosystems, 107 Tg C yr⁻¹ was allocated to offsite products (HWP), and 13 Tg C yr⁻¹ 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⁻¹ 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⁻¹ [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⁻¹ 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⁻¹ and mortality of 10.5 Tg C yr⁻¹ 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₂ yr⁻¹ (57 Tg C yr⁻¹) 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 long-term 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.

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⁻¹ (above ground live carbon only) between 2003 and 2012. However, annual emissions were estimated at only 2 Tg C yr⁻¹, 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⁻¹ of live biomass during the 20^th century while Dahal et al [66] estimated a loss of 18.2 Tg C yr⁻¹ 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⁻¹) which resulted in the combined annual loss of 58 Tg C yr⁻¹ 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⁻¹, compared with our estimate of 26 Tg C yr⁻¹. Furthermore, Butman et al estimated an additional 85.3 Tg C yr⁻¹ of emissions from water bodies which was offset by 20.6 Tg C yr⁻¹ 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⁻¹, 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₂ 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 fully-coupled 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₂ 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.