Characterizing Spatial Neighborhoods of Refugia Following Large Fires in Northern New Mexico USA

Abstract

The spatial patterns resulting from large fires include refugial habitats that support surviving legacies and promote ecosystem recovery. To better understand the diverse ecological functions of refugia on burn mosaics, we used remotely sensed data to quantify neighborhood patterns of areas relatively unchanged following the 2011 Las Conchas fire. Spatial patterns of refugia measured within 10-ha moving windows varied across a gradient from areas of high density, clustered in space, to sparsely populated neighborhoods that occurred in the background matrix. The scaling of these patterns was related to the underlying structure of topography measured by slope, aspect and potential soil wetness, and spatially varying climate. Using a nonmetric multidimensional scaling analysis of species cover data collected post-Las Conchas, we found that trees and forest associates were present across the refugial gradient, but communities also exhibited a range of species compositions and potential functions. Spatial patterns of refugia quantified for three previous burns (La Mesa 1977, Dome 1996, Cerro Grande 2000) were dynamic between fire events, but most refugia persisted through at least two fires. Efforts to maintain burn heterogeneity and its ecological functions can begin with identifying where refugia are likely to occur, using terrain-based microclimate models, burn severity models and available field data.

1. Introduction

Large fires are known to result in a diversity of habitats along a spectrum of change relative to the pre-disturbed landscape [1,2,3,4]. Within this heterogeneity, places less affected by fire, including unburnt islands and peninsulas surrounded by the burnt matrix, play a key role in the ongoing function of ecosystems [5,6]. These places serve as refugia during the fire event, by fostering survival of roots and rhizomes, soil and canopy seedbanks and a wide array of species and communities. Following fire, surviving legacies allow repopulation of surrounding areas and ecosystem recovery. Over evolutionary time scales, the diverse ecological functions implicit in refugia hold promise for adaptation to future environmental change [3,7,8].

Various processes underlie the patterns of refugia on burn mosaics. Many studies have focused on identifying refugia through establishing a link to topographic heterogeneity (e.g., [9,10,11]). Specifically, refugia may be found where fire interacts with stable landscape features that influence microclimate, reduce flammability or create topographic or hydrologic barriers to fire spread [4,12,13,14,15]. In many cases, refugia occur where topography and fuels are heterogeneous, but in more homogeneous settings, fuel characteristics are more likely to influence patterns of refugia [4]. Stochastic factors including weather during a fire event are also increasingly recognized as an important driver of pattern in post-fire refugia [4,11,16]. The relationship between refugia and environment is also mediated by legacies from past fires [17].

A prior disturbance event can drive interactions with subsequent events by altering the inherent resistance or resilience of ecosystem components [17,18]. For example, structural properties of biological legacies from past fires can provide a mechanism for resistance to change in a subsequent fire event (e.g., fine fuel reduction [19]). Moreover, legacy species and communities from one fire event influence recovery in the surrounding landscape [20], altering fuel properties, including spatial connectivity, in the interval between past events and those that follow. Interactions between persistent legacies of past fires, including spatial patterns of refugial habitats, are important drivers of disturbance heterogeneity and ecological response [17,21].

Patterns of refugia following fire have also been defined by specific ecological characteristics and functions. For example, refugial habitats can include old growth forest remnants protected from fire [9,22] that maintain community processes and biodiversity [23]. Refugia may also support species compositions and habitat specialists not found elsewhere on the burn mosaic [24,25], including species that require post-fire seedling recruitment, that is, obligate seeders [26].

Change in remotely sensed indices through time has been used to map patterns of refugia; these indices include a green spot index derived using a time series of Fraction of Absorbed Photosynthetically Active Radiation imagery (MODIS fPAR; [25,27] and the Normalized Burn Ratio derived from spectral bands representative of biomass and moisture content of soil and vegetation, differenced from pre- to post-fire (Landsat TM dNBR; [11,28,29,30]). Remotely sensed data provide the opportunity to quantify refugial spatial patterns and explore the relationship of those patterns to both underlying processes and ecological functions.

Spatial statistics are commonly used in landscape ecology to visualize spatial trends and quantify spatial pattern and covariance in ecological data [31]. For example, kernel density estimates can be used to convert point data to continuous surfaces showing event density or intensity at varying scales defined by neighborhood size [32]. Kernel surfaces enable identification of neighborhoods with high densities, separated by areas of low density [33]. On burn mosaics, kernel density patterns of refugia represent habitat mosaics with abrupt structural changes from disturbed to undisturbed vegetation. Such structural transitions are important in determining rates of species turnover, providing various spatial configurations of seed sources for re-establishment of vegetation within severely disturbed areas [34]. In addition, the range of kernel densities allows the evaluation of a full spectrum of heterogeneity, from neighborhoods entirely protected from severe fire, to those where higher severity was predominant.

Fire regimes in the Jemez Mountains of northern New Mexico have been the focus of extensive research, due to the unique perspectives gained from the long history of interactions between climate, fire and human activity on this landscape [35,36,37,38]. Recent studies have documented rapid and ongoing changes in the landscape ecology of the region following several large fires [18,20,39]. The latest of the fires, Las Conchas, burned approximately 63,000 ha in the summer of 2011, including areas burned in several, previous fires [18]. To advance a landscape ecological perspective on these fire events, we pursued the following objectives: (1) quantify spatial patterns of refugia on the Las Conchas burn mosaic using kernel density surfaces; (2) evaluate underlying physical processes that generate those patterns; and (3) develop a link to ecological function by describing the plant communities associated with the spatial patterns. Finally, we explored the spatial dynamics of refugia by examining changes in refugial neighborhoods where the Las Conchas burn mosaic included a past fire event.

2. Materials and Methods

2.1. Study Area

The study area is located in northern New Mexico, USA, in the Northwestern Forested Mountains of the Western Cordillera, Southern Rockies ecosystem (Figure 1a,b) [40]. This ecoregion occupies the portion of the Rocky Mountains extending from southern Wyoming, through Colorado and into northern New Mexico. Climate patterns are characterized by fluctuations at a decadal scale influenced by the El Niño Southern Oscillation [41]. Annually, frequent, strong thunderstorms occur during July through early September. Storms in winter (December–March) also bring moisture, but in lesser amounts than received during the summer monsoon. Snow accumulates in winter at elevations >1500 m; below freezing, overnight lows occur throughout winter and are possible in any season.

The Las Conchas fire in 2011 occurred at the southern end of the ecoregion, in the eastern portion of the Jemez Mountains and on the Pajarito Plateau (longitude −106.4155°, latitude 35.84923°; Figure 1a,b). Spatial patterns of three older burns were also examined in the present study, in places where these areas were reburned by Las Conchas: La Mesa 1977, Dome 1996 and Cerro Grande 2000 (Figure 1c). The perimeters of the three older burns overlapped in some places. Dome reburned the southern edge of La Mesa, and Cerro Grande reburned La Mesa’s northern end. The Las Conchas fire overlapped substantially with the areas burned previously by La Mesa, Dome and Cerro Grande, as well as additional areas outside the older burn perimeters. All four fires occurred during extreme drought, and in each case, drought was preceded by a wetter El Niño period [42]. During the Las Conchas fire, a convergence of weather and fuel conditions led to an extremely rapid rate of spread during the first 13 h of burning [43].

Topographic variability is dramatic at the edges of steep canyons (Frijoles, Alamo and Capulin) and broad mesas that define the Pajarito Plateau. Ash-flow tuffs, erupted from the Jemez Mountains, define the plateau; its alternating broad mesas and steep canyons drain eastward to White Rock Canyon of the Rio Grande [44]. The San Miguel Mountains, a sub-range of the Jemez, overlap the Dome fire perimeter. Topography within the Cerro Grande fire perimeter is less heterogeneous at lower- to mid-elevations, increasing at mesa-canyon interfaces and with elevation; higher reaches include the Cerro Grande and Pajarito mountains. The area burned during Las Conchas that was outside of the older perimeters included a series of dissected canyons, higher elevation stretches of the Jemez Mountains and a portion of the Valles Caldera to the north and west.

Forest communities generally correspond with changes in elevation. Pinyon-juniper (Pinus edulis-Juniperus spp.) woodlands occur at the lowest elevations; ponderosa pine (P. ponderosa) forests appear as elevation increases. Continuing along the elevational gradient, there is a transition from ponderosa pine to mixed forests composed of two or more tree species, including ponderosa pine, Douglas-fir (Pseudotsuga menziesii), occasional Colorado blue spruce (Picea pungens), white fir (Abies concolor) and some Southwestern white pine (P. strobiformis). Quaking aspen (Populus tremuloides) can be found intermixed or in relatively pure stands. Engelmann spruce (P. engelmannii) occurs at higher elevations. The historical fire regime across these forest and woodland types was composed of a wide range of frequencies and severities that varied with fuels and climate [45,46]. Following the three older fires, the elevational gradient of vegetation communities became more of a patchwork of forested and non-forested vegetation types. For example, 28 years after the 1977 La Mesa fire, severely burned areas were repopulated by assemblages of shrub, grasslands and young trees surrounded by surviving forests [20,39,47]. The region’s ecology has evolved with human influence including the use of fire and cultivation of plants over thousands of years [37]. Diverse cultures of Native American Pueblos persist in many places, alongside private communities and lands managed by Bandelier National Monument and the Dome Wilderness, the Santa Fe National Forest and the Los Alamos National Laboratory, a United States Department of Energy facility.

2.2. Las Conchas RG

Within the Las Conchas fire perimeter, we constructed a spatial Refugial Gradient (RG) representing low- to high-density of locations classified as refugia (Figure 2). The classification was based on areas relatively unchanged using a satellite change-detection image (dNBR). The Normalized Burn Ratio (NBR) was calculated using Bands 4 and 7 from pre-fire (02/09/2010) and post-fire (05/09/2011) Landsat Thematic Mapper images (Path/Row 33/35); the difference was calculated to obtain dNBR (Eidenshink, et al. [48]). Fire perimeters for Las Conchas (Day 1 of burning: 17,654 ha; final perimeter: 61,057 ha) were obtained from the Geosciences and Environmental Change Science Center Outgoing Datasets (GeoMAC website) [49]. Using the change-detection image, we applied a threshold (−200 ≤ dNBR ≤ +200) to classify areas of little or no change from pre-fire as refugia (1); all other dNBR values were classified as non-refugia (0). Refugia, defined by the threshold, include low severity in vegetation types that historically had frequent fire, as well as low severity where fire has been less frequent with mixed severity. Then, we applied a Gaussian kernel estimator to the classified map in a moving window approximately 10 ha in size (11 $\times$ 11 window; 30-m cells). Specifically, a Gaussian weights matrix defined by its standard deviation (sd = 50; sum of values = 1) was applied to the 0, 1 values in each window, and the sum of the products was assigned to the center cell. The size of the neighborhood over which values are weighted and summed controls the amount of smoothing; small moving windows highlight finer-scale patterns, while larger moving windows reveal broader-scale patterns [31]. We considered a range of neighborhood sizes and settled on the 10-ha scale, because it produced a wide range in patterns, from neighborhoods highly populated (RG values close to or equal to 1) to those where refugia were rare (RG values close to 0). We used the statistical software R for all computations and analyses, including the analysis of spatial data [50]. Archived data and R code used in the analysis are available (https://0-doi-org.brum.beds.ac.uk/10.2737/RDS-2017-0005).

2.3. Refugial Environments

We modeled the post-Las Conchas refugial environment using terrain metrics derived from a digital elevation model at 30-m resolution and spatial climate data available for North America [51] (

Table 1. Spatial variables used to characterize the refugial environment, including terrain metrics derived from a 30-m DEM (calculations were done using R raster [52] and RSAGA [54]) and climate variables mapped at 1-km resolution [51].

Terrain Metric	Description
Local morphometry	Method: 2nd degree polynomial
Slope	In radians
Aspect	In radians
Overall curvature	Positive values: convex
Profile curvature	Positive values: convex
Plan curvature	Positive values: divergent
Characteristics of the Catchment Basin	Method: multiple flow direction
Slope	Average slope of the contributing area in radians
Size	Size of the contributing area in m2
Height	Altitudinal difference (relief)
Length	Length of flow path (m)
Aspect	Direction of slope (radians)
Other Indices
Topographic Convergence Index	Cold air pooling: log(size of the contributing area/local slope)
Topographic Wetness Index	Potential soil wetness based on area and slope of the contributing area
Relative position	Position of a raster cell compared to its surrounding cells (500-m radius window)
Topographic Wind Shelter	Maximum gradient in topographic relief from eight upwind directions (500-m radius window)
Climate Variable (Bioclim)
MAT	Mean Annual Temperature (°C)
MAP	Mean Annual Precipitation (mm)
NFFD	The Number Of Frost-Free Days

). The terrain metrics included catchment and local morphometrics, relative position, topographic convergence index, topographic wetness index and wind shelter. All metrics were calculated using raster [52] or RSAGA packages [53,54]. We tested three climate variables (1-km resolution) in the models: mean annual temperature, mean annual precipitation and number of frost-free days [51].

The goal of the modeling was to associate physical gradients with refugial habitats, measured at the kernel scale (~10 ha), and ultimately to represent the environment in relation to a multivariate community ordination. Thus, the models were meant to describe observed patterns, with the scope of inference limited to the study area. To that end, we sampled the terrain metrics, climate variables and the RG at random locations where the RG value was >0 (0.1 sample point per ha; n = 3611). First, we compared multiple regression Generalized Additive Models with no interaction terms (GAMs; R mgcv [55,56]) using four alternative combinations of the terrain metrics: (1) local morphometrics only; (2) catchment metrics only; (3) wetness, convergence and relative position only; and (4) two models that combined significant variables from Models 1 and 2. For local and catchment aspect predictors, we used a circular smoother available in the mgcv::gam function. Wind shelter variables were added one at a time to the best combined model, based on deviance explained, to check for improvement (∆AIC > 2; [57]). The climate variables were highly correlated (Pearson’s r > 0.8), so we developed three models, each containing only one climate predictor. The RG values followed a continuous distribution with a finite range, consistent with the beta family. We included a factor for random effects as a penalized regression term accounting for markedly different (i.e., more extreme) fire weather during Day 1 of the Las Conchas fire progression. Finally, we tested the usefulness of our final models by comparing statistics (AIC, p-values, deviance explained and adjusted r²) with statistics for null (intercept only) models and models with randomized predictor values.

2.4. Modeling Plant Communities

To characterize refugial plant communities, we used field data compiled in a recent study of upland vegetation patterns on the Las Conchas landscape [18]. Field data collection was done using sampling protocols in Muldavin et al. [39]; cover was recorded in 20 $\times$ 20-m plots by species and strata (tree, woody plants ≥5 m; shrub, woody plants <5 m, including tree seedlings and saplings; graminoid; and forb). We used two criteria to select field plots for our analysis: first, plots must have been measured post-fire (in either 2013 or 2014), and second, the RG value at that point had to be >0. These criteria resulted in 57 plots (

Table 2. Spatial data sources, processing methods and field sample count for the four burns. The Las Conchas dNBR image was provided by Steve Howard, US Geological Survey. For Dome and Cerro Grande, dNBR images were downloaded from the Monitoring Trends in Burn Severity website (http://mtbs.gov/). Details for La Mesa image processing (change in Normalized Difference Vegetation Index from pre- to post-fire (dNDVI)) are described in [18] and from Lisa Holsinger, US Forest Service, pers. com. The sample size of field plots within older burn perimeters is given (note: two plots burned in both La Mesa and Cerro Grande); seven plots did not burn in any of the previous fires.

Fire Name (Start Date)	Fire Size (ha)	Data Source (Resolution)	Pre- and Post-Fire Image Dates	Processing Algorithm	Field Plots (n)
La Mesa (16/06/1977)	5777	Multispectral Scanner Path/Row 35/36 (60-m)	29/05/1977; 18/05/1978	$dNBR=\left(1.467 \times dNDVI\right) + 24.86$	15
Dome (25/04/1996)	6387	Thematic Mapper Path/Row 33/35 (30-m)	23/07/1995; 14/09/1997	$dNBR=preNBR-postNBR$	9
Cerro Grande (05/05/2000)	17919	Thematic Mapper Path/Row 33/35 (30-m)	14/04/2000; 25/04/2001	$dNBR=preNBR-postNBR$	28
Las Conchas (26/06/2011)	61057	Thematic Mapper Path/Row 33/35 (30-m)	02/09/2010; 05/09/2011	$dNBR=preNBR-postNBR$	7

Date format: dd/mm/yyyy.

).

We conducted a Nonmetric Multidimensional Scaling Analysis (NMDS) using a matrix of species cover data (57 sites $\times$ 159 species). The process included Wisconsin double standardization with square root transformation of the data; the Bray–Curtis dissimilarity index was used to create the distance matrix (R vegan; [58]). We tracked stress over a range of dimensions (1–6) using the default maximum of 20 random starts, and compared results over several independent runs in the search of a stable solution. Then, we fitted a smooth surface to the ordination scores in two dimensions using GAMs, with Las Conchas refugial gradient values as the response and ordination scores from the chosen solution as predictors. We displayed the fitted surface in ordination space to aid interpretation of plant community distributions across refugial environments.

2.5. Gradient Dynamics

In order to assess how the refugial neighborhoods had changed from one fire event to another, we constructed spatial refugial gradients for each of the three older burns. Dome and Cerro Grande RGs were derived using the same methods as Las Conchas, described above (Table 2). Because La Mesa pre-dated Landsat Thematic Mapper, dNBR was predicted from Normalized Difference Vegetation Index (dNDVI) values derived using Multispectral Scanner Bands 3 and 4 [18]. We sampled the RGs within the older burn perimeters randomly, and at the 57 field plot locations and used this information to interpret the ordination results relative to past burn event. In addition, we identified where RG values at the field plots had undergone a dramatic shift from one fire to the next, defined by difference in RG ≥ |0.5| and described the plant communities in places that had shifted toward more or less densely populated neighborhoods.

3. Results

3.1. Las Conchas RG

The Las Conchas RG included several areas with high-density neighborhoods (RG ~1; Figure 2b). Moderate-density neighborhoods (0.4 < RG < 0.7) were scattered throughout the Las Conchas footprint and also typically occurred in zones situated between low and high density neighborhoods. Sparsely inhabited neighborhoods (RG < 0.4) were also widely distributed, but were most predominant in the northwestern portion of the Day 1 perimeter (refer to Figure 1b), which coincides with a high elevation area of the Jemez Mountains.

3.2. Refugial Environments

The RG varied in relation to particular topographic and bioclimatic environments, based on a random sample of locations (n = 3611). During model development, we determined that effects of wind shelter were important contributors to combination models (∆AIC > 2) in some cases, but not in others, and the direction of influence also varied. Other terrain predictors were either never significant (i.e., catchment height and area; all curvature metrics; relative position) or lost significance in combined models (i.e., catchment flow path length; local aspect and slope; topographic convergence index).

We selected a final terrain model containing catchment slope, catchment aspect and the topographic wetness index. These three predictors were significant (p < 0.05) in all multiple regression models that included these variables (i.e., Model 2: catchment metrics only; Model 3: wetness, convergence and relative position only; and Model 4: two models that combined significant variables from Models 1 and 2). The random effects factor representing Day 1 versus other days of fire progression was significant in all model trials. The statistics for the final model with three terrain predictors plus the random effects ‘day’ factor were: deviance explained = 13% and adjusted r² = 0.10. Climate variables were significant predictors of RG when added to the terrain variables (p < 0.001). Adding a climate variable to the final terrain model increased the significance of the terrain predictors (p < 0.001) and improved overall model fit. Adding either mean annual temperature or mean annual precipitation to the terrain model resulted in deviance explained = 23% and adjusted r² = 0.18; adding number of frost-free days as a predictor resulted in deviance explained = 21% and adjusted r² = 0.17. The final models, with three terrain variables and one climate predictor, represented significant improvement over the null (intercept only) model (AIC_final > AIC_null). Likewise, final models represented a significant improvement over models with randomized predictors. None of the randomized variables were significant at p < 0.05; deviance explained and adjusted r² were less than 1.0; the randomized models consistently contained less information (AIC_final > AIC_randomized).

Based on the GAM smooth plots for the models with three terrain and one climate predictor, high density neighborhoods occurred on gentle slopes, but low density neighborhoods were more likely on moderate to steep slopes (Figure 3a). Neighborhood density was greater on east-southeasterly aspects (Figure 3b) and in places with potentially higher soil moisture (Figure 3c). The field plots were located at the lower and mid-range of the catchment slope, but included a fairly wide range of variability in catchment aspect and topographic wetness (top rug on plots in Figure 3a–c).

Refugia increased in density under varied spatial climates (Figure 3d–f). A positive response occurred at low and high ends of the range in temperature, precipitation and length of growing season. The mid-range of the climate variables was more conducive to low density neighborhoods. The field plot locations fell in the mid-range of temperature and precipitation, but only represented the upper range of variability in number of frost-free days, compared to the random sample.

3.3. Community Models

The ordination revealed a wide distribution of plant assemblages across Las Conchas refugial environments (Figure 4). The two- and three-dimensional NMDS solutions (stress = 0.21 and 0.19, respectively) showed similar patterns when results were displayed on the modeled surface (GAM fit to NMDS 1 and 2: adjusted r² = 0.12; deviance explained = 0.38; and estimated degrees of freedom = 3.86). Most of the field plots were located within the perimeter of one or more past fires (Figure 4a). Those previously burned in Cerro Grande (Figure 4a; purple circles) were located in the left-hand portion of the ordination space, apart from other sites previously burned in La Mesa and Dome in the right-hand portion of the space (orange and violet circles). Among the sample sites, seven occurred in locations that were outside of the older burn perimeters (black circles); these sites were mainly clustered in one region of the ordination space, along the low to mid-range of the RG.

Communities exhibited a range of species compositions, including those with few or no trees and shrubby understory, grassy meadows, open canopy forest and rocky substrates with sparse grasses and forbs. To enhance the interpretation of the results, please refer to the table of the 159 species recorded at the 57 field plots (Appendix A). Data include scientific and common names, observed frequency across the field plots, species classifications [59], forest understory designation, and high, medium and low RG classes based on mean value.

Many species were widely distributed across the RG, including species that were abundant, as well as those less frequently observed (Appendix A). For example, beardlip penstemon (Penstemon barbatus), a characteristic understory species of open forests, was present across the full range of refugial neighborhood densities. Ponderosa pine also occurred widely, as did Gambel’s oak (Quercus gambelii), New Mexico locust (Robinea neomexicana) and Thurber’s fescue (Festuca thurberi). Species linked to forest habitats that were more limited in their distribution included creeping barberry (Mahonia repens), thimbleberry (Rubus parviflorus) and pineywoods geranium (Geranium caespitosum). Greater canopy cover of legacy trees occurred in the medium to high range of the RG (Figure 4b); species included ponderosa pine, Douglas-fir, white fir and aspen. Forest understory species (Appendix A, denoted by asterisk) were widely distributed across the refugial environments, including some sites where trees were absent (Figure 4c).

3.4. Gradient Dynamics

The Las Conchas RG values at the field plots had changed in comparison to gradient values of past fires, but the magnitude and direction of change varied (Figure 5). Out of the 57 sample plots, 88 percent were located within the fire perimeters of La Mesa (n = 14; orange circles), Dome (n = 9; violet circles) or Cerro Grande (n = 27; purple circles). For the 15 plots that previously burned in La Mesa (Figure 5a), the maximum position shift was dramatic in either direction (median change = −0.11; min = −0.95, max = 0.83). For six of these plots, gradient values after Las Conchas were higher (i.e., greater spatial density of refugia). All nine field plots that previously burned in the Dome fire (Figure 5b) decreased in gradient value after Las Conchas, with dramatic shifts from dense to sparse in a few cases (median change = −0.50; min = −0.99, max = −0.08). The 28 plots that previously burned in Cerro Grande (Figure 5c) were also split in direction of change, and similar to La Mesa, there were plots that saw a shift from one end of the gradient to the other (median change = 0.01; min = −0.96, max = 0.98); 17 of these 28 field plots moved to higher positions on the RG after Las Conchas. Only two plots burned at higher severity in past fires (i.e., RG = 0); one at La Mesa and one at Cerro Grande. The patterns of direction of change we observed at the field plots were consistent with a random sample of points within each of the older burns that reburned in Las Conchas (La Mesa, n = 390; Dome, n = 641; Cerro Grande, n = 815; results not shown). Moreover, the gradient values of random points pre-Las Conchas were poorly correlated with values post-Las Conchas (Pearson’s r < 0.33 in all cases).

Some field plots where the RG value increased substantially after Las Conchas (difference ≥0.5) were inhabited by trees (i.e., canopy >0: La Mesa: 3 out of 3; Dome 0 out of 0; Cerro Grande 3 out of 5). Several of the plots where the value decreased (difference ≤−0.5) were also still forested after Las Conchas (La Mesa 2 out of 6; Dome 2 out of 5; Cerro Grande 0 out of 2). Most species were present at one or more plots where the gradient value exhibited a dramatic change after Las Conchas (difference ≥|0.5|; Appendix A). Many species were found at plots that experienced a negative change (sparser neighborhood) at some places and a positive change (denser neighborhood) at others, including ponderosa pine, quaking aspen, Rocky Mountain iris (Iris missouriensis) and New Mexico locust. Kingcup cactus (Echinocereus triglochidiatus), pine dropseed (Blepharoneuron tricholepis) and mountain ninebark (Physocarpus monogynus) were examples of species located at one or more field plots that saw a negative shift only; examples of species found at one or more field plots that saw a positive shift only were Douglas-fir, wavyleaf oak (Quercus x pauciloba), mountain mahogany (Cercocarpus montanus) and Colorado barberry (Berberis fendleri).

4. Discussion

Landscape ecological studies have made significant contributions to understanding the patterns of biological legacies following large disturbance events, notably after the Yellowstone fires of 1988 [60,61]. Here, we qualified patterns of legacies, or refugia, by describing their neighborhood and environmental setting. On the study landscape, we found that refugia occurred in differing spatial densities at the selected scale of measurement; high density neighborhoods were clustered in space, and sparsely populated neighborhoods, including unpopulated areas, occurred in the background matrix (Figure 2).

The scaling of neighborhood patterns of refugia following the Las Conchas fire was correlated with three topographic predictors: slope inclination, aspect and topographic wetness index (Figure 3). Each of these variables may interact with fire behavior in ways that support refugia formation. Refugia were denser at relatively flat sites and places with higher potential soil moisture. In our study region, these settings are generally associated with shifts from upper to lower slope positions, with corresponding increases in soil moisture, finer soil texture, colder temperatures, greater herbaceous cover and reduced tree cover and needle litter [62]. Transitions to wetter and cooler conditions are expected to reduce surface fire line intensity, flame length and crown fire initiation. Likewise, low tree density in these settings apparently did not support active crown fire. Moreover, higher spatial density of refugia often occurred where valley-bottom meadows bordered low-density stands of trees.

Refugia also occurred on moderate to steep slopes, but at lower densities. These locations often support scattered trees above cliffs, talus slopes and on rocky unproductive substrates (J.D. Coop, personal observation), all of which could reduce surface fire intensity, crown fire initiation and spread. Finally, refugia at the high end of the RG on south and east aspects may be related to both grassy openings and low-density stands on south-facing slopes in the study area. The predominant direction of fire spread from west-to-east during Day 1, during which time fire moved uphill (at higher speed and intensity) on west-facing slopes, may have resulted in greater refugia occurrence on south and east slopes, as well. Our observations of refugial environments complement recent findings regarding the importance of topographically-mediated fuel moisture on burn severity in dry forests elsewhere [63,64].

In addition, spatial climate influenced patterns of refugia on the post-Las Conchas landscape. Dense configurations (i.e., high RG) were associated with more extreme climates (Figure 3). High and low temperatures generally occur at the extremes of elevation; dense refugia at low elevation included areas that burned previously, in the easterly side of the burn (Figure 1 and Figure 2). Refugia in more scattered configurations (i.e., low RG) were common within the Day 1 perimeter, leading to a logical association of these patterns with extreme fire weather (Figure 2). However, our models indicate that given random effects associated with Day 1, spatial climate characterized by moderate temperature and precipitation and growing season of intermediate length played a significant role in generating these sparse refugial patterns (Figure 3).

In the context of other research, refugia found in wetter or cooler places conform to common definitions of climate refugia under a warming climate [65], but our observation that spatial patterns of refugia were associated with varied environments motivates the consideration of how refugia function in diverse neighborhoods. The differing spatial density of legacy species and communities could afford both advantages and disadvantages for near-term recovery of severely burnt areas. At the high density end of the RG, post-fire legacies are likely buffered from extreme conditions on the surrounding burn mosaic by other survivors. These buffered habitats may provide favorable conditions for continued survival through plant-soil interactions (e.g., nitrogen availability [66]). Dense neighborhoods may also provide additional opportunities for increased diversity through seed dispersal by wind or animals [67] at local, within-neighborhood scales.

In contrast, widely scattered refugia (i.e., low density neighborhoods) may be more exposed to drying from wind and erosion of soils by wind or water in their immediate neighborhood. Legacies that survive in these more isolated, potentially extreme conditions, however, might be particularly important in advancing reforestation rates. For example, these individuals occupy an advantageous position for the dispersal of seeds into places where soil nutrients have been released by fire [68,69]. In the mid-range of the gradient, the ecotones present in mixed neighborhoods (i.e., with a varying proportion of low- to high-severity context) have been found to represent a diversity of resources and habitats that support important ecological functions [69,70,71,72].

The presence of trees and forest associates even at the low end of the RG appears encouraging for forest recovery (Figure 4; Appendix A), but studies indicate that repopulation of the surrounding landscape in the near term may be limited due to drought conditions [73] and competition from early post-fire resprouting shrubs [74]. The long-term protection of seedbanks, through topographically defined microclimates [65], as well as active management may be necessary to maintain these refugial populations until conditions in the surrounding landscape become more favorable for forest species.

Although forest species are of great interest, refugia comprised of non-forest species deserve equal attention for their role in ecosystem function [75]. The species we observed included both forest and non-forest species that hold important traditional uses among Pueblo peoples, including Colorado barberry, Carruth’s sagewort (Artemisia carruthii), ragleaf bahia (Bahia dissecta) and many others [37]. Grass and shrub species, present across the RG (Appendix A), are often found in the surrounding burned matrix and tend to persist through multiple fire events by rebounding quickly [18]. Where meadows and shrublands are protected from severe fire, these communities may harbor legacies, including seedbanks of obligate seeders (including the conifer trees found in the region), and may provide habitat for animal seed dispersers [34,76]. The montane meadows we observed can be rich in biodiversity and serve as indicators of environmental change [77,78]. Furthermore, early successional grassy openings created by fire have been observed to transition to young forests in some places [20].

Our field sampling excluded riparian areas, but the map of the refugial gradient reveals that refugia occurred in riparian environments on the Las Conchas landscape along creeks such as El Rito de Los Frijoles. In other studies, riparian areas served as refugia for P. sabiniana in the southern Sierra Nevada, California, USA, enabling pines to recolonize severely-burned areas during fire-free intervals [79]. In frequent-fire environments of southeastern Australia, gullies were found to play a role in preserving heterogeneity in forest structure [15]. These wetter, cooler areas may be critical climate and fire refugia in rapidly-changing environments [80].

Refugial habitats, overall, are expected to change through time, in both form and function [80]. At Las Conchas, forest remnants experienced both increases and decreases in neighborhood refugial density in comparison to past fires. It is possible that some of these increases may be due to post-fire tree regeneration that occurred between fire events, particularly where increases occurred between the 1977 La Mesa fire and 2011. Recolonization of these landscapes has been found to result from the migration of regenerating trees from the edge of surviving trees in a wave-form succession [47]; with sufficient time since fire and favorable climate conditions for germination and growth, young trees may have reached sufficient stature to survive after the Las Conchas fire and thus increased the size of refugia. Decreasing density of refugia may occur in places where mechanisms such as environmental setting (e.g., cooler or wetter microclimates; topographic position) are insufficient to maintain patterns of refugia observed following earlier fire events. Furthermore, protection afforded by terrain can be lessened under severe fire weather [11].

Most of the refugial neighborhoods after Las Conchas coincided with refugia after a previous fire (Figure 5), indicating that these refugia persisted through at least two fires. Environments and associated communities may differ within particular burn events and with fire intervals, as evidenced by the separation of sites in ordination space that burned in Cerro Grande from those burned in La Mesa and Dome (Figure 4). Likely mechanisms for persistence of refugia include topographic setting (Figure 3), such as wetter environments and steep, grassy slopes with scattered trees on Cerro Grande, Pajarito and other mountains. Burn patterns from the 2000 Cerro Grande fire were important contributors to successful containment efforts during the 2011 Las Conchas fire [81]; knowledge of these interactions can be used to promote a persistent template of refugia defined by particular physical environments and species compositions.

Our results indicate that when a neighborhood increases or decreases in density of refugia, habitats and their inhabitants may also change. A wide variety of species were represented at locations that became more or less populated by refugia. The implications of this observation need further investigation, to understand each species’ independent response to its changing neighborhood environment and the resulting complex dynamics [82]. Species interactions and dispersal patterns, as well as habitat parameters (e.g., temperature and water availability) likely change with neighborhood context.

For many of the species we observed, including several tree species, our study area at the southern end of the Rocky Mountains is located at the rear edge, or the lower latitudinal margin of their geographic distribution. As climate change creates shifts in species distributions, the persistence of rear edge populations may depend on areas of greater climatic stability that occur in regions of diverse topography [80]. Stressors including drought [83] and fire [73] limit the distribution of vulnerable species to refugial environments [80,84]. Identifying the regional context of refugia within burn mosaics is important to gaining a broader view of places protected from changing climate and disturbance regimes [76,85].

5. Conclusions

The landscape ecological perspective we adopted illustrates one approach to understanding refugia and their role in the ecological functioning of post-fire landscapes. Alternative perspectives, definitions and models will continue to emerge, based on studies conducted in different locations and at varying spatial scales. Moreover, examining additional factors that interact to define refugial environments remains a key priority and could serve to improve the limited explanatory power of the models we developed specifically for our study landscape. These factors include the temporal progress of burning and associated fire weather that have been shown to shape refugial environments [11]. Future development and tests of new indices derived from remotely sensed data could also result in improved understanding of refugial characteristics following fire events.

Refugia on burn mosaics are increasingly acknowledged as a key resource for recovery and adaptation to change that needs to be protected and maintained [12,23,72,76,86]. However, widespread efforts to exclude, mitigate or contain severe fire through forest thinning, suppression and fire-fighting activities (including back-burning), as well as efforts to speed recovery from severe fire do not currently incorporate knowledge of existing or potential refugia on managed landscapes [87]. Large areas of tree mortality can easily dominate the early post-fire perspective in the absence of a quantitative assessment of burn heterogeneity and associated ecological functions.

Recognizing where refugia are likely to occur following fire could aid in maintaining burn heterogeneity and prevent severe degradation of refugial resources. We recommend starting with a coarse-filter based on the physical template. Terrain metrics that indicate microclimates where refugia tend to form in mountains [65] are easily computed from digital elevation models using open source software and code, such as R [50]. These spatial metrics can be viewed in conjunction with burn severity models; by applying a range of low-change thresholds, potential refugia can be identified.

Additionally, surveys of old growth forest structure [88] and species requiring specialized habitats protected from severe fire (e.g., Jemez Mountain salamander [89]) can also be taken as indicators of refugia on disturbed or undisturbed landscapes. We concur with recommended conservation measures considering the potential impacts of climate change that include maintenance of the greatest possible number of local refugia as habitat networks [25]. In this way, the potential contribution of refugia to regional diversity and landscape resilience under different future scenarios could be maximized [80,90].