ROCKY MOUNTAIN LANDSCAPE SIMULATOR (RMLANDS)
Historic Range of Variability Analyses
This section provides a description of the HRV simulation and the methods employed to summarize the results of the simulation.
Simulation Study Design
For purposes of this application, we established a single disturbance scenario representative of the historic reference period (referred to as the HRV scenario). This scenario was intended to allow us to quantify the HRV in landscape structure and wildlife habitat and to serve as the basis for comparison with future management scenarios (to be completed in the next phase of this project). The parameterization of the HRV scenario represented our best estimate of how disturbance and succession processes operated during the reference period (see RMLANDS Model Parameterization). The HRV simulation was conducted with the following constraints:
• Simulation Length - The simulation was run for an 800-year period. Note, the simulation length was not related to the length of the reference period. Instead, the simulation was run for a period equal to roughly twice the longest return interval in any cover type of the dominant disturbance process (wildfire), after accounting for a model equilibration period of roughly 100 years (see below). In this case, the longest return interval was determined to be roughly 350 years in spruce-fir forest. This simulation length (i.e., twice the rotation period of the dominant disturbance) was deemed necessary to adequately describe the dynamic changes in vegetation. This does not mean that we believe the disturbance regime was stationary (i.e., unchanging) for this length of time. Unless otherwise noted, we excluded the first 100 years of the simulation from all results. Because we were interested in characterizing the disturbance and vegetation dynamics under reference period conditions, it was necessary to first transform the current landscape (i.e., starting point of our simulations) into one in balance with the reference period disturbance regime. Preliminary analyses suggested that most measured landscape attributes equilibrated within 100 years (i.e, 10 time steps), and many did so in less time. Thus, for our purposes we simply ignored the first 100 years of the simulation for reporting the HRV in landscape structure and wildlife habitat.
• Time Step -The simulation employed 10-year time steps. Thus, a single 800-year simulation produced 80 output landscapes (i.e., snapshots) at 10- year intervals. The landscape snapshots associated with each time step were summarized to characterize the range of variation in landscape structure.
• Number of Simulations - We ran 5 replicate simulations to characterize the HRV scenario. Note, each replicate simulation was a stochastic outcome of the same model parameterization. Therefore, while the model parameters remained the same under each run, the outcome of each run (i.e., landscape trajectory) was a stochastic realization of that model. In addition, while the exact landscape trajectories varied among runs, the basic patterns of variation in landscape structure remained remarkably consistent. Thus, it was unnecessary to replicate further.
• Spatial Extent - The simulation was run for the full UPL landscape, including a 10-km buffer around the Forest. The buffer was included to avoid boundary effects in the implementation of disturbance processes (see RMLANDS Spatial Data).
Characterization of Disturbance Regime
For each disturbance process, we quantified the following overall temporal and spatial characteristics of the disturbance regime:
• Disturbance Frequency & Extent - To examine the dynamics in the frequency and extent of each disturbance process, we plotted (as a line chart) number of disturbance initiations and total disturbed area (as a percentage of the landscape) pooled across cover types against time. We chose a single representative 800-year simulation run for these plots. In addition, to evaluate the interval between decades experiencing various thresholds of disturbance, we plotted (as a line chart) total disturbed area (as a percentage of the landscape) pooled across cover types against recurrence interval (in years). Here, we defined ‘recurrence interval’ as the number of years between disturbances exceeding a particular threshold in total disturbed area. We combined data across simulation runs (i.e., we ignored the separate runs); hence, each chart represented 350 observations (5 runs x 70 snapshots at 10-year intervals). Lastly, to illustrate the spatial pattern of disturbance produced under ‘extreme’ conditions, we selected a single decade representing the maximum or near maximum total disturbed area and depicted a map of the landscape showing the disturbed cells by mortality level (i.e, low versus high) against all eligible cells as background.
• Climate Effect - To examine the relationship between climate (as represented by the climate modifier parameter) and disturbance frequency and extent, we plotted (as a scatter plot) number of disturbance initiations and total disturbed area (as a percentage of the landscape) pooled across cover types against time. We combined data across simulation runs (i.e., we ignored the separate runs); hence, each scatter plot represented 350 observations (5 runs x 70 snapshots at 10-year intervals).
• Return Interval & Rotation Period - To characterize the overall disturbance regime, we computed the cell-specific mean return interval (i.e., average number of years between disturbances at a single 25-m cell) for each eligible cell and depicted it as a continuous surface map, where the color intensity reflected the mean return interval at that location. For this map we chose a single representative 800-year simulation run, inclusive of the first 100 years (i.e., equilibration period). This map was used to display the geographic variation in cell-specific return intervals. We summarized the cell-specific mean return interval information in a histogram that depicted the percentage of eligible cells that experienced each possible mean return interval. For this purpose, we combined data across simulation runs. This histogram was used to emphasize the variation in mean return intervals. We computed the average mean return interval across all eligible cells (i.e., the mean of the return interval distribution) and referred to this as the population mean return interval. Note, this method of calculating mean return interval is completely unbiased because it is based on the full population of cells and incorporates the complete disturbance history of each cell.
It is important to recognize that return intervals are often calculated differently from our approach. In particular, field studies designed to reconstruct fire histories in areas with low-severity fire regimes often employ dendrochronological methods. This field method involves sampling (i.e. not exhaustively) recorder trees (which may exhibit bias in their propensity to record fires) and observing the intervals between recorded fires on each tree. Here, each recorded interval may be treated as an independent observation. The mean return interval calculated in this manner is likely to be somewhat less than the population mean return interval described above for a number reasons (e.g., the exclusion of long intervals at locations that did not burn, exclusion of the interval between stand origin and the first burn [i.e., origin-to-scar interval], etc.). We refer to the mean return interval derived using this approach (i.e., treating each recorded interval as a separate observation) as the sample mean return interval. In addition, it is quite common to calculate a composite mean return interval, which represents the frequency in which disturbance occurs anywhere within a sampling area of fixed size - usually corresponding to a forest “stand”. The composite mean return interval is profoundly influenced by the size of the sampled area: a larger area will encompass more disturbances and will therefore have a shorter composite mean return interval. In addition, disturbances often only disturb a portion of the specified sampling area, which can substantially inflate the estimate of the mean return interval to any particular location. To help resolve this uncertainty, and to help reduce the effect of small, localized disturbances on the composite mean return interval, some investigators also compute a filtered composite mean return interval, based only on disturbances that are recorded on, say, >25% of recorder trees. These are the disturbances that presumably affect a greater proportion of the stand. The filtered composite mean return intervals are often twice as long as the unfiltered, but filtering does not resolve the fundamental issue of what fraction of a sampling area actually was affected by each disturbance event. In this report, we report all results using the population (i.e., cell-specific) mean return interval approach; readers should exercise extreme caution when comparing our results to published return intervals.
Lastly, we summarized the cell-specific population mean return interval as a rotation period (i.e., number of years required to disturb an area equivalent to the total eligible area). Note, rotation period is equivalent to the cell-specific population mean return interval for the entire landscape - including only eligible cells - and is a nonspatial representation of the disturbance regime because it does not depend on the explicit spatial distribution of disturbances; rather, it depends only on the total area disturbed each time step in relation to the total eligible area. We computed the rotation period separately for low-mortality disturbances, high-mortality disturbances, and any-mortality disturbances (low and high combined). Furthermore, we computed the rotation period (by mortality level) separately for each simulation run and reported the mean, minimum and maximum rotation period across runs. We decided to treat simulation runs as independent observations to highlight the stochastic variation in disturbance regimes.
• Susceptibility Dynamics - To examine the changes in susceptibility of the landscape to disturbance over time, we computed the susceptibility of each eligible cell at the beginning of each time step. We summarized the susceptibility of the landscape at each time step by calculating the proportion of eligible area in 10 equal-interval susceptibility classes (0.1, 0.2, ..., 1.0), as well as the overall average susceptibility across all eligible cells; this was a non-spatial summary of landscape susceptibility. We summarized the range of variation in susceptibility across all 5 simulation runs as follows. We treated each time step as a single independent observation of the susceptibility distribution. Based on the combined data set across runs (350 observations: 5 runs x 70 snapshots at 10-year intervals), we computed the 0th, 5th, 25th, 50th, 75th, 95th and 100th percentiles of the distribution. We compared the current landscape condition (i.e., proportion of eligible area in each susceptibility class, and the overall average cell susceptibility) to this simulated historic range of variability to determine whether the current landscape deviates, and to what degree, from the HRV.
In addition, for each cover type, we quantified the following specific temporal and spatial characteristics of the disturbance regime:
• Disturbance Frequency & Extent - Similar to above, to examine the dynamics in the frequency and extent of each disturbance process, we plotted (as a line chart) number of disturbance initiations and total disturbed area (as a percentage of the cover type) against time. We chose a single representative 800-year simulation run for these plots. In addition, to evaluate the interval between decades experiencing various thresholds of disturbance, we plotted (as a line chart) total disturbed area (as a percentage of the cover type) against recurrence interval (in years). We combined data across simulation runs (i.e., we ignored the separate runs); hence, each chart represented 350 observations (5 runs x 70 snapshots at 10-year intervals).
• Return Interval & Rotation Period - Similar to above, to characterize the overall disturbance regime, we computed the cell-specific mean return interval for each cell of the corresponding cover type and depicted it as a continuous surface map that depicted the geographic variation in return intervals and as a histogram that depicted the percentage of cells of the corresponding cover type that experienced each possible mean return interval. Lastly, we summarized the cell-specific mean return interval information as a rotation period, and did so separately for low mortality disturbances, high mortality disturbances, and low and high mortality disturbances combined. As before, we computed the rotation period (by mortality level) separately for each simulation run and reported the mean, minimum and maximum rotation period across runs.
Characterization of Vegetation Dynamics
For each cover type that exhibited distinct successional stages and changes (see RMLANDS Model Parameterization), we quantified the following temporal and spatial vegetation dynamics:
• Survivorship - To examine the age structure and dynamics of each cover type, we plotted (as a line chart) the survivorship distribution (i.e., the cumulative proportion of the cover type surviving longer than time t; which represents the probability of surviving without a stand-replacing disturbance longer than time t). The data were pooled across simulation runs. We plotted the mean probability of survivorship and the 5th and 95th percentiles of the distribution. The range of variation in survivorship depicted by the 5-95% envelope reflected the variation over time in the instantaneous age distribution.
• Seral Stage Distribution & Dynamics - To examine the distribution of area among stand conditions (i.e., seral stages) and its changes over time, we created a 100% stacked area chart in which the total area of the corresponding cover type was apportioned among the condition classes (summing to 100% of the cover type area) at each time step. We chose a single representative 800-year simulation run and portrayed the first 100 years in the chart to illustrate the equilibration process. In addition, we summarized the range of variation in the seral stage distribution across all 5 simulation runs as follows. We treated each time step as a single independent observation of the seral stage distribution and computed the proportion of the corresponding cover type in each condition class. Based on the combined data set across runs (350 observations: 5 runs x 70 snapshots at 10-year intervals), we computed the mean, median, minimum, maximum, 5th, 25th, 75th and 95th percentiles of the distribution. We compared the current landscape condition (i.e., proportion of cover type in each condition class) to this simulated historic range of variability to determine whether the current landscape deviates, and to what degree, from the HRV. In some cases due to limitations in the spatial data for the current condition it was necessary to combine condition classes in order to meaningfully compare the simulated HRV with the current landscape condition. Note, we also conducted a much more thorough analysis of HRV departure for each cover type (and for the landscape as a whole) by examining the range of variation in several metrics (computed using FRAGSTATS) describing the extent and spatial configuration of each seral stage (see Methods).