Effects of Bark Beetle Outbreaks on Forest Landscape Pattern in the Southern Rocky Mountains, U.S.A.

: Since the late 1990s, extensive outbreaks of native bark beetles (Curculionidae: Scolytinae) have affected coniferous forests throughout Europe and North America, driving changes in carbon storage, wildlife habitat, nutrient cycling, and water resource provisioning. Remote sensing is a crucial tool for quantifying the effects of these disturbances across broad landscapes. In particular, Landsat time series (LTS) are increasingly used to characterize outbreak dynamics, including the presence and severity of bark beetle-caused tree mortality, though broad-scale LTS-based maps are rarely informed by detailed ﬁeld validation. Here we used spatial and temporal information from LTS products, in combination with extensive ﬁeld data and Random Forest (RF) models, to develop 30-m maps of the presence (i.e., any occurrence) and severity (i.e., cumulative percent basal area mortality) of beetle-caused tree mortality 1997–2019 in subalpine forests throughout the Southern Rocky Mountains, USA. Using resultant maps, we also quantiﬁed spatial patterns of cumulative tree mortality throughout the region, an important yet poorly understood concept in beetle-affected forests. RF models using LTS products to predict presence and severity performed well, with 80.3% correctly classiﬁed (Kappa = 0.61) and R 2 = 0.68 (RMSE = 17.3), respectively. We found that ≥ 10,256 km 2 of subalpine forest area (39.5% of the study area) was affected by bark beetles and 19.3% of the study area experienced ≥ 70% tree mortality over the twenty-three year period. Variograms indicated that severity was autocorrelated at scales < 250 km. Interestingly, cumulative patch-size distributions showed that areas with a near-total loss of the overstory canopy (i.e., ≥ 90% mortality) were relatively small (<0.24 km 2 ) and isolated throughout the study area. Our ﬁndings help to inform an understanding of the variable effects of bark beetle outbreaks across complex forested regions and provide insight into patterns of disturbance legacies, landscape connectivity, and susceptibility to future disturbance. (1) develop (2) use maps of beetle-caused mortality at a regional extent and 30-m resolution, and (3) quantify the spatial patterns of bark beetle attack. Our study advances prior research because of the inclusion of a robust ﬁeld validation dataset (n = 239 plots), the focus on the cumulative effects of multiple bark beetle species rather than a single mortality agent, and the characterization of patch-size distributions and spatial autocorrelation of beetle-caused tree mortality across a broad, heterogeneous region.


Introduction
Native bark beetles (Curculionidae: Scolytinae) are key drivers of ecosystem structure and function in Earth's temperate and boreal forests [1][2][3]. Bark beetles can cause extensive tree mortality through coordinated mass attacks, altering important ecosystem services such as carbon storage, wildlife habitat, nutrient cycling, and water resource provisioning [4][5][6]. For instance, in the western U.S., c. 3.8 × 10 9 trees were killed by bark beetles since the late 1990s [7], nearly twice that killed by wildfire during the same period [8]. Yet despite the ubiquity of these insects and their importance for forest dynamics, characterizing even basic attributes such as the affected area and severity of effects remains challenging, particularly at relatively fine (c. 30 m) spatial grains and across broad extents. While the spatial patterns of ecological disturbances are known to influence vegetation dynamics and a wide range of ecosystem processes [9][10][11], the spatial patterns of bark beetle outbreaks have rarely been explored, particularly across regions with irruptions of multiple beetle species. Improved mapping efforts that combine detailed field data and remotely sensed products are critical to refining understanding of the extent, severity, and spatial patterns of outbreaks.
Most bark beetles are specialist herbivores that target only a single host-tree species or genus [12][13][14]. Typically, bark beetles persist at low population levels and attack damaged or weakened trees (i.e., endemic conditions); outbreaks are initiated by population irruptions that require abundant susceptible hosts as well as suitable climate conditions [3,15,16]. Beetles in the genera Dendroctonus and Dryocoetes, which are responsible for much of the recent tree mortality in U.S. Rocky Mountain forests, preferentially target denser stands with larger trees, higher stem density, and greater host-tree abundance [17][18][19][20]. Given a suitable landscape, outbreaks may be incited by warm temperatures and drought conditions which increase tree susceptibility to attack and enhance larval development and overwinter survival [20][21][22][23]. Further climate warming is expected to amplify bark beetle activity [24,25], but with the potential for negative feedbacks as beetle attacks and climatic limitations deplete suitable host trees [26][27][28]. Climate, tree species composition, and forest structure vary dynamically across mountainous regions, and these factors establish the template for spatially heterogeneous patterns of bark beetle outbreaks.
Remote sensing techniques are the most effective means of mapping insect-caused tree mortality across large areas [29,30]. One such approach is the US Forest Service Aerial Detection Survey (ADS), an airborne monitoring program in which trained observers manually identify areas with tree damage, defoliation, or mortality [31,32]. Though ADS data are widely used to monitor coarse-scale patterns of insect activity [7,33], locational error, uncertain locations of tree mortality within mapped perimeters, and temporal gaps in data acquisition limit their utility for detailed fine-grain mapping of bark beetle outbreaks. To address this need, a number of semi-automated remote sensing approaches have been used to quantify the presence and severity of bark beetle attack. Landsat imagery has emerged as the most widely used data for these purposes because of their extended temporal record and relatively fine spatial grain [29]. In particular, Landsat image time series (LTS), which combine data from successive time steps, are increasingly used to characterize the effects of bark beetles and defoliators [34][35][36]. Annual LTS products can be used to effectively quantify the progression of outbreaks [22,37]. Temporal segmentation or post-classification correction of LTS products can also reduce image irregularities in individual scenes and help to isolate the effects of specific mortality agents [35,36,38]. While past LTS-based approaches have proven effective for mapping bark beetle outbreaks, the use of multispectral ensembles [39], inclusion of spatial context surrounding individual pixels [40,41], and high-quality field validation have not been widely used in past studies, but are likely to improve mapping efforts.
To quantify the extent, severity, and spatial patterns of recent (c. 1997-2019) bark beetle outbreaks throughout the Southern Rocky Mountains, USA (SRM; EPA Level III Ecoregion 21; Figure 1), we used extensive field data, LTS-based mapping products, and Random Forest models to map bark beetle-caused tree mortality at a regional scale. We focused on Remote Sens. 2021, 13, 1089 3 of 18 subalpine forests in the SRM, where bark beetle-caused mortality has been the most severe and widespread [7,21,22]. Our specific objectives were to: (1) develop and evaluate models to predict the presence and severity of bark beetle-caused tree mortality from spectral characteristics, (2) use these models to create maps of beetle-caused mortality at a regional extent and 30-m resolution, and (3) quantify the spatial patterns of bark beetle attack. Our study advances prior research because of the inclusion of a robust field validation dataset (n = 239 plots), the focus on the cumulative effects of multiple bark beetle species rather than a single mortality agent, and the characterization of patch-size distributions and spatial autocorrelation of beetle-caused tree mortality across a broad, heterogeneous region. focused on subalpine forests in the SRM, where bark beetle-caused mortality has been the most severe and widespread [7,21,22]. Our specific objectives were to: (1) develop and evaluate models to predict the presence and severity of bark beetle-caused tree mortality from spectral characteristics, (2) use these models to create maps of beetle-caused mortality at a regional extent and 30-m resolution, and (3) quantify the spatial patterns of bark beetle attack. Our study advances prior research because of the inclusion of a robust field validation dataset (n = 239 plots), the focus on the cumulative effects of multiple bark beetle species rather than a single mortality agent, and the characterization of patch-size distributions and spatial autocorrelation of beetle-caused tree mortality across a broad, heterogeneous region.

Study Area, Bark Beetles, and Host Tree Species
This study focuses on the effects of tree-killing bark beetles in subalpine forests of the SRM, a mountainous area with complex topography (elevations 1450 to 4400 m) spanning 145,000 km 2 in southern Wyoming, Colorado, and northern New Mexico, USA ( Figure 1)
To restrict analyses to subalpine forests that may have been affected by bark beetles, we defined the study area using the following criteria: (1) within 500 m of ADS-mapped tree mortality attributed to bark beetles (subalpine tree species hosts only) from 1997 to 2019 [48], (2) forested areas in the 1992 National Land Cover Dataset (NLCD) [49] and any vegetation type in 2016 NLCD [44], (3) at least 1 m 2 ha −1 combined basal area of the dominant subalpine conifer species [42], (4)

Data Sources-Field Data and Image Interpretation
To inform remotely sensed estimates of the extent and severity of bark beetle attack throughout the study area, we compiled data from 239 field plots throughout the SRM (Figures 1 and 2, Table S3). Field data were collected from 2010 to 2019 in previous studies of tree mortality, forest structure, and tree regeneration following bark beetle outbreaks in the subalpine zone [20,38,[54][55][56][57][58]. Because all field plots had evidence of bark beetle presence (i.e., 0.7-99.7% basal area mortality), we also supplemented the 239 field plots with an equal number of 'control' points, in which we used interpretation of aerial imagery, Landsat image time series (LTS), and ancillary spatial data describing the locations of fire and timber harvest to identify forested sites with the absence of visible tree mortality from bark beetles or other causes from 1996 to 2019. Using field data, we quantified the severity of tree mortality within each plot as the percentage of pre-disturbance live basal area for all tree species killed c. 1990s-2010s. Because of incomplete records of the causal agents of tree mortality in a subset of field plots (Table S3), we included a small amount of mortality unrelated to bark beetles (e.g., competition, drought stress). However, where causes were identified, beetles were the primary tree mortality agent, driving >90% of total basal area mortality in plots with severe outbreaks [57,59]. Thus, we assumed that the majority of the tree mortality identified in the field plots and final maps could be attributed to bark beetles.

Data Sources-Landsat Time Series
We used LTS to quantify spectral changes indicative of tree mortality in subalpine forests throughout the SRM ( Figure 2). First, we extracted and processed all Landsat Tier 1 Surface Reflectance scenes from 1996-2019 (i.e., the 1997-2019 study period and 1996 to capture initial conditions) using Google Earth Engine (GEE) [60]. These data included imagery from the Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operation Land Imager (OLI); ETM+ images collected 2003-2019 were excluded to remove striping artifacts in final maps ( Figure S2). To account for wavelength differences among sensors, OLI-derived images were harmonized to Figure 2. Flowchart of image processing and statistical modeling steps involved in developing regional maps of the presence and severity of bark beetle attack in the SRM. Steps performed in Google Earth Engine [60] are represented with dashed boxes, and steps performed in the R environment [61] are represented with solid boxes. LandTrendr (Landsat-based detection of trends in disturbance and recovery) is a pixel-based temporal segmentation algorithm used to identify homogeneous periods of spectral increase, stability, and decline [62,63].

Data Sources-Landsat Time Series
We used LTS to quantify spectral changes indicative of tree mortality in subalpine forests throughout the SRM ( Figure 2). First, we extracted and processed all Landsat Tier 1 Surface Reflectance scenes from 1996-2019 (i.e., the 1997-2019 study period and 1996 to capture initial conditions) using Google Earth Engine (GEE) [60]. These data included imagery from the Landsat 5 Thematic Mapper (TM), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 8 Operation Land Imager (OLI); ETM+ images collected 2003-2019 were excluded to remove striping artifacts in final maps ( Figure S2). To account for wavelength differences among sensors, OLI-derived images were harmonized to TM/ETM+ equivalents using coefficients from [64]. In each scene, we used the CFMask-derived quality assurance band to remove pixels obstructed by clouds, shadows, snow, and surface water [65]. We then developed annual image composites for the summer growing season (June 1-September 30) using the multi-dimensional median of all unmasked pixels [66]. From annual image composites, we extracted values from TM-equivalent bands 1-5, and 7, and calculated nine spectral indices sensitive to plant photosynthesis and foliar moisture content ( Figure 2, Table 1). Because spectral indices from the winter season may help to isolate the signals of mortality and growth of evergreen conifers [67][68][69], we also developed annual composites of the Normalized Difference Vegetation Index (NDVI) from December 1 focal yr to April 1 focal yr+1 (Table 1; Supplementary Materials).  To develop predictors of bark beetle activity from yearly time series of the selected spectral bands and indices, we used the GEE implementation of LandTrendr (Landsat-based detection of trends in disturbance and recovery), a temporal segmentation algorithm that partitions LTS into homogeneous periods of spectral increase, stability, and decline [62,63] ( Figure 2). For each band and index, we used LandTrendr to identify spectral decline events (i.e., disturbance segments) from 1996 to 2019, and calculated the total magnitude of all declines, consistent with our focus on cumulative tree mortality throughout the study period. To limit the effects of tree mortality unrelated to bark beetles, we used a maximum slope filter following [36] to exclude spectral declines that occurred with a greater rate of change than is typical of bark beetle attack (Table S5). Because LandTrendr is a pixel-based algorithm that does not incorporate spatial context from the surrounding area, we developed three additional predictors from each band and index that incorporated neighborhood information ( Figure 2). Specifically, we used Simple Non-Iterative Clustering [70] at three different spatial scales (5-, 10-, and 20-pixel spacing of seed locations) to perform an objectbased smoothing of annual band/index values prior to constructing yearly time series and using LandTrendr ( Figure S1). In total, we developed four LTS predictors from each of the 16 spectral bands and indices (for a total of 64 predictors), one using raw values from annual maps, and three using spatial smoothing at a range of scales prior to temporal segmentation using LandTrendr. For comparison with field data, we extracted values of each LTS predictor at each plot center location. Though the footprints of our field plots sometimes intersected multiple pixels, this approach maintains raw values in the LTS data.

Objectives 1 and 2-Predicting Presence and Severity of Bark Beetle-Caused Tree Mortality
We developed two Random Forest (RF) models to link field data to LTS predictors ( Figure 2). The response variables in RF models were the presence (i.e., any tree mortality due to bark beetles) and severity of bark beetle-caused tree mortality (i.e., the percent of cumulative tree basal area loss c. 1990s to 2010s). While presence models are useful for identifying areas most likely to have been affected, severity models provide insight into the ranges of potential ecological effects. We used Variable Selection Using Random Forests (VSURF) [71] to identify a parsimonious subset of c. 10 LTS predictors for inclusion in each RF model and tested for multicollinearity using the 'rfUtilities' package [72]. With the VSURF-selected predictors, we fit final RF models using the 'ranger' package [73] and used 10-fold cross-validation in the 'caret' package [74] to optimize hyperparameters. To compare the relative influences of different bands and indices in final models, we calculated variable importance using the permutation-based mean decrease in accuracy statistic. Because RF regression can reduce predicted values toward the mean of the response, biasing predictions of extreme values [75,76], we corrected RF predictions of severity using the following equation:ŷ whereŷ unc andŷ c are the predicted pixel values at a given field plot, before and after bias correction, respectively. µ rf and µ obs are the means and σ rf and σ obs are the standard deviations of RF-predicted and observed pixel values across all field plots in 10-fold crossvalidation (Supplementary Materials). Finally, we assessed the influences of potential confounding variables not included in RF models such as pre-outbreak forest density, data contributor, bark beetle species, and co-occurring outbreaks of defoliating insects (i.e., non-lethal defoliation of conifers), by comparing RF model residuals with field data. Following model fitting and evaluation, we used final RF models and LTS products to make predictive maps of the presence and severity of bark beetle attack at a 30-m grain size throughout the study area. All variable selection, model fitting, validation, and spatial predictions were performed in R v. 3.6.0 [61].

Objective 3-Quantifying Spatial Patterns of Beetle-Caused Tree Mortality
We used the newly developed maps of bark beetle presence and severity to quantify spatial patterns of bark beetle-induced tree mortality throughout the study area. To describe the spatial relationships of severity with distance, we created empirical variograms in the 'gstat' package [77] in R. For computational feasibility in variogram calculations, we used a sample of 100,000 individual pixels (mean [min] nearest neighbor distance = 301 [30] m). We used an omnidirectional variogram with a lag width of 150 m (a minimum of 12,096 point pairs in the smallest bin) and a maximum distance cutoff of 300 km (approximately 1/3 of the maximum diagonal of the bounding box) [78]. To characterize patch-size distributions of more severely affected forests, we calculated the cumulative proportion of total area with ≥50%, ≥60%, ≥70%, ≥80%, and ≥90% basal area mortality in patches of different sizes (using eight-neighbor connectivity).

Additional Information
Additional information regarding research methods is provided in the Supplementary Materials. Additionally, all field data, spatial data, statistical code, and model outputs are available through Dryad Digital Repository [79].

Objective 1-Models of the Presence and Severity of Bark Beetle-Caused Tree Mortality
Using predictors derived from LTS, RF models were able to predict field-derived measures of the presence and severity of tree mortality due to bark beetle attack. Based on 10-fold cross-validation, the RF presence model had a classification accuracy of 80.3% and Cohen's Kappa of 0.61; the RF severity model had an R 2 value of 0.68 and a root-meansquared-error (RMSE) of 17.3. LTS predictors in final models included shortwave infrared bands (TM-equivalent B5 and B7), spectral indices sensitive to foliar moisture (NBR, NDMI), Tasseled Cap indices (TCA, TCB, TCW), and winter NDVI (Figure 3). LTS variables without spatial segmentation were the highest-ranked predictors, though predictors that included spatial segmentation of annual maps were also retained in each RF model (Figure 3). Importantly, there were no substantial biases in RF predictions based on pre-outbreak forest density, data contributor, or dominant bark beetle species (Supplementary Materials). Non-lethal defoliation by other insect species led to modest overestimates of the severity of beetle-caused tree mortality (Supplementary Materials), but these effects were only notable at the highest defoliation intensities, and defoliation attributed to any agent (e.g., Choristoneura freemani Razowski, Malacosoma spp.) was present in just 15.5% of the study area [48]. Overall, RF presence and severity models appeared robust for use in predictions throughout the SRM.

Additional Information
Additional information regarding research methods is provided in the Supplementary Materials. Additionally, all field data, spatial data, statistical code, and model outputs are available through Dryad Digital Repository [79].

Objective 1-Models of the Presence and Severity of Bark Beetle-Caused Tree Mortality
Using predictors derived from LTS, RF models were able to predict field-derived measures of the presence and severity of tree mortality due to bark beetle attack. Based on 10-fold cross-validation, the RF presence model had a classification accuracy of 80.3% and Cohen's Kappa of 0.61; the RF severity model had an R 2 value of 0.68 and a root-meansquared-error (RMSE) of 17.3. LTS predictors in final models included shortwave infrared bands (TM-equivalent B5 and B7), spectral indices sensitive to foliar moisture (NBR, NDMI), Tasseled Cap indices (TCA, TCB, TCW), and winter NDVI (Figure 3). LTS variables without spatial segmentation were the highest-ranked predictors, though predictors that included spatial segmentation of annual maps were also retained in each RF model ( Figure 3). Importantly, there were no substantial biases in RF predictions based on preoutbreak forest density, data contributor, or dominant bark beetle species (Supplementary Materials). Non-lethal defoliation by other insect species led to modest overestimates of the severity of beetle-caused tree mortality (Supplementary Materials), but these effects were only notable at the highest defoliation intensities, and defoliation attributed to any agent (e.g., Choristoneura freemani Razowski, Malacosoma spp.) was present in just 15.5% of the study area [48]. Overall, RF presence and severity models appeared robust for use in predictions throughout the SRM.  Table 1. Subscripts refer to the scale of spatial segmentation (i.e., 5-, 10-, and 20-pixel seed spacing in Simple Non-Iterative Clustering) of annual maps prior to temporal segmentation, where higher values refer to a coarser-scale segmentation with  Table 1. Subscripts refer to the scale of spatial segmentation (i.e., 5-, 10-, and 20-pixel seed spacing in Simple Non-Iterative Clustering) of annual maps prior to temporal segmentation, where higher values refer to a coarser-scale segmentation with larger image objects. Variables without subscripts were not spatially segmented, and thus represent the total disturbance in a 30-m voxel without accounting for neighborhood context.

Objective 2-Mapping Beetle-Caused Tree Mortality across the SRM
Maps derived from RF models and LTS products identified coarse-scale patterns in bark beetle-caused tree mortality throughout the SRM (Figure 4a,c) but also provided Remote Sens. 2021, 13, 1089 9 of 18 insight into subregional and local variation in the effects of attack (Figure 4e-g). Beetlecaused tree mortality was most widespread in northern Colorado, southern Wyoming, and southwestern Colorado, and less common in central Colorado (Figure 4a,c). The map of bark beetle presence identified 39.5% (10,256 km 2 ) of the total study area as affected by bark beetles (Figure 4a,c). However, detection of lower-severity attack was difficult using LTS. For example, omission error was just 9.6% in field plots with at least 40% basal area mortality, as compared to 70.6% omission for plots below 40% basal area mortality. Thus, many places classified as bark beetle absence may have low or moderate levels of tree mortality. The map of bark beetle severity indicates substantial variation in tree mortality throughout the SRM (Figure 4b,d). Lower levels (0-30%), moderate levels (30-70%), and higher levels of tree mortality (70-100%) were predicted across 29.6% (7674 km 2 ), 51.1% (13,245 km 2 ), and 19.3% (5017 km 2 ) of the study area, respectively (Figure 4b,d). Therefore, low and moderate levels of tree mortality were widespread throughout the region, but higher levels were more limited, primarily constrained to forests affected by mountain pine beetle and spruce beetle in northern Colorado and southern Wyoming, and forests affected by spruce beetle in southwestern Colorado. sent the total disturbance in a 30-m voxel without accounting for neighborhood context.

Objective 2-Mapping Beetle-Caused Tree Mortality across the SRM
Maps derived from RF models and LTS products identified coarse-scale patterns in bark beetle-caused tree mortality throughout the SRM (Figure 4a,c) but also provided insight into subregional and local variation in the effects of attack (Figure 4e-g). Beetlecaused tree mortality was most widespread in northern Colorado, southern Wyoming, and southwestern Colorado, and less common in central Colorado (Figure 4a,c). The map of bark beetle presence identified 39.5% (10,256 km 2 ) of the total study area as affected by bark beetles (Figure 4a,c). However, detection of lower-severity attack was difficult using LTS. For example, omission error was just 9.6% in field plots with at least 40% basal area mortality, as compared to 70.6% omission for plots below 40% basal area mortality. Thus, many places classified as bark beetle absence may have low or moderate levels of tree mortality. The map of bark beetle severity indicates substantial variation in tree mortality throughout the SRM (Figure 4b,d). Lower levels (0-30%), moderate levels (30-70%), and higher levels of tree mortality (70-100%) were predicted across 29.6% (7674 km 2 ), 51.1% (13,245 km 2 ), and 19.3% (5017 km 2 ) of the study area, respectively (Figure 4b,d). Therefore, low and moderate levels of tree mortality were widespread throughout the region, but higher levels were more limited, primarily constrained to forests affected by mountain pine beetle and spruce beetle in northern Colorado and southern Wyoming, and forests affected by spruce beetle in southwestern Colorado.

Objective 3-Spatial Patterns of Beetle-Caused Tree Mortality
Spatial patterns of bark beetle-caused tree mortality throughout subalpine forests in the SRM reflect the complex and heterogeneous nature of these biotic disturbances ( Figure 5). The severity of beetle-caused tree mortality was most strongly autocorrelated at lag distances of 0-5 km. However, 30-m pixels were increasingly dissimilar up to distances of 250 km, indicating spatial dependence in patterns of tree mortality across extremely broad extents (Figure 5a). Cumulative patch-size distributions varied markedly for forests exceeding different mortality thresholds in the SRM (Figure 5b). For example, while half of the area exceeding 50% tree mortality was contained in patches larger than 7.74 km 2 , half of the area exceeding 90% tree mortality was in patches smaller than 0.24 km 2 . Thus, areas exceeding 50% tree mortality were relatively large, widespread, and contiguous, but areas exceeding 90% mortality (i.e., with a near-total loss of the overstory tree canopy) tended to be small and disconnected throughout the region.
summarize maps (a,b), respectively. (e-g) Inset maps give an example of an area that experienced a high-severity spruce beetle outbreak c. 2010 (37.75 N and 106.91 W). (e) A false-color composite of 1-m aerial photography from the National Agriculture Imagery Program from 2019 alongside maps of (f) presence and (g) severity in the same area, shows the ability of maps to identify within-stand variation in beetle-caused tree mortality. In (e) the false-color composite, red represents live vegetation while grey represents standing dead trees or a lack of vegetation.

Objective 3-Spatial Patterns of Beetle-Caused Tree Mortality
Spatial patterns of bark beetle-caused tree mortality throughout subalpine forests in the SRM reflect the complex and heterogeneous nature of these biotic disturbances ( Figure  5). The severity of beetle-caused tree mortality was most strongly autocorrelated at lag distances of 0-5 km. However, 30-m pixels were increasingly dissimilar up to distances of 250 km, indicating spatial dependence in patterns of tree mortality across extremely broad extents (Figure 5a). Cumulative patch-size distributions varied markedly for forests exceeding different mortality thresholds in the SRM (Figure 5b). For example, while half of the area exceeding 50% tree mortality was contained in patches larger than 7.74 km 2 , half of the area exceeding 90% tree mortality was in patches smaller than 0.24 km 2 . Thus, areas exceeding 50% tree mortality were relatively large, widespread, and contiguous, but areas exceeding 90% mortality (i.e., with a near-total loss of the overstory tree canopy) tended to be small and disconnected throughout the region. (a) Sample variogram shows the relationship of tree mortality (i.e., percent basal area loss) among sites at distances from 0 to 300 km, where higher values on the y-axis indicate greater dissimilarity among sites at a given distance. Note that the x-axis in (a) is broken with different scales for greater interpretability from 0 to 50 km. (b) The cumulative patch-size distribution gives the proportion of the total area exceeding different mortality thresholds contained within patches of different sizes. In (b), a convex line shape indicates that the majority of the area is contained (a) Sample variogram shows the relationship of tree mortality (i.e., percent basal area loss) among sites at distances from 0 to 300 km, where higher values on the y-axis indicate greater dissimilarity among sites at a given distance. Note that the x-axis in (a) is broken with different scales for greater interpretability from 0 to 50 km. (b) The cumulative patch-size distribution gives the proportion of the total area exceeding different mortality thresholds contained within patches of different sizes. In (b), a convex line shape indicates that the majority of the area is contained within small patches (e.g., ≥90% basal area loss), while a concave shape indicates that large patches are more common (e.g., ≥50% basal area loss).

Discussion
We integrated field data and remotely sensed products to develop seamless maps of the presence and severity of bark beetle outbreaks in the subalpine zone of the SRM, an area that has experienced extensive tree mortality since the late 1990s. The methods used in the present study provide a framework that can be extended to other regions by leveraging freely available spatial data (e.g., Landsat imagery, land cover maps, insect and disease surveys) and existing field datasets to quantify the effects of biotic disturbances on forest structure and landscape pattern. Herein, we noted four key findings: (1) the presence and severity of bark beetle outbreaks of multiple species could be reliably mapped across a complex region using remotely sensed data, (2) at least 10,256 km 2 of subalpine forests (39.5% of the study area) experienced bark beetle-caused tree mortality from 1997 to 2019, (3) low to moderate levels of tree mortality (i.e., <70% basal area loss) were the most common in the SRM, while the most heavily affected areas were primarily in northern and southwestern Colorado, and (4) severity was autocorrelated at broad distances of up to 250 km, and areas with ≥ 90% tree mortality were relatively small (<0.24 km 2 ) and isolated.

Implications for Remotely Sensed Detection of Tree Mortality
Previous studies using LTS have effectively quantified the presence and severity of bark beetle-caused tree mortality in a range of contexts, with reported accuracies ranging 75-90% (i.e., classification accuracy in presence/absence models) and R 2 values ranging from 0.6 to 0.8 (i.e., observed vs. predicted values in severity models) [34][35][36][37][38]. Accuracies for our models were within the reported ranges of these prior studies (i.e., presence accuracy c. 80%; severity R 2 c. 0.7), which is notable considering the large, heterogeneous study area and our focus on multiple insect species with differing spectral signals. Still, separating low-severity tree mortality from unaffected areas remains a difficult challenge when using LTS and other remotely sensed data [34,80,81]. For example, the timing and occurrence of spruce beetle outbreaks are difficult to identify when <35% of a Landsat pixel is affected [82]. Similarly, here, we found that field plots with <40% basal area mortality were erroneously omitted at much higher rates in presence maps. Our use of widespread field data in model fitting and validation also allowed us to evaluate potential confounding factors that are not typically assessed in broad-scale mapping efforts, including pre-outbreak stand density, bark beetle species, and co-occurring outbreaks of defoliators (Supplementary Materials). RF presence and severity models appeared robust to many of these factors, though we demonstrated that severe defoliation can lead to overestimates of beetle-caused tree mortality, aligning with prior work that has demonstrated confusion among insect agents [35]. Thus, non-lethal defoliation events should be considered in the future development of maps characterizing tree mortality.
This study also provides insight into methods used in future regional and national disturbance mapping efforts. We found that shortwave-infrared bands (i.e., B5 and B7) and spectral indices that incorporate them (e.g., NBR, NDMI, TCW) were particularly sensitive to tree mortality, in agreement with prior research of bark beetles and other agents [29,39,[83][84][85][86]. We also found that incorporating spatial context in the development of LTS predictors led to modest improvements in final RF models. Though we used spatial segmentation of annual maps to reduce within-stand spectral noise in LTS, object-based approaches may also be helpful for change attribution because they can be used to quantify the size and shape of events [87,88]. Still, a combined spatial and temporal segmentation approach may be less effective in bark beetle-affected forests where tree mortality can unfold gradually over several years [37,82], than in forests affected by discrete events such as wildfire and timber harvest. Finally, NDVI during the winter period was retained as a predictor in each of the final RF models, as opposed to summer NDVI, which had low predictive ability. Thus, winter imagery is a valuable resource for mapping the mortality and regrowth of evergreen conifers in temperate forests [67][68][69].

Patterns of Bark Beetle-Caused Tree Mortality in the SRM
Recent bark beetle outbreaks have altered forest ecosystems in the SRM and similar ecosystems throughout the Northern Hemisphere [2,21,22,89,90]. By developing a mapping approach that links field data with LTS products, we found that over 10,000 km 2 of subalpine forest area has been affected by three important bark beetle species-mountain pine beetle, spruce beetle, and western balsam bark beetle-in just two decades. Furthermore, the total affected area is probably much higher because of omission error in stands with <40% mortality. The total extent of this bark beetle-caused tree mortality is striking, yet regional maps of severity also illustrate notable variability in the effects of outbreaks throughout the region, likely attributed to past disturbances, tree species composition, stand structure, landscape connectivity, topography, and climate [21,22,27,91,92].
The spatial patterns of cumulative tree mortality reflect the heterogeneous nature of bark beetle attack throughout the SRM, as shaped by factors influencing outbreak initiation and spread. Using maps of the severity of bark beetle attack, we found that individual 30-m pixels were autocorrelated at broad distances of up to 250 km, but were most closely related at distances of 0-5 km. For individual beetle species, regional drivers such as drought and warm temperatures can lead to the synchronous development of outbreaks that extend far beyond the dispersal limits of individual populations [93,94], leading to correlated outbreak dynamics at distances of 80-900 km or more [21,22,95,96]. Yet, stand-scale variation in the effects of bark beetle outbreaks is also promoted by the abundance of host trees and landscape connectivity [21,89]. These same drivers, of broad-scale climate, and local-scale structure and composition, may also shape patterns of the cumulative effects of multiple bark beetle species in the SRM. Interestingly, we found that areas of ≥90% basal area mortality were relatively small and isolated throughout the region, with more than 50% of that area contained in patches smaller than 0.24 km 2 . In contrast, c. 50% of the forest area burned at high severity (>90-95% basal area mortality) in recent fires throughout the Northern Rocky Mountains, USA, was in patch sizes larger than 10 km 2 [97], a difference of over an order of magnitude when compared with similarly severe patches caused by outbreaks presented here. In addition, even in the areas most severely affected by beetle-caused tree mortality, seedlings and smaller trees typically survive beetle outbreaks, resulting in a divergent post-disturbance response when compared to stand-replacing fires [57,[98][99][100][101]. Fire and bark beetles, the two dominant disturbance types in subalpine forests throughout the SRM [102], have fundamentally different impacts on these important ecosystems.
The extensive bark beetle-caused tree mortality in forests throughout the SRM raises questions about trajectories of these ecosystems in a warmer future. Following bark beetle outbreak, abundant advance regeneration and new seedling establishment suggest that many forest stands have the potential to return to a similar structural state within a few decades [99,103,104]. Still, the composition and density of post-outbreak stands are highly variable due to outbreak severity, topography, weather, local microclimatic conditions, and prior disturbance by fire or blowdown that influence regeneration dynamics [57,101,[105][106][107][108]. Many of these same factors can also cause differences in individuallevel growth [57], leading to broad-scale variation in rates of post-disturbance recovery that may influence susceptibility to future disturbance [27,109]. While warming climate conditions have the potential to amplify bark beetle activity in the future, these effects will also be moderated by the depletion of suitable host trees, shifts in species composition, and variable rates of recovery across the landscape [26,28]. Maps of the effects of recent bark beetle outbreaks, such as those developed in the present study, can help to inform expectations of future forest dynamics in the SRM and similar regions, but future work is needed to tie this information to records of post-outbreak successional trajectories, tree growth, and the potential for subsequent disturbances.

Study Limitations
While the maps of the presence and severity of bark beetle outbreaks developed here provide a useful resource for understanding the effects of recent biotic disturbances throughout the SRM, important limitations of these data should also be recognized. Unlike active remote sensing products (e.g., lidar, radar) and other data capable of capturing three-dimensional forest structure (e.g., structure-from-motion), Landsat imagery characterizes the spectral reflectance of the visible overstory canopy and is limited in its ability to monitor subcanopy effects [110,111]. Thus, it is likely that some mortality of subcanopy trees has been omitted in our maps, though subcanopy trees are not commonly targeted by aggressive bark beetles [12,112] and represent just a small percentage of stand basal area. Restricting analyses to US Forest Service lands without recent fires or timber man-agement activities was necessary to focus on the effects of bark beetle attack. However, these restrictions may have led to slight underestimates of patch sizes due to edge effects in spatial pattern analyses. While we infer that bark beetles were the primary mortality agent in regional maps, we could not fully exclude factors such as competition and drought stress, sudden aspen decline, white pine blister rust, and root disease that are important in parts of the SRM [58,113,114]. Forested landscapes with a diverse array of structural conditions and tree species compositions are likely to be influenced by a wide range of mortality agents, which complicates change attribution.

Conclusions
A warming climate is accelerating tree mortality due to bark beetle outbreaks, wildfire, and other drought-mediated forest disturbances. Because disturbances play an important role in the structure, composition, and spatial patterns of vegetation, an understanding of their effects is critical to understanding future trajectories of forest ecosystems. To this end, remote sensing and field inventories are two of the most important tools for understanding ecosystem dynamics at broad spatial scales. In the present study, we combined extensive field surveys of tree mortality in stands attacked by bark beetles c. 1997-2019 and LTS products to predict the presence and severity of bark beetle attack at a c. 30-m resolution throughout the SRM, USA. Bark beetles have driven extensive tree mortality in these ecosystems, but their effects were highly variable across the study area. Still, the legacies of recent bark beetle outbreaks may influence forest dynamics for decades to centuries in subalpine forests, with crucial implications for wildlife habitat, carbon budgets, and other important ecosystem services.
Supplementary Materials: The following are available online at https://www.mdpi.com/2072-4 292/13/6/1089/s1: Additional descriptions of research methods, Figure S1: Example of image segmentation effects on annual maps, Figure S2: A comparison of Random Forest model results when including and excluding Landsat ETM+ imagery, Figure S3: Example of maximum slope thresholds for beetle-caused tree mortality, fire, and timber harvest, Figure S4: Observed vs. predicted values from the Random Forest model of severity with and without bias correction, Table S1: Summary of total mapped area from US Forest Service Aerial Detection Surveys from 1997 to 2019, Table S2: US Forest Service Aerial Detection Survey damage agent and host combinations used to restrict analyses, Table S3: Summary of field datasets, data contributor, and collection protocols, Table S4: Summary of forest structure variables and Random Forest model accuracies for the different field datasets, Table S5: Maximum slope thresholds used to exclude disturbances that occurred with a greater rate of change than beetle-caused tree mortality.  Data Availability Statement: All field data, spatial data, statistical code, and model outputs are available through Dryad Digital Repository at https://doi.org/10.5061/dryad.1rn8pk0sn (accessed on 1 February 2021) [79].