Landscape Topoedaphic Features Create Refugia from Drought and Insect Disturbance in a Lodgepole and Whitebark Pine Forest

: Droughts and insect outbreaks are primary disturbance processes linking climate change to tree mortality in western North America. Refugia from these disturbances—locations where impacts are less severe relative to the surrounding landscape—may be priorities for conservation, restoration, and monitoring. In this study, hypotheses concerning physical and biological processes supporting refugia were investigated by modelling the landscape controls on disturbance refugia that were identiﬁed using remotely sensed vegetation indicators. Refugia were identiﬁed at 30-m resolution using anomalies of Landsat-derived Normalized Difference Moisture Index in lodgepole and whitebark pine forests in southern Oregon, USA, in 2001 (a single-year drought with no insect outbreak) and 2009 (during a multi-year drought and severe outbreak of mountain pine beetle). Landscape controls on refugia (topographic, soil, and forest characteristics) were modeled using boosted regression trees. Landscape characteristics better explained and predicted refugia locations in 2009, when forest impacts were greater, than in 2001. Refugia in lodgepole and whitebark pine forests were generally associated with topographically shaded slopes, convergent environments such as valleys, areas of relatively low soil bulk density, and in thinner forest stands. In whitebark pine forest, refugia were associated with riparian areas along headwater streams. Spatial patterns in evapotranspiration, snowmelt dynamics, soil water storage, and drought-tolerance and insect-resistance abilities may help create refugia from drought and mountain pine beetle. Identiﬁcation of the landscape characteristics supporting refugia can help forest managers target conservation resources in an era of climate-change exacerbation of droughts and insect outbreaks.


Introduction
Droughts and insect outbreaks are primary disturbance processes by which climate change has been implicated in tree mortality and forest change globally [1,2]. Droughts can trigger outbreaks of cambium-feeding insects by exerting physiological stress that weakens host tree defenses [3][4][5]. Although many forest ecosystems are adapted to episodic droughts and insect outbreaks, intensification of these disturbances (i.e., increased frequency or severity) is a growing conservation concern in the context of projected warmer and drier climates in much of western North America [6][7][8]. Combined effects from drought and insect outbreaks can be especially damaging in forests that are already affected by other stressors. For example, whitebark pine (Pinus albicaulis Engelm., WP) population declines across western North America are linked to interactions between drought, mountain pine beetle (Dendroctonus ponderosae, MPB), white pine blister rust (causal agent Cronartium ribicola), importance as a keystone species, WP was recommended for listing under the U.S. Endangered Species Act [7,11].
One approach to help conserve threatened natural communities in the face of climate change involves identification of refugia, or localized areas that are relatively buffered from environmental change over time [12]. Refugia identification has primarily emphasized physical climate attributes, such as cool or moist microclimates to buffer against long-term regional warming and drying [13][14][15]. Changing disturbance regimes represent an important manifestation of climate change in forest ecosystems [2], indicating a need to locate and investigate possible refugia from disturbances such as drought and fire [16,17]. Drought refugia can be conceptualized as locations where vegetation productivity is maintained relative to the surrounding landscape during drought [18] and can be expanded to include drought-induced disturbances such as insect outbreaks. Although refugia identification in a given area can be locally helpful to natural resource managers, application of refugia concepts to regional conservation requires that mechanistic linkages be established in time and space between refugia and landscape drivers such as topography, soil characteristics, and vegetation patterns [12,15,16,19,20].
A conceptual model can guide identification and analysis of refugia from droughts and insect outbreaks ( Figure 1). "Top down" drivers of disturbance dynamics at landscape scales interact with "bottom up" controls at the scale of individual forest stands [5] to influence the spatial distribution of potential refugia. These finer-scale controls include landscape topoedaphic characteristics that shape local hydrologic processes such as snowpack accumulation and melt timing, runoff and infiltration, growing-season evaporative demand, transpiration, soil water storage, and shallow groundwater availability [13,[21][22][23][24][25][26][27]. Interactions between these hydrologic processes can produce hydrologic microrefugia in the form of localized buffering of soil moisture reserves from the effects of regional drought [15,28].

Figure 1.
Conceptual diagram of potential factors affecting the presence of refugia in forests during drought-insect disturbances. Orange boxes are "top down" drivers of disturbance dynamics operating at regional to landscape scales. Gray boxes are finer-scale landscape characteristics that produce spatial heterogeneity in soil moisture and ecological conditions including locations of refugia (blue and green boxes). Ecological refugia from disturbance may be supported by interactions between hydrologic refugia and spatial patterns in forest structure, including traits that confer drought and insect resistance. Conceptual diagram of potential factors affecting the presence of refugia in forests during drought-insect disturbances. Orange boxes are "top down" drivers of disturbance dynamics operating at regional to landscape scales. Gray boxes are finer-scale landscape characteristics that produce spatial heterogeneity in soil moisture and ecological conditions including locations of refugia (blue and green boxes). Ecological refugia from disturbance may be supported by interactions between hydrologic refugia and spatial patterns in forest structure, including traits that confer drought and insect resistance. Fine-scale biological controls are also important, especially forest structure characteristics that can produce spatial variability in drought tolerance and insect resistance at the stand level, including distribution of species, growth rates, age and size classes, and tree density [9,[29][30][31][32]. Spatial patterns in tree functional traits may also redistribute soil moisture and potentially help shape hydrologic refugia through processes such as hydraulic lift, canopy interception, and routing of precipitation [15]. While ecological refugia (locations where trees are relatively buffered from physiological stress and mortality) may be supported in part by hydrologic refugia during droughts, they will likely be shaped by complex interactions between soil hydrology and stand-level biological traits that influence forest resistance and resilience to disturbances.
Existing knowledge concerning landscape controls on drought and MPB severity provides a reasonable basis to form hypotheses concerning spatial determinants of refugia (Table 1). Empirical evaluation of these hypotheses has potential both to advance understanding of disturbance processes and to clarify mechanisms supporting disturbance refugia. For example, existing susceptibility indices for MPB are based largely on stand characteristics such as age, species composition, and density [29] that do not account for topographic and soil controls on soil moisture availability during droughts. Topographic controls have been identified for refugia from fire [16] and temperature increases [14], and could plausibly operate similarly for drought-insect disturbance refugia. This is an unresolved question because some studies have reported topographic buffering from drought and insect disturbance [33], while others have not [9,30]. Table 1. Hypothesized relationships between landscape characteristics (topographic, soil, and forest variables) and likelihood of refugia from drought and mountain pine beetle.

Elevation (m) positive
More refugia are expected at higher elevations. Higher elevations are generally associated with greater precipitation [34] and lower evaporative demand [35]. Greater snowpack and later onset of spring snowmelt at higher elevations may help maintain soil moisture, especially on leeward slopes and in topographically shaded areas. Cooler temperatures at higher elevations are less conducive to MPB brood survival [36,37].

Slope (percent rise) unknown
Depending on which physical mechanisms dominate, refugia might be associated with flatter areas or alternatively with steeper slopes. Steeper slopes promote faster runoff and less infiltration of rainfall and snowmelt, potentially creating drier soils [33]; however, steep leeward slopes near ridgelines where snowdrifts accumulate may have deeper snowpack that persists later into the growing season [38][39][40]. Snowmelt may be further delayed and growing season temperatures and evapotranspiration (ET) may be reduced on steep slopes that are topographically shaded [22,24]. Steeper slopes may also support lower density stands with fewer large-diameter mature trees, resulting in less competition for soil water during droughts and less favorable conditions for MPB [41].
Topographic position index (TPI) calculated with 300-m radius negative More refugia are expected in topographically concave areas. Low TPI indicates landscape concavity, e.g., coves or valley bottoms where cold-air pooling (CAP) may occur. CAP may enable dew formation and/or higher humidity, moderating the effects of drought [15,42]. Cooler temperatures and higher humidity associated with CAP might suppress MPB infestation. More refugia are expected in topographically shaded areas on poleward-facing slopes). In low-HLI (topographically shaded) areas, reduced evaporative demand could maintain greater soil moisture during droughts [30,35,42]. Spring snowmelt may occur later in these areas, allowing soil moisture to last longer into the growing season [22,24]. Topographically shaded slopes have been shown to have greater soil water retention [23]. MPB generally favors warmer, south-facing slopes (high HLI areas) that are more favorable to brood survival [43].
Compound topographic index (CTI) Higher CTI associated with streams & riparian areas, lower values are ridgetops positive More refugia are expected in riparian areas. Topographic convergence predicts groundwater expression in streams and riparian water tables [15] and acts as a steady-state wetness index [44]. CTI is positively related to soil depth, silt and clay content, and water holding capacity [41]. Areas of high CTI (stream channels, riparian areas) are expected to maintain greater soil moisture during droughts than low CTI-areas (ridgetops).
Soil bulk density (SBD; kg/m 3 ) SBD 0 cm = SBD at the soil surface; SBD 100 cm = SBD at 1-m depth negative More refugia are expected where soils are less dense/compacted. Lower soil bulk density allows greater infiltration of rainfall and snowmelt and is associated with higher porosity and thus greater water-holding capacity of soil [23,45]. Also, lower bulk density facilitates greater root elongation and density [45,46]. In turn, trees with more extensive root systems may be more resilient to drought.
Total basal area (m 2 /ha) Higher basal area indicates greater forest density negative More refugia are expected in areas with low forest density. Densely stocked forests may have increased vulnerability to drought mortality [30]. Forest density is positively associated with severity of MPB damage and mortality, particularly for large-diameter trees [9,31]. Remote-sensing approaches, particularly analysis of multi-temporal Landsat imagery, are commonly used to quantify the effects of disturbances such as droughts and MPB outbreaks in conifer forests [47,48]. Trees killed by MPB progress from a green-attack stage (retaining green foliage) through red-attack (a color shift that typically occurs 6 to 12 months after initial attack), to the gray-attack stage once foliage is lost [43]. Landsat-derived spectral indices such as the Normalized Difference Vegetation Index (NDVI), Normalized Differenced Moisture Index (NDMI), and the Tasselled Cap Transformation (TCT) have demonstrated accuracies ranging from 70 to 90% in classification of the red-attack stage [41,43,49,50]. Multi-temporal Landsat imagery has also demonstrated relatively high accuracy in quantifying ecophysiological effects of drought, including reductions in leaf-area index and leaf chlorophyll due to stomatal closure, and mortality [51,52].
The objectives of this study were to (1) identify refugia from drought and insect disturbance using multi-temporal Landsat data representing spatial and temporal patterns of vegetation disturbance severity, (2) model the landscape controls (topographic, soil, and forest characteristics) on locations of identified refugia, and (3) evaluate hypotheses concerning the biophysical mechanisms supporting refugia ( Table 1). The study was conducted in an area with complex topography that experienced a single-year drought in 2001 without any insect disturbance as well as a multi-year drought from 2007 through 2010 that coincided with a severe MPB outbreak. Because these two droughts differed in their magnitudes and disturbance interactions, a quantitative comparison was performed of the ability of landscape characteristics to explain and predict refugia locations. In addition, because this study area includes stands of WP as well as lodgepole pine (Pinus contorta Dougl. ex. Loud.; LP), comparisons were made between different stand types regarding landscape controls on refugia.

Study Area
This study was conducted in the Gearhart Mountain Wilderness, an approximately 95-km 2 designated wilderness area in southern Oregon ( Figure 2). The study area is located in a dry conifer forest transition zone between mesic conifer forests of the Cascade Mountains to the west and the semi-arid sage steppe of the Great Basin to the east. Elevation ranges from approximately 1750 m near the study area perimeter to 2530 m at the summit of Gearhart Mountain. Surficial geology consists primarily of andesite and basalts, overlain discontinuously by glacial and fluvial deposits [53]. The topography of the study area has been shaped by glacial erosion and includes several prominent cirques (U-shaped glacial valleys) separated by ridges. Upper slopes of cirques commonly have slopes exceeding 50% and, in some places, exceeding 100%.
At lower elevations, mean monthly precipitation (rain and snow-water equivalent) is approximately 10 cm in winter and 3 cm in summer; mean monthly minimum and maximum temperatures, respectively, are 5 • C and 21 • C in summer and −5 • C and 5 • C in winter [54]. At higher elevations, mean monthly precipitation is approximately 15 cm in winter and 4 cm in summer; mean monthly minimum and maximum temperatures are 2 • C and 17 • C in summer and −8 • C and 0 • C in winter. Lower-elevation forests within valleys of the study area are predominantly a mixture of LP, ponderosa pine (Pinus ponderosa Douglas ex P. Lawson & C. Lawson), and white fir (Abies concolor (Gord. & Glend.) Lindl. ex Hildebr.). At approximately 2100 m, there is an ecotone from this mixed forest below to relatively homogenous stands of LP above. At elevations above 2200 m are stands of WP and mixed lodgepole-whitebark pine (LWP) stands. Interspersed throughout the forest at a range of elevations are meadows consisting of herbaceous vegetation with small stands of aspen (Populus tremuloides Michx.). ponderosa pine (Pinus ponderosa Douglas ex P. Lawson & C. Lawson), and white fir (Abies concolor (Gord. & Glend.) Lindl. ex Hildebr.). At approximately 2100 m, there is an ecotone from this mixed forest below to relatively homogenous stands of LP above. At elevations above 2200 m are stands of WP and mixed lodgepole-whitebark pine (LWP) stands. Interspersed throughout the forest at a range of elevations are meadows consisting of herbaceous vegetation with small stands of aspen (Populus tremuloides Michx.).

Vegetation Classification
Forests within the study area were classified into eight categories using modeled 30-m gridded basal area estimates for major tree species [55] and a hierarchical rule-based decision tree (Appendix A). Grid cells with total basal area ≥5 m 2 /ha were defined as forest and classified as dominated by LP, WP, ponderosa pine, or fir species if the corresponding species basal area percentage exceeded 60%. Otherwise, cells were classified as LP-WP co-dominated (LWP) or LP-fir co-dominated if the corresponding combined basal area percentages exceeded 60%. Forest cells that could not be classified according to these rules were classified as No Dominant. Grid cells with total basal area ≤5 m 2 /ha were classified as non-forest. Small treeless areas such as meadows and scree slopes that were not adequately captured by this approach were hand-digitized based on high-resolution imagery from the National Agriculture Imagery Program (NAIP), rasterized at 30-m resolution to match the forest classification grid, and classified as non-forest.

Disturbance Histories and Drought Identification
The disturbance history of the study area since 1985 includes a single-year drought in 2001 and a multi-year drought from 2007 through 2010 based on time series of the Standardized Precipitation Evapotranspiration Index (SPEI) (Figure 3a, [56]). Although MPB is endemic in the study area, U.S. Forest Service annual aerial surveys indicated no major outbreaks of MPB between 1985 and 2000 ( Figure 3b,c, [57]). A severe MPB outbreak began in the early 2000s, affecting over 80% of the study area at its peak in 2007 and 2008. No fire activity has been recorded in the study area since 1985 [58,59], and human disturbances such as logging, mining, development, and mineral and energy extraction are prohibited. Additional details on data sources used to evaluate the disturbance history of the study area are provided in Appendix B. , the mortality spike in the mid-1980s was attributable to Modoc budworm (Choristoneura viridis); the spike in the mid-1990s was attributable to fir engraver (Scolytus ventralis).

Landsat Images
Vegetation conditions were assessed using NDMI, a metric of canopy water content that compares a near-infrared (NIR) spectral band to a mid-infrared (MIR) band: NDMI differs from the more widely used NDVI in its reliance on the MIR band, for which absorbance is primarily determined by vegetation water content [60]. Other remotely-sensed vegetation indices have also been used to monitor drought effects in forests, including NDVI, TCT, and the Enhanced Vegetation Index (EVI) [49,61,62]. Among these indices, however, NDMI had the strongest relationship to field measurements of leaf-area index and canopy-gap fraction in a semiarid forest [51], two forest metrics that are responsive to drought-induced physiological stress and mortality. NDMI also has demonstrated success in quantifying mortality and canopy moisture loss resulting from MPB infestation [50,60].
Gridded 30-m NDMI datasets derived from Tier-1 Landsat Thematic Mapper (TM) 5 images were obtained from the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on-demand interface (Appendix C). Only NDMI data from Landsat 5 TM were used, which limited the time period of analysis to 1985 to 2011, because changes in sensor bandwidth between various Landsat missions may affect interpretation of spectral indices of vegetation [63] and because quantifying and correcting for these changes was beyond the scope of this study. All Landsat images were from August, to limit phenology effects on vegetation spectral characteristics, and had cloud cover ≤50%. ESPA surface reflectance products from which NDMI was derived were radiometrically and atmospherically corrected using the USGS Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) software [64]. Automated masking of clouds, cloud shadows, water, and snow was performed using the CFmask file included with the NDMI datasets. Cells with saturated NDMI values (>1.0) were also removed from subsequent analysis. Based on an integrated assessment of the disturbance history of the study area, two drought years were selected for analysis and identification of possible refugia: 2001 (a single-year drought with no insect disturbance) and 2009 (the third year of a multi-year drought following peak MPB outbreak). The year 2009 was selected to represent the multi-year drought and MPB outbreak because it was 1-2 years after the peak of the MPB outbreak based on aerial survey data, suggesting that most MPB-infested trees in the study area in 2009 would be in the red-attack or gray-attack stage rather than in the earlier green-attack stage, allowing MPB-induced mortality to be adequately captured in NDMI values [47]. The period 1985-2000 was selected as a reference period (a 15-year period without any severe, pervasive disturbance events) against which to compare remotely sensed vegetation characteristics during the drought years 2001 and 2009. Climatic moisture conditions averaged over this reference period approximated the long-term average for the study area and insect disturbance impacts were minimal during the reference period.
The study area was well-suited for the identification of refugia from drought and insect disturbance, as it provided an opportunity to examine and compare two droughts that differed in their magnitudes and disturbance interactions. Because the study area contained both LP and WP stands, it provided an opportunity to assess differences and similarities in how these species responded to drought and insect disturbance, including how refugia for both species related to landscape characteristics and particular landforms. Furthermore, unlike much of the surrounding landscape in the Pacific Northwest of the USA, the study area was free of other disturbance types (fire, development, and natural resource extraction) during the period of analysis that might have otherwise confounded potential identification of refugia from drought and insect disturbance.

Landsat Images
Vegetation conditions were assessed using NDMI, a metric of canopy water content that compares a near-infrared (NIR) spectral band to a mid-infrared (MIR) band: NDMI differs from the more widely used NDVI in its reliance on the MIR band, for which absorbance is primarily determined by vegetation water content [60]. Other remotely-sensed vegetation indices have also been used to monitor drought effects in forests, including NDVI, TCT, and the Enhanced Vegetation Index (EVI) [49,61,62]. Among these indices, however, NDMI had the strongest relationship to field measurements of leaf-area index and canopy-gap fraction in a semi-arid forest [51], two forest metrics that are responsive to drought-induced physiological stress and mortality. NDMI also has demonstrated success in quantifying mortality and canopy moisture loss resulting from MPB infestation [50,60].
Gridded 30-m NDMI datasets derived from Tier-1 Landsat Thematic Mapper (TM) 5 images were obtained from the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on-demand interface (Appendix C). Only NDMI data from Landsat 5 TM were used, which limited the time period of analysis to 1985 to 2011, because changes in sensor bandwidth between various Landsat missions may affect interpretation of spectral indices of vegetation [63] and because quantifying and correcting for these changes was beyond the scope of this study. All Landsat images were from August, to limit phenology effects on vegetation spectral characteristics, and had cloud cover ≤50%. ESPA surface reflectance products from which NDMI was derived were radiometrically and atmospherically corrected using the USGS Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) software [64]. Automated masking of clouds, cloud shadows, water, and snow was performed using the CFmask file included with the NDMI datasets. Cells with saturated NDMI values (>1.0) were also removed from subsequent analysis.

Identification of Refugia
Refugia were identified in drought years 2001 and 2009 within an area consisting of 25,412 30-m grid cells classified as LP, WP, and LWP on Gearhart Mountain in the center of the study area. Refugia were conceptualized as locations within each forest type for which forest responses to drought and MPB infestation were substantially less than for the surrounding landscape. Forest responses to these disturbances were quantified using anomalies from the reference period 1985-2000, following similar approaches in previous studies [49,62]. For the drought years 2001 and 2009, anomalies were calculated as: where A i(t) is the NDMI anomaly for grid cell i in August of year t (t = 2001 or 2009), NDMI i(t) is mean August NDMI for grid cell i in year t (N = 2 Landsat images in each year), and R i is reference NDMI for grid cell i calculated as median August NDMI for grid cell i for the reference period 1985 to 2000 (N = 27 Landsat images). Thus, negative A i(t) values (anomalies) indicated grid cells for which NDMI during drought years was lower relative to the reference (baseline) values.
The study area contained forest types that are MPB hosts (LP and WP) as well as fir-dominated stands that are relatively resistant to MPB infestation. Forest responses to drought even in the absence of insect infestation may vary by stand composition. A common approach is to normalize anomalies by dividing by reference (baseline) standard deviation [62] to enable analysis of anomaly values across spatial gradients, including across different vegetation types. In this study, however, use of absolute (not normalized) anomalies enabled explicit comparisons between different forest types, such as between MPB host versus non-host forest types. Because of canopy-type differences in reference NDMI and NDMI anomalies, refugia identification was stratified by forest type. To identify refugia in each drought year, 10% of cells within each forest type (LP, WP, and LWP) were identified as having the highest (least negative) NDMI anomalies for that forest type. These grid cells indicated areas with the least spectral deviation from reference conditions during the drought in 2001 or combined drought and MPB attack in 2009. For each drought year, a composite map of refugia locations was produced by overlaying refugia for each forest type. In addition, to explore the effects of identifying single, isolated grid cells as refugia, a subset of identified refugia was created representing multi-cell refugia, defined as clumps of two or more spatially contiguous refugial grid cells. All raster-based geospatial analysis was performed in the R statistical environment [65] using the Raster package [66]. All R code and resulting data are available from [67].

Landscape Controls on Refugia
Hypotheses regarding the ecohydrologic determinants of refugia from drought and drought-MPB disturbance ( Table 1) were evaluated by assessing the relationships between refugia locations and landscape characteristics (topographic, soil, and forest composition variables). Because of the size of the study area and the resolution of the analysis (30 m), spatial patterns in climate were best represented by elevation and other topographic variables known to produce distinct microclimates [13,24] rather than by coarser-scale regional macroclimate variables. Gridded 30-m elevation data from the National Elevation Dataset (NED) [68] were used to calculate slope (% rise) and topographic position index (TPI) using a 10-cell (300-m) radius. Topographic heat load index (HLI) and compound topographic index (CTI) at 30-m resolution were obtained from [44]. Modeled 250-m gridded estimates of soil bulk density (kg/m 3 ) at the soil surface (SBD 0cm ) and at 1-m depth (SBD 100cm ) from SoilGrids250 [69] were resampled to 30-m resolution using nearest-neighbor interpolation. To represent forest density, total basal area (m 2 /ha) was calculated for each grid cell as the sum of basal area estimates for LP, WP, ponderosa pine, and fir species. For cells classified as dominated by LP, WP, or co-dominated LWP, the distance to the ecotone with fir forest was calculated as the Euclidean distance to the nearest cell that was classified as fir-dominated or co-dominated LP-fir.
The spatial relationships between refugia locations and landscape characteristics were investigated using boosted regression tree (BRT) models, a machine-learning approach that is well suited to modeling complex nonlinear relationships [70,71]. Explanatory variables were checked for collinearity, which does not bias BRT model results but can make interpretation of the individual effects of collinear predictors more difficult [72,73]. The maximum variance inflation factor (VIF) was 5.2 for elevation; all other VIFs were <3.5. Only one variable pair (elevation and distance to the fir ecotone) had a Spearman's rank correlation (ρ) ≥0.65. Because no extreme collinearity was detected in the explanatory variables, and because of the potential for elevation and distance to the fir ecotone to exert mechanistically distinct influences on forest disturbance and presence of refugia (Table 1), all explanatory variables were used in BRT modeling. Cross-validated BRT models were developed and calibrated independently for each combination of drought year (2001 and 2009) and forest type (LP, WP, and mixed LWP) with refugium occurrence as a binary response. To assess the effects on model performance of removing single-cell refugia from analysis, two versions of each model were independently calibrated, and the results compared: one with all identified refugia as the response and one with only multi-cell refugia as the response.
To aid in inference, final BRT models of multi-cell refugia were bootstrapped using 20 iterations (independent model runs) with a random subsample of 2000 grid cells for each model. Models allowed for natural variability in landscape characteristics within multi-cell refugia by using 30-m grid cells as units of analysis and not averaging landscape characteristics within refugia. For all BRT model runs, the tree complexity was 5 and bag fraction was 0.5. The gbm.step procedure was used to test a range of learning rates. This procedure uses 10-fold cross-validation to identify an optimal number of trees by minimizing predictive deviance as described in [71]. Model calibration involved selecting an optimal learning rate based on (a) maximized predictive accuracy (area under the curve of the Receiver Operating Characteristic, AUC-ROC), (b) minimized predictive deviance, (c) a steep increase in predictive deviance after the minimum was reached, to avoid overfitting, and (d) a number of trees between 1000 and 5000 [71,74,75]. All BRT modeling was performed in the R statistical environment [65] using the gbm [76] and dismo [77] packages.
Model predictive ability was evaluated using AUC-ROC, with higher values indicating better ability to predict both presence and absence of refugia [16]. Values of AUC-ROC can range from 0.5 (predictive ability no better than random) to 1.0 (perfect prediction). Model explanatory ability was assessed based on percent variance of the response variable (refugium occurrence or non-occurrence in a grid cell) that was explained by the model. Because BRT models with AUC-ROC <0.75 are considered insufficiently accurate in their predictive performance [75], only models with mean AUC-ROC ≥0.75 were retained for interpretation of landscape controls on refugia locations to avoid misleading conclusions based on poorly performing models. The relationships between landscape characteristics and refugia locations were evaluated based on relative influence and on the shapes of partial-dependence (marginal response) curves [70]. Relative influence is a metric of each explanatory variable's importance in determining occurrence of refugia. Partial-dependence plots use fitted BRT models to depict the marginal effect of each explanatory variable on the predicted prevalence of refugia.

Temporal and Spatial Patterns of Drought and Insect Disturbance
During the reference period from 1985 to 2000, NDMI for all forest types was fairly stable (Figure 3). Reference NDMI showed considerable spatial variation ( Figure 4a) and was generally higher in fir and LP-fir forest and lower in areas of LP and WP forest. Across forest types, NDMI declined somewhat during the 2001 drought (a single-year drought without insect disturbance) and then quickly rebounded (from 2002 to 2005) to levels similar to those observed during the reference period ( Figure 3). NDMI anomalies in 2001 were lowest in fir and LP-fir stands, but were relatively modest for all forest types with mean NDMI anomalies >−0.1 (Figure 4).
By contrast, NDMI in 2009 (the third year of a multi-year drought that coincided with a severe MPB outbreak) showed dramatic declines from the reference period that varied across forest types. In the LP and WP forest types dominated by MPB hosts, NDMI declined sharply beginning around 2007 ( Figure 3) and mean anomalies in 2009 were <−0.1 ( Figure 4). In fir forest, NDMI declines were observed at the lower end of the NDMI gradient, but at the higher end NDMI remained fairly constant and similar to the reference period (see 25th and 75th percentile lines in Figure 3), resulting in less extreme NDMI anomalies in 2009 ( Figure 4).  Figure 4). In fir forest, NDMI declines were observed at the lower end of the NDMI gradient, but at the higher end NDMI remained fairly constant and similar to the reference period (see 25th and 75th percentile lines in Figure 3), resulting in less extreme NDMI anomalies in 2009 ( Figure 4).

Landscape Controls on Refugia
In each year, 2544 grid cells (10% of the model area grid cells, by definition) were identified as refugia (Appendix D). The spatial overlap between years was 28%, i.e., of grid cells identified as

Landscape Controls on Refugia
In each year, 2544 grid cells (10% of the model area grid cells, by definition) were identified as refugia (Appendix D). The spatial overlap between years was 28%, i.e., of grid cells identified as refugia in 2001, 28% were also refugia in 2009. Approximately 76% and 91% of refugial grid cells belonged to multi-cell refugia in 2001 and 2009, respectively. Of grid cells belonging to multi-cell refugia in 2001, 30% also belonged to multi-cell refugia in 2009. For all forest types in both years, BRT models of landscape controls on only multi-cell refugia showed improved predictive and explanatory performance over models that also included single-cell refugia, as indicated by AUC-ROC and percent deviance explained, respectively (Appendix E). Consequently, final BRT models used only the multi-cell refugia to evaluate the relative influence of landscape characteristics and the shapes of the relationships between landscape characteristics and prevalence of refugia.
For all forest types, BRT model explanatory and predictive abilities were greater for refugia in 2009 than in 2001 (Table 2). BRT models explained from 53% to 66% of deviance in occurrence of refugia in 2009 but only 33% to 40% in 2001. Mean AUC-ROC values for 2009 models ranged from 0.85 to 0.88, indicated very good ability to predict presence and absence of refugia [16], whereas mean AUC-ROC values for 2001 models ranged from 0.71 to 0.77. Two models from 2001 had mean AUC-ROC <0.75 (LWP and WP; AUC-ROC = 0.71 and 0.74, respectively), and thus were considered insufficiently accurate in their predictive performance [75]; these models were not retained for interpretation. The other four models (LP 2001, LP 2009, LWP 2009, and WP 2009) were retained and used to interpret landscape controls on refugia locations. The LP 2001 model had mean AUC-ROC = 0.77, which was marginally above the threshold for model retention. All models for 2009 were considered highly accurate in their predictions of refugia (all AUC-ROC ≥0.85; [16]). The four retained models explained from approximately 40% to 66% of deviance in occurrence of refugia and demonstrated cross-validated correlations ranging from 0.3 to 0.51 (Table 2). Table 2. Boosted regression tree model performance statistics and variable relative influence relating occurrence of multi-cell refugia to landscape characteristics. Landscape characteristics with greatest relative influence included total basal area, SBD 100cm , elevation, slope, and heat load index ( Table 2). SBD 0cm , topographic position index, compound topographic index, and distance to the fir ecotone had somewhat lower relative influence. Percent fir had minimal relative influence across all models.

BRT Model Parameters
The relationships between landscape characteristics and probability of refugium occurrence were depicted in partial-dependence plots overlain from each of the four retained BRT models ( Figure 5; see Appendix F for plots from individual bootstrapped model runs). Modeled relationships were typically non-monotonic, however, general patterns could be discerned by examining ranges of each landscape variable for which probability of identifying refugia was relatively high or low (indicated by fitted function values on the vertical axes of Figure 5). Distributions of landscape characteristics within the study area (indicated by gray histograms in Figure 5) were used to guide interpretation of partial dependence plots, to avoid interpretation of outlier effects represented in the horizontal ranges of partial-dependence plots that were much greater than or much less than the central distribution of each landscape characteristic. each landscape variable for which probability of identifying refugia was relatively high or low (indicated by fitted function values on the vertical axes of Figure 5). Distributions of landscape characteristics within the study area (indicated by gray histograms in Figure 5) were used to guide interpretation of partial dependence plots, to avoid interpretation of outlier effects represented in the horizontal ranges of partial-dependence plots that were much greater than or much less than the central distribution of each landscape characteristic.

Total Basal Area
In all models, refugium occurrence generally decreased with increasing stand density up to approximately 40 m 2 /ha (Figure 5a) in agreement with the hypothesized negative relationship between stand density and refugia (Table 1). Although partial-dependence plots also showed

Total Basal Area
In all models, refugium occurrence generally decreased with increasing stand density up to approximately 40 m 2 /ha (Figure 5a) in agreement with the hypothesized negative relationship between stand density and refugia (Table 1). Although partial-dependence plots also showed relatively high refugium occurrence at high stand densities (i.e., >60 m 2 /ha), those portions of plots should be interpreted with caution due to relatively few grid cells with such high stand densities (see histogram in Figure 5a).

Soil Bulk Density
Despite some disagreement across models, partial-dependence plots generally suggested an overall negative relationship between refugium prevalence and soil bulk density, in support of the hypothesized relationship (Table 1). In the LP and LWP forest models, SBD 100cm was an important control on refugium occurrence (Table 2), with refugia occurring preferentially in areas of lower SBD 100cm (Figure 5b). In the WP 2009 model, however, SBD 100cm demonstrated lower relative influence and a more complex non-monotonic relationship with refugium occurrence. The relative influence of SBD 0cm was less than that of SBD 100cm for the LP and LWP models ( Table 2). The hypothesized negative relationship between refugia and SBD 0cm was observed in the LP 2009 and WP 2009 models, but not the LP 2001 or LWP 2009 models (Figure 5f).

Slope and Heat Load Index.
Although partial-dependence plots for all 2009 models showed refugia occurring preferentially on the steepest slopes (>50% rise), only approximately 5% of grid cells in the study area were this steep (Figure 5c). For slopes between 0 and 30% rise (representing approximately 80% of the study area grid cells), the shapes of the relationships between slope and refugia differed across models making overall conclusions more difficult. Across all models, refugium prevalence generally decreased with HLI across the primary range of the HLI distribution in the study area (0.5 ≤ HLI ≤ 0.9), indicating that refugia were most likely to occur in the most topographically shaded areas (Figure 5e), as hypothesized (Table 1). Although slope is used in the calculation of HLI, slope and HLI in this study area were weakly correlated (ρ = 0.08) indicating that slope and HLI played non-redundant roles in relation to refugia.

Elevation
Elevation was a moderately important control on refugium occurrence in both LP and LWP forests (Table 2), however, the shapes of relationships with elevation differed somewhat across models (Figure 5d). In all models from 2009, refugia were generally found at lower elevations, e.g., more common around 2200 m than around 2400 m, contradicting the hypothesized relationship (Table 1). However, the opposite pattern was observed in the LP 2001 model, where refugium probability increased with elevation.

Topographic Position Index
The hypothesized occurrence of refugia in topographically concave areas such as valley bottoms (low TPI grid cells) was observed for all forest types in 2009 across most of the distribution of TPI values in the study area, i.e., from approximately −20 to 40 (Figure 5g). Conversely, in LP forest in 2001, refugia occurred preferentially in topographically more convex (higher TPI) areas, such as along ridgetops.

Compound Topographic Index
Refugia were generally most likely to occur in grid cells with the highest and lowest CTI values (Figure 5h). Refugia associated with relatively high CTI values, indicating areas of topographic convergence, agreed with the hypothesized relationship based on surface water routing and soil characteristics (Table 1). Notably in this study area, CTI values >8 were visibly aligned with headwater streamlines from the high-resolution (1:24,000) National Hydrography Dataset (NHD). Models for refugia in WP and LWP in 2009, but not the LP models, showed thresholds near CTI = 8, possibly suggesting refugia associated with riparian areas along headwater streams. However, all models also showed relatively high probability of refugium occurrence at CTI <5, corresponding to upper slopes and ridgetops, in contrast to the hypothesized relationship.

Distance to the Fir Ecotone
For all models, refugium prevalence generally declined with increasing distance to the fir ecotone, with refugia most likely to be found within~100 m of the ecotone (Figure 5i), however, this geographic characteristic had generally low relative influence ( Table 2).

Landforms Associated with Refugia
Visual examination of multi-cell refugia showed that some were associated with particular landforms (Figure 6). Prominent landforms associated with refugia from the 2009 drought and MPB outbreak included steep north-and northeast-facing slopes (Figure 6b through 6d), including upper slopes at the heads of cirques (Figure 6c) and adjacent to riparian areas (Figure 6d). Other aggregations of refugial cells were located in upper valley positions in the headwaters of streams (upstream ends of high-resolution NHD flowlines). A relatively large aggregation of refugial grid cells was observed on a moderately sloping, northwest-facing hillslope in LP forest in the northwest of the area in which refugia were identified and modeled (dashed red box in Figure 6a). These refugia were not associated with any prominent topographic features, but were associated with low values of SBD 100cm . Within these refugia, SBD 100cm ranged from approximately 1350 to 1425 kg/m 3 and was significantly lower than the distribution SBD 100cm for the entire modeling area (Kolmogorov-Smirnov two-sample test: p < 0.001).

Discussion
In this study of a dry conifer forest, refugia from drought and insect disturbance identified from multi-year Landsat imagery were associated with a suite of topographic, soil, and vegetation characteristics. However, forest drought effects such as tree physiological stress, tree mortality, triggering of insect outbreaks may differ depending on drought magnitude and severity [1,30], as illustrated in this study by differences between 2001 and 2009 in forest responses (Figures 3 and 4) and model performance ( Table 2). During a multi-year drought from 2007 to 2010-which Figure 6. Examples of landforms associated with refugia. Refugia locations mapped in (a) included, framed in solid red boxes from north to south, (b) a steep north-facing slope in lodgepole pine forest, (c) a steep and north-facing slope in lodgepole-whitebark pine forest at the head of a cirque, and (d) a steep and northeast-facing slope in whitebark and lodgepole-whitebark pine forest adjacent to a riparian area. Refugia locations in the dashed red box were not associated with any discernible topographic features but were associated with low soil bulk density at 1-meter depth. All imagery is from the National Agriculture Imagery Program (NAIP), acquired 29 June, 2016. Streamlines are from the National Hydrography Dataset (NHD).

Discussion
In this study of a dry conifer forest, refugia from drought and insect disturbance identified from multi-year Landsat imagery were associated with a suite of topographic, soil, and vegetation characteristics. However, forest drought effects such as tree physiological stress, tree mortality, triggering of insect outbreaks may differ depending on drought magnitude and severity [1,30], as illustrated in this study by differences between 2001 and 2009 in forest responses (Figures 3 and 4) and model performance ( Table 2). During a multi-year drought from 2007 to 2010-which contributed to a severe coincident MPB outbreak in LP and WP forest-refugia from combined drought-insect disturbance were well-predicted and moderately well-explained by landscape characteristics. By contrast, refugia during a single-year drought in 2001 that did not trigger any insect outbreak were not as well explained or predicted by landscape characteristics, and poor model performance for WP and LWP refugia in 2001 rendered these models insufficient for interpretation. A possible explanation is that WP forest impacts from a single-year drought without insect disturbance were not sufficiently severe to manifest the full expression of topographically determined refugia, at least not refugia discernible by remotely sensed NDMI.

Landscape Controls on Refugia from Drought and Mountain Pine Beetle
Modeled relationships generally supported the hypotheses that refugia from drought and MPB impacts would be associated with topographically shaded slopes, convergent areas such as valleys, areas of lower soil bulk density, thinner forest stands (areas of lower total basal area), and areas closer to the ecotone with fir forest (Table 1, Figure 5). However, some modeled relationships did not conform to expectations, such as prevalence of refugia at lower elevations in 2009 models and in areas of low CTI in all models. In some cases, modeled relationships differed somewhat among forest types and included complex, multimodal responses ( Figure 5), suggesting that there may be multiple interacting physical and biological processes promoting refugia that may be sensitive to forest stand characteristics such as dominant species. These findings are consistent with a conceptual model (Figure 1) in which the spatial distribution of refugia is influenced not only by physical characteristics of the landscape but also by forest structure and biological traits of individuals, populations, and species.
The association between refugia and topographically shaded slopes in this study is consistent with previous studies showing greater drought-induced tree mortality and bark-beetle attack rates on warmer south-facing slopes than cooler north-facing slopes [33,50,78,79]. On shaded slopes, delayed snowmelt and reduced evaporative demand may support hydrologic refugia in the form of growing-season soil moisture reserves [15,22,24,35,42,80]. Soil water retention curves have demonstrated that topographically shaded slopes retain up to 25% more soil water at a given soil water pressure than nearby poleward-facing slopes [23]. Cooler temperatures in topographically shaded microenvironments might also decrease MPB brood survival [43].
The hypothesis that refugia from drought and MPB attack would be found primarily at higher elevations was not supported by this study; indeed, the WP and LWP 2009 models suggested the opposite. Traditionally, high elevations have been considered to impose thermal limitations on MPB brood survival and make completion of univoltine (single-year) life cycles more difficult [36]. However, thermal limitations may be reduced by regional warming, contributing to concerns that high elevations alone may be insufficient to protect WP from MPB-induced mortality. Indeed, univoltine populations have been observed as high as 3000 m [6,43,50], which is almost 500 m higher than the highest elevations in this study area.
Previous studies have found MPB attack rates and tree mortality to be positively associated with slope [33], negatively associated [41], or to have no relationship to slope [9]. Models from this study indicated the greatest probability of refugia on very steep slopes, but these were relatively rare in the study area. In some cases, very steep, topographically shaded slopes immediately below ridgelines ( Figure 6) might accumulate snow drifts in winter and maintain snowpack and soil moisture longer into the growing season compared to other topographic positions [22,24,38,39]. Across the range of more subdued slope angles present in most of the study area, refugia models differed across forest types, producing no compelling evidence for overall positive or negative relationship between refugia and slope angle.
Modeled relationships for the 2009 drought-MPB disturbance suggested prevalence of refugia in topographically concave areas such as valley bottoms for all forest types and in riparian areas along headwater streams in WP and LWP stands. Similarly, topographically sheltered coves and riparian areas were recently shown to provide cool, moist refugia for limber pine (Pinus flexilis James) [42]. Cold air pooling in coves, canyons, and valleys may provide pockets of high relative humidity and promote formation of dew [15,24], which can help relieve plant water stress during droughts [81,82]. Suppression of evapotranspiration by high relative humidity can have substantial ameliorating effects on drought mortality [83]. In addition, wind can redistribute precipitation in rugged terrain preferentially to valleys [84,85]. Lower nocturnal minimum temperatures associated with cold-air pools [86] might negatively affect MPB habitat suitability and contribute to refugia from insect damage. Although cooler temperatures might provide more favorable conditions for mature trees, it should be noted that survival of WP seedlings may be depressed in frost pockets affected by cold air drainage and in moist sites where vegetative competition and animal disturbance may be higher [87]. Notably, modeled relationships from this study also suggested that refugia could be associated with topographically convex environments, such as along upper slopes and ridgelines.
Soil bulk density is inversely related to porosity and an important determinant of water-holding capacity [45,88]. Trees growing in low bulk-density soils may also have greater drought tolerance owing to more extensively developed root networks [45,46]. Relative influence ( Table 2) and partial-dependence plots (Figure 5b,f) suggest that in LP and LWP forests, low soil bulk density at 1-m depth might have a greater role in supporting refugia than low soil bulk density at the surface, with refugia preferentially identified in areas with SBD 100cm < 1450 kg/m 3 . In WP forest, refugia occurred preferentially in areas of low soil bulk density at the soil surface (Figure 5f) but the shape of the relationship with SBD 100cm appeared multimodal (Figure 5b), suggesting a need for improved understanding of the role of soil physical characteristics in shaping refugia in WP ecosystems.
In all models from this study, refugia were associated with areas of low total basal area, supporting the notion that overall tree vigor leading up to major disturbances may be suppressed in dense stands due to increased competition, and that low vigor may predispose trees to drought mortality and MPB attack [30,89]. For example, [9] reported thresholds of 10 and 18 m 2 /ha below which WP and LP stands, respectively, were not attacked by MPB. Fire suppression can result in densely stocked stands that may contain more mature larger-diameter host trees, making them more susceptible to MPB attack [29]. Thinner stands also increase wind penetration, helping to disperse beetle pheromones and disrupt chemical communications needed to coordinate attacks [33]. Current management recommendations in LP and ponderosa pine stands include thinning to improve tree vigor and decrease habitat suitability for bark beetles [90]. Even in the absence of insect disturbance, increased competition for limited soil water during multi-year droughts (e.g., in 2009) could make denser stands more vulnerable to drought mortality [51] and thus less likely to support drought refugia. Evidence from thinning treatments in LP stands suggests the potential to increase snowpack accumulation, possibly contributing to greater growing-season soil moisture reserves, although these effects may depend on the spatial configuration of residual trees [91]. Although models from this study suggested that refugia were more likely to be located within 200 m of the ecotone with fir forest, this effect was small based on relative influence values and might be affected by reduced accuracy of the forest-type map used to classify grid cells near this ecotone. For WP especially, areas near the fir ecotone might not be optimal locations for refugia because of competing tree species and abundance of seed predators in fir forests [92].

Directions for Future Research
Confidence in disturbance refugia identified by remote sensing could be improved by field validation. In particular, misclassification of forest types could produce "false positive" refugia, whereby local errors in the tree-species maps used as inputs to analysis could cause areas dominated by non-host trees to be mistakenly identified as refugia from MPB outbreaks. For example, examination of high-resolution aerial imagery suggests that the forest stand in Figure 6b may differ in species composition from its surroundings, calling into question the accuracy of the refugium identified at this site and highlighting a need for field validation. Additionally, field investigations could discern whether isolated single-cell refugia grid cells represented actual microrefugia (very small patches of muted forest response to disturbance) or were instead artifacts of Landsat data processing. In some cases, field assessments might clarify the biophysical bases supporting refugia. For example, field observations of snow accumulations along topographically shaded upper slopes identified as LP refugia would lend credence to the speculation that topographic controls on snowmelt and soil moisture dynamics helped account for these refugia. This study used a 10% threshold on NDMI anomaly to delineate refugia. Threshold choice might be improved by a sensitivity analysis across a range of delineation thresholds and a variety of vegetation indices (e.g., NDVI, TCT, and EVI), calibrated to field observations of disturbance refugia.
Field observations could also provide critically important information to contextualize refugia identification within larger patterns and processes of forest ecology. Such information could include (1) local prevalence of disease agents such as blister rust and western gall rust (causal agent Endocronartium harknessii) and parasites such as dwarf mistletoe (Arceuthobium species) within refugia and within the larger forest system; (2) population genetics with emphasis on the traits affecting long-term population viability such as genetic diversity and disease resistance; and (3) understory population demographics within identified refugia, including prevalence of seedlings and saplings to support regeneration. Factors controlling drought and insect-induced mortality are less understood in WP forests than LP forests and may not be directly transferable from the LP context [31]. For example, although WP is generally more drought-tolerant than LP [93], its defense capacities to MPB may be less [94] and there may be different relationships between precipitation patterns and tree mortality from those observed in LP forests [7]. In addition, because WP is dependent on Clark's nutcracker (Nucifraga columbiana) for dispersal, stand characteristics associated with nutcracker seed caching-e.g., cone production capacity and pre-dispersal cone predation by American red squirrel (Tamiasciurus hudsonicus)-are important attributes that may determine whether putative refugia are capable of sustaining WP populations long term [10,92].
Future investigations could attempt to quantify spatial change or stability in locations of refugia between disturbance events or as multi-year disturbances (droughts and insect outbreaks) progress. This study identified refugia in 2009 to represent a multi-year disturbance event (drought and MPB outbreak from 2007 to 2010). However, the relationships between MPB attack rates and topography can change over the course of an outbreak [5,50]. Refugia that persist through multi-year disturbances or multiple disturbance events may be of special interest to natural resource managers. Thus, spatial analysis of refugia locations across diverse disturbance types would be helpful, including comparison of mild versus severe droughts, and of droughts with and without insect outbreaks. In this study, the modest (~30%) spatial overlap of refugia in 2001 and 2009 was unsurprising given the differences in drought magnitudes and disturbance interactions between these two drought years. This finding highlights the importance of differentiating between droughts and insect outbreaks in refugia identification-although these disturbances are often causally linked [3,4], they can and do occur independent of each other and can produce different ecological consequences. The relationships between landscape characteristics and refugia locations in LP forest agreed between the 2001 and 2009 models for some variables (i.e., refugia in both years were associated with thinner forest stands, lower SBD 100 , and lower HLI), but not others (e.g., topographic position index). Such comparisons were not possible for WP and LWP forest due to inadequate performance of 2001 models. Importantly, comparison of only two disturbance years is insufficient to assess the long-term spatial stability of disturbance refugia. In general, robust analysis of refugia stability using these methods would be difficult in most landscapes due to the limited temporal extent of remote-sensing archives relative to disturbance recurrence intervals [95].
Another important area of investigation is how disturbance refugia function in heavily managed forests, including areas subject to multiple interacting types of disturbance such as fire, timber harvest, droughts, and insect outbreaks. This study was intentionally conducted in an area free from recent burns or harvests, however, these forest disturbances are widespread in the western USA, including in the national forest land immediately surrounding the study area [59]. Disturbance legacy effects can influence ecological responses to subsequent disturbances with potential for complex, cascading effects [96]. Increasing attention to fire refugia has shown the importance of unburned areas within large wildfires to post-fire regeneration and biodiversity [97]. Further research is needed to understand (1) how refugia from one disturbance type (e.g., fire) shape forest responses to subsequent disturbance (e.g., insect outbreaks), (2) whether certain landscape characteristics could support refugia from multiple forms of disturbance, and (3) how forest management practices-such as harvest scheduling, intensity, and spatial patterns-affect disturbance refugia.

Conclusions
In LP and WP forest in a mountainous study area that experienced a single year drought without insect disturbance (2001) and a multi-year drought combined with a severe MPB outbreak (2007-2010), refugia from drought and insect disturbance were identified using a simple anomaly-based approach and time-series of remotely sensed NDMI. BRT models of landscape characteristics (topography, soil, and forest stand variables) performed better in predicting and explaining the locations of refugia in LP and WP stands in 2009 than 2001, possibly because the 2009 disturbance was much more severe in these forest types. In general, refugia from drought and MPB impacts were found on topographically shaded slopes, in convergent areas such as valleys, in areas of low soil bulk density and in thinner forest stands. These findings suggest a variety of physical and biological processes that may interact to create refugia from drought and insect disturbance-possibly including spatial variations in evaporative demand, snow accumulation and melt timing, soil water storage capacity, and drought-tolerance abilities at the stand level-and that may function differently in WP than in LP forests. Future efforts to identify and characterize refugia from drought and insect outbreaks could benefit from the integration of field observations with remotely sensed vegetation indicators and from tracking refugia over time as multi-year disturbances progress.
Improved understanding of the landscape controls on refugia from drought and insect disturbance can help support effective management of threatened ecosystems-such as high-elevation five-needle pine forests-that may face increased stress in an era of rapid climate change. In the western USA and Canada, WP conservation and restoration initiatives range in scales from individual trees to large landscapes, employing approaches such as thinning, prescribed fire, cutting of shade-tolerant competitors, and targeted planting of rust-resistant seedlings [10]. Such efforts rely on information including disturbance histories, population demographics and genetics, and site characteristics including accessibility and desirability for recreation. In planning for future droughts and insect outbreaks, the likelihood of a given site to function as a WP refugium from these disturbances could be an additional consideration in WP restoration planning and prioritization.  tree species in the study area: lodgepole pine (LP), whitebark pine (WP), ponderosa pine, and fir species. Species-level basal area estimates were not available for individual fir species, so an aggregate fir-species dataset was used. Although other tree species such as aspen and mountain hemlock are present within the study area, their limited spatial extents are unlikely to have substantially affected forest-type classification for most of the study area, as indicated by a coarser-resolution (250 m) forest type map [98].
For each grid cell, total basal area (m 2 / ha) was calculated as the sum of the four tree-specific basal area estimates and was used to calculate the percent basal area for each tree species. Cells were classified into one of eight forest types according to a hierarchical rule-based decision tree with three criteria (Table A1). Cells with total basal area ≤5 m 2 /ha were classified as non-forest. This threshold was chosen based on visual comparison of total basal area estimates with high-resolution imagery from the National Agriculture Imagery Program (NAIP). All other cells were classified as dominated by a single tree species if that species basal area percentage exceeded 60%. Cells that could not be classified according to these rules were then classified as LP-WP co-dominated (LWP) if the combined LP and WP basal area percentages exceeded 60%, LP-fir co-dominated if the combined LP and fir basal area percentages exceeded 60%, or else as No Dominant.
In addition, small meadows and other treeless areas that were not adequately captured by the basal-area approach to forest classification were hand-digitized based on high-resolution NAIP imagery. These treeless areas were rasterized at 30-m resolution and classified as non-forest. 12-month antecedent conditions) was computed for the study area from 1985 through 2011. These SPEI-12 values were averaged across June, July, and August to produce a summer SPEI-12 value for each year. Summer SPEI-12 fluctuated over the period 1985-2000, with a relatively dry period followed by a wetter period, without any major droughts, i.e., all SPEI-12 >−1.25 (Figure 3a). Mean summer SPEI-12 over 1985-2000 was 0.07 (i.e., near zero), indicating that on average over this period, moisture conditions approximated the long-term average for this study area. Summer SPEI-12 values <−1.

Appendix B.2. Insect and Disease Outbreaks
Annual aerial insect and disease survey data were obtained from the USDA Forest Service, Pacific Northwest Region, Forest Health Protection program [57]. Each yearly dataset includes polygons of mapped insect and disease disturbance with an accompanying agent code representing the most likely insect or disease agent and an estimate of severity in trees per acre. Each yearly polygon dataset was rasterized at 30-m resolution and the areal percent of the study area affected by mapped disturbance agents was calculated along with the mean severity (in trees per hectare) of affected areas. These annual values were compiled to represent time-series of mortality extent and severity in the study area (Figure 3b,c, respectively).
The insect and disease disturbance history of the study area was characterized by three outbreaks: (1) a mild outbreak of Modoc budworm (Choristoneura viridis) in the mid-1980s affecting~50% of the study area at fewer than 0.5 trees per hectare, (2) a mild outbreak of fir engraver (Scolytus ventralis) in the mid-1990s affecting~25% of the study area at fewer than 2 trees per hectare, and (3) a severe outbreak of Mountain pine beetle (Dendroctonus ponderosae; MPB) beginning in the early 2000s which at its peak in 2007 and 2008 affected approximately 80% of the study area at greater than 12 trees per hectare. MPB is endemic to the study area, however, from 1979 to 2004 MPB outbreaks never affected more than 10% of the study area and were generally at mild levels of severity (typically affecting fewer than 5 trees per hectare).
The single-year drought (2001) occurred in the absence of any insect outbreaks, whereas the multi-year drought (2007-2010) coincided with a severe and widespread MPB outbreak. The 15-year reference period selected for this study area (1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000) contained two insect outbreaks, however, both outbreaks were mild in severity compared to the MPB outbreak in the late 2000s and neither outbreak was associated with a dramatic decline in NDMI as was the case for the MPB outbreak. this period were of mild severity in comparison to the later MPB outbreak in the late 2000s (Figure 3b,c), and (c) NDMI trajectories through time were relatively stable (Figure 3d through Figure 3f).        To designate the type of refugia used in each model, All indicates that the response variable was the occurrence of refugia using all identified refugial grid cells (both single-cell and multi-cell refugia); Multi indicates that the response variable was the occurrence of multi-cell refugia only. Each model was calibrated independently to determine the optimal learning rate (see methods in main text). AUC-ROC is the area under the curve of the receiver operator characteristic. For the purpose of comparing single-cell and multi-cell refugia models, bootstrapping was not conducted; instead, all available grid cells were used in each model (for this reason, model parameters vary slightly from those reported in Table 2).