Prototyping a Methodology for Long-Term (1680–2100) Historical-to-Future Landscape Modeling for the Conterminous United States

Abstract

Land system change has been identified as one of four major Earth system processes where change has passed a destabilizing threshold. A historical record of landscape change is required to understand the impacts change has had on human and natural systems, while scenarios of future landscape change are required to facilitate planning and mitigation efforts. A methodology for modeling long-term historical and future landscape change was applied in the Delaware River Basin of the United States. A parcel-based modeling framework was used to reconstruct historical landscapes back to 1680, parameterized with a variety of spatial and nonspatial historical datasets. Similarly, scenarios of future landscape change were modeled for multiple scenarios out to 2100. Results demonstrate the ability to represent historical land cover proportions and general patterns at broad spatial scales and model multiple potential future landscape trajectories. The resulting land cover collection provides consistent data from 1680 through 2100, at a 30-m spatial resolution, 10-year intervals, and high thematic resolution. The data are consistent with the spatial and thematic characteristics of widely used national-scale land cover datasets, facilitating use within existing land management and research workflows. The methodology demonstrated in the Delaware River Basin is extensible and scalable, with potential applications at national scales for the United States.

1. Introduction

A variety of interdependent natural and socioeconomic processes drive Earth system change, with destabilization in these processes potentially having detrimental impacts on society. Land system change, defined as the amount and pattern of change in terrestrial biomes, and biogeophysical processes in land systems that regulate climate, has been identified as one of four major processes where change has passed a destabilizing threshold (along with climate change, change in biosphere integrity, and biogeochemical flows) [1]. Anthropogenic influences on the Earth’s land system have been occurring for millennia, yet have been amplified in recent decades. Factors that have contributed to anthropogenic-based land cover change include growth in human populations and societal and technological advances that have increased our ability to rapidly influence the landscape [2,3]. In order to understand landscape change and the impacts it has on natural and anthropogenic systems, a unified Earth System Science (ESS) framework has emerged around three interrelated foci: (1) observation and characterization of ESS, (2) computer simulations of future system dynamics, and (3) assessments and synthesis of interrelationships of ESS components [4].

One key element of a unified ESS framework is the characterization of the Earth’s landscape via land cover. Land cover data capture anthropogenic activity on the physical landscape, thereby representing a unique connection point between physical systems (e.g., hydrologic, biogeochemical, ecological) and human behavior. The effects of policy decisions, land management, and energy usage are expressed in quantifiable ways on the landscape. To understand landscape change and the impacts it has on human and natural systems, a historical record of land cover is required, along with corresponding data on the impacts of that change [5,6]. Similarly, scenarios of future landscape change are required to facilitate landscape planning and potentially mitigate negative impacts [7,8].

Remote sensing provides an ideal platform for the mapping, monitoring, and characterization of landscape change over broad geographic regions. The Landsat platform, for example, provides the most consistent, continuous coverage across the globe from 1972 to the present, and offers the best means to map both current and historical change. For the conterminous United States (US), Landsat data are the basis of multiple national-scale mapping efforts. The National Land Cover Database (NLCD) is the most widely used land cover dataset for the United States, with the 1992 NLCD being the first such product for the Nation, and subsequent iterations (2001, 2006, 2011, 2016) continually improving product accuracy and content [9,10]. The LANDFIRE project similarly uses Landsat data to produce national-scale landscape information for the US, with a focus on the mapping of landscape variables that inform pre- and post-fire management [11]. The Land Change Monitoring, Assessment, and Projection (LCMAP) project is the latest national-scale effort, using Landsat analysis-ready data to map land cover change on an annual basis from 1986 to 2016 [12,13].

However, there is a disconnect between the availability of remote sensing data upon which landscape monitoring can be based (1972–present) and the stakeholder requirements for long-term time series of land change data. Landsat data only provides consistent 30-m coverage back to 1984. An adequate exploration of the historical impacts of landscape change on human and natural systems often requires a longer time period and reconstruction of historical time periods [14,15], while anticipating and planning for the impacts of future landscape change requires the development of robust future scenarios [16,17]. A growing array of landscape models offer the capability to both reconstruct historical landscapes, as well as model scenarios of future landscape change. For example, landscape modeling is part of the interdisciplinary, integrated assessment modeling (IAM) done in support of Intergovernmental Panel on Climate Change (IPCC) global activities, and these IAMs can provide long-term, consistent, historical-to-future data [18,19]. However, IAM modeling applications are typically very coarse spatially (often 0.5-degree grid cells) and are not directly comparable to US landscape monitoring projects such as NLCD, LCMAP, or LANDFIRE. Radeloff et al. [20] produced a 30-m set of landscape projections for the US, but they have a simplified classification scheme compared to NLCD and a limited temporal span from 2001 to 2051.

To facilitate analyses of the long-term impacts of landscape change and support planning and mitigation efforts related to future change in the US, a long-term time series dataset is needed that provides:

• Long-term historical landscape reconstructions for the years prior to the availability of remote sensing data;
• Similar spatial and thematic resolution to the most widely used land use and land cover (LULC) datasets for the US such as NLCD;
• Parcel-based modeling approach, ensuring the representation of realistic landscape patterns through the use of real land ownership and land management parcels;
• Multiscenario future projections, offering the capability to explore uncertainties related to future landscape trajectories.

The Forecasting Scenarios of land use (FORE-SCE) model has been previously used to model historical “backcasts” and future projections for the conterminous US from 1938 to 2100 [21,22], with a similar thematic resolution to national-scale NLCD products but at a substantially coarser spatial resolution of 250-m. FORE-SCE uses historical datasets and a spatial framework of real land ownership and land management parcels to reconstruct long-term historical landscape change, as well as long-term, scenario-based future projections. The result is an unprecedented long-term time series of landscape data that matches the high spatial and thematic resolution of widely used national-scale landscape data, facilitating seamless integration into existing research and land management workflows. Here, we describe a methodology to produce a long-term (1680 through 2100) LULC dataset for the Delaware River Basin (DRB), with an assessment of the results and next steps towards completing a similar product for the conterminous United States.

2. Materials and Methods

The FORE-SCE model can produce both historical and projected land cover simulations. Study region, time period, and scenarios of a given application are typically dictated by project resources and stakeholder interests. The most recent broad-scale application of the FORE-SCE model provided parcel-based modeling for a number of future climate and socioeconomic scenarios in the US Great Plains from 2014 to 2100 [13,23]. Our methodology here builds from that work, with substantial adaptations to: (1) accommodate the challenge of modeling long-term historical landscape reconstruction at a moderately high spatial and thematic resolution, (2) revamp urban modeling including the first FORE-SCE application with multiple urban density classes, and (3) provide additional improvements in model function and efficiency. The overarching goal of this study was to:

• Provide a long-term 1680 to 2100 database of LULC with consistency throughout the time period, at 10-year intervals;
• Provide spatial and thematic consistency with existing large-scale US LULC datasets such as NLCD and the US Department of Agriculture (USDA) Cropland Data Layer (CDL; [24]), facilitating seamless use of these data in existing land management and research applications;
• Develop a methodology that is extensible and scalable, to support longer-term goals of a national-scale historical-to-future high-resolution landscape database for the US.

2.1. Study Area

The DRB is an active area of focus for US Geological Survey (USGS) research activities, making it an ideal test bed for the development of a long-term landscape time series (Figure 1). The spatial coverage of the study area was defined from the watersheds defined in the Watershed Boundary Dataset component of the USGS National Hydrography Dataset (Hydrologic Unit Codes 020401 and 020402, [25]). The DRB is small yet ecologically diverse and provides drinking water to ~15 million people in the watershed and surrounding areas. Based on the most recent NLCD product (NLCD 2016), ~20% of the region is currently developed, with another ~15% used for agricultural purposes [10]. Societal concerns in the region include the historic impacts of drought, including events such as the 1960s “drought of the century” that resulted in low water flows and saltwater intrusion in the estuary that threatened drinking water supplies for Philadelphia, Pa; New York City, NY; Camden, NJ; and other urban centers [26,27]. The timing of this drought precludes the use of direct Landsat observations for characterizing historical LULC. Historical change in anthropogenic land use has had strong impacts on water resources in the region, while anticipated future land use and climate change are likely to further stress those resources [28]. The long-term landscape time series developed here will be used to explore historical relationships between land use, water use, and drought, as well as facilitate future water resource planning in the region.

2.2. Initial Conditions and Modeling Parameters

Here, we lay the foundation for potential broader-scale application at the national scale by prototyping DRB-specific land cover backcasting datasets to the year 1680 (determined to be “pre-settlement” as described in Section 2.5 below) and forecasting datasets to the year 2100. All land cover modeling was completed at 10-year intervals, at a 30-m spatial resolution, and for 20 thematic land cover categories (

Table 1. Thematic land cover classes modeled for the Delaware River Basin (DRB). Classes represent a majority of National Land Cover Database (NLCD) 2016 classes, with major crop types from the US Department of Agriculture (USDA) Cropland Data Layer (CDL). Change type refers to spatial allocation method enacting landscape change.

Thematic Class	Description	Change Type
Hay/Alfalfa 1	Areas of grasses, legumes, or grass-legume mixtures planted for livestock grazing or the production of seed or hay crops, typically on a perennial cycle. Hay/alfalfa vegetation accounts for greater than 20% of total vegetation.	Parcel
Corn 1	Cropland areas comprised primarily of corn.	Parcel
Soybeans 1	Cropland areas comprised primarily of soybeans.	Parcel
Fallow/Idle 1	Cropland areas comprised primarily of fallow and idle lands.	Parcel
Wheat 1	Cropland areas comprised primarily of wheat.	Parcel
Fruits/Vegetables 1	Cropland areas comprised primarily of fruit and vegetable crops. Represents an aggregate of fruit and vegetable classes from the Cropland Data Layer.	Parcel
Tobacco 1	Cropland areas comprised primarily of tobacco.	Parcel
Other Crops 1	Cropland areas comprised of crop types not represented in classes 1–7 above.	Parcel
Open Water 2	Areas of open water, generally with less than 25% cover of vegetation or soil.	Pre-defined Patch
Developed Low Intensity 2	Representing NLCD classes “Developed-Open Space” and “Developed-Low Intensity”, areas where impervious surfaces account for 0 to 49% of total cover. Includes single family housing units, parks, golf courses, and other recreational urban settings.	Parcel and Random Patch
Developed Med-High Intensity 2	Representing NLCD classes “Developed-Medium Intensity” and “Developed-High Intensity”, areas where impervious surfaces account for 50–100% of total cover. Includes dense single-family housing units, apartment complexes, and commercial/industrial areas.	Parcel and Random Patch
Barren 2	Areas of bedrock, desert pavement, scarps, talus, slides, volcanic material, glacial debris, sand dunes, strip mines, gravel pits, and other accumulations of earthen material. Generally, vegetation accounts for less than 15% of total cover.	Pixel
Deciduous Forest 2	Areas dominated by trees generally greater than 5 m tall, and greater than 20% of total vegetation cover. More than 75% of the tree species shed foliage simultaneously in response to seasonal change.	Pixel
Evergreen Forest 2	Areas dominated by trees generally greater than 5 m tall, and greater than 20% of total vegetation cover. More than 75% of the tree species maintain their leaves all year. Canopy is never without green foliage.	Pixel
Mixed Forest 2	Areas dominated by trees generally greater than 5 m tall, and greater than 20% of total vegetation cover. Neither deciduous nor evergreen species are greater than 75% of total tree cover.	Pixel
Shrubland 2	Areas dominated by shrubs; less than 5 m tall with shrub canopy typically greater than 20% of total vegetation. This class includes true shrubs, young trees in an early successional stage, or trees stunted from environmental conditions.	Pixel
Grassland/Pasture 2	Areas dominated by graminoid or herbaceous vegetation, generally greater than 80% of total vegetation. These areas are not subject to intensive management such as tilling but can be utilized for grazing.	Parcel
Woody Wetlands 2	Areas where forest or shrubland vegetation accounts for greater than 20% of vegetative cover and the soil or substrate is periodically saturated with or covered with water.	Pixel
Herbaceous Wetlands 2	Areas where perennial herbaceous vegetation accounts for greater than 80% of vegetative cover and the soil or substrate is periodically saturated with or covered with water.	Pixel
Perennial Grasses 3	Areas where perennial grasses used for the production of cellulosic-based biofuel represent greater than 80% of total vegetation.	Parcel

¹ Derived from the US Department of Agriculture Cropland Data Layer [24]. ² Derived from the National Land Cover Database [10]. ³ Represents a new land use class form only present in future biofuel scenarios.

). FORE-SCE can run at any temporal time step, but 10-year intervals were selected due to computational constraints imposed by higher temporal resolutions. The 10-year analysis interval is also consistent with recent applications of the FORE-SCE model in the Great Plains [13,23]. A 30-m spatial resolution is a substantial improvement over past national-scale applications of FORE-SCE which were modeled at 250-m [21,22]. The selected 30-m resolution also provides consistency with widely used, national-scale datasets such as NLCD and CDL. Thematic classes represent most classes from the 2016 NLCD [10], with NLCD’s “cultivated cropland” class further differentiated into major crop types in the region from CDL data [24]. To maintain thematic consistency with land cover modeling work underway in adjacent regions 2018-era CDL data were used and the crop category tobacco, though not present in DRB scenarios, remains in the classification schema. The four NLCD classes related to urban development density were collapsed into two classes, a “low-intensity” class and a “medium/high intensity” class. Forecasting was conducted under seven socioeconomic scenarios and three climate projections, resulting in 21 total future landscape projections (Section 2.6.1).

2.3. FORE-SCE Structure

The structure of the FORE-SCE version used here mimics the overall “demand” and “spatial allocation” paradigm from recent FORE-SCE applications [13,23], with substantial model updates as we note below. The primary driver in the FORE-SCE workflow is a temporal series of quantitative change prescriptions for each LULC category in a specific geographic region, known as demand. Spatial allocation algorithms ingest demand at the regional scale and place pixels and parcels of change on the landscape until demand for a given scenario and timestep is met. The same model paradigm exists for both historical landscape reconstruction, and for future projections, with a requirement for regional-level demand (with a region consisting of any geographic framework, such as by state, ecoregion, or watershed), and subsequent spatial allocation to create spatially explicit LULC maps at requisite time intervals. However, given the unique nature of temporal LULC change characteristics when reconstructing historical landscapes versus projecting future landscapes, model function and parameterization differ depending upon temporal trajectory, as noted in Section 2.5 and Section 2.6 below, respectively. The primary objective of FORE-SCE is to represent any scenario of prescribed landscape change with the highest degrees of spatial and thematic resolutions and pattern fidelity as possible.

A key element of the spatial allocation algorithms are suitability surfaces which quantify, for a given pixel, the probability of each LULC class occurrence. Suitability surfaces play a key role in driving where change occurs on the landscape. Suitability surfaces can also be dynamic through time, with past applications updating suitability surfaces in lockstep with projected changes in climate or water availability [13]. The spatial allocation of prescribed change is implemented at either the pixel- or patch-level. Patch-level change refers to land cover transitions that occur in discrete, contiguous groups of pixels, often corresponding with anthropogenic activity. Pixel-level change typically corresponds with natural vegetation land covers and does not employ discretely defined units.

Patch-level change is further divided into two styles: randomly generated patches and predefined ownership-based parcels. Allocation style is defined for each LULC class and may be modified dependent on region and scenario (Table 1). The use of ownership-based parcels ensures realistic representation of landscape pattern changes at fine spatial scales, but the parcel data themselves are unavailable to FORE-SCE data users due to privacy and distribution restrictions. Parcels for the DRB are derived from three sources, as no one parcel dataset provided adequate ownership and management information for all land use sectors:

• USDA Common Land Unit (CLU) data were used as the primary dataset to represent agricultural land management parcels [29];
• Agricultural parcels extracted from Landsat data by Yan and Roy [30] were used to augment CLU data and fill in gaps where CLU data were unavailable;
• Parcels from the Homeland Infrastructure Foundational-Level Data provided ownership parcels for urban and other non-agricultural lands [31].

FORE-SCE accounts for the complete suite of LULC classes found on the landscape. The driver of change for “demand” is primarily assumed to be anthropogenic change, with natural land types responding to changes in anthropogenic land use related to urban development, agriculture, forestry, mining, and reservoir construction. These are also the land use and land cover classes for which both historical data and future scenarios are much more widely available, with information on trends and spatial distribution of natural land covers much more difficult to acquire. FORE-SCE thus maintains a flexible status for land cover classes without explicit prescriptions for change. For example, in a scenario calling for agricultural expansion, natural vegetation categories (e.g., grassland, shrubland, forest) in locations that satisfy suitability requirements for cropland will convert to agricultural categories. For scenarios exploring a contraction of agriculture, natural vegetation classes will increase and flexibly displace crop categories as needed. Flexibility is defined at the class level and dependent on region and scenario. For instance, grassland/pasture was supplied with explicitly prescribed change for historical backcasting but allowed to respond flexibly to other classes for future scenarios. Classes with explicitly prescribed change are further addressed in Section 3.4.

Natural vegetation can also change over time even in the absence of anthropogenic driving forces. Dynamic suitability surfaces drive shifts in the distributions and proportions of natural vegetation classes including grassland, shrubland, and forest classes. Tied to projected changes in climate input into FORE-SCE, pixels of these classes that are characterized as “unsuitable” to support a given cover class under new climate conditions may shift to more suitable classes, as described by Sohl et al. in 2019 [13].

Additional details about the overall FORE-SCE paradigm and structure can be found in prior publications [13,22,23]. What follows are new FORE-SCE developments applied to the DRB (Section 2.4), specific modeling methodologies for both historical landscape reconstruction (Section 2.5) and future projections (Section 2.6) in the DRB.

2.4. FORE-SCE Improvements for Delaware Basin Application

The general FORE-SCE structure and paradigm above was suitable for both historical landscape reconstruction and future scenario modeling in the DRB. However, a number of enhancements and new modeling elements were developed to improve model performance in the basin, while other modifications were made to improve model flexibility and automation, as we aim to apply this framework across the conterminous US.

2.4.1. Urban Module Advances

The FORE-SCE urban module operates similarly to a cellular automata model, with urban growth patterns empirically dependent upon the existing urban pattern and growth elements that can mimic both edge and spontaneous urban center growth. A priority is placed on edge growth from existing urban pixels (including road networks), with urban growth prioritized for existing high-density urban areas, but model parameters also create opportunities for smaller urban areas to coalesce or for development to leapfrog from larger urban areas to smaller, more rural areas. Prior work has demonstrated that this approach produced results more consistent with theoretical, long-term patterns of urban growth [13].

Additional advances in FORE-SCE’s urban module used in the DRB include the following:

• Development of capabilities for the urban change module to operate over descending time steps, supporting historical backcasting;
• Flexibility in parameterization, permitting users to supply spatially explicit likelihood datasets that influence urban change;
• Support of land cover category schemas containing multiple urban classes (past FORE-SCE applications modeled a single urban class).

These modifications facilitate historical backcasting reconstruction and improve the functionality and flexibility of modeling future scenarios. For example, the addition of ancillary data that influence the spatial allocation of change allows for the use of historical data to guide the removal of urban features as the model operates in reverse, replicating actual historical landscape patterns. Support for multiple classes enables transitions between urban intensity levels over time.

2.4.2. Improving Automation and Stochasticity

Suitability surfaces are a key element in converting quantitative demand into spatially explicit land cover realizations [13]. These surfaces are derived by mapping the likelihood-of-presence resulting from stepwise logistic regression using biophysical and anthropogenic covariates and initial land cover conditions (i.e., 2018-CDL), [22] and are utilized by the model to establish “growth zones”, spatial zones where change in a given land cover class is most likely to occur [13]. Prior to this application, the establishment of those class-by-class zones was often an iterative process, requiring interpreter intervention and analysis, potentially introducing subjectivity to the process. The need for a more objective process necessitated the development of new, more automated methodologies for generating suitability surfaces and establishing spatial zones of potential change. This methodology is expected to scale up towards the national scale.

A new methodology was developed that: (1) reduces interpreter intervention, (2) automates establishment of spatial zones of change, and (3) enhances model stochasticity by using initial landscape conditions and weighted random selection to automatically calibrate and establish spatial-growth zones for each class. Figure 2 demonstrates this new process for the land cover category “corn” (Table 1) in a 54 km × 54 km subregion northeast of Pennsylvania (Figure 1, bottom left). Initial land cover class presence/absence and corresponding suitability values (Figure 2a,b, respectively) are used to generate the total operating characteristic (TOC) curve (Figure 2c) for each land cover category [32]. The TOC curve explicitly reveals the quantities of correctly predicted presence (true positives) (TP), correctly predicted absence (true negatives) (TN), incorrectly predicted presence (false positives) (FP), and incorrectly predicted absence (false negatives) (FN) for all possible suitability thresholds of a given category. The maximum and minimum boundaries of the plot define the space of diagnostic potential. Sensitivity (Equation (1)) and specificity (Equation (2)) are computed by dividing TP by the number of actual pixels in the initial landscape with category present (P) and TN by the number of actual pixels with category absent (N).

$$\mathrm{Sensitivity} \left(\mathrm{True} \mathrm{Positive} \mathrm{rate}\right)=\frac{\mathrm{TP}}{\mathrm{P}}$$

(1)

$$\mathrm{Specificity} \left(\mathrm{True} \mathrm{Negative} \mathrm{rate}\right)=\frac{\mathrm{TN}}{\mathrm{N}}$$

(2)

The threshold maximizing the sum of sensitivity and specificity (MSPS) for the category is then computed. In Figure 2d the category’s underlying suitability values for its initial condition of presence is plotted as a histogram along with MSPS. For land covers prescribed to increase on the landscape, an iterative process of random selection centered on MSPS converges on a selection of pixels with a suitability-value distribution that is similar in shape to the category’s starting condition. The selection is from pixels that do not correspond with class presence and therefore represents the spatial extent of potential expansion for the given class. In Figure 2e the histogram for all eligible expansion pixels (those where the given class is absent) is displayed in green (note: pixel counts for very low-suitability values exceed y-axis plotting limits). Pixels for expansion selected from all eligible pixels under the (1) manual percentile-cutoff method and (2) the new automated method are plotted (pink and brown, respectively). Note: suitability values are relative to each land cover category and are not used as raw prediction values when establishing growth zones.

2.5. Historical Reconstruction: 1680–2018

FORE-SCE was used in a “backcasting” mode to iteratively reconstruct landscape change from our starting 2018 date backward to 1680, at 10-year intervals. We selected 1680 to represent “pre-settlement” vegetation in the region, as the first European settlement in the DRB did not occur until 1631, and impact on pre-settlement vegetation was deemed to be negligible prior to 1680 [33,34]. There is also a lack of useful population, agricultural land use, or other historical data sources related to landscape change prior to 1680.

As with past FORE-SCE applications, long-term historical reconstruction is based on quantifying regional demand for modeled timesteps with spatial allocation to produce spatially explicit maps. Regional demand, in turn, is based on the acquisition and harmonization of historical data sources quantifying regional LULC change. With a goal to develop a methodology that was extensible and scalable to other regions of the US, we relied on data sources for historical landscape reconstruction that were consistent and nationally available.

2.5.1. Establishing Pre-Settlement Vegetation

The LANDFIRE Biophysical Settings (BPS) dataset [11] represents vegetation most likely to have been prevalent prior to Euro-American settlement (i.e., 1680) and provides consistent data at our 30-m spatial resolution. We recoded the LANDFIRE BPS dataset for the DRB to match our categorical schema (Table 1) and used the resulting dataset as a pre-settlement target for historical landscape reconstruction. As FORE-SCE models backward in time towards a historical end year before Euro-American influence, anthropogenic landscape features (i.e., cropland, urban, reservoirs) eventually transition to pre-settlement vegetation types. The rate, pattern, and spatial distribution of change, as “current” (2018) landscape conditions revert to pre-settlement conditions, are driven by a variety of historical data sources for the region, as explained in Section 2.5.2 below.

2.5.2. Establishing Historical Demand

Agricultural and urban land uses generally have the most information available historically and serve as the backbone of the historical reconstruction process. These and select other anthropogenic land uses are actively modeled within FORE-SCE, while non-anthropogenic, “natural” vegetative covers such as grassland, shrubland, and forest are passively modeled as described in Section 2.3. Below, historical demands for FORE-SCE applications are described in terms of agricultural and urban demand.

Agricultural Demand

The demand for historical change for agricultural classes was derived using Agricultural Census records [35] and historical population and labor force information [36,37,38]. The Agricultural Census provides county-level estimates for crop production beginning in 1840. Crop acreage estimates, however, were not introduced until 1880. From 1840–1870 Agricultural Census estimates consisted of crop yields. County-level acreages provided by the census were aggregated to match modeled thematic classes and converted to pixel-level prescriptions of change. Interpolation was used to estimate acreage for years in which crops of interest were absent from the census. Extrapolation was used to estimate acreage for census years from 1840–1870. Population and labor force data corresponding with Agricultural Census acreages were used to compute acres-per-crop-laborer. This relationship could then be used to derive acreage using estimated labor population for the years prior to the Agricultural Census (1680–1830).

Urban Demand

Developing the historical demand for urbanization required harmonization of multiple historical data sources, with substantial disagreement on rates of historical urban change. The Historical Settlement Data Compilation for the United States (HISDAC-US) is derived from historical property records and provides temporal (five-year increments) gridded (250-m resolution) settlement data between 1810 and 2015 for the conterminous United States [39]. The HISDAC-US First Built-up Year layer provides the first recorded year a structure was built within a given cell extent. Another dataset on historical urban change is the US Census American Community Survey (ACS) with five-year estimates that include census block group-level median year-built estimates for building structures [40]. While neither of these two historical datasets directly indicate when a given spatial unit would have been classified as “urban”, they can both be used to generate independent estimates for historical urbanization trends, as follows.

Year-built pixel counts were based on all urban land cover pixels present in the model’s initial (2018) land cover and then converted to cumulative urban land cover pixel counts where each year’s cumulative value is the sum of all urban pixels with a year-built less than or equal to the given year. Plotting cumulative urban pixel counts against historical source year-built values (median for ACS and first year for HISDAC-US) reveals curves that are well fit by logistic growth functions and estimate urbanization trends for the region (Figure 3a). ACS cumulative urban pixel counts suggest a more abrupt historical growth rate for the region than do the HISDAC-US counts. This observed trend is likely due to the influence of urban redevelopment within some ACS block groups. Block groups experiencing large amounts of redevelopment of existing urban cover are represented with a temporal shift forward of median year-built and original structure-built dates are lost as redevelopment occurs. Conversely, the HISDAC-US may underestimate the timing of urbanization as “first built-up year” may precede the year at which a pixel would be considered “urban” in our thematic classification system.

Cumulative population [36,37] plotted on the same time interval (using the left vertical axis) reveals a logistic growth curve more similar in shape to HISDAC-US than to ACS according to the comparison of their logistic growth steepness factors (Figure 3a). To further assess both historical data sources, urban areas from the USGS Historical Topographic Map Collection were manually digitized to estimate urban cover for the year 1900 [41]. The estimation of total urban cover for 1900 from historical topographic maps lies between the cumulative growth curve fits for HISDAC-US and ACS (Figure 3b, red square). We also compared HISDAC-US and ACS data with urban land cover quantities as mapped by contemporary remote-sensing-based land cover classifications. Remotely-sensed urban pixel counts were derived from the Land Change Monitoring, Assessment, and Projection (LCMAP) time series for years 1985 and 2000 [12] and from the NLCD for 2016 [10] (Figure 3b).

Ultimately, we selected a harmonization approach to incorporate both HISDAC-US and ACS in the derivation of urban land cover estimates for backcast land cover modeling. A third logistic growth model of cumulative urban land cover quantity was fit using both sets of historical time series data (Figure 3b, green curve). The fit is weighted more heavily to HISDAC-US as that dataset provides a curve more consistent to overall population trends and provides a broader temporal range of data points. The final, harmonized logistic growth curve was used to estimate the total urban land cover area for each backcasting 10-year interval. To separate historical urban change into the low- and high-intensity categories, we assumed a simple, constant rate of change between categories. The rate of conversion between urban types was then solved to arrive at known proportions in the model’s initial land cover conditions, assuming pre-settlement conditions of no urban land cover. The conversion rate was then applied to calculate the amount one urban category is expected to transition to another for any modeled time interval. Demand used by the model to spatially allocate prescribed change for a given category is then simply the difference in category quantity estimates between subsequent backcasting time steps.

Other Demand

Agricultural and urban demand are the two largest anthropogenic land use impacts in the DRB, but demand for other land cover classes was also included. Temporal and spatially explicit information on wetlands is lacking in the DRB, but coarser-scale information provided guidance for aggregate demand across the watershed. We used drained land information from the US Census of Agriculture (1920–1960) [42,43,44,45] and a US Fish and Wildlife Service report to Congress on wetland loss (1780–1980) [46] to characterize historical change. With LANDFIRE BPS serving as the pre-settlement historical endpoint, wetland quantities for historical reconstruction years without reported estimates could be interpolated.

The establishment of reservoirs over the historical 1680 to 2018 period also served as a form of “demand”, with the need to remove reservoirs as the model iterated backward in time. A separate methodology for reservoir removal (and subsequent replacement with other land covers) was established in the spatial allocation element of FORE-SCE, as described in Section 2.5.3 below. Natural land cover classes (e.g., forest, shrubland, and grassland classes) are not provided an explicitly prescribed demand within FORE-SCE and instead are treated as “passively” modeled land covers that respond to prescribed changes in anthropogenic land use activity, as described in Section 2.5.3 below.

Table A1. Datasets used in FORE-SCE’s historical reconstruction.

Dataset	Description	Time Period	Resolution
Historical Settlement Data Compilation for the United States (HISDAC-US) *	Historical gridded settlement data derived from property records	1810–2015	250 m
American Community Survey (ACS) *	Supplement to US Census; includes median year-built estimates for structures	1939–present	block group
Historical agricultural land use time series *	Time series maps represent history of agricultural development	1850–1997	county
LANDFIRE Biophysical Settings (BPS)	Gridded data represent the vegetation system that may have been dominant on the landscape prior to Euro-American settlement	variable	30-m
United States Census of Agriculture **	Historical census of agricultural land usage in the US	1840–present	county
United States Census of Agriculture: Drained Lands	Survey of lands drained for agricultural use in the US	1920–1969	county
United States Census: Population **	Historical census of population in the US	1790–present	county
United States Census: Labor Force	Historical survey of people employed in the agricultural and non-agricultural sectors	1820–1940	national
Population Census prior to formation of United States	Census population estimates of colonial states	1600–1780	state
Wetlands Losses in the United States	US Fish and Wildlife Service report to Congress on the status of wetlands resources in the US	1780–1980	state

* Accuracy reported. ** Accuracy reported for some years.

of Appendix A provides a summary of the historical datasets used in the application, including time period, resolution, and internal accuracy assessment.

2.5.3. Historical Spatial Allocation

Data to support historical backcasting are often not abundant and are generally spatially coarse, but our spatial allocation process is able to incorporate that information to produce spatially explicit land cover maps consistent with the starting 2018 land cover data. For historical backcasting, the landscape is altered for 10-year intervals from 2018 backward to 1680, with a target endpoint of “pre-settlement” vegetation as defined by the LANDFIRE BPS data. Demand for anthropogenic land covers was developed based on multiple historical data sources, with land covers such as the two urban classes and eight agricultural classes eventually disappearing from the landscape as we meet pre-settlement conditions (Section 2.5.2). Natural land cover classes are “passively” modeled, replacing “lost” anthropogenic land cover classes as the model iterates backward in time, eventually converging on pre-settlement vegetation (Section 2.3).

Suitability surfaces used to guide the spatial allocation of change can be dynamic over time (Section 2.3). While suitability surfaces for future scenario-based projections have typically been parameterized to respond to corresponding projected changes in climate or water availability, suitability surfaces for historical backcasting were parameterized to effectively calibrate to existing historical data. For example, the finest-scale long-term historical data available for agricultural classes were county-level Agricultural Census data. Those data were used to not only construct “demand” (Section 2.5.2) but also to modify our starting point, contemporary suitability surfaces to create historically dynamic suitability surfaces for agricultural classes. The modified, 1850–1997 Agricultural Census data developed by Waisanen and Bliss [47] were used to modify historical suitability surfaces for agricultural classes, improving the ability of the model to match historical distributions of individual crop types. The time series of agricultural development provides county-level proportions of harmonized total “improved farmland”. A contemporary land-in-agriculture ratio was computed and included among covariates in the initial regression model. Following the same process used to derive dynamic suitability surfaces for forecasting scenarios described by Sohl et al. [48], dynamic backcasting suitability surfaces were derived by substituting the historical ratios for the contemporary ratios before mapping the regression parameters over the covariates to produce spatially explicit suitability values.

Similarly, coarse-scale historical data can be used to help guide the spatial allocation of urban lands as the model iterates backward in time. The HISDAC-US Built-up Intensity (BUI) product provides indoor-building gross-area sums (square meters) per grid cell (250 m) at five-year increments from 1810 to 2015 [39]. For backcasting the years in which HISDAC-US data were available, normalized moving-window sums of historical BUI were used to weigh the selection of pixels most likely to experience transition. This increases the likelihood that modeled urban decline is most likely to occur where historical BUI was lowest.

Another anthropogenic land use accounted for in the historical modeling is reservoir construction. We developed a semiautomated process that utilizes the National Inventory of Dams [49] and extracts reservoirs from the initial land cover, labeling them with the construction year of the associated dam. We added to the FORE-SCE model the capability to use labeled reservoir datasets as change patches to transition reservoirs to dominant local natural vegetation. We employed the new capability to remove reservoirs from the landscape during historical backcasting, replacing “water” land cover associated with a reservoir with suitable upland land cover classes at the first 10-year interval prior to the recorded dam construction date. The reservoir replacement class is determined by the majority natural cover within a variable buffer around the feature (we used 90 m for the DRB). The methodology described here was used for historical landscape reconstruction but could also support forecasting scenarios that include reservoir removal.

2.6. Future Projections: 2018–2100

Future projections are more straightforward than historical reconstruction, with the basic FORE-SCE structure (Section 2.3) and the recent model improvements (Section 2.4) capable of addressing modeling needs in the region. To facilitate long-term analyses related to the impacts of landscape change, model runs were done from current (2018) to 2100, with nominal 10-year intervals that match the historical backcasting intervals. A scenario-based approach was used to develop and visualize a variety of future landscape futures, facilitating the exploration of uncertainties related to future landscape change.

2.6.1. Scenario-Based Demand

As with historical backcasting, FORE-SCE requires quantitative demand for each modeled land cover class for each future modeled time interval. Several different approaches were used to construct demand for past FORE-SCE applications including: (1) straightforward “business-as-usual” (BAU) scenarios extrapolated recent trends (typically derived from remote-sensing-based LULC databases such as NLCD), and (2) downscaling of existing, coarse resolution scenarios, or custom econometric scenarios. Demand may also be constructed in partnership with regional stakeholders by exploring “what-if” scenarios that facilitate visualization of landscape outcomes for various policies.

DRB scenario development followed a similar process to that employed for recent model applications of FORE-SCE in the Great Plains [13,23].

Table 2. Forecasting scenarios developed for the DRB.

Forecasting Scenario	Acronym	Description
Business-as-Usual	BAU	Extrapolation of existing trends in 2001–2016 NLCD and 2008–2018 CDL.
Billion Ton Update, $40 Farmgate	BTU $40	US Department of Energy (DOE) scenario for national biomass production >= 1 billion tons assuming farmgate price of USD 40/ton for dry biomass.
Billion Ton Update, $60 Farmgate	BTU $60	US DOE scenario for national biomass production >= 1 billion tons assuming farmgate price of USD 60/ton for dry biomass.
Billion Ton Update, $80 Farmgate	BTU $80	US DOE scenario for national biomass production >= 1 billion tons assuming farmgate price of USD 80/ton for dry biomass.
Global Change Assessment, no mitigation	GCAM Ref	Global Change and Assessment Model land cover projections assuming no climate change mitigation.
Global Change Assessment, mitigation pathway 2.6	GCAM 2.6	Global Change and Assessment Model land cover projections assuming mitigation in line with Representative Concentration Pathway (RCP) 2.6 framework.
Global Change Assessment, mitigation pathway 4.5	GCAM 4.5	Global Change and Assessment Model land cover projections assuming mitigation in line with RCP 4.5 framework.

describes the seven individual scenarios that were modeled, based on three different sources. The BAU scenario extends recent trends from two different sources, with trends in individual crop categories based on recent trends between the 2008 and 2018 dates of the USDA CDL data, with all other land cover classes based on extrapolated trends between the 2001 and 2016 NLCD dates. Three scenarios were based on the US Department of Energy Billion Ton Update (BTU) [50], which provides a county-level exploration of the impacts of various economic scenarios on the agricultural landscape. The BTU scenarios are based on variations in modeled farmgate prices for biomass for biofuel production. The final three scenarios were derived from the Global Change Assessment Model (GCAM), with GCAM annual land cover percentages at 0.5-degree resolution downscaled to finer resolutions to support regional and national applications in the United States [51]. West et al. provide three GCAM scenarios to explore a range of human effort to mitigate climate change—a reference scenario provides landscape conditions with the assumption of no targeted climate-change mitigation efforts, while two additional scenarios assume mitigation efforts necessary to achieve total radiative forcing levels of 2.6 and 4.5 W/m² by the year 2100, respectively. Sohl et al. [13] provide additional information on development of BAU, BTU, and GCAM scenarios.

2.6.2. Projected Spatial Allocation

Spatial allocation of the scenarios described above largely follows the paradigm used for the recent Great Plains application [13] and is described in Section 2.3, with the addition of new modeling elements described in Section 2.4. The new urban model was applied for the two urban density classes described in Table 1. As with past applications [13,23,48], projected climate information was incorporated into future projection modeling. Downscaled temperature and precipitation projections [52] averaged from multiple Coupled Model Intercomparison Project Phase 5 (CMIP5) models (

Table A2. Coupled Model Intercomparison Project Phase 5 (CMIP5) models providing unrestricted data used to compute projected climate variables from downscaled, gridded, monthly temperature, and precipitation projections.

Modeling Center	Model	Institution
CSIRO-BOM	ACCESS1.0	Commonwealth Scientific and Industrial Research Organisation, Australia and Bureau of Meteorology, Australia
CSIRO-QCCCE	CSIRO-Mk3.6.0	Commonwealth Scientific and Industrial Research Organisation, Australia and Queensland Climate Change Centre of Excellence, Australia
INM	INM-CM4	Institute for Numerical Mathematics
IPSL	IPSL-CM5B-LR	Institut Pierre-Simon Laplace
MIROC	MIROC-ESM	Japan Agency for Marine-Earth Science and Technology, Atmosphere and Ocean Research Institute (The University of Tokyo), and National Institute for Environmental Studies

) were used to generate dynamic suitability surfaces sensitive to multiple projected climate change pathways. Pixel-level suitability for a given LULC class thus evolves over time, with LULC classes expanding, contracting, or shifting in space in response to projected changes in climate.

2.7. Model Assessment

FORE-SCE model performance was evaluated at multiple levels:

• Inter-scenario model comparison—We examine differences between modeled land cover over the same time intervals under different scenarios.
• Intra-scenario model performance—We assess FORE-SCE’s internal abilities to represent prescribed change (demand) in a modeled land cover realization for a given future scenario and to spatially match observed change between two reference LULC datasets, NLCD 2001 and NLCD 2016.
• Historical reference—We compare modeled land cover with spatially-explicit historical references.
• Cross-model scenario comparison—We compare FORE-SCE’s simulation of land cover change to the simulations produced by a suite of global IAMs.

As with past applications, the concepts of quantity and allocation disagreement were used to evaluate model results [53]. Quantity disagreement is the amount of difference between model results and a reference dataset arising from variances in overall proportions of landscape categories. Allocation disagreement results from an imperfect match in the spatial arrangement of landscape elements. Allocation disagreement can be further divided into the categories exchange and shift. Exchange spatial disagreement is the difference between two categories in which one location corresponds with the reverse in another location. Shift disagreement entails spatial differences between more than two categories [54].

To track the model’s ability to represent prescribed change (demand) on the landscape, we computed the absolute deviation (both positive and negative) of modeled change from the prescribed change at every timestep, for every land cover category that has explicitly supplied demand. This procedure was conducted for all scenarios and summarized across modeled time intervals.

Applying the FORE-SCE model to historical land cover reconstruction requires additional parameterization and the flexibility to use what is available from limited historical reference data. As described above, we made use of historical datasets [39,47] to guide the spatial allocation of agricultural and urban land cover categories. To assess the impact of incorporating these historical datasets on resulting patterns we employed tandem model runs. These runs are identically parameterized other than the historical parameter of interest and serve to isolate the effect of the varying parameter. We quantify and examine the effects below.

Hurtt et al. [19] have produced Land Use Harmonization (LUH2) datasets to smoothly connect historical data from the History of the Global Environment database (850–2015) and multiple alternative future scenarios represented by IAMs (2015–2100). We compare the trends of aggregated categories modeled by FORE-SCE for historical reconstruction and multiple future scenarios to their counterparts as represented by LUH2.

3. Results

The overarching goal was to produce consistent, long-term, historical, current, and future representations of landscape change, with one historical landscape reconstruction that strove to match reference data as closely as possible, and multiple future scenarios to represent the considerable uncertainties.

3.1. Modeled Trends for Aggregated Categories

Figure 4 depicts temporal trends of modeled results for the DRB from 1680 through 2100, depicting trends in major, aggregated land cover classes, including the two biggest anthropogenic driving forces of change in the region (agriculture and urban), and the largest “passively modeled” natural vegetation class (forest). Agricultural land covers increased steadily for the first two centuries from the starting 1680 pre-settlement land cover, plateauing in the latter half of the 19th century and peaking near the turn of the century, before declining in extent from about 1920 onwards. As agriculture started to decline, the rate of urbanization dramatically rose, rising steadily until the present day. Forested lands decreased sharply in the 19th century in response to agricultural expansion, with reforestation occurring until about 1970 as marginal agricultural lands were then abandoned. However, forests again began to lose area after 1970, as agricultural abandonment abated and urbanization continued a strong increase. The U-shaped pattern of forest decline and subsequent rebound from 1680 to 1970 is consistent with forest transition theory [55]. The renewed declines in forest area after the 1970s fit the pattern of a “contemporary regressive stage” of forest transition found for the Eastern US as a whole during this time period, driven by land use intensification, urban expansion, and diminution in processes of agricultural abandonment and reforestation [56]. The aggregated categories simplify the visualization of trends and together comprise the majority of the DRB’s initial landscape (83%).

Projected land cover trends from 2018 to 2100 show substantial variability among scenarios, yet none of the scenarios anticipate volatility approaching that of 20th-century landscape change in the DRB. The BAU scenario is the most dynamic, with extrapolations of recent trends resulting in offsetting trends of strong declines in agricultural land, and continued expansion of urbanized lands, with relatively stable forested lands. Other scenarios show similar trends in agriculture and urban but with lower rates of urban growth, with the exception of the GCAM Reference scenario where agricultural lands increase. Overall, as described in Sohl et al. [57], there tends to be more similarity within a scenario family than between them. For example, changes in agriculture for the three GCAM scenarios are all clustered more closely to “no-change”, while the three BTU scenarios all show higher agricultural loss than any GCAM scenario. Similarly, BTU scenarios tend to show higher forest increases than most of the GCAM scenarios. However, across all scenarios, substantial variability exists in modeled landscape trends (as described in the next section), facilitating analyses of the uncertainties associated with landscape change impacts on human and natural systems.

3.2. Examining Scenarios in a Subregion

Figure 5 displays initial land cover conditions, all forecast scenarios at the year 2100, and three backcasting dates for the same subregion in Figure 1. In each image, we observe the interplay of agricultural and urban anthropogenic classes with the dominant natural vegetation, forest classes. Historical landscape reconstructions spatially depict the patterns described in Section 3.1. The historical landscape at “pre-settlement” in 1680 is dominated by deciduous and mixed forest, while a substantial proportion of forest was cut by the time the region as a whole reached peak agricultural extent near the start of the 20th century. From 1900 to 1970, marginal agricultural lands throughout much of the region were abandoned and allowed to regenerate back into forests, leaving much smaller pockets of higher-value cropland within the growing forest mosaic. Urban lands increased dramatically during the 20th century as shown in the northern part of the image.

Most forecasting scenarios generally agree on trends in urban expansion and agricultural decline. The displacement of agriculture by forest is most prevalent in the BTU $60 and $80 scenarios which, despite prescribing modest growth for perennial grass (associated with assumptions on cellulosic biofuel production in those scenarios), forecast general agricultural decline. As noted above, the GCAM Ref scenario stands alone in its prescription for agricultural gains, as it is a scenario without any assumptions on climate mitigation and thus less likely to convert marginal agricultural lands to forests. The BAU scenario maintains urban expansion at rates which the other forecasting scenarios taper.

3.3. Assessing Inter-Scenario Variability with Difference Metrics

Figure 6 shows the components of quantity difference and allocation difference (including “exchange” and “shift” components) for every future scenario combination across the forecast period. We expect to see quantity difference between scenarios as, by definition, a scenario contains prescriptions of quantitative change that differentiate it from other scenarios. We also expect to see allocation difference. Quantity difference without any allocation difference between different scenarios would indicate that the modeled projections result from a completely deterministic, non-stochastic process of spatial allocation. Whereas, high levels of allocation difference in relation to prescribed change may suggest that the model is poorly calibrated or highly random.

Not surprisingly, we observe consistent, increasing trends in difference for all pairs of scenarios as they are modeled into the future. The difference present by the year 2100 is largely cumulative given that most changing land cover categories trend in one direction under forecast scenarios. Figure 7a charts all components of difference for the final year of projection (2100), when scenarios have fully diverged, and the difference is at its highest. Figure 7b isolates the spatial components of difference described above and Figure 7c separates allocation difference for the final forecast year into its land cover category components. Note: considering allocation difference separated by category results in each pixel in the difference being counted twice—once for each different category in the pair—the actual total allocation difference across the region is thereby half of the allocation difference summed across categories. The bars in Figure 7c are thereby twice the height of the actual allocation difference bars in Figure 7b.

For the projected year 2100, the total allocation difference across all category combinations ranges from 5.4% to 8.5% of all the region pixels (Figure 7b). On average, “exchange” represents a majority of allocation difference when total allocation difference is lower, while “shift” on average grows as a proportion of total allocation difference when total difference itself increases. The land cover categories contributing most to allocation difference are both low- and high-intensity developed classes and deciduous forest (Figure 7c, pink, red, and green, respectively). This is consistent across all scenario comparisons and reflects the high degree of potential spatial variability present in the urban-to-non-urban interface.

Figure 7d shows total absolute demand (both prescribed gain and loss) separated by land cover category and plotted by scenario as a proportion of the region. BAU and GCAM 4.5 prescribed the highest amounts of change, respectively. The contribution of agriculture land cover categories to allocation differences is relatively minor in comparison to the contribution of urban categories when prescribed absolute change is considered. As described above, both agriculture and urban change are parcel-based and guided by suitability surfaces. However, the urban module enacts change (be it growth or decline) at the urban-to-non-urban interface and here finds more opportunities for divergence between scenarios than does agricultural change. This results in distinct urban patterns between scenarios with a less overall agreement.

3.4. Modeled vs. Prescribed Change (Demand)

3.4.1. Backcasting Assessment

The absolute deviation of modeled quantities from prescribed demand summed over all years of historical backcasting is approximately 0.94% (

Table 3. Gross prescribed change (both positive and negative) and absolute deviation from prescribed change aggregated over all backcast years, separated into categories with explicitly supplied demand.

	Backcasting
Total Gross Demand (ha)	3,579,796
Total Gross Deviation (%)	0.94
Hay/Alfalfa	0.16
Corn	0.05
Soybeans	0.00
Fallow/Idle	0.04
Wheat	0.07
Fruits/Vegetables	0.08
Other Crops	0.13
Grassland/Pasture	0.39
Herbaceous Wetlands	<0.01
Perennial Grass	<0.01

). In Figure 8, each land cover category’s modeled outcomes are plotted against the corresponding demand for each year of backcasting. The passive, flexible land cover categories have no explicit demand and instead respond to prescribed demand in the anthropogenic land cover classes. Demand for these classes is thus plotted as a horizontal line representing the initial quantity.

The region’s historical trends in agriculture see a peak near the year 1900. Passive, natural vegetation classes experience a corresponding trough in quantity but grow to dominate the landscape by pre-settlement. Urban categories are reduced to zero by pre-settlement. Notably, prescribed demand for urban categories representing low- and high-intensity development does not include urban land cover representing roads. The deviation between prescribed demand and modeled change that is apparent in the plot for urban/developed categories is due to the aforementioned removal of roads, which the model conducts under prescriptions of urban decline but which is not explicitly represented in the prescribed demand. For this reason, urban categories are omitted in the overall assessment of deviation from demand (

Table 4. Gross prescribed change (both positive and negative) and absolute deviation from prescribed change aggregated over all forecast years, separated into categories with explicitly supplied demand.

	BAU	BTU $40	BTU $60	BTU $80	GCAM Ref	GCAM 2.6	GCAM 4.5
Total Gross Demand (ha)	419,457	249,332	345,830	394,145	222,392	263,244	183,015
Total Gross Deviation (%)	0.16	0.56	0.36	0.26	1.12	0.55	0.61
Hay/Alfalfa	0.01	0.02	0.02	0.01	0.35	0.03	0.03
Corn	0.02	0.03	0.02	0.02	0.18	0.04	0.04
Soybeans	0.03	0.03	0.05	0.02	0.12	0.08	0.03
Fallow/Idle	<0.01	0.11	0.02	0.01	0.06	<0.01	0.01
Wheat	0.01	0.08	0.02	0.01	0.04	0.04	0.04
Fruits/Vegetables	0.02	0.08	0.01	0.01	0.06	0.03	0.05
Other Crops	0.01	0.04	0.01	<0.01	0.04	0.01	0.07
Low Density Developed	0.01	0.03	0.01	0.01	0.01	0.01	0.01
High Density Developed	<0.01	0.01	0.01	<0.01	0.01	0.02	0.03
Woody Wetlands	0.01	0.01	<0.01	0.02	0.01	<0.01	<0.01
Herbaceous Wetlands	0.05	0.01	<0.01	0.01	0.02	0.01	<0.01
Perennial Grass	<0.01	0.11	0.20	0.15	0.24	0.28	0.32

).

Woody Wetlands are present in the pre-settlement natural vegetation dataset based on LANDFIRE BPS. However, Herbaceous Wetlands are not. Woody Wetlands were left flexible and allowed to displace anthropogenic classes as necessary, whereas prescribed change was supplied for Herbaceous Wetlands to reach pre-settlement levels.

3.4.2. Forecasting Assessment

Table 5. Agreement and components of disagreement between reference NLCD 2016 and FORE-SCE modeled output for the same year and land cover categories for the DRB.

Measure	Percent of Region
Agreement	95.00
Quantity	0.03
Exchange	2.52
Shift	2.45

displays deviation from demand for the seven forecasting scenarios under RCP-average climate assumptions. Overall, for six of the scenarios, the total modeled change deviates from prescribed demand by less than 1%. The GCAM Ref scenario stands out with ~1.12% deviation from the prescribed demand. Notably, this scenario is the only forecasting scenario prescribing net gains for aggregated crop classes (Figure 4). Overall, FORE-SCE excels at matching prescribed demand, and thus quantity difference for prescribed scenario land cover.

The FORE-SCE model can use pre-processed ownership parcels to represent anthropogenic change on the landscape, which results in more realistic patterns at higher spatial resolutions [13]. However, when a category expands into a new area, it may do so via parcels that have not historically been of the expanding land cover type and struggle to replicate historical patterns. We observe that the categories contributing most to deviation from demand for the GCAM Ref scenario are agriculture categories that miss prescribed growth (both overshooting and undershooting) due to reliance on coarse, overly-large ownership parcels. This is specifically observable for the Perennial Grass category across BTU and GCAM scenarios, which all call for growth and primarily see expansion among large ownership parcels, as these are the only “suitable” parcels for the class to occupy. FORE-SCE currently does not account for the potential land ownership parcels to evolve over time (e.g., subdivision of a large agricultural lot as urban centers expand). Automating the development of ownership boundaries that vary in size and shape under different scenarios is an active area of research for the project.

3.5. Modeled vs. Contemporary Reference

Using NLCD 2001 as the starting point, we modeled land cover quantity changes from 2001 to 2016 in a single 15-year time interval for the DRB. Observed quantity changes for all categories between 2001 and 2016 in the NLCD datasets served as FORE-SCE’s demand input. The model’s simulated 2016 land cover dataset was then compared to the reference NLCD 2016 land cover to assess FORE-SCE’s ability to maintain spatial allocation integrity between reference points inside the remote sensing period of record.

FORE-SCE’s simulation of 2016 NLCD change finds total pixel-level agreement/disagreement at 95.00%/5.00% against the 2016 NLCD reference dataset (Table 5). As observed in Section 3.4.1 and Section 3.4.2 the model closely matches prescribed changes at the category quantity level. Here quantity represents 0.65% of all disagreement. Allocation disagreement represents almost all the difference between the modeled and reference 2016 datasets and is split approximately evenly between the spatial components exchange and shift. We consider this result an appropriate reflection of the stochasticity inherent in the interface between anthropogenic and natural land covers and validation of FORE-SCE’s ability to model with high spatial fidelity in the region.

3.6. Modeled vs. Historical References

3.6.1. Modeled vs. Historical Agriculture

Section 3.4.1 and Section 3.4.2 demonstrate FORE-SCE’s ability to meet overall demand for a region, but options for assessing the model’s ability to accurately represent historical landscape patterns at finer scales are limited, given the lack of spatially explicit, historical land cover data (the reason this modeling work is necessary). We can do some comparisons with our finest-scale available data, which for many historical datasets is county-based. We modified suitability surfaces using historical county-level agricultural time series (described above) with the goal of improving the spatial allocation of agricultural classes in backcasting. Below we compare the resulting county-level modeled area-in-agriculture to the historical time series dataset for two backcasting model runs: model run A without historically dynamic suitability surfaces and model run B with historically dynamic surfaces.

For comparison purposes, the modeled crop categories were aggregated to match the Agricultural Census data’s generalized class “improved farmland” [47]. The absolute difference in area (10³ hectares) between modeled and historical land in agriculture was computed for the counties with the majority area inside the study region (Figure 9a) and then summed across counties for each year of comparison (Figure 9b). Using the absolute difference of modeled vs. historical avoids the effect of counties with positive differences canceling those with negative differences when all county-level differences are summed across the region.

We observe that the incorporation of historically dynamic suitability surfaces reduces the amount of difference between the historical dataset and modeled results. Taking the year 1880 as an example, we see that the absolute difference between modeled and historical cropland when aggregated at the county level is approximately 507,000 ha—about 34% of all historical cropland for that year—when the model does not incorporate historically dynamic suitability surfaces (model A). When historically dynamic suitability surfaces are incorporated to guide historical spatial change (model B) the absolute difference at the year 1880 across counties falls to approximately 320,000 ha (~22% of all historical cropland).

In Figure 9c, we separate disagreement with the historical reference by county. Most counties show very modest differences between modeled and reference data, with an expected but modest increase in disagreement as the model iterates backward in time. In Figure 9d–f, we examine the three counties contributing the most to total disagreement with the historical reference, where the overall trend in agriculture is modeled correctly, but disparities exist in the quantity of change. The common feature of each of these three counties is dramatic patterns of agricultural change throughout the historical period. In the adjacent counties of Wayne County, Pennsylvania, and Delaware County, New York, agriculture peaked around 1900, and then sharply declined as a majority of marginal agricultural land in the region was abandoned and allowed to regenerate into forests. In Montgomery County, Pennsylvania, agricultural extent peaked earlier in the 1800s, and then continuously declined to very low levels today as urbanization supplanted most agricultural land. Even with the dynamically altered suitability surfaces, the model can have trouble reconstructing landscapes where conditions defining suitability for agriculture are dramatically different than the present day.

In combination with the results shown in Section 3.4.1 and Section 3.4.2, we thus show (1) FORE-SCE can very precisely match prescribed total proportions of landscape change at the aggregate (DRB) scale, (2) the use of dynamic suitability surfaces dramatically improves the spatial allocation of change within the DRB, and (3) at the county level, representation of agricultural patterns is very good.

3.6.2. Modeled vs. Historical Urban

To assess the impact of historical HISDAC-US BUI [39] in the urban change module (described in Section 2.5.3) we examined the relationship between BUI and modeled urban over time for two backcasting model runs—model run C did not incorporate historical BUI whereas model run D did. To simplify the assessment, we aggregated the two modeled urban intensity categories.

Figure 9a–c demonstrates the pattern differences resulting from incorporation/non-incorporation of available historic BUI in the urban module for a 6400 km² area centered on Philadelphia, Pennsylvania. To quantify the differences between models C and D the percent of total historical BUI for the region that is coincident with the modeled urban footprint was calculated for each year of backcasting. Incorporating historical BUI in the urban module results in modeled urban land cover that captures a consistently higher percentage of historical BUI for the region (Figure 10d).

The percent of total BUI on the landscape captured by the modeled urban footprint drops notably for the years prior to 1925. Initially, this appears to indicate model shortcomings. However, further analysis reveals this is an appropriate outcome, reflecting historical population dynamics that increasingly resist discrete urban categorization as we move back in time.

Figure 10e displays a time series of Global Moran’s I computed using all BUI values for the region. Here we set spatial weights for non-zero BUI pixels to include the four nearest neighbors using a variable-bandwidth Gaussian kernel function [11]. Global Moran’s I ranges from −1 (negative autocorrelation) to +1 (positive autocorrelation) [58]. Though the index remains positive (indicating clustering), it declines moving back in time with a notable drop near 1925, corresponding with the drop in the modeled urban’s capture of BUI.

US Census data for the same time period and region fleshes out the narrative, revealing a steady trend in shifting population proportions from rural to urban from 1840 onward (Figure 10f) [36,37]. During this period, the region becomes increasingly urbanized with much of the labor force transitioning from farming to industry [59,60]. Moving forward in time, as populations shift from rural to urban areas, built-up intensities concentrate in densities sufficient to be categorically represented in land cover classification. However, moving backward in time, we see BUI become more evenly distributed as a variable with a spatial expression that is increasingly rural and less likely to be captured by the urban land cover footprint. For instance, in 1950, 80% of the total BUI for the region is spatially concentrated in the top 14% of all non-zero BUI pixels, which correspond with urban centers. In 1810, 62% of the highest-value BUI pixels are required to reach 80% of cumulative BUI for the region, indicating a significant lowering in concentrations and a more dispersed population. This historical trend results in lower densities and higher fragmentation of built-up intensity and makes discrete urban cover increasingly unlikely in the landscape as the model moves back in time. The drop-off in the modeled urban capture of historical BUI demonstrated in Figure 10d is the result, as build-up intensities fall below levels of discrete representation.

3.6.3. Modeled Pre-Settlement (1680) vs. LANDFIRE BPS

Anthropogenic landscape features with historical change prescriptions (cropland, urban, reservoirs) revert to pre-settlement defined by LANDFIRE BPS during modeled backcasting (Section 2.5.1).

Table 6. Modeled land cover at 1680 vs. LANDFIRE BPS pre-settlement dataset percent-of-region.

Category Name	BPS	Modeled 1680
Hay/Alfalfa	0.00%	0.00%
Corn	0.00%	0.00%
Soybeans	0.00%	0.00%
Fallow/Idle	0.00%	0.00%
Wheat	0.00%	0.00%
Fruits/Vegetables	0.00%	0.00%
Tobacco	0.00%	0.00%
Other Crops	0.00%	0.00%
Open Water	2.90%	1.96%
Developed Low Intensity	0.00%	0.00%
Developed High Intensity	0.00%	0.00%
Barren	0.78%	0.81%
Deciduous Forest	58.65%	65.65%
Evergreen Forest	3.77%	3.61%
Mixed Forest	22.31%	15.02%
Shrubland	0.06%	0.91%
Grassland/Pasture	0.61%	0.27%
Woody Wetlands	10.92%	8.22%
Herbaceous Wetlands	0.00%	3.54%
Perennial Grass	0.00%	0.00%

compares land cover category histograms for the modeled result at 1680 and the BPS dataset. Pontius difference metrics [54] for the two datasets find total agreement at 94.68%. Disagreement is comprised of 3.27% quantity, 1.45% exchange, and 0.59% shift. Overall, FORE-SCE’s end landscape reconstruction for 1680 very closely matches the LANDFIRE BPS data. Some of the differences can be attributed to the BPS dataset containing no counterpart to the historically modeled category Herbaceous Wetlands. BPS also represents Open Water reservoirs that would not have existed prior to settlement. Finally, FORE-SCE does not currently attempt to convert contemporary non-anthropogenic, natural vegetation categories to pre-settlement types (as represented in the BPS data) given the dearth of historical information on such transitions.

3.7. Cross-Model Scenario Comparison

In Figure 11 we compare trends for the region’s majority categories (described in Section 3.1) between FORE-SCE and LUH2 over multiple scenarios. The y-axes represent the percent of the DRB occupied by a given category. We observe considerable agreement for historical trends, a validation of our historical demand methodology. Excepting the LUH2 scenario RCP3.0 SSP3, variation among projected scenarios is modest. The future scenarios selected for FORE-SCE generally trend upward in quantity for forest classes and downward for agriculture whereas the opposite generally occurs across the LUH2 scenarios with forest decreasing and agriculture increasing. We do not attempt to validate one scenario against another or produce a “best” scenario given that the value of the scenario approach lies in plausible diversity [61]. We do note that the level of cross-model agreement we observe here is typically not the norm across models and scenarios [57].

4. Discussion

We have developed a methodology capable of generating long-term temporal time series of LULC, with spatial and thematic consistency with the most widely used, national-scale LULC databases for the US. The pilot application for the DRB demonstrated a capability to accurately reconstruct historical landscape conditions based on spatially coarse or non-spatial historical data and quickly generate a wide range of future LULC scenarios (Figure 4) that facilitate the exploration of uncertain future landscape conditions on a host of socioeconomic and biophysical processes.

Data harmonization is a major challenge for both historical landscape reconstruction and for future projections. For historical reconstruction, inconsistencies among multiple historical datasets are frequent, leading to challenging decisions on how to harmonize data differences. The differences between urban extent and trends from HISDAC-US, ACS, population data, LCMAP, and NLCD noted in Section 2.5.2 are typical, while similar challenges exist in attempting to harmonize data sources such as agricultural census data and LULC derived from remotely sensed imagery. For future projections, our use of existing scenario frameworks such as the Billion Ton Report or IPCC scenarios is confounded by initial starting conditions (i.e., starting LULC proportions) that are inconsistent with our starting LULC (based on NLCD or CDL). Unless there is one, universally accepted “right” dataset to represent landscape change, either historically or for future scenarios, data harmonization will continue to be an issue for FORE-SCE applications.

We strove to develop a methodology that was extensible and scalable, with the capability to produce historical landscape reconstructions and future projections for any geography or scale in the US, including potential national-scale application. We thus relied on data sources that were consistently available for the conterminous US. Given the reliance of FORE-SCE on land ownership and land management parcels, consistently available parcel data remain the biggest data challenge. Privacy concerns and policy limit the availability of data such as USDA’s Common Land Unit data, by far the best parcel data for our agricultural modeling. Land ownership data are increasingly available even at national scales through private groups that have assembled thousands of individually acquired county and local data sources, but access costs are a challenge, particularly for a potential national-scale application. Another long-term need is a module for the dynamic evolution of parcels, with, for example, a capability to sub-divide large agricultural parcels at the urban fringe as cities expand.

While FORE-SCE is fully capable of replicating remote-sensing-based LULC databases such as NLCD or CDL, we are still improving the representation of some processes governing landscape change. For example, we are aiming to quantify interactions between LULC and hydrology, including both water use and surface water/groundwater availability. FORE-SCE/hydrologic model coupling efforts are underway at USGS. These efforts involve loose coupling with input/output data from FORE-SCE and USGS water balance, energy-balance, and hydrologic models. To address these challenges, open-source cloud computing ecosystems are expected to facilitate processing and exploit these linkages amongst land use, water use/availability, and climatology.

For extended timelines, estimating probabilities for different scenarios of change becomes increasingly difficult as well as irrelevant. Scenario planning frameworks benefiting from an array of land cover projections, as we have produced for the DRB, assume a level of irreducible future uncertainty and attempt to establish a range of plausible futures. Policy-making and risk-mitigation strategies based on scenario planning that includes worst-case, boundary projections will be more resilient than those aligned with projected averages or focused on short-term scenario agreement [61]. We are continuing to improve scenario linkages with FORE-SCE, including current efforts to better link with IAMs produced as part of IPCC activities.