1. Introduction

The old adage that one must learn from the past to understand the present and plan for the future finds its way into land management through the concept of the historical range of variability (HRV). The HRV can be defined as the variation in ecological processes (e.g. disturbance regimes such as fire) or characteristics (e.g. forest composition and structure) for a specific time and place [1]. The concept has been used to inform ecosystem restoration [2, 3], conservation action [4] and land management planning [5, 6, 7]. It offers both a benchmark against which current ecosystem conditions can be measured [8] and provides an understanding of the processes that drive ecosystem change over time which can inform land management activities [9].

Despite their utility, data to help us understand the HRV are often lacking [10, 11]. Detailed information about ecosystem characteristics from tree rings or sediments, such as age distribution, structure characteristics and disturbance regimes, are available for very few ecosystems, expensive to collect and, even when available, may be limited in geographical and temporal extent [12]. In the absence of data, models can be used heuristically to explore a range of possible past conditions and improve our knowledge of how ecosystems respond to change over time [1]. The state-and-transition simulation modeling framework is one approach that has been employed for estimating the HRV [2, 4, 13, 14].

State-and-transition simulation models (STSMs) divide ecosystems into their component state classes and define rates and pathways for probabilistic and deterministic transitions between them [15]. In the LANDFIRE STSMs, deterministic transitions are used to represent growth or ageing-related changes between classes, and probabilistic transitions are used to represent the effects of disturbance processes. Software tools such as the Vegetation Dynamics Development Tool (VDDT) [16], and its successor ST-Sim [17] can be used to assign probabilities to the transitions and to stochastically simulate multiple iterations of the model [15]. Simulation results include the percent of the landscape in each state class over time, which has been used as the basis for HRV estimates [18, 19].

A challenge with this type of modeling approach is to develop reasonable estimates of the range of possible reference conditions for an ecosystem. Many ecosystem assessment efforts in the United States have implemented the Fire Regime Condition Class (FRCC) method which uses STSMs to estimate an average reference condition against which the current landscape condition can be compared [19, 20, 21, 22]. The FRCC reference conditions are typically represented as the mean relative abundance of each state class calculated by averaging the results of multiple Monte Carlo simulations [4, 23, 24]. Lack of historical data presents a significant challenge to simulating a range of possible reference conditions. Moreover, the expense and time associated with gathering the necessary historical data may be prohibitive [25].

Various approaches have been tried to estimate the range of variability around mean reference conditions of various ecosystems. Haugo and others [13] estimated HRV as the stochastic variation present in multiple iterations of an STSM, an approach that captures state class variability due to model stochasticity but not due to uncertainty in transition probability estimates or variability in disturbance regimes over time. Similar non-spatial modeling efforts have adjusted the disturbance probabilities using a time series of multipliers in an STSM to estimate state class variability over multiple Monte Carlo simulations [2, 4, 14]. The time series are generally derived using some combination of data and expert judgment. These approaches approximate a more ecologically meaningful HRV, but the data required to generate the time series are frequently unavailable for the time period and landscapes where HRV estimates are required. Other researchers have used more computationally and data-intensive spatially explicit simulation models to estimate the HRV [26, 27]. The data needed to parameterize these models are extensive and frequently lacking [12], and the models themselves may be so complex that they are difficult to use [1]. Despite the data limitations, we do understand that there is not a single disturbance return interval over time and we do have some knowledge of the range of possible disturbance frequencies and/or estimates of the uncertainty in these ranges. Our method fills a gap in the existing approaches by incorporating information on disturbance ranges and/or the uncertainty in the range estimates for broad areas where detailed time series data, spatial data or spatial modeling expertise are not available.

Our goal was to develop an STSM approach for estimating a reasonable and ecologically relevant range of state class variation that 1) relied primarily on existing information and 2) did not require extensive modeling expertise or input data. We chose to use LANDFIRE’s STSMs which describe the vegetation dynamics for more than 1,000 ecosystems, called Biophysical Setting (BpS), in the United States [19]. A BpS represents a vegetation community that may have existed prior to Euro-American settlement based on the current biophysical environment and an approximation of the historical disturbance regime [19]. Each BpS has both an STSM and an associated description document that describes the vegetation concept and its model.

Every disturbance in a LANDFIRE STSM is assigned an average probability of occurrence. We assumed that if we could sample from a reasonable probability distribution around that average for the modeled disturbances, that we could estimate the bounds of uncertainty for the HRV in the state classes over time. There is little information available about the historical frequency distribution of various disturbances, but the LANDFIRE STSM description documents sometimes include information about the estimated range of fire frequencies that we were able to use. The results documented here demonstrate the range of state class variation resulting from varying fire frequencies across model iterations. We only varied fire frequencies because it was simpler to demonstrate the method using a single disturbance, and fire was the only modeled disturbance that contained information on the minimum and maximum return interval. The same approach could be used to incorporate the variability in other disturbance regimes if data were available. This paper offers a proof of concept for this approach by demonstrating its application to three BpS representing different fire regimes.

2. Materials and Method

We tested our approach on three LANDFIRE STSMs: Northern Rocky Mountain Dry-Mesic Montane Mixed Conifer Forest (Mixed Conifer; BpS number 10450) Great Basin Pinyon-Juniper Woodland (Pinyon-Juniper; BpS number 10190) and Boreal White Spruce-Fir-Hardwood Forest-Coastal (Spruce-Fir-Hardwood; BpS number 13652; Figures 1 and 2) [28]. These BpS were chosen because they represented a range of biophysical environments and disturbance regimes and each contained all the information required for our technique - an estimate of the minimum, mean and maximum fire return interval. The model documentation for each BpS includes a mean fire return interval (FRI) for each fire severity type (LANDFIRE models use surface, mixed and replacement fire severity types) in the model, which was the result of a 1,000 timestep simulation in VDDT. These particular models also included estimates of the minimum and maximum FRI for each fire severity type based on the modelers’ knowledge of the literature and experience with the BpS (Table 1).

Figure 1. STSM box and arrow diagrams: A) Mixed Conifer, B) Pinyon-Juniper and C) Spruce-Fir-Hardwood. Each box represents a state class which is defined by its cover type and structural state (labeled in the upper left of the box) and its age range (labeled in the lower left of the box). Fire pathways are shown with their transition probability.

DownLoad: Full-Size ImgPowerPoint

Figure 2. Distribution map showing the spatial extent to which the Mixed Conifer, Pinyon-Juniper and Spruce-Fir-Hardwood STSMs were applied [26].

DownLoad: Full-Size ImgPowerPoint

Table 1. FRIs by severity type for each BpS. The all fire return interval is the aggregate FRI for all fire severity types. The mean and all FRI values were the result of the LANDFIRE STSMs. The minimum and maximum FRIs were estimated by the LANDFIRE modelers.

Biophysical Setting	Fire Severity Type	Fire Return Interval (years)
		Minimum	Mean	Maximum	All Fire
Mixed Conifer	Surface	20	30	35	20
	Mixed	70	110	175
	Replacement	70	135	200
Pinyon-Juniper	Surface	5	715	1,000	166
	Mixed	10	370	1,000
	Replacement	10	525	1,000
Spruce-Fir-Hardwood	Surface	NA	NA	NA	526
	Mixed	300	1,250	1,500
	Replacement	300	909	1,500

 | Show Table

DownLoad: CSV

For each fire probability, we used a beta distribution with parameters that satisfied the minimum, mean and maximum FRI values identified in the LANDFIRE STSMs. The beta distribution was chosen because it can be used to represent the variability within a fixed range [29] such as the minimum and maximum FRI. This allowed us to restrict our distribution to a range within the bounds specified. However, because the beta distribution can assume a variety of shapes depending on the value of the alpha and beta parameters that define it, we had to make assumptions about the shape. Lacking information about the historical FRI probability distributions, we heuristically varied the variance in probability distribution in order to select a value that resulted in a bell-shaped FRI distribution, with values close to the mean being the most probable and those towards the specified minimum and maximum being the least probable (for more information see section 2.1.3). We ran simulations in ST-Sim using the transition multiplier function to sample from the beta distributions we defined for each fire severity type (for more information see section 2.1).

The Mixed Conifer and Spruce-Fir-Hardwood BpS models were simulated in ST-Sim for 100 iterations, over 1,000 timesteps, using 10,000 simulation cells. In the absence of information about initial conditions, we initialized the models with an equal proportion of cells in each state class. For each iteration, one FRI was sampled for each fire severity type from the fire frequency distributions. To describe the range of variation in the simulations, we summarized the percent of cells in each state class and the minimum and maximum values for the last 800 timesteps for each BpS. By graphing the results we determined that within first 200 years the model results stabilized and overcame the influence of the initial state class distribution.

The Pinyon-Juniper STSM took longer to reach equilibrium so we adjusted it to run for 1,400 timesteps and summarized the results for the last 800 timesteps. For this BpS we pooled the Late1 Open and Late2 Open state class results into a single Late Open class because we were aware, based on our contribution to the LANDFIRE mapping effort, that the LANDFIRE current condition spatial data layer was unlikely to distinguish these classes from each other accurately. LANDFIRE maps current conditions primarily based on vegetation cover and height [19], but this particular model specified that the late classes should be distinguished based on tree size class which was not a LANDFIRE current condition mapping criterion. In all other regards, the Pinyon-Juniper STSM was parameterized and summarized similarly to the other BpS.

To illustrate the application of the simulated range of state class variation in a management context, we calculated the current percent of the BpS in each state class by combining the LANDFIRE 2010 BpS [30] and Succession Class [31] spatial data layers in ArcGIS. The current condition data layers include the percent of the landscape in each of the state classes and in an uncharacteristic class. The uncharacteristic class includes both exotic species composition and structures (typically either cover or height) that would not have occurred under the reference condition.

Calculating the simulation parameters

Transition multipliers were used to vary the simulation fire probabilities. The parameters required for this function are the mean probability multiplier value, the distribution type, the standard deviation and the minimum and maximum multipliers (Table 2).

Table 2. Transition multipliers were used to vary the simulation fire probabilities for each fire severity type in the model.

BpS	Transition Group	Standard Deviation	Minimum Multiplier	Maximum Multiplier
Mixed Conifer	Surface Fire	0.0700	0.8571	1.5000
	Mixed Fire	0.1700	0.6286	1.5714
	Replacement Fire	0.1700	0.6750	1.9286
Pinyon-Juniper	Surface Fire	0.1800	0.7150	143.0000
	Mixed Fire	0.3500	0.3700	37.0000
	Replacement Fire	0.3100	0.5250	52.5000
Spruce-Fir-Hardwood	Surface Fire	--	--	--
	Mixed Fire	0.1300	0.8333	4.1667
	Replacement Fire	0.2200	0.6060	3.0300

 | Show Table

DownLoad: CSV

2.1. Multiplier

The multiplier is used to adjust the transition probabilities by the specified factor. In all cases we set the mean of the multiplier distribution to 1.000 so that on average across all iterations the mean FRI was equivalent to the value originally specified in the source LANDFIRE model.

2.2. Distribution type

ST-Sim allows users to specify either a normal or a beta distribution. We used the beta distribution because it can be bounded between the minimum and maximum multipliers specified.

2.3. Standard deviation

We selected a standard deviation of the beta distribution for each fire severity type heuristically by visualizing the resulting changes in the frequency distribution shape and selecting the final simulation standard deviation that best approximated a bell-shaped distribution while maintaining the widest range of possible values (i.e. the minimum and maximum FRI; Figure 3).

Figure 3. The fire frequency distribution was selected heuristically to approximate a bell shape while maintaining the widest possible range of values (i.e. the minimum and maximum FRI). This example shows the shape of the FRI distributions for surface, mixed and replacement severity fires for the Mixed Conifer BpS.

DownLoad: Full-Size ImgPowerPoint

2.4. Minimum and maximum multipliers

The minimum and maximum multipliers are used to bound the distribution to satisfy the minimum and maximum FRI values specified in each LANDFIRE model. We used the inverse of the estimated minimum and maximum relative to the inverse of the mean FRI values for each fire severity as reported in the model documentation. By taking the inverse of the minimum, maximum and mean, the FRIs are converted to annual probabilities. The minimum and maximum probabilities are then relativized to the mean probability so that across all model iterations the mean FRI is equal to the mean FRI specified in the LANDFIRE model. The minimum and maximum multipliers were calculated using the formulae:

3. Results

Varying the fire frequency did result in a range of values for each state class, and some state classes varied more than others (Figure 4). The Mixed Conifer BpS varied least in the Mid1 Closed class and most in the Late1 Open class. The current distribution for the Early1 All and Mid1 Open classes fell within the simulated range, but the Mid1 Closed and Late1 Closed classes were approximately 12 percentage points above the maximum simulated value and the Late1 Open was 25 percentage points below the minimum simulated value (Figure 4A).

Figure 4. Results of HRV simulations showing the state class mean, range and variability for the (A) Mixed Conifer (B) Pinyon-Juniper and (C) Spruce-Fir-Hardwood BpS.

DownLoad: Full-Size ImgPowerPoint

The Pinyon-Juniper BpS varied least in the Early1 Open class and most in the Late Open. The current amount was outside of the simulated range for all classes except Early1. Although the Mid1 Open and Mid2 Open current amounts were outside of the simulated range, both were within eight percentage points of the simulated mean. The Late Open class had a 20% spread between the minimum and maximum values, the largest spread of any of the state classes for the BpS presented here (Figure 4B).

The Spruce-Fir-Hardwood BpS exhibited the least variation in the Early1 Closed class and the most in the Late1 Closed class. The current amount was outside of the simulated range for all classes except the Early1, and the Uncharacteristic class accounted for 40% of the current condition. The current amount of the Late1 Closed class was 66 percentage points below the minimum simulated value for the class, the widest gap between the current amount and the simulated range for the three BpS. This BpS exhibited the narrowest range of variation across classes and had the greatest proportion of area mapped into the Uncharacteristic class of the models we tested (Figure 4C).

4. Discussion

Using our method, the simulated range of state class variability is largely a function of the interaction between the state class age ranges and the fire frequency distribution. One might anticipate that the greater the spread in the FRI range (as translated into minimum and maximum transition multipliers), the greater the resulting state class range should be. If this were the case, then the Spruce-Fir-Hardwood model, which had the largest FRI spread, would have the greatest range of state class variability, but this was not the result. In fact, the relationship between the model inputs and the results was more complex and reflected the fire regime characteristics of each BpS.

The dynamics of the Mixed Conifer model were largely controlled by its most frequently modeled disturbances, surface and mixed severity fires. Surface fire that maintained the open classes occurred on average every 22‒25 years, and mixed fire that transitioned closed to open classes occurred on average every 33 years. These disturbances mediated the dynamic between the open and closed classes in the model, and they occurred frequently enough to affect the Mid and Late classes that had start ages of 30 years and above. Although replacement fire occurred less frequently than the other fire severities, it occurred often enough to increase the variability in the Early 1 class.

The Pinyon-Juniper BpS was most sensitive to mixed severity fire, its most frequently modeled disturbance. Despite the wide mixed FRI range (10 to 1,000 years), the average overall mixed FRI was still relatively infrequent at 370 years. On average, less than one fire would have occurred for any given simulation cell in this BpS before it reached the Late2 Open class, which started at 300 years. In the Late2 Open class the wide mixed FRI range had the greatest impact and resulted in the greatest range of state class variation in the Late Open class (which pooled results for Late1 Open and Late2 Open) for the BpS.

The fire dynamics of the Spruce-Fir-Hardwood model were dominated by infrequent replacement severity fires occurring about every 909 years. On average, replacement fire would only occur 1.1 times for any given cell during a 1,000 timestep simulation, and given the low disturbance rate, most cells would reach the Late1 Closed class (with a starting age of 71 years) prior to any replacement fire event occurring. Therefore, the Late1 Closed class exhibited the widest state class range for the Spruce-Fir-Hardwood BpS, but it had relatively narrow state classes ranges compared to the other BpS we tested. It should also be noted that this BpS was the only one we tested where fire was not the dominant disturbance, and therefore, other disturbances played an important role in the dynamics of the model.

4.1. Management implications

Simulating the potential HRV in state class distribution can provide for a more nuanced assessment of current vegetation conditions. For example, the current amount of the Pinyon-Juniper classes, with the exception of the early state, fell nearly within the extent of the simulated range. The HRV results illustrate that these classes, while departed from the mean, are not too far outside of the possible simulated range. Similarly, the Mixed Conifer Mid1 Open class could be considered departed based on determining departure from the mean, even though its current amount was well within the simulated state class range. In contrast, the Mixed Conifer Mid1 Closed, Late1 Open and Late1 Closed classes were not only departed from the mean, but they fell well outside of the simulated range. While beyond the scope of this study to explain the difference between the current amount and the simulated HRV, we noted that the Mixed Conifer results were consistent with patterns found in other studies that have documented increasing closed canopy forest conditions in the 20^th century, possibly as a result of fire suppression and other land use changes [32, 33, 34].

Restoration activities that incorporate information about the possible range of state class values could allow for more ecological variation on the landscape. For example, the Mixed Conifer Early1, Mid1 Open and Late1 classes and the Pinyon-Juniper Late Open class had wider simulated ranges of variability than the other state classes within each BpS. Restoration actions based on the HRV could provide managers with greater flexibility in their management actions and allow for the possibility of more ecological variability where Mixed Conifer and Pinyon-Juniper are present.

The Spruce-Fir-Hardwood results depict a situation where a substantial portion of the ecosystem has shifted to an uncharacteristic state that was not present in the reference condition model. We speculate that conversion to northern hardwoods and low canopy cover relative to the reference condition, possibly as a result of past forest management activities, could account for these uncharacteristic conditions. The LANDFIRE model description also notes that spruce budworm (Choristoneura fumiferana) can cause changes in the BpS that would be outside of the HRV [28]. While there may be silvicultural options for shifting the uncharacteristic state closer to HRV, one might question whether that is possible or desirable. In this case, the value of the HRV analysis may be to provide an objective assessment of the current situation and help determine whether these conditions are attainable under current ecological and social conditions.

4.2. Impacts of scale

Our model results for Spruce-Fir-Hardwood BpS show a relatively narrow range of state class variation partially because fire occurred relatively infrequently in the STSM. At the scale of the entire BpS this result may be appropriate, but at finer scales, for example, at the scale of individual watersheds, the replacement severity fire disturbance events that characterize this BpS could cause major fluctuations in state class distribution. These fluctuations cannot be resolved using our method because we did not tie our simulations to an actual geography or spatial extent and our simulations were non-spatial. We also assumed that the mean FRI for each model iteration did not vary temporally. For these reasons, we submit that the technique presented here is best suited to broad scale analyses at least an order of magnitude larger than the upper end of the size distribution of fire events typically observed on a particular landscape. Studies have shown that landscape variability can increase as the size of the modeled assessment area decreases relative to the size of disturbance events [35, 36, 37]. More realistic spatial and temporal variability patterns could be simulated by modifying our method to include the use of a spatially explicit STSM and the addition of temporally varying fire (or other disturbance) probabilities and a fire size distribution.

4.3. Utility and limitations of the STSM approach

An advantage of our approach is that it relies primarily on existing data and software tools to develop an ecologically meaningful estimate of HRV that does not require extensive spatial modeling expertise. The demand for model parameterization using our method is relatively low compared to spatially explicit HRV modeling approaches, and expert judgment can supplement empirical data, if needed. The approach is suitable for the hundreds of existing LANDIFRE models as well as other STSMs that are in or can be transferred into the ST-Sim modeling framework.

In addition to the STSM, our method requires an estimate of the FRI mean and range and the standard deviation of the fire frequency distribution. About 50% of the LANDFIRE models in the conterminous United States and Hawaii specify a FRI range for at least one of the modeled fire severity types; however the level of information supporting the existing estimates varies between BpS [25, 38]. Literature review, field studies and/or expert estimate could be used to create new and refine existing disturbance range estimates. Methods from recent studies on systematically eliciting expert opinion for STSMs [39, 40] could be applied to LANDFIRE models where FRI ranges have not been provided.

We chose to use a bell-shaped beta distribution, but the method is flexible with regard to the type and shape of the distribution. For systems where fire is not the dominant disturbance, other probabilities could be altered to improve the HRV estimates. The models presented here included non-fire disturbances that could be varied using the technique we outlined for fire if data were available, but LANDFIRE models only estimate the mean frequency of non-fire disturbances and do not include estimates of the minimum and maximum return interval.

4.4. HRV in a changing climate

A criticism of management approaches that use the HRV to set management targets is that historical conditions may no longer be relevant under a climate system that is without historical precedent [41, 42]. Management targets based on historical conditions may be increasingly costly to maintain and create ecosystems that are poorly adapted to current conditions [41]. This criticism does not necessarily discount the value of understanding the HRV [43, 44]. In fact, a changing climate may make understanding the historical context of an ecosystem more important for addressing current and future management challenges [8, 44, 45]. Higgs and others [8] suggest that restoration-based ecosystem management activities can use historical information not as a target but as a guide to increase knowledge of how ecosystems responded to change in the past that can provide clues for how they might respond to change in the future.

Information about historical conditions can be used as a guide to plan for anticipated climate conditions and inform the development of forest management strategies consistent with that future climate. For example, in The Rogue Basin Action Plan for Resilient Watersheds and Forests in a Changing Climate, the Southern Oregon Forest Restoration Collaborative and partners targeted treatments that promoted the development of historically open forest structures by removing small-diameter trees and fire-sensitive tree species [46]. These treatments, informed by historical stand structure and composition, should reduce wildfire risk and may over time increase stable carbon stocks [46, 47]. It is expected that these changes will help forests in the Rogue Basin to better adapt to the predicted hotter, drier and more fire-prone climate conditions anticipated in the future [46].

Given the lack of certainty in climate projections, Keane and others [1] suggest that the use of HRV to inform land management may be a prudent approach in the near-term. However, HRV analysis is probably most useful when combined with approaches that account for changing climate conditions and the projected future ranges of variability that define the characteristics of probable future ecosystems [1, 48]. With modifications, the method for determining a range of state class variability presented here could be adapted for exploring future or desired ranges of variability. Modeling future scenarios using this approach may identify policies and management strategies that are sustainable and increase ecosystem resilience as well as identify ecosystems where restoring or maintaining the HRV is unlikely regardless of available resources.

5. Conclusion

Providing objective and reasonable estimates for the HRV for a particular ecosystem is a vexing challenge given the limited historical data available. The approach presented here offers one way to estimate the HRV using readily available data and free software tools that we believe could provide additional and useful information for land management planning and decision-making. For broad scale assessments, this method could provide a first approximation of the HRV that could be refined with local information. The technique is flexible and the assumptions we made about the type and shape of the fire frequency distribution could be changed. Future work in this area could examine the impact of these assumptions on the HRV results. In light of the ongoing debate about the relevance of HRV, we suggest that if, as the old adage reminds us, change is constant, then perhaps a better understanding of past ecosystem conditions will provide useful insights to guide future land management activities.

Acknowledgments

We thank two anonymous reviewers for helpful comments on the manuscript. This work was partially funded bythe LANDFIRE Programthrough Cooperative Agreement USFS 14-CA-11132543-043.

Conflict of Interest

Authors declare no conflicts of interest in this paper.