Next Article in Journal
Converting Larch Plantations to Larch-Walnut Mixed Stands: Effects of Spatial Distribution Pattern of Larch Plantations on the Rodent-Mediated Seed Dispersal of Juglans mandshurica
Next Article in Special Issue
Whitebark and Foxtail Pine in Yosemite, Sequoia, and Kings Canyon National Parks: Initial Assessment of Stand Structure and Condition
Previous Article in Journal
A General Leaf Area Geometric Formula Exists for Plants—Evidence from the Simplified Gielis Equation
Previous Article in Special Issue
Managing Wildfire for Whitebark Pine Ecosystem Restoration in western North America
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Jennifer Cartwright
U.S. Geological Survey, Lower Mississippi-Gulf Water Science Center, 640 Grassmere Park, Suite 100, Nashville, TN 37211, USA
Forests 2018, 9(11), 715; https://doi.org/10.3390/f9110715
Submission received: 24 October 2018 / Revised: 14 November 2018 / Accepted: 16 November 2018 / Published: 18 November 2018
(This article belongs to the Special Issue Ecology and Restoration of Whitebark Pine)

Abstract

:
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 identified using remotely sensed vegetation indicators. Refugia were identified 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. Identification 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.

1. 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), and successional replacement due to fire suppression [9,10]. Based on these threats and its 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].
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].
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.

2. Materials and Methods

2.1. Study Area

This study was conducted in the Gearhart Mountain Wilderness, an approximately 95-km2 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.).

2.2. 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 m2/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 m2/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.

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

2.4. 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 = (NIR − MIR) / (NIR + MIR)
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.

2.5. 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:
Ai(t) = NDMIi(t) − Ri
where Ai(t) is the NDMI anomaly for grid cell i in August of year t (t = 2001 or 2009), NDMIi(t) is mean August NDMI for grid cell i in year t (N = 2 Landsat images in each year), and Ri 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 Ai(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].

2.6. 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/m3) at the soil surface (SBD0cm) and at 1-m depth (SBD100cm) from SoilGrids250 [69] were resampled to 30-m resolution using nearest-neighbor interpolation. To represent forest density, total basal area (m2/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.

3. Results

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

3.2. 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).
Landscape characteristics with greatest relative influence included total basal area, SBD100cm, elevation, slope, and heat load index (Table 2). SBD0cm, 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.
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.

3.2.1. Total Basal Area

In all models, refugium occurrence generally decreased with increasing stand density up to approximately 40 m2/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 m2/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).

3.2.2. 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, SBD100cm was an important control on refugium occurrence (Table 2), with refugia occurring preferentially in areas of lower SBD100cm (Figure 5b). In the WP 2009 model, however, SBD100cm demonstrated lower relative influence and a more complex non-monotonic relationship with refugium occurrence. The relative influence of SBD0cm was less than that of SBD100cm for the LP and LWP models (Table 2). The hypothesized negative relationship between refugia and SBD0cm was observed in the LP 2009 and WP 2009 models, but not the LP 2001 or LWP 2009 models (Figure 5f).

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

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

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

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

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

3.3. 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 SBD100cm. Within these refugia, SBD100cm ranged from approximately 1350 to 1425 kg/m3 and was significantly lower than the distribution SBD100cm for the entire modeling area (Kolmogorov–Smirnov two-sample test: p < 0.001).

4. 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 (Figure 3 and Figure 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.

4.1. 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 SBD100cm < 1450 kg/m3. 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 SBD100cm 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 m2/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].

4.2. 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 SBD100, 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.

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

Funding

This research was funded by the U.S. Geological Survey and the Department of the Interior Northwest Climate Adaptation Science Center.

Acknowledgments

The author wishes to thank Timothy Assal and Scott Worland (U.S. Geological Survey), and Meg Krawchuk and Garrett Meigs (Oregon State University) for their helpful feedback on the study design. This manuscript was improved based on reviews by Jeffrey Hicke (University of Idaho) and anonymous reviewers. Arjan Meddens (University of Idaho) and Mary Frances Mahalovich (USDA Forest Service) provided helpful comments on an early draft. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest

The author declares no conflict of interest.

Appendix A. Forest-Type Classification

Modeled 30-m gridded basal area estimates were obtained from the U.S. Department of Agrigulture (USDA) Forest Service, Forest Health Technology Enterprise Team [55] for the predominant 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 (m2/ 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 m2/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.
Table A1. Criteria for forest-type classification based on 30-m gridded basal area estimates for lodgepole pine, ponderosa pine, whitebark pine, and fir species.
Table A1. Criteria for forest-type classification based on 30-m gridded basal area estimates for lodgepole pine, ponderosa pine, whitebark pine, and fir species.
Criterion 1Criterion 2Criterion 3Forest-type Classification
Total basal area ≤5 m2/ha--Non-forest
Total basal area >5 m2/ha% Lodgepole >60%-Lodgepole-dominated
% Ponderosa >60%-Ponderosa-dominated
% Whitebark >60%-Whitebark-dominated
% Firs >60%-Fir-dominated
No species with >60%Lodgepole + Whitebark >60%Lodgepole–Whitebark co-dominated
Lodgepole + Firs >60%Lodgepole-fir co-dominated
Does not meet any of above criteriaNo dominant

Appendix B. Disturbance Histories of the Study Area

The disturbance histories in the Gearhart Mountain Wilderness study area (droughts, insect and disease outbreaks, and fires) were assessed using multiple data sources. Because the time period of Landsat data availability was from 1985 to 2011 (see methods section), assessment of disturbance histories focused on this time period. Because the study area has been a designated wilderness area since 1964, human disturbances such as logging, mining, development, and energy and mineral extraction are prohibited and were assumed not to have occurred. An integrated assessment of disturbance history was used to select two drought years for the identification of refugia (2001 and 2009) and a 15-year reference period (1985–2000) against which to compare the selected drought years.

Appendix B.1. Droughts

Drought history was assessed using 0.5-degree gridded Standardized Precipitation Evapotranspiration Index (SPEI) [56], which accounts for both precipitation and potential evapotranspiration in determining drought intensity. SPEI is similar to the commonly used standardized precipitation index (SPI) except that it also incorporates the effects of increased temperatures on water demand. On a monthly time step, median SPEI-12 (SPEI calculated with 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.25 were observed in 2001, 2007, and 2009. Preliminary analysis of the SPEI-12 time series suggested that the drought in 2001, and in 4 consecutive years from 2007 to 2010 could be useful for the identification of refugia. Because ecological response to drought depends in part on drought duration and severity, this study area presented an opportunity to compare the effects of a single year drought (2001) with those of a multi-year drought (2007–2010).

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

Appendix B.3. Fires

Fire history data were obtained from two sources: (1) yearly fire extent and severity from 1999–2011 [59] and (2) fire perimeter and burn severity data from 1984 to 2014 from the Monitoring Trends in Burn Severity (MTBS) consortium [58]. Both data sources show lack of fire activity within the study area during the time period of this study (1985–2011).

Appendix B.4. Integrated Assessment of Disturbance History

Lack of disturbance from fire and human activities (logging, mining, development) suggests that the vegetation responses to disturbance in the study area have been primarily driven by droughts and insect outbreaks, making this study area well-suited to the identification of refugia from these disturbances. The period from 1985 to 2000 was chosen as a reference (i.e., control or baseline) period because (a) climatic moisture conditions were relatively well balanced between wet and dry periods based on SPEI-12 with a long-term average SPEI-12 near zero (Figure 3a), (b) insect outbreaks during 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).

Appendix C. Landsat Data

Table A2 lists the Landsat scenes for which Normalized Differenced Moisture Index (NDMI) datasets were obtained from the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center Science Processing Architecture (ESPA) on-demand interface. All scenes are Tier-1 Landsat Thematic Mapper (TM) 5 images for Path 044, Row 030, from August of the years 1985 through 2011 and had cloud cover ≤50%. These Landsat scenes are available for download from http://earthexplorer.usgs.gov.
Table A2. Landsat scenes for which Normalized Differenced Moisture Index (NDMI) was analyzed in this study.
Table A2. Landsat scenes for which Normalized Differenced Moisture Index (NDMI) was analyzed in this study.
Landsat Product IdentifierAcquisition DateCloud Cover
LT05_L1TP_044030_19850809_20161004_01_T18/9/19850
LT05_L1TP_044030_19850825_20161004_01_T18/25/19850
LT05_L1TP_044030_19860812_20161004_01_T18/12/19860
LT05_L1TP_044030_19860828_20161003_01_T18/28/198621
LT05_L1TP_044030_19870815_20161003_01_T18/15/198723
LT05_L1TP_044030_19870831_20161002_01_T18/31/198715
LT05_L1TP_044030_19880801_20161002_01_T18/1/19880
LT05_L1TP_044030_19890804_20161002_01_T18/4/19890
LT05_L1TP_044030_19890820_20161002_01_T18/20/19894
LT05_L1TP_044030_19900807_20161002_01_T18/7/199017
LT05_L1TP_044030_19900823_20161002_01_T18/23/199011
LT05_L1TP_044030_19910810_20161001_01_T18/10/19910
LT05_L1TP_044030_19910826_20161001_01_T18/26/19911
LT05_L1TP_044030_19920828_20160929_01_T18/28/19920
LT05_L1TP_044030_19930831_20160927_01_T18/31/19930
LT05_L1TP_044030_19940802_20160927_01_T18/2/19944
LT05_L1TP_044030_19940818_20160926_01_T18/18/19940
LT05_L1TP_044030_19950805_20160927_01_T18/5/19950
LT05_L1TP_044030_19950821_20160926_01_T18/21/19952
LT05_L1TP_044030_19960807_20160924_01_T18/7/19960
LT05_L1TP_044030_19960823_20160925_01_T18/23/19960
LT05_L1TP_044030_19970810_20160923_01_T18/10/199749
LT05_L1TP_044030_19980813_20160923_01_T18/13/19980
LT05_L1TP_044030_19980829_20160923_01_T18/29/19980
LT05_L1TP_044030_19990816_20160919_01_T18/16/19990
LT05_L1TP_044030_20000802_20160918_01_T18/2/20000
LT05_L1TP_044030_20000818_20160918_01_T18/18/20000
LT05_L1TP_044030_20010805_20160917_01_T18/5/20010
LT05_L1TP_044030_20010821_20160917_01_T18/21/20017
LT05_L1TP_044030_20020808_20160916_01_T18/8/20020
LT05_L1TP_044030_20020824_20160916_01_T18/24/20025
LT05_L1TP_044030_20030811_20160915_01_T18/11/20039
LT05_L1TP_044030_20030827_20160915_01_T18/27/20036
LT05_L1TP_044030_20040813_20160913_01_T18/13/200433
LT05_L1TP_044030_20040829_20160913_01_T18/29/20040
LT05_L1TP_044030_20050816_20160912_01_T18/16/20050
LT05_L1TP_044030_20060803_20160911_01_T18/3/20060
LT05_L1TP_044030_20060819_20160909_01_T18/19/20060
LT05_L1TP_044030_20070806_20160907_01_T18/6/20078
LT05_L1TP_044030_20070822_20160910_01_T18/22/20072
LT05_L1TP_044030_20080808_20160905_01_T18/8/20083
LT05_L1TP_044030_20080824_20160906_01_T18/24/20086
LT05_L1TP_044030_20090811_20160903_01_T18/11/20090
LT05_L1TP_044030_20090827_20160905_01_T18/27/20090
LT05_L1TP_044030_20100814_20160901_01_T18/14/20100
LT05_L1TP_044030_20110801_20160831_01_T18/1/20115
LT05_L1TP_044030_20110817_20160831_01_T18/17/20112

Appendix D. Refugia Identified in 2001 and 2009

Single cell and multi-cell refugia were mapped in the study area for 2001 (Figure A1) and 2009 (Figure A2). Spatial overlap was assessed for multi-cell refugia in both years (Figure A3).
Figure A1. Refugia identified in 2001 (10% of grid cells, stratified by forest type, that had the highest NDMI anomalies in 2001) and multi-cell refugia identified in 2001 (clusters of two or more contiguous grid cells identified as refugia in 2001).
Figure A1. Refugia identified in 2001 (10% of grid cells, stratified by forest type, that had the highest NDMI anomalies in 2001) and multi-cell refugia identified in 2001 (clusters of two or more contiguous grid cells identified as refugia in 2001).
Forests 09 00715 g0a1
Figure A2. Refugia identified in 2009 (10% of grid cells, stratified by forest type, that had the highest NDMI anomalies in 2009) and multi-cell refugia identified in 2009 (clusters of two or more contiguous grid cells identified as refugia in 2009).
Figure A2. Refugia identified in 2009 (10% of grid cells, stratified by forest type, that had the highest NDMI anomalies in 2009) and multi-cell refugia identified in 2009 (clusters of two or more contiguous grid cells identified as refugia in 2009).
Forests 09 00715 g0a2
Figure A3. Spatial overlap between multi-cell refugia identified in 2001 and 2009. Approximately 30% of multi-cell refugia grid cells in 2001 were also identified as multi-cell refugia in 2009.
Figure A3. Spatial overlap between multi-cell refugia identified in 2001 and 2009. Approximately 30% of multi-cell refugia grid cells in 2001 were also identified as multi-cell refugia in 2009.
Forests 09 00715 g0a3

Appendix E. Comparison of Boosted Regression Tree (BRT) Model Performance for Models of All Refugia (Single-cell and Multi-cell) Versus Models only for Multi-cell Refugia

Table A3. Boosted regression tree (BRT) model results relating probability of refugia to landscape characteristics, comparing models of all refugia (single-cell and multi-cell) to models only for multi-cell refugia.
Table A3. Boosted regression tree (BRT) model results relating probability of refugia to landscape characteristics, comparing models of all refugia (single-cell and multi-cell) to models only for multi-cell refugia.
BRT Model ParametersForest Type
LodgepoleLodgepole–whitebarkWhitebark
Year200120092001200920012009
Type of refugiaAllMultiAllMultiAllMultiAllMultiAllMultiAllMulti
Learning rate0.020.030.040.0550.0010.00250.0030.0040.00250.00350.0060.007
Number of trees250040003500300030002500300025002000250025002500
AUC-ROC0.770.850.910.920.670.720.840.850.690.760.860.87
Percent deviance explained0.320.590.670.730.180.340.470.510.250.40.570.62
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).

Appendix F. Partial-Dependence Plots from Boosted Regression Tree Models across 20 Bootstrapped Model Runs

For each of the four retained BRT models (LP 2001, LP 2009, LWP 2009, WP 2009; see Results section), partial-dependence plots were overlain for the 20 bootstrapped model runs and were averaged to produce an overall relationship for each combination of landscape characteristic and model (Figure A4, Figure A5, Figure A6 and Figure A7). Thick red lines (model averages) plotted in Figure A4, Figure A5, Figure A6 and Figure A7 correspond to plots in Figure 5 in Results.
Figure A4. Partial-dependence plots from boosted regression tree models of drought refugia in lodgepole pine areas in 2001. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Figure A4. Partial-dependence plots from boosted regression tree models of drought refugia in lodgepole pine areas in 2001. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Forests 09 00715 g0a4
Figure A5. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in lodgepole pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Figure A5. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in lodgepole pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Forests 09 00715 g0a5
Figure A6. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in lodgepole–whitebark pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Figure A6. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in lodgepole–whitebark pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Forests 09 00715 g0a6
Figure A7. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in whitebark pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Figure A7. Partial-dependence plots from boosted regression tree models of refugia from drought and mountain pine beetle impacts in whitebark pine areas in 2009. Black lines represent loess-smoothed partial-dependence relationships from 20 independent bootstrapped model runs using a random subsample of n = 2000 grid cells for each model run. Red lines indicate means across 20 bootstrapped model runs. In all plots, higher values on the vertical axis represent greater prevalence of refugia.
Forests 09 00715 g0a7

References

  1. Allen, C.D.; Macalady, A.K.; Chenchouni, H.; Bachelet, D.; McDowell, N.; Vennetier, M.; Kitzberger, T.; Rigling, A.; Breshears, D.D.; Hogg, E.T.; et al. A global overview of drought and heat-induced tree mortality reveals emerging climate change risks for forests. For. Ecol. Manag. 2010, 259, 660–684. [Google Scholar] [CrossRef]
  2. Dale, V.H.; Joyce, L.; Mcnulty, S.; Neilson, R.P.; Ayres, M.P.; Flannigan, M.D.; Hanson, P.J.; Irland, L.C.; Lugo, A.E.; Peterson, C.J.; et al. Climate change and forest disturbances. Bioscience 2013, 51, 723–734. [Google Scholar] [CrossRef]
  3. Mattson, W.J.; Haack, R. Role of drought in outbreaks of plant-eating insects. Bioscience 1987, 37, 110–118. [Google Scholar] [CrossRef]
  4. Raffa, K.F.; Aukema, B.H.; Bentz, B.J.; Carroll, A.L.; Hicke, J.; Turner, M.G.; Romme, W.H. Cross-scale drivers of natural disturbances prone to anthropogenic amplification: The dynamics of bark beetle eruptions. Bioscience 2008, 58, 501–517. [Google Scholar] [CrossRef]
  5. Chapman, T.B.; Veblen, T.T.; Schoennagel, T. Spatiotemporal patterns of mountain pine beetle activity in the southern Rocky Mountains. Ecology 2012, 93, 2175–2185. [Google Scholar] [CrossRef] [PubMed]
  6. Logan, J.; Powell, J. Ghost forests, global warming, and the mountain pine beetle (Coleoptera: Scolytidae). Am. Entomol. 2001, 47, 160–173. [Google Scholar] [CrossRef]
  7. Buotte, P.C.; Hicke, J.; Preisler, H.K.; Abatzoglou, J.T.; Raffa, K.F.; Logan, J. Climate influences on whitebark pine mortality from mountain pine beetle in the Greater Yellowstone Ecosystem. Ecol. Appl. 2016, 26, 2507–2524. [Google Scholar] [CrossRef] [PubMed]
  8. Bentz, B.; Regniere, J.; Fettig, C.; Hansen, M.; Hayes, J.; Hicke, J.; Kelsey, R.; Negron, J.; Seybold, S. Climate change and bark beetles of the western United States and Canada: Direct and indirect effects. Bioscience 2010, 60, 602–613. [Google Scholar] [CrossRef]
  9. Perkins, D.L.; Roberts, D.W. Predictive models of whitebark pine mortality from mountain pine beetle. For. Ecol. Manag. 2003, 174, 495–510. [Google Scholar] [CrossRef]
  10. Keane, R.E.; Tomback, D.F.; Aubry, C.A.; Bower, A.D.; Campbell, E.M.; Cripps, C.L.; Jenkins, M.B.; Mahalovich, M.F.; Manning, M.; McKinney, S.T.; et al. A Range-Wide Restoration Strategy for Whitebark Pine (Pinus albicaulis); General Technical Report RMRS-GTR-279; Department of Agriculture, Forest Service, Rocky Mountain Research Station: Fort Collins, CO, USA, 2012. [Google Scholar]
  11. U.S. Fish and Wildlife Service. Endangered and threatened wildlife and plants; 12-month finding on a petition to list Pinus albicaulis as endangered or threatened with critical habitat. Fed. Regist. 2011, 76, 42631–42654. [Google Scholar]
  12. Morelli, T.; Maher, S.; Nydick, K.; Monahan, W.; Ebersole, J.; Daly, C.; Dobrowski, S.; Dulen, D.; Jackson, S.; Lundquist, J.; et al. Managing climate change refugia for climate adaptation. PLoS ONE 2016, 11, e0159909. [Google Scholar] [CrossRef] [PubMed]
  13. Dobrowski, S. A climatic basis for microrefugia: The influence of terrain on climate. Glob. Chang. Biol. 2011, 17, 1022–1035. [Google Scholar] [CrossRef]
  14. Lenoir, J.; Hattab, T.; Pierre, G. Climatic microrefugia under anthropogenic climate change: Implications for species redistribution. Ecography 2016, 39, 1–14. [Google Scholar] [CrossRef]
  15. McLaughlin, B.; Ackerly, D.D.; Klos, P.Z.; Natali, J.; Dawson, T.E.; Thompson, S.E. Hydrologic refugia, plants and climate change. Glob. Chang. Biol. 2017, 23, 1–21. [Google Scholar] [CrossRef] [PubMed]
  16. Krawchuk, M.; Haire, S.L.; Coop, J.D.; Parisien, M.-A.; Whitman, E.; Chong, G.W.; Miller, C. Topographic and fire weather controls of fire refugia in forested ecosystems of northwestern North America. Ecosphere 2016, 7, e01632. [Google Scholar] [CrossRef]
  17. Mackey, B.; Berry, S.; Hugh, S.; Ferrier, S.; Harwood, T.D.; Williams, K.J.; Mackey, B.; Berry, S.; Hugh, S.; Ferrier, S.; et al. Ecosystem greenspots: Identifying potential drought, fire, and climate-change micro-refuges. Ecol. Appl. 2012, 22, 1852–1864. [Google Scholar] [CrossRef] [PubMed]
  18. Gould, S.F.; Hugh, S.; Porfirio, L.L.; Mackey, B. Ecosystem greenspots pass the first test. Landsc. Ecol. 2015, 30, 141–151. [Google Scholar] [CrossRef]
  19. Ashcroft, M.B. Identifying refugia from climate change. J. Biogeogr. 2010, 37, 1407–1413. [Google Scholar] [CrossRef]
  20. Keppel, G.; Mokany, K.; Wardell-Johnson, G.W.; Phillips, B.L.; Welbergen, J.A.; Reside, A.E. The capacity of refugia for conservation planning under climate change. Front. Ecol. Environ. 2015, 13, 106–112. [Google Scholar] [CrossRef] [Green Version]
  21. Dunne, T.; Zhang, W.; Aubry, B.F. Effects of rainfall, vegetation, and microtopography on infiltration and runoff. Water Resour. Res. 1991, 27, 2271–2285. [Google Scholar] [CrossRef]
  22. Lundquist, J.D.; Flint, A.L. Onset of snowmelt and streamflow in 2004 in the western United States: How shading may affect spring streamflow timing in a warmer world. J. Hydrometeorol. 2006, 7, 1199–1217. [Google Scholar] [CrossRef]
  23. Geroy, I.J.; Gribb, M.M.; Marshall, H.P.; Chandler, D.G.; Benner, S.G.; Mcnamara, J.P. Aspect influences on soil water retention and storage. Hydrol. Process. 2011, 25, 3836–3842. [Google Scholar] [CrossRef]
  24. Holden, Z.; Jolly, W.M. Modeling topographic influences on fuel moisture and fire danger in complex terrain to improve wildland fire management decision support. For. Ecol. Manag. 2011, 262, 2133–2141. [Google Scholar] [CrossRef]
  25. Jencso, K.G.; McGlynn, B.L. Hierarchical controls on runoff generation: Topographically driven hydrologic connectivity, geology, and vegetation. Water Resour. Res. 2011, 47, 1–16. [Google Scholar] [CrossRef]
  26. Flint, L.E.; Flint, A.L. Downscaling future climate scenarios to fine scales for hydrologic and ecological modeling and analysis. Ecol. Process. 2012, 1, 2. [Google Scholar] [CrossRef]
  27. Gutiérrez-Jurado, H.; Vivoni, E.R.; Cikoski, C.; Harrison, B.J.; Bras, R.L.; Istanbulluoglu, E. On the observed ecohydrologic dynamics of a semiarid basin with aspect-delimited ecosystems. Water Resour. Res. 2013, 49, 8263–8284. [Google Scholar] [CrossRef] [Green Version]
  28. Godfree, R.; Lepschi, B.; Reside, A.; Bolger, T.; Robertson, B.; Marshall, D.; Carnegie, M. Multiscale topoedaphic heterogeneity increases resilience and resistance of a dominant grassland species to extreme drought and climate change. Glob. Chang. Biol. 2011, 17, 943–958. [Google Scholar] [CrossRef]
  29. Shore, T.; Safranyik, L.; Lemieux, J. Susceptibility of lodgepole pine strands to the mountain pine beetle: Testing of a rating system. Can. J. For. Res. 2000, 30, 44–49. [Google Scholar] [CrossRef]
  30. Guarín, A.; Taylor, A.H. Drought triggered tree mortality in mixed conifer forests in Yosemite National Park, California, USA. For. Ecol. Manag. 2005, 218, 229–244. [Google Scholar] [CrossRef]
  31. Gibson, K.; Skov, K.; Kegley, S.; Jorgensen, C.; Smith, S.; Witcosky, J. Mountain Pine Beetle Impacts in High-Elevation Five-Needle Pines: Current Trends and Challenges; R1-08-020; USDA Forest Service: Missoula, MT, USA, 2008. [Google Scholar]
  32. Cooper, L.A.; Reed, C.C.; Ballantyne, A.P. Mountain pine beetle attack faster growing lodgepole pine at low elevations in western Montana, USA. For. Ecol. Manag. 2018, 427, 200–207. [Google Scholar] [CrossRef]
  33. Kaiser, K.E.; Mcglynn, B.L.; Emanuel, R.E. Ecohydrology of an outbreak: Mountain pine beetle impacts trees in drier landscape positions first. Ecohydrology 2013, 6, 444–454. [Google Scholar] [CrossRef]
  34. Daly, C.; Neilson, R.; Phillips, D. A statistical-topographic model for mapping climatological precipitation over mountainous terrain. J. Appl. Meteorol. 1994, 33, 140–158. [Google Scholar] [CrossRef]
  35. Stephenson, N.L. Actual evapotranspiration and deficit: Biologically meaningful correlates of vegetation distribution across spatial scales. J. Biogeogr. 1998, 25, 855–870. [Google Scholar] [CrossRef]
  36. Amman, G.D. Population changes of the mountain pine beetle in relation to elevation. Environ. Entomol. 1973, 4, 541–547. [Google Scholar] [CrossRef]
  37. Régnière, J.; Bentz, B. Modeling cold tolerance in the mountain pine beetle, Dendroctonus ponderosae. J. Insect Physiol. 2007, 53, 559–572. [Google Scholar] [CrossRef] [PubMed]
  38. Essery, R.; Li, L.; Pomeroy, J. A distributed model of blowing snow over complex terrain. Hydrol. Process. 1999, 13, 2423–2438. [Google Scholar] [CrossRef]
  39. Billings, W.; Bliss, L. An alpine snowbank environment and its effects on vegetation, plant development, and productivity. Ecology 1959, 40, 388–397. [Google Scholar] [CrossRef]
  40. Sturm, M.; McFadden, J.P.; Liston, G.E.; Chapin III, F.S.; Racine, C.H.; Holmgren, J.; Stuart Chapin, F.; Racine, C.H.; Holmgren, J. Snow—Shrub interactions in arctic tundra: A hypothesis with climatic implications. J. Clim. 2001, 14, 336–344. [Google Scholar] [CrossRef]
  41. Coops, N.C.; Wulder, M.; White, J.C. Integrating remotely sensed and ancillary data sources to characterize a mountain pine beetle infestation. Remote Sens. Environ. 2006, 105, 83–97. [Google Scholar] [CrossRef]
  42. Millar, C.I.; Charlet, D.A.; Westfall, R.D.; King, J.C.; Delany, D.L.; Flint, A.L.; Flint, L.E. Do low-elevation ravines provide climate refugia for subalpine limber pine (Pinus flexilis) in the Great Basin, USA? Can. J. For. Res. 2018, 48, 663–671. [Google Scholar] [CrossRef]
  43. Wulder, M.; White, J.C.; Bentz, B.; Alvarez, M.F.; Coops, N.C. Estimating the probability of mountain pine beetle red-attack damage. Remote Sens. Environ. 2006, 101, 150–166. [Google Scholar] [CrossRef]
  44. Buttrick, S.; Popper, K.; Schindel, M.; Mcrae, B.; Unnasch, B.; Jones, A.; Platt, J. Conserving Nature’s Stage: Identifying Resilient Terrestrial Landscapes in the Pacific Northwest; The Nature Conservancy: Portland, OR, USA, 2015. [Google Scholar]
  45. Basso, A.S.; Miguez, F.E.; Laird, D.; Horton, R.; Westgate, M. Assessing potential of biochar for increasing water-holding capacity of sandy soils. GCB Bioenergy 2013, 5, 132–143. [Google Scholar] [CrossRef]
  46. Siegel-Issem, C.M.; Burger, J.; Powers, R.; Ponder, F.; Patterson, S. Seedling root growth as a function of soil density and water content. Soil Sci. Soc. Am. J. 2005, 69, 215–226. [Google Scholar] [CrossRef]
  47. Wulder, M.; Dymond, C.C.; White, J.C.; Leckie, D.G.; Carroll, A.L. Surveying mountain pine beetle damage of forests: A review of remote sensing opportunities. For. Ecol. Manag. 2006, 221, 27–41. [Google Scholar] [CrossRef]
  48. Vogelmann, J.E.; Xian, G.; Homer, C.; Tolk, B. Monitoring gradual ecosystem change using Landsat time series analyses: Case studies in selected forest and rangeland ecosystems. Remote Sens. Environ. 2012, 122, 92–105. [Google Scholar] [CrossRef]
  49. Meddens, A.J.H.; Hicke, J.A.; Vierling, L.A.; Hudak, A.T. Evaluating methods to detect bark beetle-caused tree mortality using single-date and multi-date Landsat imagery. Remote Sens. Environ. 2013, 132, 49–58. [Google Scholar] [CrossRef]
  50. Walter, J.; Platt, R.V. Multi-temporal analysis reveals that predictors of mountain pine beetle infestation change during outbreak cycles. For. Ecol. Manag. 2013, 302, 308–318. [Google Scholar] [CrossRef]
  51. Assal, T.J.; Anderson, P.J.; Sibold, J. Spatial and temporal trends of drought effects in a heterogeneous semi-arid forest ecosystem. For. Ecol. Manag. 2016, 365, 137–151. [Google Scholar] [CrossRef]
  52. Schwantes, A.M.; Swenson, J.J.; Jackson, R.B. Quantifying drought-induced tree mortality in the open canopy woodlands of central Texas. Remote Sens. Environ. 2016, 181, 54–64. [Google Scholar] [CrossRef] [Green Version]
  53. Walker, G.; Ridenour, J. Mineral. Resource Potential of the Gearhart Mountain Wilderness and Roadless Area, Lake and Klamath Counties, Oregon; U.S. Geological Survey: Reston, VA, USA, 1982. [Google Scholar]
  54. PRISM Climate Group. PRISM Data Explorer: Time Series Values for Individual Locations. Available online: http://prism.oregonstate.edu/explorer/ (accessed on 1 December 2016).
  55. Ellenwood, J.R.; Krist, F.J.J.; Romero, S.A. National Individual Tree Species Atlas; U.S. Forest Service, Forest Health Technology Enterprise Team: Fort Collins, CO, USA, 2015. [Google Scholar]
  56. Beguería, S.; Vicente-Serrano, S. Global SPEI Database, SPEIbase Version 2.4. Available online: http://spei.csic.es/database.html (accessed on 1 May 2017).
  57. U.S. Forest Service Aerial Insect and Disease Survey GIS Data for Oregon and Washington 1947-Present. Available online: https://www.fs.usda.gov/detail/r6/forest-grasslandhealth/insects-diseases/?cid=stelprd3791643 (accessed on 1 December 2016).
  58. Eidenshink, J.; Schwind, B.; Brewer, K.; Zhu, Z.; Quayle, B.; Howard, S. A project for monitoring trends in burn severity. Fire Ecol. 2007, 3, 3–21. [Google Scholar] [CrossRef]
  59. LandFire Program LandFire Data Distribution Site. Available online: https://www.landfire.gov/viewer/ (accessed on 1 November 2016).
  60. Goodwin, N.R.; Coops, N.C.; Wulder, M.; Gillanders, S.; Schroeder, T.; Nelson, T. Estimation of insect infestation dynamics using a temporal sequence of Landsat data. Remote Sens. Environ. 2008, 112, 3680–3689. [Google Scholar] [CrossRef]
  61. Czerwinski, C.J.; King, D.J.; Mitchell, S.W. Mapping forest growth and decline in a temperate mixed forest using temporal trend analysis of Landsat imagery, 1987–2010. Remote Sens. Environ. 2014, 141, 188–200. [Google Scholar] [CrossRef]
  62. Hope, A.; Fouad, G.; Granovskaya, Y. Evaluating drought response of Southern Cape Indigenous Forests, South Africa, using MODIS data. Int. J. Remote Sens. 2014, 35, 4852–4864. [Google Scholar] [CrossRef]
  63. Huntington, J.; Mcgwire, K.; Morton, C.; Snyder, K.; Peterson, S.; Erickson, T.; Niswonger, R.; Carroll, R.; Smith, G.; Allen, R. Assessing the role of climate and resource management on groundwater dependent ecosystem changes in arid environments with the Landsat archive. Remote Sens. Environ. 2016, 185, 186–197. [Google Scholar] [CrossRef]
  64. U.S. Geological Survey. Product Guide: Landsat Surface Reflectance-Derived Spectral Indices; USGS Earth Resources Observation and Science (EROS): Sioux Falls, SD, USA, 2016.
  65. R Core Team. R: A Language and Environment for Statistical Computing. Available online: https://www.r-project.org (accessed on 1 February 2017).
  66. Hijmans, R.J. Raster: Geographic Data Analysis and Modeling, R Package Version 2.5-8. Available online: https://cran.r-project.org/web/packages/raster/raster.pdf (accessed on 1 May 2017).
  67. Cartwright, J. Analysis of Remotely-sensed Vegetation Conditions during Droughts and a Mountain Pine Beetle Outbreak, Gearhart Mountain Wilderness, Oregon. U.S. Geological Survey data release. Available online: https://doi.org/10.5066/f74q7swx (accessed on 16 November 2018).
  68. U.S. Geological Survey National Elevation Dataset (NED). Available online: https://lta.cr.usgs.gov/NED (accessed on 23 February 2017).
  69. Hengl, T.; de Jesus, J.M.; Heuvelink, G.B.M.; Gonzalez, M.R.; Kilibarda, M.; Blagotí, A.; Shangguan, W.; Wright, M.N.; Geng, X.; Bauer-Marschallinger, B.; et al. SoilGrids250m: Global gridded soil information based on machine learning. PLoS ONE 2017, 12, e0169748. [Google Scholar] [CrossRef] [PubMed]
  70. Friedman, J.; Hastie, T.; Tibshirani, R. Additive logistic regression: A statistical view of boosting. Ann. Stat. 2000, 28, 337–407. [Google Scholar] [CrossRef]
  71. Elith, J.; Leathwick, J.R.; Hastie, T. A working guide to boosted regression trees. J. Anim. Ecol. 2008, 77, 802–813. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  72. Parisien, M.; Mortiz, M. Environmental controls on the distribution of wildfire at multiple spatial scales. Ecol. Monogr. 2009, 79, 127–154. [Google Scholar] [CrossRef]
  73. Perry, G.; Wilmshurst, J.; Mcglone, M.; Napier, A. Reconstructing spatial vulnerability to forest loss by fire in pre-historic New Zealand. Glob. Ecol. Biogeogr. 2012, 21, 1029–1041. [Google Scholar] [CrossRef]
  74. Swets, J. Measuring the accuracy of diagnostic systems. Science 1988, 240, 1285–1293. [Google Scholar] [CrossRef] [PubMed]
  75. Trumbo, D.R.; Burgett, A.; Hopkins, R.L.; Biro, E.G.; Chase, J.M.; Knouft, J.H. Integrating local breeding pond, landcover, and climate factors in predicting amphibian distributions. Landsc. Ecol. 2012, 27, 1183–1196. [Google Scholar] [CrossRef]
  76. Ridgeway, G. Gbm: Generalized Boosted Regression Models, R Package Version 2.1.3. Available online: https://cran.r-project.org/web/packages/gbm/gbm.pdf (accessed on 1 May 2017).
  77. Hijmans, R.; Phillips, S.; Leathwick, J.; Elith, J. Dismo: Species Distribution Modeling, R Package Version 1.1-4. Available online: https://cran.r-project.org/web/packages/dismo/dismo.pdf (accessed on 1 May 2017).
  78. Powers, J.S.; Sollins, P.; Harmon, M.E.; Jones, J. Plant-pest interactions in time and space: A Douglas-fir bark beetle outbreak as a case study. Landsc. Ecol. 1999, 14, 105–120. [Google Scholar] [CrossRef]
  79. Huang, C.Y.; Anderegg, W.R.L. Large drought-induced aboveground live biomass losses in southern Rocky Mountain aspen forests. Glob. Chang. Biol. 2012, 18, 1016–1027. [Google Scholar] [CrossRef]
  80. Adams, H.R.; Barnard, H.R.; Loomis, A.K. Topography alters tree growth–climate relationships in a semi-arid forested catchment. Ecosphere 2014, 5, art148. [Google Scholar] [CrossRef]
  81. Matimati, I. The Relevance of Fog and Dew Precipitation to Succulent Plant Hydrology in an Arid South African Ecosystem. Ph.D. Thesis, University of the Western Cape, Bellville, South Africa, May 2009. [Google Scholar]
  82. Hill, A.J.; Dawson, T.E.; Shelef, O.; Rachmilevitch, S. The role of dew in Negev Desert plants. Oecologia 2015, 178, 317–327. [Google Scholar] [CrossRef] [PubMed]
  83. Fischer, D.T.; Still, C.J.; Williams, A.P. Significance of summer fog and overcast for drought stress and ecological functioning of coastal California endemic plant species. J. Biogeogr. 2009, 36, 783–799. [Google Scholar] [CrossRef]
  84. Sevruk, B. Regional dependency of precipitation-altitude relationship in the swiss alps. Clim. Chang. 1997, 36, 355–369. [Google Scholar] [CrossRef]
  85. Watson, F.G.R.; Anderson, T.N.; Kramer, M.; Detka, J.; Masek, T.; Cornish, S.S.; Moore, S.W. Chapter 5: Effects of wind, terrain, and vegetation on snow pack. In The Ecology of Large Mammals in Central Yellowstone; Garrott, R., White, P., Watson, F., Eds.; Elsevier Inc.: Boston, MA, USA, 2009; Volume 3, pp. 68–84. [Google Scholar]
  86. Lundquist, J.D.; Pepin, N.; Rochford, C. Automated algorithm for mapping regions of cold-air pooling in complex terrain. J. Geophys. Res. 2008, 113, D22. [Google Scholar] [CrossRef]
  87. Scott, G.L.; McCaughey, W.W. Whitebark pine guidelines for planting prescriptions. In National Proceedings: Forest and Conservation Nursery Associations—2005; Proc. RMRS-P-43; Riley, L.E., Dumroese, R.K., Landis, T.D., Eds.; U.S. Department of Agriculture, Forest Service, Rocky Mountain Research Station: Fort Collins, CO, USA, 2006; pp. 84–90. [Google Scholar]
  88. Miller, D.A.; White, R.A. A conterminous United States multilayer soil characteristics dataset for regional climate and hydrology modeling. Earth Interact. 1998, 2, 1–26. [Google Scholar] [CrossRef]
  89. Larsson, S.; Oren, R.; Waring, R.H.; Barrett, J.W. Attacks of mountain pine beetle as related to tree vigor of ponderosa pine. For. Sci. 1983, 29, 395–402. [Google Scholar]
  90. Whitehead, R.; Safranyik, L.; Shore, T. Preventative management. In The Mountain Pine Beetle: A Synthesis of Biology, Management, and Impacts on Lodgepole Pine; Safranyik, L., Wilson, B., Eds.; Natural Resources Canada, Canadian Forest Service: Victoria, BC, Canada, 2006; pp. 173–193. [Google Scholar]
  91. Woods, S.W.; Ahl, R.; Sappington, J.; Mccaughey, W. Snow accumulation in thinned lodgepole pine stands, Montana, USA. For. Ecol. Manag. 2006, 235, 202–211. [Google Scholar] [CrossRef]
  92. McKinney, S.T.; Fiedler, C.E. Tree squirrel habitat selection and predispersal seed predation in a declining subalpine conifer. Oecologia 2010, 162, 697–707. [Google Scholar] [CrossRef] [PubMed]
  93. Weaver, T. Whitebark pine and its environment. In Whitebark Pine Communities: Ecology and Restoration; Tomback, D.F., Arno, S.F., Keane, R.E., Eds.; Island Press: Washington, DC, USA, 2001; pp. 42–73. [Google Scholar]
  94. Raffa, K.F.; Mason, C.J.; Bonello, P.; Cook, S.; Erbilgin, N.; Keefover-Ring, K.; Klutsch, J.G.; Villari, C.; Townsend, P.A. Defence syndromes in lodgepole—Whitebark pine ecosystems relate to degree of historical exposure to mountain pine beetles. Plant Cell. Environ. 2017, 40, 1791–1806. [Google Scholar] [CrossRef] [PubMed]
  95. Keppel, G.; van Niel, K.; Wardell-Johnson, G.; Yates, C.J.; Byrne, M.; Mucina, L.; Schut, A.G.T.; Hopper, S.D.; Franklin, S.E. Refugia: Identifying and understanding safe havens for biodiversity under climate change. Glob. Ecol. Biogeogr. 2012, 21, 393–404. [Google Scholar] [CrossRef]
  96. Buma, B. Disturbance interactions: Characterization, prediction, and the potential for cascading effects. Ecosphere 2015, 6, 1–15. [Google Scholar] [CrossRef]
  97. Meddens, A.J.H.; Kolden, C.A.; Lutz, J.A.; Smith, A.M.S.; Cansler, C.A.; Abatzoglou, J.T.; Meigs, G.W.; Downing, W.M.; Krawchuk, M.A. Fire refugia: What are they, and why do they matter for global change? Bioscience 2018, 1–11. [Google Scholar] [CrossRef]
  98. Ruefenacht, B.; Finco, M.V.; Nelson, M.D.; Czaplewski, R.; Helmer, E.H.; Blackard, J.A.; Holden, G.R.; Lister, A.J.; Salajanu, D.; Weyermann, D.; et al. Conterminous, U.S. and Alaska forest type mapping using Forest Inventory and Analysis data. Photogramm. Eng. Remote Sens. 2008, 74, 1379–1388. [Google Scholar] [CrossRef]
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.
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.
Forests 09 00715 g001
Figure 2. Location of the Gearhart Mountain Wilderness study area in southern Oregon, USA. Refugia from drought and a mountain pine beetle outbreak were identified within lodgepole pine, lodgepole pine–whitebark pine, and whitebark pine forests.
Figure 2. Location of the Gearhart Mountain Wilderness study area in southern Oregon, USA. Refugia from drought and a mountain pine beetle outbreak were identified within lodgepole pine, lodgepole pine–whitebark pine, and whitebark pine forests.
Forests 09 00715 g002
Figure 3. Time series of (a) Standardized Precipitation Evapotranspiration Index calculated with 12-month antecedent conditions (SPEI-12); (b) areal percent of the study area affected by tree mortality from mapped insect and disease agents; (c) severity (in mean trees per hectare, TPH) for areas affected by mortality; (df) August values of Normalized Difference Moisture Index (NDMI) by forest type, including median value (bold line) and interquartile range (thinner lines). In all plots, the reference period 1985–2000 is shaded gray and the two drought years selected for analysis (2001 and 2009) are marked with dashed vertical lines. In (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).
Figure 3. Time series of (a) Standardized Precipitation Evapotranspiration Index calculated with 12-month antecedent conditions (SPEI-12); (b) areal percent of the study area affected by tree mortality from mapped insect and disease agents; (c) severity (in mean trees per hectare, TPH) for areas affected by mortality; (df) August values of Normalized Difference Moisture Index (NDMI) by forest type, including median value (bold line) and interquartile range (thinner lines). In all plots, the reference period 1985–2000 is shaded gray and the two drought years selected for analysis (2001 and 2009) are marked with dashed vertical lines. In (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).
Forests 09 00715 g003
Figure 4. Maps of median August values for normalized difference moisture index (NDMI) over (a) the reference period from 1985 to 2000, (b) 2001, a year of moderate drought, (c) 2009, a year of severe drought and mountain pine beetle outbreak. NDMI anomalies for (d) 2001 and (e) 2009 were calculated by subtracting reference NDMI from drought-year NDMI. The area for refugia identification is outlined in black in (a) through (e). Distributions of NDMI anomalies for 2001 and 2009 by forest type are presented in (f).
Figure 4. Maps of median August values for normalized difference moisture index (NDMI) over (a) the reference period from 1985 to 2000, (b) 2001, a year of moderate drought, (c) 2009, a year of severe drought and mountain pine beetle outbreak. NDMI anomalies for (d) 2001 and (e) 2009 were calculated by subtracting reference NDMI from drought-year NDMI. The area for refugia identification is outlined in black in (a) through (e). Distributions of NDMI anomalies for 2001 and 2009 by forest type are presented in (f).
Forests 09 00715 g004
Figure 5. Partial-dependence plots from boosted regression tree models depicting the relationships between landscape characteristics and probability of refugium occurrence in 2001 (dashed lines) and 2009 (solid lines), in lodgepole pine forest (blue), lodgepole–whitebark pine forest (purple), and whitebark pine forest (red). Plots represent loess-smoothed averages across 20 independent bootstrapped model runs using a random subsample of 2000 grid cells for each forest type in each model run (Appendix F). Higher values on the vertical axes indicate greater probability of a grid cell being identified as a refugium. Histograms at the bottom of each plot indicate distributions of landscape characteristics across all grid cells in the study area used for modeling.
Figure 5. Partial-dependence plots from boosted regression tree models depicting the relationships between landscape characteristics and probability of refugium occurrence in 2001 (dashed lines) and 2009 (solid lines), in lodgepole pine forest (blue), lodgepole–whitebark pine forest (purple), and whitebark pine forest (red). Plots represent loess-smoothed averages across 20 independent bootstrapped model runs using a random subsample of 2000 grid cells for each forest type in each model run (Appendix F). Higher values on the vertical axes indicate greater probability of a grid cell being identified as a refugium. Histograms at the bottom of each plot indicate distributions of landscape characteristics across all grid cells in the study area used for modeling.
Forests 09 00715 g005
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).
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).
Forests 09 00715 g006
Table 1. Hypothesized relationships between landscape characteristics (topographic, soil, and forest variables) and likelihood of refugia from drought and mountain pine beetle.
Table 1. Hypothesized relationships between landscape characteristics (topographic, soil, and forest variables) and likelihood of refugia from drought and mountain pine beetle.
VariableHypothesized RelationshipMechanism
Elevation (m)positiveMore 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)unknownDepending 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 radiusnegativeMore 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.
Topographic heat load index (HLI) Low HLI indicates topographic shadingnegativeMore 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 ridgetopspositiveMore 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/m3) SBD0 cm = SBD at the soil surface; SBD100 cm = SBD at 1-m depthnegativeMore 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 (m2/ha) Higher basal area indicates greater forest densitynegativeMore 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].
Distance to ecotone with fir forest (m)negativeMore refugia (in LP and WP stands) are expected closer to the ecotone with fir forest. MPB host tree species (LP and WP) may be buffered from MPB outbreak severity near forest ecotones to MPB-resistant fir species.
Percent fir (%)positiveMore refugia (in LP and WP stands) are expected in grid cells with a greater percentage of fir trees. MPB host species may be buffered from MPB outbreak severity if they are in less-homogenous forest stands or embedded in a matrix of MPB-resistant fir species.
MPB = mountain pine beetle; LP = lodgepole pine, WP = whitebark pine.
Table 2. Boosted regression tree model performance statistics and variable relative influence relating occurrence of multi-cell refugia to landscape characteristics.
Table 2. Boosted regression tree model performance statistics and variable relative influence relating occurrence of multi-cell refugia to landscape characteristics.
BRT Model ParametersCanopy Type
LodgepoleLodgepole–WhitebarkWhitebark
Year2001 *2009 *20012009 *20012009 *
Learning rate0.00150.00250.00150.0030.0020.0045
Number of trees3750 (769)4250 (596)3475 (617)3125 (535)3525 (786)3725 (658)
AUC-ROC0.77 (0.02)0.86 (0.02)0.71 (0.01)0.85 (0.01)0.74 (0.01)0.88 (0.01)
Cross-validated correlation0.3 (0.03)0.44 (0.05)0.19 (0.02)0.42 (0.02)0.27 (0.03)0.51 (0.02)
Percent deviance explained39.4 (4.21)56.75 (4.08)33.4 (3.45)52.65 (3.53)39.35 (4.52)66 (3.71)
Relative influence
Total basal area (15.9)14.29 (2.11)16.97 (2.08) 20.69 (1.17) 11.72 (0.77)
SBD100cm (13.4)16.35 (2.49)15.26 (2.48) 12.51 (0.86) 9.28 (0.79)
Slope (12.5)8.36 (1.63)10.16 (0.93) 12.27 (0.57) 19.22 (1.18)
Elevation (11.2)10.86 (1.78)11.54 (1.51) 13.9 (0.71) 8.38 (0.88)
HLI (10.9)9.96 (1.97)9.26 (1.43) 11.81 (0.83) 12.5 (1.09)
SBD0cm (10.3)9.2 (2.54)14.06 (2.78) 6.78 (0.72) 10.98 (1.09)
TPI (9.0)11.43 (2.14)6.67 (1.04) 7.34 (0.7) 10.54 (0.95)
CTI (8.1)8.8 (1.38)7.11 (1.05) 8.33 (0.51) 8 (0.92)
Distance to fir (7.3)7.68 (1.97)6.28 (1.23) 5.7 (0.45) 9.41 (0.61)
Percent fir (1.6)3.04 (1.18)2.69 (0.76) 0.66 (0.18) 0 (0)
AUC-ROC = area under the curve of the receiver operator characteristic, SBD100cm = soil bulk density at 100 cm soil depth, HLI = topographic heat load index, TPI = topographic position index, CTI = compound topographic index, SBD0cm = soil bulk density at the soil surface; units for each variable are presented in Table 1. All reported values represent means (with standard deviations in parentheses) across 20 bootstrapped model runs. Four models (denoted with *) were retained for interpretation based on AUC-ROC ≥ 0.75. For each retained model, the three landscape characteristics with highest relative influence values are in bold. Values in parentheses next to variable names represent mean relative influence across the four retained models.

Share and Cite

MDPI and ACS Style

Cartwright, J. Landscape Topoedaphic Features Create Refugia from Drought and Insect Disturbance in a Lodgepole and Whitebark Pine Forest. Forests 2018, 9, 715. https://doi.org/10.3390/f9110715

AMA Style

Cartwright J. Landscape Topoedaphic Features Create Refugia from Drought and Insect Disturbance in a Lodgepole and Whitebark Pine Forest. Forests. 2018; 9(11):715. https://doi.org/10.3390/f9110715

Chicago/Turabian Style

Cartwright, Jennifer. 2018. "Landscape Topoedaphic Features Create Refugia from Drought and Insect Disturbance in a Lodgepole and Whitebark Pine Forest" Forests 9, no. 11: 715. https://doi.org/10.3390/f9110715

APA Style

Cartwright, J. (2018). Landscape Topoedaphic Features Create Refugia from Drought and Insect Disturbance in a Lodgepole and Whitebark Pine Forest. Forests, 9(11), 715. https://doi.org/10.3390/f9110715

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop