Use of landscape simulation modeling to quantify resilience for ecological applications

Goals of fostering ecological resilience are increasingly used to guide U.S. public land management in the context of anthropogenic climate change and increasing landscape disturbances. There are, however, few operational means of assessing the resilience of a landscape or ecosystem. We present a method to evaluate resilience using simulation modeling. In this method, we use historical conditions (e.g., in North America, prior to European settlement), quantified using simulation modeling, to provide a comparative reference for contemporary conditions, where substantial departures indicate loss of resilience. Contemporary ecological conditions are compared statistically to the historical time series to create a resilience index, which can be used to prioritize landscapes for treatment and inform possible treatments. However, managing for resilience based on historical conditions is tenuous in the Anthropocene, which is characterized by rapid climate change, extensive human land use, altered disturbance regimes, and exotic species introductions. To account for the future variability of ecosystems resulting from climate and disturbance regime shifts, we augment historical simulations with simulations of ecosystem dynamics under projected climate and land use changes to assess the degree of departure from benchmark historical conditions. We use a mechanistic landscape model (FireBGCv2) applied to a large landscape in western Montana, USA, to illustrate the methods presented in this paper. Spatially explicit ecosystem modeling provides the vehicle to generate the historical and future time series needed to quantify potential resilience conditions associated with past and potential future conditions. Our methods show that given selection of a useful set of metrics, managers could use simulations like ours to evaluate potential future management


INTRODUCTION
Management for ecological resilience, mandated by some U.S. public land policies, is intended to guide land stewardship in a context of profound environmental challenges caused by complex and potentially novel interactions of anthropogenic climate changes, shifting fire regimes, exotic plant, insect, and pathogen invasions, and industrial, agricultural, and urban development (Moritz and Agudo 2013, Joyce et al. 2014, Bone et al. 2016, Kolb et al. 2016, Smith et al. 2016, Schoennagel et al. 2017. For example, the National Cohesive Wildland Fire Management Strategy identifies creating resilient landscapes as the first of its three major goals (USDOI and USDA 2014) and the U.S. Forest Service is mandated to restore natural resources to be "more resilient to climate change" as stated in the 2012 National Forest Planning Rule (USFS 2012). However, operational protocols for real-world implementation are not readily available to evaluate or quantify resilience for most regions (Angeler andAllen 2016, Falk 2016). To manage for ecological resilience, land managers must have a standardized and scientifically credible method of quantifying resilience based on concrete concepts that can be translated into multiple use analysis and planning (Bone et al. 2016, Colavito 2017, Dale et al. 2018. Quantifying ecological resilience for operational use confronts many challenges. A first challenge is adoption of a working definition of resilience that is appropriate for terrestrial, disturbance-prone ecosystems. For example, Gunderson's (2000) definition, "the capacity of a system to absorb impacts before a threshold is reached where the system changes into a different state," fails to recognize alternative states that are part of the system, which is similar to Holling's (1973) definition of ecological resilience as "a measure of the persistence of systems and of their ability to absorb change and disturbance and still maintain the same relationships between populations or state variable." Other authors define resilience as the time required for a system to return to a steady state following a perturbation (Ives 1995), assuming that the system has a single or global equilibrium condition (Mittelbach et al. 1995, Neubert andCaswell 1997). These definitions are limited in an operational context because they are focused on single disturbance events and do not explicitly recognize (1) potential for multiple states; (2) spatial and temporal domains of ecological dynamics; (3) inherent properties of the perturbed system; (4) disturbance regimes (type, duration, frequency, and intensity); and (5) potential for multiple interacting disturbances with synergistic or additive effects (Loehman et al. 2014, Fisichelli et al. 2016. For example, for systems with inherent, fire-adapted characteristics, a lack of disturbance (e.g., fire exclusion) over a long time period may create a new ecological state with low resilience as defined above (Merschel et al. 2014). In dry forests of interior western North America, where historically frequent fire often maintained open forests dominated by large, fire-tolerant trees (Veblen et al. 2000, Savage and Mast 2005, Hessburg et al. 2015, limited competition increased tree vigor, which in turn increased the ability of the fire-adapted trees to withstand insects, diseases, and fire, that is, conferring resilience to other disturbances (Bentz et al. 2010, van Mantgem et al. 2016. The exclusion of fire over the past century has heightened competitive stress among trees and increased surface, canopy, and ladder fuel loads, causing higher surface fire intensities, greater crown fire spread, and higher tree mortality than occurred historically. These interacting disturbances increase the potential for conversion from forests to persistent open savannas, shrublands, or grasslands after fire (Coop et al. 2016, Walker et al. 2018, resulting in a loss of resilience (Savage andMast 2005, Collins et al. 2016).
A second challenge to operationalizing models of ecosystem resilience is the selection of the number and type of variables to use to represent resilience; to be used effectively, these variables must align with multiple management objectives and alternatives (Keane 2012, Bone et al. 2016. Choosing only one or two assessment variables may ignore other important ecosystem attributes that contribute greatly to both resilience and management objectives, but including too many variables in an analysis may overcomplicate an already complex process (Moritz et al. 2005). Calculating resilience from only species composition variables, for example, may fail to capture key components of demography or development, such as structure and multiple canopy layers, or keystone ecosystem processes (Zedler 1996). Ideally, the variables used to evaluate resilience must be readily measurable in the field and represent attributes that are used commonly in land management planning (Bone et al. 2016).
A third challenge is scale: Over what scales of time and space should resilience be evaluated? An analysis area that is too small (e.g., a stand), for example, might ignore the influence of adjacent lands on resilience, whereas an area that is too large may be inappropriate for planning management actions (Karau andKeane 2007, Keane 2012). Suitable timescales will vary by ecosystem and productivity gradients and should match species life history traits. For example, success of tree regeneration after a fire is often used as a measure of resilience (Stevens-Rumann et al. 2018), but successful tree regeneration may take decades in some forest types, such as after burns in high-elevation settings (Agee andSmith 1984, O'Connor et al. 2014). Additionally, bounding the ecological context for resilience is important. A parking lot, strip mine, or an area dominated by exotic species, for example, may be resilient in that they are difficult to alter to a different state, but not ecologically healthy, socially desirable, or a viable management target.
We propose that an operational resilience metric must first have a descriptive foundation that matches the resolution and details of the analytical procedure (Dale et al. 2018), but be simple enough to be understood and computed by managers and scientists (Colavito 2017). The resilience metric uses a benchmark or reference condition to form the basis of comparison across appropriate time and space scales and ecological context (Cissel et al. 1994, Laughlin et al. 2004. As advocated by Seidl et al. (2016), reference ranges of variability can be used to "locate and delineate the basins of attraction of a system. . . [to]. . .measure resilience." Contemporary and projected future conditions can then be compared to this reference to determine departures, and to calculate a statistical index of departure of contemporary or future conditions to the historical reference period (Hessburg et al. 1999, Swetnam et al. 1999. As in Seidl et al. (2016), we propose to use historical reference conditions as a benchmark for characterizing and quantifying ecological resilience for land management. Historical conditions have been used extensively and successfully as references, benchmarks, or targets by researchers and land managers over the last two decades (Hessburg et al. 1999, Dickinson 2014. The concept behind historical range and variation (HRV) is that historical ecosystem characteristics, described by management-relevant variables such as burned area, species composition, or patch size distribution, represent the broad envelope of responses possible for a persistent (resilient) ecosystem under natural perturbations of climate, competitive stress, disturbances, and other stressors (Fig. 1). The HRV concept captures the fact that ecosystems are dynamic and highly variable, and their potential responses to changing conditions are consequently best represented by past variability of ecological characteristics (Veblen 2003).
A next critical step in the use of HRV for evaluating resilience is the quantification of historical conditions at the appropriate scale, context, and level of realism (Keane 2012). HRV has been described for many ecosystems across the western United States using a wide variety of techniques (see Keane et al. 2009). Time series HRV estimates can be generated for any scale of analysis, especially if they are produced by simulation models (Keane et al. 2015). The degree of departure of contemporary or future conditions from historical analogs can then be used as a resilience metric for prioritizing ecosystems or landscapes for treatment and to inform potential treatments (Morgan et al. 1994, Landres et al. 1999. In this paper, we describe the use of simulation modeling as a means for generating the reference time series for quantifying HRV that is then used to calculate resilience metrics to inform land management. The use of HRV as a benchmark to guide future management has been challenged, primarily because climate and other environmental conditions are changing so rapidly that the future may be significantly different than analogs in the ecologically relevant past (Millar and Woolfenden 1999, Millar et al. 2007, Shive et al. 2014. As climate projections suggest, future environments may become so distinct from the past that, at some point, managing to maintain strict HRV may be impractical if not irrelevant Woolfenden 1999, Millar et al. 2007). While HRV may reflect the range of conditions under which ❖ www.esajournals.org the current flora and fauna have evolved, future, no-analog climate regimes may make it difficult to manage landscapes within historical envelopes (Falk and Millar 2016). On the other hand, ecosystems possess powerful stabilizing legacies that can buffer the impacts of change, including ❖ www.esajournals.org 5 September 2018 ❖ Volume 9(9) ❖ Article e02414 topoclimatic refugia, genetic adaptations to emerging conditions, and unexpressed phenotypic variation (Johnstone et al. 2016). On balance, we posit that HRV remains a valuable tool to bracket potential ecosystem responses and frame resilient management goals (Alagona et al. 2012, Seidl et al. 2016). We do not suggest that managers ignore future climates; instead, we present a way to map HRV onto multiple future ranges of variability (FRV) to manage for resilience. In this paper, we propose a method (called hereafter "the HRV-resilience method") for quantifying resilience for ecological applications, based on a resilience index that evaluates statistically significant departures from HRV. We use a spatially explicit ecosystem simulation model (FireBGCv2; Keane et al. 2011) to generate HRV time series-the resilience envelope for a particular landscape (Seidl et al. 2016)-and potential FRVs. The comparison of FRVs with HRV allowed us to explore potential effects of changing climate on ecosystem resilience, and to identify any overlap between HRV and FRV time series that might indicate possible future target areas for resilience management (Hansen et al. 2014). We use simulated data for a landscape in western Montana, USA, to illustrate this method. Our approach is intended for scientists and managers who are attempting to integrate measures of ecological resilience into research and land management planning.

Background
The ball-and-cup model is the most widely recognized concept of ecological resilience (Holling 1973;Fig. 2). An ecosystem (the ball) remains within a set of bounding conditions (sides of the cup), representing a resilient state because when perturbed, recovery trajectories are convergent to the bottom of the cup (by gravitational analogy). Perturbations-such as wildfires or extreme climate episodes-can act on the ecosystem with sufficient force to move the ball across a threshold (cup lip) into a different resilient state (new cup). As a heuristic, the ball-and-cup model is intuitive (Nolting and Abbott 2016), but is limited in its suitability for translation to operational metrics (Desjardins et al. 2015), particularly for complex, multiscale, and highly dynamic disturbance-adapted ecosystems. For example, in the simple heuristic, the ball returns to the cup bottom following perturbation; this implies that ecosystems always return to the same stable, steady state after disturbance, and remain in this state until perturbed again. In nature, ecosystems are highly dynamic over time and space and rarely return to the same state after disturbance, or by the same trajectory (Suding et al. 2004). Recovery pathways are often hysteretic, and a range of forces are always acting on the position and trajectory of the ecosystem (ball), influencing pre-and post-disturbance structure, composition, and distribution of the biota White 1985, Pickett et al. 1987). Further, the balland-cup model does not portray ecosystems' spatial heterogeneity (i.e., large, flat cup bottom) or temporal variability (i.e., the shape of the cup can change over time; Desjardins et al. 2015), attributes that influence biotic diversity, resistance to disturbance, and overall ecosystem function (Hansen et al. 1991, Lehman andTilman 2000). The slope and height of the sides of the cup represent the magnitude of the perturbation needed to move the ball across the lip of the cup. However, the magnitude is always relative to the force acting on the ball and the position of the ball in the cup; some forces, such as exotic invasions or decadal climate change, move the ball to a new state by accumulated pressure (Stevens and Falk 2009). Other forces, such as large punctuated disturbances, rapid human development, or extreme climate episodes, can shift ecosystem states abruptly. Lastly, the ball-and-cup model is implicitly univariate, whereas community and ecosystem resilience are inherently multivariate problems with multiple driver and response variables that may not act in concert (Gunderson 2000, Desjardins et al. 2015, Palmer et al. 2016. Important in the context of climatic changes and other disturbances, the ball-and-cup model can reflect both discrete (pulse) disturbances and more long-lasting (press) stressors (Bender et al. 1984). The height of the sides of the cup represents the magnitude of either a single perturbation or the cumulative effects of persistent stress, such as drought or multiple disturbance events. Ecosystems can be perturbed by multiple interacting forces, some of which require little force to move the ball to a new state because they operate across spatiotemporal domains (Creutzburg et al. 2015). For example, invasion of sagebrushsteppe communities in the U.S. Great Basin by the non-native annual cheatgrass (Bromus tectorum) stresses native communities both via direct resource competition and indirectly through changes in the fire regime (Whisenant 1990, Billings 1992) that can convert broad areas rapidly to persistent cheatgrass monocultures. Other ecosystems may experience apparently severe disturbances that rarely drive the ball out of the cup; the widespread North American interior lodgepole pine (Pinus contorta var. contorta) ecosystem, for example, commonly experiences high-intensity, stand-replacing fires that kill most trees, but the species is adapted to populate postfire environments via cone serotiny (Turner et al. 1999, Schoennagel et al. 2008). Yet, a low-intensity fire that burns young lodgepole pine stands before saplings have produced serotinous cones may shift the ecosystem from a forest to a grassland or shrubland that may require centuries for pine to re-invade (Pierce andTaylor 2011, Clark et al. 2017). This delay in the recovery time does not mean that the system is not resilient, but that the magnitude of disturbance may not necessarily predict the magnitude or duration of ecological response.

Our HRV-resilience model
Our HRV-resilience model can be visualized as the movement of a marble (ecosystem or landscape) in two-dimensional phase space (flat plane or table; Fig. 3). The dimensions of the table define the biophysical envelope or fundamental niche space of key system variables, such as vegetation composition and structure, that are important indicators of ecological pattern and process and are also relevant for management. In this sense, the table is analogous to the basin in the ball-and-cup model, with the difference that by being planar, the table does not imply the existence of a single, stable steady state. A system's fundamental setting is finite, because some forces, such as meteor impacts, human development, and volcanic eruptions, can push the ecosystem (marble) off the table and onto another domain (another biophysical envelope). The marble is constantly acted on by exogenous forces, causing it to move (range) through various states on the table, each an expression of realized niche space. We illustrate four major forces (arrows) in this paper-climate, disturbance, ecological succession, and human activity-but other forces could be used to represent the leading influence on specific ecosystems, or to account for evolutionary processes. The direction of the arrows can vary over time and may or may not be consistent and equal among forces.
Because exogenous (climate, disturbance, humans) and endogenous (succession, evolution) forces act continually on the ecosystem, it is constantly in a state of flux as indicated by the path of the marble over time (Fig. 3A). However, during periods of more rapid environmental change, such as the Pleistocene-Holocene transition or the current emergence of the Anthropocene, the path of the system may respond to secular trends in the climate or disturbance signal, developing a non-zero net trajectory in phase space ( Fig. 3B; Bond et al. 2005). Paleoecological records indicate major shifts in species composition, hydrology, and fire regimes in many regions of North America over the last 4000 yr of the Holocene (Gavin et al. 2007, Whitlock et al. 2011), although there is considerable variation as represented by the long trajectory (Fig. 3B). In general, our knowledge of ecosystem paths is confined to the time domain of existing data for the landscape, as there are few detailed data sets of sufficient temporal depth, spatial extent, and appropriate resolution, to evaluate millennialscale historical dynamics (Keane et al. 2009).
Our modeling HRV-resilience method requires two types of data: values for variables used to describe the current state of the landscape or ecosystem and its representative HRV time series (Fig. 3C). Our starting assumption in using HRV ❖ www.esajournals.org as a reference for resilience is that ecosystems operating within HRV for key variables are resilient ipso facto because their behavior falls within the bounds of ecosystem responses (cup in Fig. 2; HRV domain in Fig. 3B). When the ecosystem (marble) is outside of HRV for key variables, we infer that ecosystems are less resilient because ecosystem structure and composition may not facilitate an expected return to HRV conditions following disturbance. Thus, the departure of current conditions from HRV can be used operationally as an index of resilience: The lower the departure value, the more resilient the landscape.
Variations of this method could be used to evaluate resilience of the contemporary landscapes or ecosystems with respect to future climate change, exotic invasions, or land management activities. Managers can augment HRV departure analysis with similar analyses of current conditions compared to simulated time series of future ranges and variations (FRVs) of a Fig. 3. The marble and table analogy to illustrate the concept of resilience. (A) The marble's position represents the state of a hypothetical forest system in two dimensions at any given time; the path varies over time depending on the spatial and temporal domain of evaluation and the environmental forces (arrows) acting on the marble. In this example, the forces are climate, disturbance, humans, and succession. (B) Realignment of drivers for the last 4000 yr produces a trend driven by warming climate and increasing rates of disturbance. (C) HRV (colored red line circle) is defined using 800-100 yr BP because of limited historical data and highest relevance to contemporary conditions. Cumulative human influence on ecosystem states and trajectories remains relatively small and non-directional during this period. (D) Recent (last 100 yr) conditions have pushed the ecosystem (current state) outside of the 800-yr HRV envelope, as climate and disturbance regimes become more aligned and stronger, forcing the realignment of successional trajectories. Human influences become dominant during this period. Projections of the ecosystem under future climate by simulation models (blue and red lines) may or may not overlap with HRV, illustrated here for two future climate change scenarios (IPCC RCP4.5-blue and RCP8.5-red). landscape or ecosystem ( Fig. 3D; , Hansen et al. 2014. With different assumptions about climate (e.g., low or high emissions) and land management (e.g., fuel treatments, fire suppression), various future scenarios can be modeled. From these projections, we can determine whether FRV envelopes overlap with the HRV envelope. Where overlaps exist, they can be evaluated as possible focal management areas for resiliency. We acknowledge that the estimation of FRVs using simulation entails uncertainty as mentioned previously. Nonetheless, it can represent quantitatively our best expectations of future conditions and thereby provide useful information informing the design and implementation of possible restoration measures (Hansen et al. 2014).

Creating HRV and FRV time series
To apply the HRV-resilience method, we used simulation modeling to derive time series representing HRV and FRVs. Simulations provide the necessary temporally deep, spatially explicit historical data that can be difficult to obtain elsewhere (Humphries and Bourgeron 2001, Keane et al. 2009, Dickinson 2014. Moreover, modeling provides a single, consistent platform for generating the required data to characterize HRV for multiple ecological attributes and for generating projections of FRV under future climates. Importantly in the current context, modeling allows the generation of multiple FRVs that include both future climates and various resource management strategies. We recognize the limitations of a simulation approach to quantify HRV and FRparticularly all models simplify reality and are subject to bias from input parameters and model mechanics. Because of this, it is essential to perform model validations to quantify uncertainty in predictions (Keane 2012, Keane et al. 2015. The example presented here was generated using the mechanistic landscape model Fire-BGCv2 ) as implemented for the 128,000-ha East Fork of the Bitterroot River (EFBR) watershed, located in the interior northern Rocky Mountains in the Bitterroot National Forest, Montana, USA. The lower elevations of the EFBR are dry, mixed-conifer ecosystems of ponderosa pine (Pinus ponderosa var. ponderosa) and Douglas-fir (Pseudotsuga menziesii var. glauca) generally with a primarily frequent, low-severity fire regime (Holsinger et al. 2014). Vegetation at montane elevations is mixed-conifer forest (primarily lodgepole pine, Douglas-fir, and subalpine fir (Abies lasiocarpa)) with mixed-severity fire regimes, while high elevations are whitebark pine (Pinus albicaulis), subalpine fir, and spruce (Picea engelmannii) forests with a long fire-free intervals and relatively high-severity fire regimes ( Table 1). The EFBR has been used in past FireBGCv2 simulation studies with its initialization, parameterization, calibration, and validation described in various papers , Holsinger et al. 2014, Clark et al. 2017. Holsinger et al. (2014) parameterized and calibrated the model for the EFBR fire histories that represent fire regimes approximately 500 yr BP. We simulated historical conditions (i.e., no fire suppression or land management, using a cycling daily weather stream from 1950 to 2005) and three future scenarios that combined future climate (CRM-C5 RCP8.5 [+5.5°C, 95% baseline precipitation]) with three wildfire management scenarios representing varying levels of fire  Table 2). Field data measured in 2009-2010(Holsinger et al. 2014) were used to quantify contemporary conditions and to initialize 500-yr simulations, and we used only the last 400 yr of output to eliminate the influence of initial conditions. We produced five replicates per scenario and produced output at ten-year intervals, for a total of 200 observations per scenario. We analyzed 14 landscape-level response variables commonly used in land management (Table 3).

Quantifying resilience using single and multiple response variables
We used box-and-whisker plots to summarize HRV for one variable (univariate) and principal components analysis (PCA) to describe HRV for multiple variables (multivariate) to quantify resilience and produce resilience indices based on the departure of current conditions from HRV. Comparison of these two approaches illustrates relative strengths and drawbacks of using multiple Idealized fire management (50% fires suppressed); CRM-C5 RCP8.5 climate (+5.5°C, 95% baseline precipitation) FRV3 Contemporary fire management (98% fires suppressed); CRM-C5 RCP8.5 climate (+5.5°C, 95% baseline precipitation) Notes: We simulated the EFBR landscape for 500 yr for the HRV and FRV scenarios. The landscape was inventoried in 2009-2010 to describe current conditions, which we used to quantify initial conditions (Holsinger et al. 2014). We generated values for seven output variables (Table 3) every 10 yr for the last 100 yr of simulation for each of five replicates for a total of 100 observations (we removed the first 100 yr to eliminate the influence of initial conditions). variables to represent HRV in the context of shifting climate and disturbance regimes.
Univariate HRV. -Our univariate method uses a single key variable to evaluate the departure of current conditions from HRV and FRV to compute resilience indexes. We used average tree basal area (BA, m 2 /ha) calculated across all stands in the EFBR landscape from FireBGCv2 output; BA is an important metric to forest managers as it is used as a measure of forest productivity and biomass. Boxes in univariate box-and-whisker diagrams show median, 25th, and 75th percentile BA (Fig. 4). We used Student's t-tests (a = 0.05) to compare current to HRV and FRV time series.
Range and variation of BA differed among the four scenarios (Fig. 4). The current BA (Present) is well within the HRV interquartile range (IQR, 25-75th percentile), indicating that this variable is not significantly departed from the modeled HRV (P < 0.001). The median BA for the FRV3 scenario also falls within the HRV IQR (i.e., is not significantly departed). The FRV3 IQR is narrower and median BA is higher (pairwise t-tests, P < 0.001) than for the FRV1 and FRV2 scenarios because the high (98%) level of fire suppression implemented in this scenario minimizes firecaused biomass loss. The FRV2 and FRV1 scenarios produce progressively lower median BA, consistent with expected lower fire suppression levels and increased fire-caused biomass loss. The zone of overlap in the IQR between each FRV scenario and the HRV scenario indicates the percent of simulation replicate-years (200) where the comparison variable (BA) resides within the HRV distribution of that variable. For FRV2, BA responses are outside of the HRV-resilience envelope for at least half of the simulation years (P < 0.001). There is almost no zone of overlap between HRV and FRV1 because tree mortality from frequent fires and likely climate stress results in persistently lower BA than the HRV reference. Fig. 4. An illustration comparing historical (HRV) and future (FRV) variability in basal area (m 2 /ha) variability compared with current conditions on the EFBR landscape (Present: the initial conditions at the start of the simulation). There appears to be a zone of overlap between HRV and FRV3, which may provide a possible reference for management. FRV1, FRV2, and FRV3 are future simulations with RCP8.5 climate with 0%, 50%, and 98% of fire ignitions suppressed, respectively. The box in this figure represents the 25th and 75th interquartiles, and the whiskers represent the range of the data.
In the univariate method, a simple percentile number can be used as a resilience index. In our example above, we calculated the percentile in which the current (Present) landscape BA resides within the HRV distribution of BA (Present was in the 64th percentile) and then used that percentile as a resilience score (64) where 50 would be high resilience and below 25 and above 75 would be low resilience. Other central tendency statistics can be used to determine where in the HRV BA probability distribution is the current value for BA and whether it is significantly different (departed) from the HRV value. For current BA in the EFBR, the probability of the current landscape condition in the HRV distribution is 0.69, which is less than our designated alpha level (P > 0.05), so this landscape could be considered resilient. We also calculated a resilience index for each FRV scenario where pairwise t-tests indicate significant departure from HRV (P > 0.05; Steele et al. 2006). In our example, the two scenarios where BA was significantly departed from HRV are FRV1 and FRV2.
Multivariate HRV. -We used PCA to assess multivariate landscape resilience for the 14 variables ( Table 3) that were used to represent an HRV and the contemporary conditions (Fig. 5). Relative component loadings show that NPP, FLOW, and three structural stages (Sap, Mat, and VLrg) have the least influence on landscape dynamics (e.g., exert the least force on the marble on the table). Multivariate approaches capture the interplay Fig. 5. Results of the principal components analysis of the output from FireBGCv2 for the HRV scenario ( Table 2) that represents the movement of the EFBR landscape in multivariate space. Each dot is a year of simulation, and each color represents a century. The condition of today's landscape is shown by the green asterisk in the bottom left of the graph, which was the initial condition of the landscape at the start of simulation. It is obvious that conditions of today's landscape are well outside of the point cloud portrayed by the HRV simulation across the multiple variables in Table 3. among multiple structuring factors that can be used to characterize complex ecosystems (Dray et al. 2012). We defined the dimensions of our HRV analysis space using the first two principal components (PC1, PC2), which together explained about 60% of the variance in the simulation variables (loadings of each of the 14 variables for both PC1 and PC2 are detailed in Table 3). The track of Fig. 6. Results of PCA of FireBGCv2 simulations for the EFBR landscape for the historical scenario (HRV; blue dots, reference) and for the three future scenarios (red dots; FRV1, FRV2, FRV3; Table 2) that are three fire suppression scenarios (no suppression, 50% ignitions suppressed, 98% ignitions suppressed) under an RCP8.5 climate (A, C, E). Also shown are simplifications of the scatterplots with circles that contain 68% of the variation in the spread of the points for the three FRV scenarios (B, D, F). The green asterisk at the lower left of graphs A, C, and E represents the condition of today's landscape. FRV1 is shown in A and B,FRV2 in C and D,and FRV3 in E and F. Variable names in B,D,and F are defined in Table 3 and indicate the important variables in the PC1 and PC2 scores. the marble (ecosystem) on the table (analysis space) is the connection of each reported year (dot) in PC1-PC2 space for all simulated years across all five replicates (Fig. 5). Unlike the univariate BA analysis, the current condition of the EFBR landscape as measured in 2010 (green asterisk) lies well outside the point cloud for HRV (i.e., departed from historical reference conditions). The large departure of years 100 to 120 from the primary HRV cloud is because it takes more than a century to eliminate initial effects from the model, and the initial EFBR conditions reflect 100 yr of fire exclusion. Also evident in the graph is the effect of drift in ecosystem condition perhaps from climatic changes (Keane et al. 2008).
Comparison of PCA results across the three FRV climate and fire management scenarios and HRV scenario provides insight into the potential impacts of changing climate and fire regimes on future landscape resilience (Fig. 6). Unlike results from the univariate analysis (BA; Fig. 5), all three fire management scenarios (0, 50, and 98% fire suppression) under RCP8.5 climate depart from HRV, especially FRV3. Moreover, the state of the contemporary landscape (green asterisk) is well outside the PC1-PC2 point clouds of the HRV and of all three FRVs, indicating that it has low resilience when multiple variables are used, regardless of climate or fire management scenario. This illustrates the value of using multiple variables when evaluating resilience. The zones of overlap among the three future fire management scenarios and the HRV scenario become smaller as suppression increases; the overlap for FRV1 and FRV2 scenarios (Fig. 6B, D, F) includes all of the HRV space, such that any treatment or wildfire that moves the landscape toward HRV should be more viable in the future. There appear to be two separate point clouds for each of the three FRV scenarios (Fig. 6A, C, E) which is the result of the slow ramping of predicted climate over the first 100 yr of simulation (Holsinger et al. 2014).
Creating a resilience metric from multivariate time series is a bit more difficult. We used a statistical approach to determine if the current condition is significantly outside the PC1-PC2 centroid, and the magnitude of that distance relative to the HRV point cloud provided us a suitable resilience index. In our example, we arbitrarily set the confidence ellipse of HRV at 0.68 (i.e., an ellipse containing 68% of HRV observations). We then computed the probability of obtaining the value for the current EFBR landscape (Present) as a random draw from the HRV point distribution. This test returned a value of 0.005, indicating substantial departure of current conditions from HRV; this probability could be used as a departure value and resilience index. The multivariate PCA HRV-resilience approach can also employ box-and-whisker diagrams using the scores of PC1 from the HRV time series to compare to current conditions similar to our use of BA above.

Managing for resilience
Enhancing resilience, especially in fireexcluded forests, will probably entail some degree of either ecosystem restoration or acceptance of change , Falk 2017. Many ecological restoration projects have used HRV of ecosystem characteristics, such as fire regimes, as a reference for designing treatments that can return stands or landscapes to historical envelopes, which are generally regarded as ipso facto resilient (White and Walker 1997, Landres et al. 1999, Dickinson 2014). Future restoration treatments should be designed to move the current landscape in the direction of HRV, while also tracking anticipated future conditions (Falk 1990). For example, the current BA in EFBR in Fig. 4 is well within HRV, but BA is lower in FRV1 and FRV2 scenarios. To keep EFBR near HRV and prepare for the uncertain future, wildland fire use (letting some wildfires burn under acceptable conditions) in remote areas or silvicultural thinnings combined with prescribed burning could reduce BA, especially the BA of fire-sensitive species, and enhance the vigor and growth of fire-adapted species in the EFBR landscape, leading to an overall increase in future resilience (Stephens et al. 2009). With this limited univariate example, it is interesting that greater fire suppression in the future (FRV3) may keep the EFBR landscape within HRV under new climates, but other important variables that are not represented in the single variable analysis might change this result (Fig. 4).
Designing treatments that enhance resilience based on multivariate analysis may be more challenging because of the reduction of many variables to a small number of synthetic axes.
However, the results presented in our multivariate example (Fig. 6, Table 2) suggest the possibility of managing for the variables with the highest factor loadings as important drivers of departure (Table 2). In this case, PCA suggests that resilience may be most positively affected by increasing seedling (Seed) stands using wildland fire use and increasing large (Lrg) diameter structural stages by thinning and prescribed burning in overly dense stands, which would ensure that some smaller diameter pole structural conditions persisted to grow into the large tree-dominated structural stages. Appropriate actions to move the EFBR toward a zone of HRV-FRV overlap might include reducing the level of fire suppression (increasing wildland fire use), because the greatest FRV-HRV overlaps occur when some wildfires are allowed to burn (0% and 50% suppression), and providing protection from uncontrolled wildfire for the large trees on the landscape via radial thinning, fuel treatments, and prescribed burning to retain those late seral structural stages.
Managers may want to aim for targets in the center part of the overlap zone or even below the zone to account for the future dynamics of the system. Such anticipatory management may increase treatment effectiveness and longevity Stephenson 2015, Bradford et al. 2018). For example, in Fig. 4 (univariate case) there is little overlap across distributions of HRV and FRV1; HRV basal areas are higher than FRV1 basal areas because of high fire occurrence in the future. Land managers, however, may de-emphasize treatments that move the landscape toward FRV1 because of uncertainty in both climate and management projections for the future (Rupp et al. 2013). Instead, managers may aim for basal areas that are closer to the low end of the HRV distribution (lower basal areas) while also falling within FRV1. In this sense, approximating the HRV remains a reasonable null criterion rather than aiming for projected targets with high climatic and ecological uncertainty (Mellert et al. 2015, Luce et al. 2016).

SUMMARY
In this paper, we demonstrate a simulation modeling method for quantifying resilience for the purpose of guiding management decisions.
The operative assumption of this HRV-resilience method is that ecosystems are resilient when they resemble conditions represented by the historical range of variability (HRV), because the HRV represents the recent conditions under which most of the biota have evolved, and under which contemporary ecological communities have formed (Bunnell 1995). At the core of our method is the use of simulation modeling to quantify the HRV time series, which are compared statistically to current conditions to compute a resilience index. Looking forward, rapid climate change, ongoing land degradation, altered and accelerated disturbance regimes, exotic species invasions, and a host of other human impacts that are occurring today and into the future demand a still broader assessment than using HRV alone (Falk 2017). To address this challenge, we integrated a companion future range of variation (FRV) evaluation into our method to help guide future land management planning. This augmentation is only possible using simulation modeling. Overlaps between HRV and FRVs may provide ideal targets for specific management-oriented environmental variables. We demonstrate variations of the method using both univariate and multivariate approaches. Inclusion of multiple variables into the assessment can be accomplished using PCAs to compress multiple variables into a smaller number of axes that define the response space of FRV and HRV. We identified several variables that are important to managers and that drive differences between current, HRV, and FRV most strongly (highest factor loadings), and then identify treatments that enhance these variables. We show examples of how to deploy this method into operational use, even in cases where there is little apparent overlap between the HRV and FRV. In these cases, managers may want to set goals that are within HRV but trending toward the FRV as a hedge, given a future that is unlikely to be similar to the past, but characterized by high uncertainty.