Wild ﬁ re activity and land use drove 20th-century changes in forest cover in the Colorado front range

. Recent shifts in global forest area highlight the importance of understanding the causes and consequences of forest change. To examine the in ﬂ uence of several potential drivers of forest cover change, we used supervised classi ﬁ cations of historical (1938 – 1940) and contemporary (2015) aerial imagery covering a 2932-km 2 study area in the northern Front Range (NFR) of Colorado and we linked observed changes in forest cover with abiotic factors, land use, and ﬁ re history. Forest cover in the NFR demonstrated broad-scale changes 1938 – 2015 and overall cover increased 7.8%, but there was notable spatial variability and many sites also experienced Forest Loss. Recent (1978 – 2015) wild ﬁ re was the largest single driver of Forest Loss, with ﬁ res burning 14.3% of the total study area. Recently burned areas showed net losses of 36.9% forest cover. Reasons for Forest Gain were more complex, with elevation, past mining density, ﬁ re history, and topographic heat load index being the strongest predictors of increases in forest cover. Historical mining activity is one of the dominant anthropogenic impacts in ecosystems in the NFR and it had a complex, non-linear relationship with 20th-century changes in forest cover. Subalpine stands originating after stand-replacing ﬁ res circa mid-1800s to early 1900s showed some of the greatest gains in forest cover, indicative of slow and continuous post- ﬁ re recovery through the 20th century. We also investigated factors such as land ownership, road density, forest management activities, and development intensity, which played detectable, but more minor roles in observed change. Twentieth-century changes in forest cover throughout the NFR are a result of ecological disturbances and anthropogenic in ﬂ uences operating at varying timescales and overlaid upon variability in the abiotic environment.


INTRODUCTION
Studies using forest inventories and remotely sensed data have identified declines in global forest area c. 1990s-2010s (Hansen 2013, Keenan et al. 2015, though longer-term trends may differ in direction and magnitude (Song et al. 2018). Forests provide countless ecosystem services to humans, are sanctums of biodiversity, and play a key role in nutrient cycling (Stocker 2013, FAO 2015, and thus, an understanding of the causes of forest change is an important research priority. Forest change is rarely a result of a single factor acting alone. The causes of land cover change vary across complex landscapes and multiple potential drivers should be considered (Black et al. 2003). The attribution of forest cover change to specific natural and anthropogenic factors, at sufficient timescales to observe both abrupt losses and gradual recovery, is likely to improve conservation planning and the understanding of global change.
Wildfire is one of the most important terrestrial disturbances (Bowman et al. 2009) and is a key determinant of the global distribution of forest vegetation (Bond et al. 2005). The effects of wildfire activity on forest cover change are complex because fire is expected to influence forest cover through both initial losses and subsequent recovery. From 1984 to 2015, warming and drying trends contributed to increases in the area burned in large wildfires across much of the western United States (Abatzoglou and Williams 2016). These increases in wildfire activity have occurred in tandem with poor post-fire recovery in many western conifer forests (Stevens-Rumann et al. 2018). Seed tree-and therefore seed availability-is one of the primary drivers of post-fire recovery, particularly in forests composed of seed obligate species without fire-adapted canopy seed banks (Bonnet et al. 2005, Haire and McGarigal 2010, Harvey et al. 2016). Due to dispersal limitations, sites in the interior of large, high-severity patches may take decades or centuries to regain forest cover (Chambers et al. 2016, Baker 2018. Furthermore, a warm and dry climate can inhibit post-fire recovery of forests (Stevens-Rumann et al. 2018). Recent wildfire activity, combined with poor postfire recovery, may be leading to forest cover declines through the western United States at the onset of the 21st century.
The longer-term history of fire activity also plays an important role in forest cover change. Slow recovery and conversion to non-forest cover types have been noted at varying spatial scales following 18th-and 19th-century fires in some coniferous forests in the Rocky Mountains (Stahelin 1943, Kaufmann et al. 2000, Huckaby et al. 2001. Thus, increases in the extent or density of some subalpine forests over the last century (Cocke et al. 2005, Zier andBaker 2006) may be partially explained by slow, yet continuous forest recovery following widespread highseverity fires in the 18th and 19th centuries. Forest cover increases have also been documented in some lower elevation montane forests where formerly frequent surface fires have declined during the 20th century (Schoennagel et al. 2004, Hessburg et al. 2005. Therefore, changes in fire occurrence during the 18th-20th centuries could account for a trend of increasing forest cover through different mechanisms (slow post-fire recovery from past severe fires vs. decreased thinning effects of surface fires) that are likely to vary in importance across the elevation gradient from lower elevation montane forests to higher elevation subalpine sites.
In addition to contemporary and historical fire activity, other factors have also played influential roles in forest cover change across the western United States over the past century. For example, moisture availability and changes in land use have been noted as important correlates with changes in dry forests of the northwest and southwest United States since the late 1800s (Merschel et al. 2014, Johnston 2017, Rodman et al. 2017. Extractive logging of large-diameter trees has altered stand structure in many coniferous forests in the western United States; in some cases, late 20th-century forest cover may still be responding to timber extraction in the 19th and early 20th centuries (Naficy et al. 2010, Merschel et al. 2014, Collins et al. 2017. Historical mining activity also led to increases in soil disturbance, logging, and fire activity in many parts of the western United States in the late 1800s and early 1900s (Gruell 1983, 2001, Veblen and Lorenz 1991, Hessburg and Agee 2003, Dethier et al. 2018). More recently, exurban housing development has occurred throughout portions of the United States (Theobald and Romme 2007, Platt et al. 2011, Radeloff et al. 2018, leading to forest fragmentation (Radeloff et al. 2005).
To effectively disentangle the potential drivers of forest cover change, detailed data covering a range of ecosystem types, ownership designations, and human land uses are needed. Though the satellite record can be used to map vegetation change since approximately 1972 (i.e., the launch of Landsat 1), these data lack both temporal depth and high spatial resolution, which is particularly important when studying long-lived tree species and areas of sparse forest cover. Aerial surveys of agricultural and forested lands began throughout the United States in the 1930s (Matthews 2005) and these photographs provide extensive records of land cover at a high (<2 m) spatial resolution (Morgan et al. 2010). These data have been used in quantitative assessments of historical forest change in the western United States (Hessburg et al. 1999, Coop and Givnish 2007, Lydersen and Collins 2018 and of the anthropogenic and biophysical factors influencing these changes (Asner et al. 2003, Black et al. 2003. We used historical and contemporary aerial photography to quantify changes in forest cover across a study area in the northern Front Range (NFR) of Colorado c. 1938-2015. The NFR (Fig. 1) is a region that has been shaped by human activities and by wildfire over the last several centuries (Veblen andDonnegan 2005, Addington 2018) and provides both a data-and disturbance-rich case study of the potential causes of forest change throughout a complex landscape. Fire history , Sherriff et al. 2014, historical stand densities (Battaglia et al. 2018), climate and tree growth (Villalba et al. 1994, Veblen et al. 2000, Lukas et al. 2014, and patterns of exurban development (Platt et al. 2011) are well documented in this region, permitting a thorough analysis of the potential drivers of vegetation change throughout the 20th century. To date, there have been no landscape-scale studies which have attempted to characterize the relative influences of land use, the abiotic environment, and wildfire to observed shifts in forest cover in the NFR. Indeed, there are few studies in similar conifer forest ecosystem in the western United States that have quantified changes in forest cover and implemented robust research designs to assess the contributions of multiple driving factors to 20th-century forest change across broad biophysical gradients. Specifically, we asked the following questions: (1) How has the extent of forest cover changed across the NFR c. 1938-2015? (2) How do elevation, topographic heat load index, land ownership, and historical land use relate to spatial patterns of increase and decrease in forest cover? and (3) How do changes in forest cover through the 20th century relate to known fire history across a range of fire regime types?

Climate and vegetation
Our study area covers 2932 km 2 on the eastern slope of the northern Front Range (NFR) of Colorado, United States (Fig. 1). Elevations range c. 1600-4200 m and climate varies substantially along Fig. 1. Overview of study area, locations of urban areas, and recent wildfires  in the northern Front Range of Colorado. Background color shows the spatial distribution of the different life zones (lower montane, 1700-2400 m; upper montane, 2400-2800 m; and subalpine, 2800-3500 m). Fig. 2. Variability in wildfire activity, tree growth, and climate 1800-2015 in the study area in the northern Front Range of Colorado. Prior to 1984, proportion burned (a) is given as the proportion of montane fire history sites recording spreading fire in a given year (min two trees and 10% of trees within a site). For 1984 burned is based on the percent surface area of the montane zone (1700-2800 m) within Monitoring Trends in Burn Severity fire perimeters in each year. Sample depth (dashed line) gives the total number of fire-scarred recorder trees (across all sites) in a given year. Variability in tree growth is derived from residual tree ring chronologies of ponderosa pine (PIPO; (b) and Engelmann spruce (PIEN; (c). Summer vapor pressure deficit (VPD; (d) is an this gradient, with lower sites being much warmer and drier than sites found at high elevations. Across elevations, minimum temperatures in January range À17 to À5°C and July maximums range 14-32°C (PRISM 2018). Average precipitation is 400-1300 mm/yr, with spring months (March-May) being the wettest (35% of annual totals; PRISM Climate Group, Oregon State University. 2018). The proportion of precipitation occurring as snowfall increases with elevation and snow depth typically peaks in April (Kittel et al. 2015). Recently, increases in temperature have been noted throughout much of Colorado (1.1°C from 1977Colorado (1.1°C from to 2006, including the NFR (Lukas et al. 2014).
Native vegetation in the study area begins with short grass prairie at the lowest sites (<1700 m), transitioning into lower montane forests (1700-2400 m) principally dominated by ponderosa pine (Pinus ponderosa) with scattered Douglas-fir (Pseudotsuga menziesii) and Rocky Mountain juniper (Juniperus scopulorum; Peet 1981, Addington 2018. Upper montane forests range 2400-2800 m and are dominated by ponderosa pine, Douglasfir, and lodgepole pine (Pinus contorta). At the highest forested sites (subalpine; 2800-3500 m), lodgepole pine, limber pine (Pinus flexilis), subalpine fir (Abies lasiocarpa), and Engelmann spruce (Picea engelmannii) are the dominant tree species, giving way to alpine vegetation at c. 3500 m. Colorado blue spruce (Picea pungens) and several broad-leafed, deciduous tree species (e.g., narrowleaf cottonwood [Populus angustifolia] and willows [Salix spp.]) are common in riparian areas. Quaking aspen (Populus tremuloides) is also present in scattered locations either as a post-fire seral component or as the persistent dominant species in some areas (Peet 1981).

Disturbance history
In low-elevation montane forests and interspersed grassland sites (below 2260 m) in the NFR, fires were relatively frequent and of low to moderate severity prior to the 20th century (Sherriff et al. 2014, Brown et al. 2015. In contrast, middle-and upper-elevation montane forests (2260-2800 m) were historically characterized by a mixed-severity fire regime, in which higherseverity patches of mortality (>70% of canopy trees; 50 ha or larger) were more common , Sherriff et al. 2014. Subalpine forests (2800-3500 m) typically burned infrequently (>200-yr intervals) and in large (>1000 ha) high-severity patches during exceptionally dry years (Buechling andBaker 2004, Sibold and. Surface fires affected only 1-3% of the forested area in the subalpine zone , although firescar records do indicate more frequent small fires at the ecotone of forest with rocky alpine sites (Sherriff et al. 2001).
The mid-late 1800s were a time of heightened fire activity across all zones in the NFR, concurrent with widespread droughts throughout the Rockies and expanding Euro-American settlement (Veblen et al. 2000, Sherriff et al. 2001, Kitzberger et al. 2007, Sherriff and Veblen 2008Fig. 2). Following this period, fire activity in the NFR was fairly limited c. 1920-2000 likely due to a combination of direct fire suppression, fuel reduction due to livestock grazing, lack of combustible fuels in higher elevation, recently burned sites, and potentially less suitable climatic conditions for fire ignition and spread (Veblen et al. 2000, Sherriff and Veblen 2008. Recent increases in fire activity have occurred since 1984 following a relatively fire quiescent period for most of the 20th century (Fig. 2).
Insects and pathogens are also important components of stand dynamics in the NFR. In this region, western spruce budworm (Choristoneura occidentallis), Douglas-fir bark beetle (Dendroctonus indicator of summer (Jun-Aug) drought stress, while winter precipitation (PPT; (e) is the sum of Nov-Dec (prior year) and January-May (focal year) precipitation across the study area. Snow water equivalent (SWE; (f) is derived from a long-term record  of snow depth and water content collected 1 May each year within the study area. All tree ring and climate values (b-f) were re-scaled using z-score transformations. Gray lines and bars give annual values, while black lines show 10-year moving averages. Regional pluvials (blue overlay) and droughts (red overlay) were defined as extended periods (five or more years) in which moving 10-year moving averages of growth were above or below (respectively) the long-term mean for both ponderosa pine and Engelmann spruce. ❖ www.esajournals.org pseudotsugae), and mountain pine beetle (MPB; Dendroctonus ponderosae) are the most important potentially lethal insect pests of the lower and upper montane zones, whereas mountain pine beetle, spruce beetle (Dendroctonus rufipennis), and western balsam fir bark beetle (Dryocoetes confusus) are the primary lethal insects of the subalpine zone. Outbreaks of many of these insect species have occurred during the 1938-2015 time period of this study Veblen 1993, Chapman et al. 2012). Although such outbreaks have killed large numbers of trees in some areas (e.g., MPB in the 1970s and the early 2000s), postoutbreak release of advance regeneration and new tree establishment typically lead to rapid forest recovery (i.e., at a decadal scale; Veblen et al. 1991, Hadley and Veblen 1993, Collins et al. 2011).

Land-use history
Human land use in the NFR has shifted over the last two centuries in ways that are likely to have influenced forest cover during the 1938-2015 study period. The mid-late 1800s and early 1900s were associated with the removal of Native American communities (and their associated land-use practices; Simmons 2000) and with increases in logging, mining, and ranching impacts from Euro-Americans (Veblen and Lorenz 1986, Veblen and Donnegan 2005, Dethier et al. 2018. The decline in Native American populations likely led to declines in intentional burning in certain areas, though historical ignitions due to lightning and humans are difficult to separate (Veblen et al. 2000). Past logging and mining have had complex impacts on forests in the NFR. Timber extraction and surface mining initially reduce forest cover. However, the longerterm effects of logging and mining depend on rates of succession and whether the site is converted to another land use such as ranching or residential development.
Although fuel reduction due to livestock grazing contributed to decreased fire spread at some sites as early as the 1860s (Veblen et al. 2000, Brown et al. 2015, widespread fire exclusion in the NFR dates from approximately 1920 when automobile roads allowed effective access to the rugged topography of the region (Veblen andLorenz 1986, Veblen andDonnegan 2005). The more direct effects of early settlers on fire activity ranged from intentionally set fires during the early settlement and mining era (e.g., 1840s-1890s) to more effective fire suppression and exclusion at permanently settled sites in lower elevations (Veblen andLorenz 1986, Veblen andDonnegan 2005). Exurban development has been widespread in the NFR since c. 1940, with home construction shifting over time from predominantly grass-covered valley bottoms to forested slopes more prone to higher-severity fires (Platt et al. 2011). The effects of this development include forest removal for home sites and changes in fire frequency-more aggressive fire suppression but also an increased risk of human-set fires.

METHODS
Interannual and interdecadal trends in climate, wildfire, and tree growth To provide important context for observed trends in forest cover throughout the study area, we compiled data describing interannual and interdecadal variability in climate, tree growth, and wildfire activity 1800-2015 (Appendix S1). We characterized climate using gridded spatial data describing summer drought stress (i.e., vapor pressure deficit) and annual winter precipitation 1896-2015 (PRISM Climate Group, Oregon State University. 2018), as well as fieldderived measurements of snow depth and density 1938-2015 (snow water equivalent from the University Camp snow course (40.03°N, À105.57°W; 3140 m; NRCS 2016)). These components of climate are important to many forest processes in the southwest United States and southern Rocky Mountains (Williams 2012, Andrus et al. 2018. We also used tree ring chronologies from ponderosa pine and Engelmann spruce (Villalba et al. 1994, Veblen et al. 2000)-two dominant conifers in the montane zone and subalpine zone, respectively-to extend the temporal coverage of the climate record to 1800. Lastly, we developed an index of fire activity in the montane zone using fieldderived fire history data (1800-1983Sherriff et al. 2014) and recent wildfire perimeters (1984Eidenshink et al. 2007). To visualize interdecadal trends in climate, tree growth, and wildfire activity, we also calculated 10-year moving averages of all annual values.

Image acquisition and digitization
The historical air photographs used in this study were captured in flights of the NFR commissioned by the U.S. Forest Service and Soil Conservation Service in 1938 and 1940 as part of broad-scale timber and resource inventories. These images have an approximate cartographic scale of 1:20,000 and over 1700 individual scenes were scanned and digitized at 600 ppi by technicians from the University of Colorado library system. The average spatial resolution of these data is 1.1 m. Of the digitized images available through the University of Colorado, we selected a subset of 308 based on the availability of pre-existing fire history and stand age data , Sherriff et al. 2014). This subset was primarily captured in 1938 and the 1940 images represent only 8.6% of the final study area (Table 1); therefore, we refer to these data as "1938 imagery" for simplicity. For comparison with historical imagery, we also acquired countylevel mosaics of NAIP (National Agriculture Imagery Program) imagery from the area that were collected in fall 2015 (Table 1). These data have a 1 m spatial resolution and three spectral bandsone each in the red, green, and blue wavelengths (United States Forest Service, National Agriculture Imagery Program 2015). To use similar methods with 1938 and 2015 imagery, we created a single panchromatic band (the mean of red, green, and blue values) from NAIP data prior to image processing and classification.

Image processing
To aid in the identification of forested areas in the 1938 and 2015 images, we used a series of image processing steps to add supplemental information describing the texture and context surrounding each c. 1 m pixel. During these steps, we identified locally dark pixels (indicative of individual trees surrounded by bright grassland) and quantified local standard deviation in brightness using moving windows at multiple scales (sensu Coburn and Roberts 2004; Fig. 3; Appendix S2). We combined these two layers with the original grayscale imagery to create three-band composite imagery (Fig. 3d). The combination of brightness, local minima, and standard deviation emphasizes differences in spectral reflectance that facilitate the separation of different forest structures (e.g., individual trees and dense stands) from non-forested areas.
We projected each of the 1938 images to NAD 83, UTM 13N (the same projection as the NAIP images) and georeferenced them using metadata that described the approximate center location of each frame, scanning resolution, and cartographic scale. Next, we mosaicked the 1938 imagery by merging overlapping areas along each N-S flight line and then by merging adjacent flight lines into a single photomosaic. Lastly, we co-registered this mosaic to 2015 NAIP imagery using 4955 tie points that were homogeneously distributed throughout the image. We assessed the accuracy of this co-registration using an independent set of tie points (n = 100). These points were located at the nearest identifiable feature to 100 randomly located points throughout the mosaic. Mean horizontal error of these points was 23.6 m (standard deviation = 17.2). Though larger than values of c. 7-16 m reported by two similar studies (e.g., Schoennagel 2009, Lydersen andCollins 2018), this alignment is reasonable considering the larger size of our study area. Still, change detection and further analyses were performed at an aggregated spatial resolution to account for this offset (see Data aggregation).

Image segmentation and classification
Following image processing and mosaicking, we performed image segmentations and developed supervised classifications of the three-band composite imagery for each date (1938 and 2015). These classifications were based on a hybrid approach of pixel-and object-based image analysis, where we used image segmentation to enhance pixel-level classification accuracy (i.e.,  Original imagery (a) was processed using an expanding window method to identify pixels darker than their surroundings (b) and to quantify local standard deviation (c). Next, these three images were merged to create a three-band composite, which we segmented to identify areas of relatively homogenous brightness, contrast, and variance (d). We then classified the segmented imagery using support vector machines to separate forest and non-forest cover types (e; forest shown in green). Following classification, we detected and removed shadows and bodies of water from the forest class by thresholding dark pixels with low local variance. enhanced pixel classification, sensu Radoux and Bogaert 2017). First, we segmented three-band composites from each date (Fig. 3d) using the segment mean shift algorithm in ArcGIS v. 10.4 (ESRI 2016). We prioritized spectral similarity over spatial scale during segmentation parameter selection, which resulted in segments ranging in size from individual trees to entire forest stands. Next, we included a hillshade model (derived from a 10 m digital elevation model and resampled to~1 m; 1938 imagery) and a red-green index (RGI; 2015 imagery) as additional bands to help account for topographic differences in illumination. RGI is a proxy for vegetation health in imagery covering the visible portion of the electromagnetic spectrum (Gartner et al. 2015); as a ratio of two bands, RGI is relatively insensitive to differences in illumination. We used a hillshade model as ancillary data in the 1938 classification because this serves as an additional way to account for topographic differences in illumination in panchromatic imagery. Finally, we assigned~1 m pixels in each set of imagery to forest or non-forest classes using support vector machines (SVMs; Fig. 3e). SVMs are a statistical learning algorithm that can identify non-linear decision boundaries between classes in multidimensional space (Hastie et al. 2009) and are commonly used in image classification.
Following classification, we post-processed classified maps using thresholds of brightness and local standard deviation to identify and remove large shadows that could be misclassified as forest (Appendix S2). We also masked lakes and reservoirs (United States Geological Survey, Geospatial Multi-Agency Coordination (GeoMAC) 2018, United States Geological Survey, National Hydrography Dataset (NHDPlus) 2018), urban areas (from spatial data available through local government offices), and elevations above 3500 m (which have very little or no forest cover) from the final classifications. Of the total study area, 93.9% (2713.2 km 2 ) was left unmasked and included in further analyses.

Data aggregation
To account for co-registration error and to permit additional analyses of the potential drivers of forest change, we aggregated the resolution of the 1 m classified maps to 250 m in each date. We selected a 250 m cell size (6.25 ha) because this extent is several times larger than our estimated co-registration error and corresponds to the coarsest-resolution dataset used in later analyses (i.e., Sohl et al. 2016). This 250 m resolution is also relevant for land managers, approximating the size of individual treatment units that are used in planning forest management activities in the NFR (Addington 2018). Henceforth, we refer to these 250-m cells as "grid cells." We performed this data aggregation by summing the areas of forest and non-forest (from the classified 1 m maps) in each grid cell; we then calculated forest cover in each grid cell as the sum of classified forest area divided by the sum of classified forest and non-forest area. This resulted in two maps of forest cover (1938 and 2015) with continuous values (0-100% in each grid cell) which were later used in change detection (see Accuracy assessment and uncertainty modeling). To limit the influence of grid cells that had little overlap with the study area, we excluded cells with more than 50% of the area masked or beyond the extent of imagery (i.e., no data). Though this criterion removed 15.4% of initial grid cells from later analysis, these cells contained <1% of the total classified area. Lastly, to provide a more detailed characterization of variability in forest cover throughout the study area, we calculated the percentage of grid cells in different forest cover classes (20% bins) in both 1938 and 2015.

Accuracy assessment and uncertainty modeling
We performed global and local accuracy assessments to quantify the classification accuracy of the SVMs, as well as the spatial variability in accuracy throughout the study area. For the global accuracy assessment (sensu Congalton 1991, Olofsson et al. 2014, we used manual photointerpretation at 2000 randomly located points throughout each image (1938 and 2015) and compared photointerpretation ("reference") to predictions from the SVM classification. We then calculated user's, producer's, and overall accuracies following the protocol of Olofsson et al. (2014). Overall accuracy was 90.2% for the 1938 imagery and 89.4% for the 2015 imagery (Appendix S3: Table S1). It should be noted, however, that the reliability of photointerpretation in the 2015 classification is generally higher (because of greater data quality) and estimates of classification accuracy in 1938 are less certain than those in 2015.
To develop local accuracy assessments, we used the 2000 manual photointerpretation points, segmented and unsegmented pixel values at each point, and SVM-predicted classes to develop binomial generalized linear models (GLMs) of correct/incorrect classification (sensu Zimmermann 2004, 2007, Appendix S3). Adjusted deviance squared values (following Guisan and Zimmermann 2000) for the final GLMs were 0.32 for 1938 imagery and 0.28 for the 2015 imagery. We then used these GLMs to predict the probability of correct classification in each 1 m pixel in the 1938 and 2015 imagery. Following GLM predictions, we calculated overall classification accuracy, forest accuracy, and nonforest accuracy as the predicted values from 1 m pixels within each 250-m grid cell. Overall accuracy was calculated as the mean GLM-predicted probability of a correct classification across all classified pixels in a grid cell, while forest and non-forest accuracies were the mean values for all pixels assigned (in final SVM classifications) to each class in each grid cell (Appendix S3: Fig. S1a-d).
Accurate results of change detection analyses are contingent upon proper thematic representations of cover types (i.e., classification accuracy) and upon spatial alignment between datasets (i.e., co-registration error; Leyk et al. 2018). To identify grid cells that had directional changes in forest cover  that exceeded the potential influence of inaccurate classification and co-registration error, we merged our local accuracy assessments with Monte Carlo simulations of error in image alignment (Appendix S3). Monte Carlo simulations consisted of 200 random shifts of the 1938 imagery to quantify the sensitivity of forest cover estimates to co-registration error. We then combined accuracy assessments and Monte Carlo simulations using the following formulas to calculate potential ranges of forest cover for each grid cell for each date: where LowerForest and UpperForest are the lower and upper estimates for forest cover at a given date ( 1938 or 2015 ), %Acc is the estimated classification accuracy for a given cover type, ForestArea and NonForestArea are the total area classified (by SVMs) as forest and non-forest, respectively, and ClassifiedArea is the total area classified as either cover type. Subscripts 6 and 195 correspond to the 6th-and 195th-ranked estimates of forest cover (95% bounds) across the Monte Carlo simulations of co-registration error for the 1938 imagery. All values in the above formulas were specific to each grid cell, thus identifying potential ranges of forest cover in each 250-m area at each date (1938 or 2015). We determined that a grid cell had experienced a significant change in forest cover 1938-2015 if the potential ranges of forest cover in 1938 and 2015 did not overlap. In other words, if a maximum estimate of 1938 forest cover was lower than the minimum estimate of 2015 forest cover, a grid cell was classified as Forest Gain. Similarly, if the maximum estimate of 2015 forest cover was lower than minimum estimate of 1938 canopy cover, a grid cell was classified as Forest Loss. If the error bounds for the two dates overlapped, a cell was then classified as Little Change (generally <15% net change in forest cover). Following this, we manually reclassified 1331 Forest Loss grid cells (3.0% of the total area) to Little Change because designation as Forest Loss in these cells was primarily due to lighting conditions in the 1938 imagery (i.e., overestimates of 1938 forest cover due to topographic shading).

Potential drivers of change in forest cover
To assess the relative influences of several known drivers of change in forest cover, we performed a Random Forest analysis to predict forest change class (i.e., Little Change, Forest Gain, and Forest Loss) in each grid cell as a function of 13 landscape variables relating to land use, land management, ownership designation, the abiotic environment, spatially modeled historical fire regime class, and recent wildfire activity (Table 2; Appendix S4). We did not include gridded climate data in this analysis, instead using elevation and heat load index (derived from slope, aspect, and latitude; McCune and Keon 2002) as proxies for average moisture availability. Because Table 2. Variables included in Random Forest analysis used to identify potential drivers of forest cover change  in the study area in the northern Front Range of Colorado. the study area spans a single mountain range, elevation and aspect are strongly tied to climate and describe much of the variability in site moisture and vegetation (Peet 1981).
We performed the analysis of these 13 landscape variables and their relationship with forest change class in each grid cell using the ranger (Wright and Ziegler 2017), mlr (Bischl et al. 2016), and caret (Kuhn 2008) packages in R (R Core Team 2018). We also tested for multicollinearity of predictors using the rfUtilities package (Evans et al. 2011). During model fitting, we tuned hyperparameters (i.e., "mtry," the number of variables available for splitting and "ntree," the total number of trees) for the final Random Forest analysis using fivefold cross-validation across a range of parameter values; selected values were "mtry" of 6 and "ntree" of 300. For ease of interpretation, permutationbased variable importance (mean decrease in accuracy) was relativized so as to sum to 1 (relative importance). We also characterized final model fit using percent classification accuracy from 100-fold cross-validation. While we present final results of these analyses at a 250 m resolution, we also quantified forest cover change, performed uncertainty thresholding, and completed an additional Random Forest analysis at a 1 km resolution using the same set of 13 predictors to assess the extent to which our results may be sensitive to spatial scale.

Comparison with fire history
In addition to the Random Forest model spanning the study area, we compared net forest change across specific montane and subalpine stands with previously collected field data describing fire history. Previous studies in the montane and subalpine zones of the NFR have compiled a combined total of 1938 samples of fire-scarred trees and 13,832 tree ages from increment cores in an effort to reconstruct past fire activity and stand-origin dates (Veblen et al. 2000, Sherriff and Veblen 2006, 2008, Gartner et al. 2012). Using these records, one of the most extensive datasets yet assembled to characterize fire history in a single region, we performed two separate analyses to quantify forest cover change within (1) montane stands historically characterized by distinctly different fire regimes (low vs. mixed severity) and (2) subalpine stands with differing time since last stand-replacing fire. For each of these analyses, we performed a sensitivity assessment as described in Accuracy assessment and uncertainty modeling to quantify potential ranges of forest cover.
Sites in the NFR with a history of frequent, low-severity fire likely experienced a greater deviation from their natural ranges of variability during the fire suppression era (c. 1920-present) than did sites with mixed-or high-severity fire regimes (Sherriff et al. 2014). For this reason, sites with a history of low-and moderate-severity fire are typically the highest priority for forest management activities such as mechanical thinning and low-intensity prescribed fire (Addington 2018). To assess differences in the magnitude of 20th-century forest cover change by historical fire regime class in fire-excluded sites throughout ( Areas with historically frequent, low-severity fire regimes are those believed to be most affected by 20th-century fire exclusion and are those in which restoration and fuel reduction objectives are believed to be most convergent Sherriff et al. (2014) the montane zone, we compared net change in forest cover using a stand-specific procedure for forests dominated by a frequent, low-severity fire regime (n = 6; 435.2 ha classified) with cover change in those dominated by a less frequent mixed-severity fire regime (n = 67; 3618.8 ha classified). Stands were classified as low or mixed severity using a combination of tree establishment data and fire-scarred samples in a previous study of the montane zone in the Colorado Front Range (Sherriff et al. 2014). The difference in sample sizes between fire regime classes reflects the relative proportions of the modern landscape characterized by the respective fire regimes as well as the accessibility to study areas not impacted by logging and exurban development (a factor that guided initial plot selection). These locations were a subset of the original fire history sites described in Sherriff et al. (2014); we selected forest stands that were the most intensively surveyed (i.e., including site-specific fire histories and/or stand age data), overlapped the area of the imagery, and did not overlap recent fire perimeters . For this analysis, sites within recent fire perimeters were removed because the focus of the montane fire history analysis was to quantify 20th-century montane forest change resulting from fire exclusion and/or longer-term post-fire recovery. In forests characterized by infrequent, highseverity fires (e.g., subalpine forests in the NFR), stand attributes such as tree density and crown cover vary with time since fire (Veblen 1986, Aplet et al. 1988, Donnegan and Rebertus 1999. Therefore, forest cover change c. 1938-2015 could be expected to differ according to stand age. To test this idea, we compared 20th-century forest cover change to stand age (time since last fire) in a 5076.2-ha area in Rocky Mountain National Park-the area of overlap between our imagery and the area studied by .  reconstructed fire history of individual forest stands greater than c. 8 ha using stand-origin methods supplemented with evidence from fire-scarred trees. The study area of this fire history reconstruction spans c. 2500-3500 m in elevation, though >75% is in the subalpine zone (>2800 m) and the vast majority is influenced by higher-severity fire (Sibold et al. 2007). For this analysis, we binned stands by age class and time since last known fire as follows: (1) prior to 1650 (stands without any recorded fire since 1650; 11.3% of the total), (2) 1650-1850 (stands originating after burning in 1650-1850 and prior to Euro-American settlement of the region; 30.4%), (3) 1850-1900 (stands originating in 1850-1900 during a period of enhanced fire activity and drought; 48.5%), (4) 1900-1938 (stands originating following fires early in the 20th century but before initial imagery; 7.0%), and (5) 1938-2015 (stands that originated following fires during the study period, 2.9%).

Changes in forest cover and contributing factors
Between 1938 and 2015, forest cover increased from 56.5% to 64.3% throughout the NFR, a net gain of 7.8% ( Fig. 4; Appendix S5: Table S1). However, there was notable spatial variability in forest cover change. Across the study area, 13.0% of 250-m grid cells were classified as Forest Loss, 39.4% were classified as Forest Gain, and 47.6% experienced Little Change (Fig. 4, Appendix S5: Table S1). Increases in cover were greatest in the subalpine zone (16%) and were more limited in the upper montane (7.4%) and lower montane zones (4.5%). The percentage of total area in different cover classes also varied among zones in 1938 and 2015 (Fig. 5). While lower montane stands had a relatively uniform distribution across cover classes, the upper montane and subalpine zones had a greater percentage of stands with high canopy cover (i.e., >60% cover; Fig. 5).
Of the total classified area (forest/non-forest) in the NFR, 14.3% overlapped with recent  fire perimeters. More than 70% of this fire activity occurred in the lower montane zone, which helps to explain lower increases in forest cover throughout this zone when compared to the upper montane and subalpine. Within recent fire perimeters, forest cover decreased from 56.6% to 19.7%, a net loss of 36.9%. So, while forest cover generally increased in the northern Front Range (NFR) of Colorado 1938-2015, a proportion of these longer-term increases in cover has been offset by recent fire activity in the lower montane zone. The years of peak fire activity between 1800 and 2015 occurred in 1859 and 1860, where 22.5% and 18.4% of fire history sites in the montane zone recorded fire, respectively (Fig. 2a). Given that these two years were not abnormally dry (Fig. 2b, c) and that the sites recording fire were primarily within the upper montane area utilized for mining (Veblen et al. 2000), the sharp rise in fire activity probably reflects increased anthropogenic ignitions-1859 marks the year of widespread increases in intensive mining activity throughout the study area. In the 1938-2015 period, 2012 was the largest  year for fire activity, with fire perimeters intersecting 9.5% of the montane zone (Fig. 2a). This year had above-average summer drought stress, below-average winter precipitation, and belowaverage spring snowpack (1 May SWE; Fig. 2). These periods of historical and contemporary fire activity had important influences on 20thcentury changes in forest cover.
The Random Forest analysis indicated that the proportion of area recently burned  was the largest single contributor to forest cover change (Fig. 6a). When at least 70% of a grid cell intersected a recent fire perimeter, the cell was more likely to belong to the Forest Loss class than to other categories (Fig. 6b). In total, 78.1% of Forest Loss grid cells overlapped 1978-2015 fire perimeters. Elevation and proportion low-severity fire-related variables given the historical prevalence of frequent, low-severity fire at lower elevations-were the second and fourth most important predictors of forest cover change, respectively. Grid cells between c. 1700-2250 m and those historically dominated by low-severity fire (according to the spatial model in Sherriff et al. 2014) showed increased probabilities of membership in the Forest Gain class (Fig. 6c, e). Grid cells at the highest and lowest elevations (i.e., upper treeline and persistent grasslands), as well as mid-upper montane sites, had lower probabilities of Forest Gain, but with some variability (Figs. 4, 6c). Grid cells with moderate to high heat load indices (HLI; i.e., SW slopes with average to above-average levels of direct insolation) showed greater probabilities of Forest Gain than did cooler, NE-facing grid cells (Appendix S5: Fig. S1). This relationship was consistent across elevations, with the exception of the grassland ecotone (<1700 m) where cool and wet sites (i.e., HLI < 0.7) were most likely to be classified as Forest Gain. At the grassland ecotone, 70% of grid cells with HLI < 0.7 were classified as Forest Gain, as compared to 25.8% of grid cells on warm/dry sites (HLI > 0.7). Mine site density was the third ranked predictor overall and the most important predictor related to land use, though effects varied nonlinearly from low to high mine densities (Fig. 6a, d). Grid cells with no recorded mine locations in the vicinity (kernel density of mines = 0/grid cell) were most likely to belong to the Little Change class. In contrast, areas with low to moderate mining density (0-0.125) had higher probabilities of Forest Gain. Sites with the highest mine densities (>0.125) showed increased probabilities of membership to the Little Change class (Fig. 6a, d).
Among the predictors with lower importance, general trends still existed ( Fig. 6a; Appendix S5: Fig. S1). For example, silvicultural treatments, primarily fuel mitigation conducted between 2006 and 2015, led to decreased probabilities of Forest Gain (Appendix S5: Fig. S1). Silvicultural treatments (c. 2006-2015) occurred on~6% of the classified area throughout the NFR and were more common in the montane zone. Surprisingly, the differences in forest cover change among current land ownership classes and 1938 land cover types were relatively minor (Fig. 6a, Appendix S5: Fig. S1). Similarly, proportion of developed land and road density were not strong predictors of forest change in our model (Fig. 6a), likely because some of these areas were developed prior to 1938 or because some of this development took place in previously non-forested areas and thus did not lead to deforestation.
Cross-validation of the Random Forest model (at 250 m resolution) indicated that overall classification accuracy was 66.3% and that accuracy was highest for the Forest Loss class (Appendix S5: Table S2). We also noted that our analyses of forest change were relatively insensitive to aggregation to a 1 km resolution. At this scale, 13.6% of 1-km cells were classified as Forest Loss, 44.9% were classified as Forest Gain, and 41.5% were classified as Little Change. The Random Forest analysis at a 1 km resolution yielded relatively similar rankings of variable importance, similar patterns in partial dependence (Appendix S5: Fig. S2), and similar classification accuracy (65.2%) to analyses at the 250 m resolution. As might be expected, fine-scale drivers such as heat load index became less important at the coarser spatial resolution of 1 km (Appendix S5: Fig. S2a) than they were at a 250 m resolution (Fig. 6a). Patterns in the partial dependence of forest cover change on mine density (Appendix S5: Fig. S2c) and on spatially modeled fire regime class (e.g., proportion of low-severity fire; Appendix S5: Fig. S2e) also became more apparent at the 1 km resolution.

Comparison of changes in forest cover with fire history
In the comparison of forest cover change and stand-specific fire history (i.e., as opposed to the spatially modeled area) in the montane zone, we documented that fire-excluded stands with a history of frequent, low-severity surface fire (prior to 1920) showed greater increases in forest cover 1938-2015 (15.7% net change) than did stands with a history of mixed-severity fire (11.7% net change), corroborating results of the Random Forest analysis with respect to trends by elevation and spatially modeled fire regime (Fig. 7, Appendix S5: Table S1). Sites within the lowseverity class were primarily at lower elevations or in close proximity to grasslands.
In higher elevation forests (primarily subalpine) in a subset of Rocky Mountain National Park with previously reconstructed fire history, forest cover change showed variation with time since fire ( Fig. 8; Appendix S5: Table S1).
Younger stands-those that burned 1850-1938showed greater increases in forest cover than did older stands (i.e., sites with the last fire prior to 1850). These results are likely due to continued but slow post-fire recovery in areas that burned in the century prior to the initial imagery c. 1938 (Fig. 8b, c). In the 1938 imagery, locations with clearly identifiable fires dating from the 1800s and early 1900s were evident elsewhere in the upper montane and subalpine zones of the NFR; many of these burned sites returned to forest cover in 2015. Taken together, these results suggest that continued post-fire recovery from highseverity events in the late 1800s and early 1900s, and a near absence of subalpine fires in the 20th century are some of the major causes for the widespread increases in forest cover at higher elevations 1938-2015. The decline in forest cover for subalpine areas that burned 1938-2015 ( Fig. 8a; À31.9% net change) was driven by a 143-ha subset of the 1978 Ouzel Fire that overlapped with our imagery and the study area of .

DISCUSSION
We quantified 20th-century forest cover change in the northern Front Range of Colorado, a complex landscape spanning broad gradients in the abiotic environment, disturbance regimes, and human influence. We then linked changes in forest cover to data describing elevation, heat load index, historical and contemporary fire history, and land use. The NFR experienced extensive changes in forest cover 1938-2015, but reasons for these changes are spatially varying and a complex combination of drivers was related to the observed changes in forest cover. However, we note three general findings:

Changes in forest cover were influenced by abiotic factors and human influences
Consistent with previous studies in the NFR and in California, United States Schoennagel 2009, Lydersen andCollins 2018), we noted that forest cover changes varied with elevation and that mesic, N-facing grid cells (with lower heat load indices) showed lower increases in cover through the 20th century than did grid cells with other aspects. In many cases, these more mesic slopes were covered by dense forest stands in the initial imagery, which limited potential increases. However, this relationship shifted at the lowest elevations. Similar to the findings of Mast et al. (1997), we found that N-facing slopes at the grassland ecotone were more likely to experience Forest Gain than were S-facing slopes in these elevations. It is logical that forest cover on S-facing sites at lower treeline is moisture-limited, and thus, forest cover remained low at these sites throughout the 20th century. Similarly, in Arizona and New Mexico, United States, increases in tree density during the fire exclusion era were more limited in xeric, moisture-stressed sites (Rodman et al. 2017), indicative of broad-scale climatic controls on tree growth and establishment. Moisture stress may be less limiting to establishment at higher elevations in the montane zone (Chambers et al. 2016, Rother and; instead, time since disturbance may be a more important driver of forest cover change at these higher sites. Fig. 7. Net change in forest cover 1938-2015 by historical fire regime class in 73 fire-excluded montane stands in the northern Front Range of Colorado with stand-specific fire history reconstructions. Stands in the low-severity class were historically characterized by frequent surface fires, while stands in the mixedseverity class experienced occasional stand-replacing patches of wildfire and tend to exist at middle elevations in the study area. Values given are area-weighted change in forest cover across all stands and error bars represent the sensitivity of forest cover change to inaccurate classification and co-registration error. Example image pair (b, c) illustrates forest expansion into grasslands, as well as recent suburban development at a site typified by frequent, low-severity surface fires prior to 1920 (39.96°N, 105.27°W).
Mining activity was the dominant historical anthropogenic impact on forest cover changes in the NFR. Mining activity beginning in the mid-1800s led to surface disturbance rates at least 50 times that of the natural rate of soil disturbance prior to the 19th century (Dethier et al. 2018). Mining brought with it the removal of large-diameter trees for railroad ties and timber, as well as increases in anthropogenic ignitions, all having important influences on forests (Veblen and Lorenz 1986, 1991, Veblen and Donnegan 2005. For example, 1859 and 1860 were the years of most widespread fire activity throughout the NFR c. 1800-2015 and 1859 marks the onset of intensive mining in the region (Veblen et al. 2000, Dethier et al. 2018. Throughout the NFR, areas with low to moderate mining density showed greater probabilities of Forest Gain through the 20th century, likely due to recovery from 1800s logging and burning. However, areas with the highest mining densities were more likely to demonstrate Little Change 1938-2015. We hypothesize that this effect could be attributed to two factors. First, areas of highest mining density were more often privately owned and had higher densities of roads in both 1938 and 2015 than many other areas in the NFR. Therefore, persistent fragmentation due to residential development and road construction is one likely mechanism for lower probabilities of Forest Gain on these intensively mined sites. Second, severe ground disturbance and seed source limitations (due to substantial removal of large-diameter timber) may have led to prolonged forest recovery in sites with high mining density, even those without persistent development or road networks. This relationship between mining density and forest cover change is complex and merits further investigation.
Like many areas throughout the western United States, forest thinning and fire hazard mitigation have been prominent goals in forests in the NFR since the early 2000s (Caggiano 2017, Addington 2018. Though silvicultural treatments were only documented in 6% of the NFR, these treatments influenced forest change in certain areas. For example, treatments were most common in the montane zone and, where present, led to decreased probabilities of Forest Gain 1938-2015. We did not perform analyses that separated silvicultural treatments by intensity or time since treatment, both of which are factors that would Fig. 8. Net change in forest cover (1938 as it relates to time since last stand-replacing wildfire in higher elevation forest stands within a~5000 ha area of Rocky Mountain National Park, CO (a). Error bars represent the sensitivity of each value to inaccurate classification and co-registration error. Areas in the 1938-2015 class are within a single wildfire that took place in 1978 (the Ouzel Fire). Example image pair (b, c) shows a subarea with a recorded fire that occurred in the late 1800s and subsequent recovery through the 20th century (40.31°N, 105.59°W). Fires in the subalpine zone typically burn infrequently (>200 yr) and at high severity. be likely to influence our results and would be useful in future studies. Similarly, we did not include data on sanitation or salvage harvests related to mountain pine beetle outbreaks in the 1970s or extractive logging prior to this time (Veblen and Donnegan 2005) and the influence of these events on forest change c. 1938-2015 is difficult to assess. However, establishment of younger cohorts has been noted following historical logging and it is likely that these mid-20th-century management activities would lead to an altered forest age structure (Veblen and Donnegan 2005) rather than a conversion of forest to non-forest c. 1938-2015.
Land ownership influences forest cover change in a multitude of ways, such as differences in patterns of development, the intensity of fire suppression activities, and differences in the degree of silvicultural treatment by land ownership category. Developed private lands make up the vast majority of the wildland-urban interface (Theobald and Romme 2007), an area where fire suppression activities are particularly intense (Gorte 2013). Forest management activities also differ across management units. For example, Easterday et al. (2018) noted substantial differences in 20th-century forest change by land ownership designation in California, likely as a result of differing long-term management strategies. In the NFR, rates of silvicultural treatment varied by ownership and were highest on city/countyowned lands (12.2% of total city/county area), intermediate on U.S. Forest Service lands (6%) and private lands (3.8%), and lowest on National Park Service lands (0.6%). Widespread and consistent implementation of silvicultural treatments has been difficult in the NFR due to the fragmented nature of land ownership and potential treatment inconsistencies across ownership designations (Calkin et al. 2014). While silvicultural treatment, land ownership, and development intensity were not the most important variables in our analyses, we do not exclude potential influences from these factors given their known influence on ecological systems.
20th-century changes in forest cover varied by life zone and historical fire regime While forest cover increased 7.8% through the entire study area from 1938 to 2015, these increases were most pronounced in the subalpine zone and unburned portions of the lower montane zone. These findings align with what is known about patterns of historical fire activity and tree establishment across life zones in the NFR over the last two centuries. We noted that lower montane stands that were historically dominated by a frequent, low-severity fire regime (<2260 m; Sherriff et al. 2014) showed greater increases in forest cover c. 1938-2015 (15.7% net gain in forest cover) than did midupper montane sites with a history of mixedseverity fire (11.7%). Stand age data indicate that abundant tree establishment occurred during the early 1900s throughout the lower montane zone Veblen 2006, Battaglia et al. 2018). Because it may take several years for tree crowns to expand to a size that could be classified in the 1938 imagery (i.e., >2 m crown diameter), some of this early 1900s establishment likely appeared to be forest cover increase during the 1938-2015 period of our study. In contrast, upper montane sites in the study area experienced abundant tree establishment in the mid-late 1800s following widespread mixed-severity fires and logging (Mast et al. 1998, Ehle and Baker 2003, Battaglia et al. 2018. Thus, some of these mid-upper montane sites may have already appeared as closed forest and relatively dense stands in 1938, limiting potential increases 1938-2015. Still, continued establishment has occurred in some areas of the upper montane zone into the 2000s Veblen 2006, Schoennagel et al. 2011), particularly areas experiencing past logging (Battaglia et al. 2018) or with low initial forest cover (Platt and Schoennagel 2009).
A previous study using historical and contemporary air photographs  in the montane zone of the NFR also identified greater increases in forest cover at low elevations in the montane zone than in the mid-upper montane (Platt and Schoennagel 2009). Yet there were important differences between our study and that of Platt and Schoennagel (2009). While this prior study identified no increases in forest cover in stands above 2432 m, we documented substantial increases in forest cover for some areas in the upper montane and subalpine zones. In addition, the 4% overall increase in forest cover identified by Platt and Schoennagel (2009) is generally lower than that presented in this study (7.8% increase). These differences could be attributed to several factors. Classification scheme, spatial extents of image coverage, and dates of contemporary imagery differed between these studies, which could have affected final outcomes. Specific percentage estimates of forest cover change using historical air photographs are imperfect due to the various limitations of these data Schoennagel 2009, Lydersen andCollins 2018). These limitations make direct comparisons between studies difficult and further highlight the importance of methods that account for potential error and uncertainty in change detection.
The pronounced increases in forest cover through the subalpine zone tell a complex story that may be linked to a combination of logging activities, widespread fires in the 1800s, and slow forest recovery from these events. In a portion of Rocky Mountain National Park (the study area of ), we noted that forest cover increases were generally higher for stands that burned 1850-1938 than within stands originating prior to 1850. This rough classification of stands by age masks some of the ecological variability inherent to subalpine forests in the NFR, where rates and trajectories of forest succession vary strongly by site and with species composition (Veblen 1986, Rebertus et al. 1991, Donnegan and Rebertus 1999, Coop et al. 2010. Still, it is notable that subalpine stands originating following abundant fire activity in the mid-late 1800s (Buechling and Baker 2004 showed much greater increases in forest cover than did older stands. Unlike sites in the lower montane zone, the direct impacts of organized fire suppression in subalpine forests are likely quite minimal given that these forests typically burn in extreme events during which suppression activities are difficult and that the period of active fire suppression is much shorter than a typical subalpine fire interval (Schoennagel et al. 2004, Naficy et al. 2016. While a small portion of the subalpine forests in the NFR had a history of frequent surface fire (e.g., 1-3%; ) and tree invasion into subalpine meadows could have resulted from fire exclusion at these sites, the thinning effects of surface fires in subalpine forests were relatively minor across the NFR landscape (Sibold et al. 2007). Thus, we conclude that forest cover in the subalpine zone in the NFR increased from 1938 to 2015 primarily as a result of continued recovery from fire, logging, and mining disturbances c. 1850-1938 rather than as a result of 20th-century fire suppression.
Contemporary fire activity drove observed declines in forest cover Recent (1978Recent ( -2015 fire activity was the largest single driver of forest change in our Random Forests analysis. Indeed, while forest cover increased throughout the study area 1938-2015, areas within fire perimeters showed substantial declines. Recent wildfires have driven rapid and broad-scale changes in landscapes of the NFR, but these changes should be considered within the context of multiple management goals. On one hand, these past fires, considered in combination with poor post-fire recovery in severely burned areas (Chambers et al. 2016, Rother and, portend a decline in forest cover into the future. Changes in fire activity and forest cover may also alter important ecosystem services (Rocca et al. 2014). However, severe wildfires can benefit many wildlife species (Hutto et al. 2016) and may moderate subsequent fire activity (Parks et al. 2014, Coop et al. 2016, thus acting as effective fuel treatments in a region with a high density of homes in fire-prone areas (Platt et al. 2011). These are important considerations as climate warming is expected to lead to continued increases in wildfire activity throughout the Rocky Mountains (Spracklen et al. 2009, Westerling et al. 2011, Liu et al. 2013).

Study Limitations
Forest cover, the focus of this study, is not necessarily indicative of forest structure and composition. For example, a 250-m grid cell with 30% forest cover could encompass a wide range of stand densities, species compositions, and finescale spatial patterns, all of which are ecologically meaningful and influence myriad ecosystem processes. Though temporal variability in climate probably played an important role in the observed changes in forest cover in the NFR through the 20th century, we did not directly quantify the influences of interannual and interdecadal climate variability on changes in forest cover. These analyses were not feasible given that we used only two sets of imagery spanning a broad time period. Instead, a review and synthesis of potential influences of climate on tree establishment, disturbance activity, and tree growth are provided in Fig. 2 and Appendix S1. We addressed spatial variability in climate by including elevation and heat load index in our analyses. These variables are strong proxies for climate in the NFR.
The seasonal timing of image acquisition has the potential to influence our results, as many quaking aspen and other broad-leafed deciduous trees at higher elevations in the study area likely had no leaves during the collection of 1938 imagery (collected primarily in late October after leaf abscission). Aspen is a fairly minor component of forests of the NFR (Peet 1981) and currently occupies only~7% of our study area in fragmented locations throughout the montane and subalpine zones (Veblen andDonnegan 2005, Rollins 2009). Still, some of the observed increases in forest cover in the upper montane and subalpine zones could reflect conversions of post-disturbance aspen cohorts in 1938 to conifer-dominated systems in 2015, rather than a conversion of non-forest to forest cover.
Shadows cast by individual trees and topography likely led to variability in our estimates of forest cover across the NFR. Recognized by previous studies Schoennagel 2009, Lydersen andCollins 2018), shadows can influence the results of classification in grayscale imagery. For example, large shadows often appear to be dense homogenous stands of forest, while shadows from individual trees appear as extensions of the tree crowns. We addressed these problems by removing large shadows and water bodies from the forest classification through post-classification correction, by using site-specific uncertainty thresholds to identify areas in which forest change exceeded classification error and uncertainty, and by manually correcting this classification when necessary. In most cases, aerial photographs in 1938 and 2015 appeared to have been collected under similar lighting conditions (mid-late morning near the spring or fall equinox), which permits a reasonable comparison through time.
The classification of standing dead trees can also be problematic in grayscale imagery without spectral information in the red and near-infrared wavelengths. Our initial accuracy assessment indicated that misclassification of standing dead trees was an important source of classification error in 2015 imagery. For this reason, we corrected 2015 forest classifications within known burn perimeters using systematic photointerpretation (Appendix S2). Though scattered tree mortality due to insects has occurred throughout much of the NFR c. 1996-2015 (United States Forest Service, R2 Aerial Detection Survey 2015), widespread insect outbreaks (i.e., near-total mortality of large stands) are primarily confined to the northwest corner of the study area. We estimate these outbreaks influenced <5% of the total study area based on manual photointerpretation and comparison with aerial detection survey data (United States Forest Service, R2 Aerial Detection Survey 2015). Insect outbreaks during the late 20th and 21st centuries have been generally much less widespread and of lower severity in our study area than in other regions in Colorado. Furthermore, forest recovery following insect mortality in the overstory can be relatively rapid due to post-outbreak release of advance regeneration as well as new tree establishment , Hadley and Veblen 1993, Collins et al. 2011; thus, conversion to non-forest is not the typical result of recent insect outbreaks in the NFR. Still, we acknowledge that the misclassification of standing dead trees in areas recently affected by insect mortality is a potential cause of error.
One additional source of error that could not be fully accounted for was the variation in image quality throughout the study area, particularly for the 1938-1940 imagery. Some scenes were well preserved, yet others had evidence of fading and creases that may have led to local error in forest cover estimates. We addressed these problems by removing the outer edges of each scene prior to mosaicking (showing the greatest likelihood of fading and distortion) and by manually removing a small number of areas due to poor image quality. Image quality is an important consideration when working with historical aerial imagery in a GIS; automatic detection and removal of poor-quality areas may lead to improved results in the future.

CONCLUSION
Our study highlights that forest cover change is rarely the result of a single causative ❖ www.esajournals.org mechanism. Observed forest cover gains in the northern Front Range 1938-2015 were related to a complex combination of effects of elevation, heat load index, historical fire activity, and 19thcentury mining activity. The effects of fire history varied along the elevation gradient, with subalpine forests showing the greatest increases in cover in part due to recovery from fire activity in 1850-1938. In the montane zone, fire-excluded stands with a history of frequent, low-severity surface fire showed a greater increase in cover than did stands with a history of mixed-severity fire (15.7% vs 11.7%). No single driver of Forest Gain applied across the entirety of the study area. In contrast, Forest Loss was primarily driven by recent fire activity in the montane zone. Late 20th and early 21st-century silvicultural treatments (e.g., forest thinning and prescribed fire), exurban development, and road construction have also played important roles in Forest Loss or in limiting Forest Gain in some areas. Changes in forest cover in the NFR c. 1938-2015 are a reflection of short-term (e.g., development, recent wildfire) and long-term (e.g., pulses of historical fire activity, fire quiescence, climatic influences on tree establishment, and 1800s mining and logging) landscape legacies which are then overlaid upon variability in the abiotic environment.

ACKNOWLEDGMENTS
This work was funded by the City of Boulder Open Space and Mountain Parks, Australian Research Council award DP170101288, and Joint Fire Science Program Graduate Research Innovation award 17-2-01-4. Publication of this article was funded by the University of Colorado Boulder Libraries Open Access Fund. Thanks to Katie Lage and Jon Jablonski from the CU Library system for help with accessing imagery and metadata. Mike Greene and Jake Zatz helped in image processing, and Tania Schoennagel contributed data. Brian Anacker, Michael Battaglia, Rosemary Sherriff, Chris Wanner, and two anonymous reviewers provided useful comments that improved the manuscript.