METHODS

Computer Modeling

Our HRV analysis employed three different computer models used in conjunction to simulate landscape changes and quantify the dynamics in landscape structure and wildlife habitat under the reference period disturbance regime. Each of these models is described in detail in the sections below.

■ RMLANDS – The landscape simulation model used to simulate changes in vegetation over time under the reference period disturbance regime.

• Model Overview – This section provides a brief overview of RMLANDS sufficient for minimally understanding its use in this application. A more complete, detailed description of RMLANDS is beyond the scope of this report, but may be necessary to fully comprehend the model parameterization.

• Spatial Data – This section provides a brief description of the spatial data base compiled for this application. Each spatial data layer is described in detail, including the data sources and processes used to derive the grid, its use in RMLANDS, and definitions of the grid values.

• Model Parameterization and Calibration – This section provides a detailed description of the model parameterization and calibration for this application. Included here is a list of values for all model parameters and the basis (i.e., justification) for their setting. Empirical data and analyses employed in the parameterization process are presented as well. Note, understanding the relevance of parameter values requires some level of comprehension of the basic structure of RMLANDS (see Model Overview).

• HRV Analyses – This section provides a description of the simulations done for this application and the methods employed to summarize the results of the simulation.

■ FRAGSTATS – The landscape pattern analysis program used to quantify the range of variability in landscape structure over the duration of the simulation.

• Model Overview – This section provides a brief overview of FRAGSTATS sufficient for minimally understanding its application in the HRV analyses. A more complete, detailed description of FRAGSTATS (provided elsewhere) is beyond the scope of this report, but may be necessary to fully comprehend the model application.

• HRV Analyses – This section provides a description of the FRAGSTATS analysis done for this application, including a description of the input data generated from the RMLANDS simulation, a description of landscape pattern metrics used, and a description of the methods employed to summarize changes in landscape pattern metrics over time.

■ HABIT@ – The wildlife habitat capability model used to quantify the range of variability in habitat capability for selected indicator species over the duration of the simulation.

• Model Overview – This section provides a brief overview of HABIT@ sufficient for minimally understanding its application in the HRV analyses. A more complete, detailed description of HABIT@ (provided elsewhere) is beyond the scope of this report, but may be necessary to fully comprehend the model application.

• Species’ Models – This section provides a description of the indicator species selected for this analysis, including the method used to select species, a general description of each species’ relevant life history and habitat associations, and the detailed habitat capability model constructed using HABIT@.

• HRV Analyses – This section provides a description of the HABIT@ analysis done for this application with reference to the species’ models described in the previous section, including a description of the input data generated from the RMLANDS simulation and a description of the methods employed to summarize changes in habitat capability over time.

Evaluation of HRV Departure (Fire Regime Condition Class)

To examine the magnitude of departure of the current landscape condition from the HRV, we modified the approach for Fire Regime Condition Class (FRCC) determination. FRCC is a categorical classification of the degree to which the current fire regime and composition and structure of a vegetation community deviates from its natural or historic range of variability (HRV) under a designated reference period (FRCC website). The FRCC approach is being developed and implemented by an Interagency and TNC (The Nature Conservancy) working group chartered and managed by the Interagency Fuels Committee and as a component under the Rocky Mountain Research Station (RMRS) Fire Effects Unit in association with the Fire Monitoring and Inventory system (FIREMON) and in parallel with the RMRS Rapid Assessment and LANDFIRE projects. FRCC has been an integral component of the development of a cohesive strategy for the restoration of fire-adapted ecosystems and for the National Fire Plan within the USDA Forest Service and USDI land management agencies. FRCC has become a driving force behind current land management activities and is widely being used as the primary (or even sole) basis for identifying and prioritizing areas for ecological restoration, including the reduction of wildfire risks associated with hazardous fuels.

Our approach, described in detail below, modified (and supplemented) the established FRCC approach in several important ways.

(1) Our approach was based on a spatially-explicit model (RMLANDS) of disturbance and succession that provides a much more realistic depiction of how disturbance processes operate in this landscape than is possible with a nonspatial approach. In particular, our model simulates the initiation and spread of disturbance across the landscape in which the spread process interacts with landscape structure (e.g., terrain, vegetation-fuels, wind direction) to create disturbance patterns. The result is that landscape context (i.e., the conditions surrounding a burning cell, for example) has an important influence on disturbance processes and ultimately has a strong influence on vegetation patterns and dynamics. This complexity cannot be captured in a nonspatial model.

(2) Our approach explicitly incorporates multiple disturbance processes. The established approach emphasizes the role of wildfire. In contrast, we model the effects of several dominant disturbance processes including wildfire in addition to several important insect/disease agents that interact to affect vegetation patterns and dynamics. Ultimately, we assess departure as the degree to which current vegetation patterns deviate from the historic range of variation in vegetation patterns, where wildfire is one of several processes driving those changes.

(3) Our approach explicitly incorporates the simulated range of variation in each parameter (see below) into the departure index, instead of calculating departure based on the mean. The established approach does assume variability about the mean, but it assumes that the variance is both symmetrical about the mean and the same for all parameters, an assumption that is highly unlikely to be met in real landscapes. In contrast, we use the actual simulated distribution of values and thus make no assumptions about symmetry and constant variance. Instead, we use percentiles of the simulated distribution to quantify the range of variability in each parameter. The use of percentiles provides a standardize measure of variation that accounts for differences in distribution shape and extent.

(4) Our approach generates a continuously-scaled departure index (ranging from 0-100), instead of a categorical index consisting of three classes. The continuous index preserves the full amount of information and does not require the specification of arbitrary threshold values for class assignment.

(5) Our approach adopts a multivariate perspective, instead of the bivariate one used in the established approach that assigns condition class on the basis of two variables: (1) vegetation-fuels departure, and (2) fire regime (frequency-severity) departure. First of all, we decided not to evaluate fire regime departure due to the difficulties in clearly establishing the current fire regime. Specifically, as there has been only one recent fire in the study area of significant size (i.e., the Missionary Ridge fire, 2002, 20,093 ha) in which fire severity determination was made (a procedure also fraught with numerous methodological problems), we did not feel confident in establishing fire severity levels across cover types. Second, the established approach for evaluating vegetation-fuels departure involves comparing the current seral-stage distribution for each cover type (i.e., percent of cover type in each seral stage) to the mean seral-stage distribution under the natural or historic reference period (obtained via a nonspatial model). Thus, vegetation-fuels departure is based on landscape composition (i.e., the abundance of each seral stage) only, without consideration of landscape configuration (i.e., the spatial pattern, distribution or arrangement of seral stages). Our approach focuses exclusively on vegetation-fuels departure, but considers both composition and configuration (evaluated using several different landscape metrics) as components of the departure index (as described below).

(6) Our approach directly addresses the issues of scale (i.e., landscape extent) and context (i.e., landscape location) in the evaluation of departure. While these issues are considered important under the established approach, and some guidance is offered for selecting the “right” scale, there is nothing in the procedure that allows for direct examination of these effects. Our HRV departure analysis was done at several different spatial scales and for different landscape contexts at each scale, as described in the previous section.

Our approach involved evaluating departure separately for each cover type and for the aggregate landscape as a whole, as follows:

*Cover Type Departure*:

Our procedure for determining departure for each cover type involved the following steps (see this figure for a schematic diagram of the method):

First, we determined the HRV distribution in seral-stage distribution for each cover type. We
used RMLANDS and the simulations described elsewhere to simulate vegetation dynamics under
the historical reference period. For each cover type, we quantified the range of variability in
seral-stage distribution pooled across simulation runs, after excluding the first 100 year
equilibration period of each run (i.e., 5 runs x 70 timesteps = 350 snapshots). The seral-stage
distribution for a single snapshot was represented as the percentage of the cover type in each
seral stage (i.e., stand condition). For each seral stage, we summarized its range of variation by
computing the 0^{th}, 5^{th}, 25^{th}, 50^{th}, 75^{th}, 95^{th} and 100^{th} percentiles of the distribution of observed
values. The 0^{th} percentile represented the lowest observed percentage of the cover type in the
corresponding seral stage. In other words, over the course of the simulations, the relative
abundance of the corresponding seral stage was never less than the 0^{th} percentile level. Similarly,
the 25^{th} percentile represented the percentage in that seral stage below which the landscape fell
25% of the time, and so on for the other percentiles.

Second, we determined the seral-stage distribution for the current landscape. For our
purposes, the “current” condition refers to the landscape in 2003 after the Missionary Ridge fire.
We simply computed the percentage of each cover type in each seral stage directly from the
initial (input) grids (i.e., timestep 0 in the simulation) and then converted these percentages to
percentiles of the HRV distribution by determining what percentage of the HRV distribution was
≥equal to the value for the corresponding seral stage. Thus, we converted the actual percentage in
each seral stage to a percentile of the corresponding HRV distribution. If the actual percentage
value was ≤ the value of the 0^{th} percentile of the HRV distribution, we assigned it the 0^{th}
percentile (since all HRV values were greater than this value). Similarly, if the actual percentage
value was ≥ the value of the 100^{th} percentile of the HRV distribution, we assigned it the 100^{th}
percentile.

Third, we compared the current seral-stage distribution to the simulated HRV distribution to
determine the degree of departure. Specifically, we compared the current condition percentile to
the percentiles of the HRV distribution. It is important to note that we compared *percentiles *not
the actual *percentages *of area in each seral stage. This was done to standardize the measure of
departure as discussed above. We computed the magnitude of departure as follows. For each
cover type and seral stage, if the current condition percentile was between 25-75^{th} percentile, we
treated it as well within the HRV distribution, and it was considered no further (i.e., it did not
contribute to the departure index, or, equivalently, it was given a weight of 0 in the departure
index). If the current condition percentile was <25^{th} percentile, we computed the absolute value
of the difference (i.e., the number of percentiles <25). Similarly, if the current condition
percentile was >75^{th} percentile, we computed the absolute value of the difference (i.e., the
number of percentiles >75). For example, if the current condition was 10 (percentiles), it received
a score of 15 (25-10); if the current condition was 95, it received a score of 20 (95-75). Thus,
each seral stage received a departure score that ranged between 0-25 based on the number of
percentiles it deviated outside the 25-75^{th} percentile range.

Fourth, for each cover type we computed the *mean departure score *across seral stages. The
mean departure score was used to account for the varying number of seral stages among cover
types and to provide an overall measure of departure for each cover type in which equal weight
was given to all seral stages. Thus, each cover type received a departure score between 0-25 that
indicated its average degree of departure from the HRV distribution.

Fifth, we rescaled the mean departure score to range between 0-100 (by dividing by 25, the
maximum score, and multiplying by 100) to create the final *Seral-Stage Departure Index*, where
a 0 indicates no departure from the 25-75^{th} percentile of variation and a 100 indicates complete
departure (i.e., current landscape is outside the HRV). While this rescaling was not necessary, it
allows for adjustments to the 25-75^{th} percentile thresholds (for considering something as
“departure”) without changing the range and interpretation of the final index. Thus, regardless of
whether the 25-75 range or any other user-specified range is considered as “within the HRV”, the
0-100 interpretation of the final index remains the same.

Sixth, we repeated the process above for several class configuration metrics (see
FRAGSTATS HRV Analyses for a complete listing and description of the metrics). Each class
configuration metric represents a different aspect of the spatial character, distribution or
arrangement of patches of the corresponding class, where the class represents a unique stand
condition or seral stage associated with a particular cover type. Briefly, each class-level
configuration metric was evaluated for the degree of departure of the current landscape from the
simulated HRV distribution. Thus, each cover type received a departure index (range 0-100) for
each metric, derived by averaging the scores across the seral stage classes. For each cover type
we computed the mean departure score across metrics to create the final *Class Configuration
Departure Index*, similar to the process for seral-stage departure.

Lastly, we combined the seral-stage and configuration departure indices into an overall *Cover
Type Departure Index *by taking the average of the two component indices. Note, the simple
mean is not the only way to combine the component indices. The final index could just as easily
be computed as from the geometric mean, weighted arithmetic mean, or maximum component
index, depending on the emphasis sought.

*Landscape Structure Departure*:

Our procedure for determining departure for full landscape involved the following steps:

First, we determined the HRV distribution in the areal extent of each unique combination of
cover type and seral-stage (i.e., stand condition) based on the RMLANDS simulations described
above. Each unique cover type and seral stage was treated as a unique patch type (or class) in this
analysis, and the mosaic of resulting patches comprised the landscape mosaic in each snapshot
(timestep) of the simulation. For each class (patch type), we summarized its range of variation by
computing the 0^{th}, 5^{th}, 25^{th}, 50^{th}, 75^{th}, 95^{th} and 100^{th} percentiles of the distribution of observed
values in the same manner as described above.

Second, we determined the areal extent of patch type for the current landscape. We simply computed the percentage of each patch type directly from the initial (input) grids (i.e., timestep 0 in the simulation) and then converted these percentages to percentiles of the HRV distribution by determining what percentage of the HRV distribution was ≤ to the value for the corresponding patch type. Thus, we converted the actual percentage in each patch type to a percentile of the corresponding HRV distribution in the same manner as described above.

Third, we compared the current landscape composition to the simulated HRV distribution to determine the degree of departure. Specifically, we compared the current condition percentile to the percentiles of the HRV distribution for each patch type in the same manner as described above.

Fourth, we computed the *mean departure score *across patch types. The resulting score
ranged between 0-25 and indicated the average degree of departure from the HRV distribution.

Fifth, we rescaled the mean departure score to range between 0-100 (by dividing by 25, the
maximum score, and multiplying by 100) to create the final *Landscape Composition Departure
Index*, where a 0 indicates no departure from the 25-75^{th} percentile of variation and a 100
indicates complete departure (i.e., current landscape is outside the HRV).

Sixth, we repeated the process above for several landscape configuration metrics (see
FRAGSTATS HRV Analyses for a complete listing and description of the metrics). Briefly, each
landscape metric was evaluated for the degree of departure of the current landscape from the
simulated HRV distribution. Thus, the landscape received a departure index (range 0-100) for
each metric. We computed the mean departure score across metrics to create the final *Landscape
Configuration Departure Index*, similar to the process for seral-stage departure.

Lastly, we combined the landscape composition and configuration departure indices into an
overall *Landscape Structure Departure Index *by taking the average of the two component
indices. Note, the simple mean is not the only way to combine the component indices. The final
index could just as easily be computed as from the geometric mean, weighted arithmetic mean,
or maximum component index, depending on the emphasis sought.

Evaluation of Scale and Context Effects

To examine the effects of scale - specifically, landscape extent - and landscape context (i.e., the unique landscape structure associated with each geographic location) on the measured dynamics in landscape structure and wildlife habitat capability (for the selected indicator species), we subdivided the full landscape into three districts (corresponding to the existing Ranger Districts) of similar size and further selected a single watershed of comparable size nested within each District (Figure-map), as follows:

Level |
Landscape |
Total Area (ha) |

1 |
San Juan National Forest |
847,638 |

2 |
Columbine District |
308,829 |

2 |
Pagosa District |
282,511 |

2 |
Dolores District |
256,298 |

3 |
Hermosa Creek Watershed |
44,103 |

3 |
Piedra River Watershed |
39,487 |

3 |
Narraguinnep Creek Watershed |
31,818 |

Note, RMLANDS was run on the full SJNF landscape only. The output landscapes were clipped to each of these sub-landscapes for purposes of landscape pattern analysis (FRAGSTATS) and wildlife habitat capability analysis (HABIT@).

To assess the effects of scale (landscape extent), we compared the range of variability in
landscape structure (using FRAGSTATS) and wildlife habitat capability (using HABIT@)
among the three scales. Our goal was not to examine a comprehensive range of scales, but rather
to examine potentially realistic scales for planning and/or analysis. With regards to landscape
structure, we combined data across simulation runs (i.e., ignored the separate runs) for a total of
350 observations (5 runs x 70 snapshots at 10-year intervals) at each scale. We computed the
90% coefficient of variation in each metric for each landscape. We defined the 90% coefficient
of variation as the difference between the 5^{th} and 95^{th} percentiles of the metric distribution,
divided by the median (50^{th} percentile) and multiplied by 100 to convert to a percentage. This
measure is similar to a conventional coefficient of variation, but instead of using the standard
deviation (which assumes a normal distribution) as a measure of variance we used the 5^{th} to 95^{th}
percentile range (which involves no assumptions about the distribution), and instead of using the
mean (which is most appropriate for a symmetrical distribution) to represent the central tendency
we used the median (which is more appropriate for non-normal distributions). This provided a
relative measure of the range of variation in each landscape metric that was comparable across
metrics and landscape extents. With regards to wildlife habitat capability, we were only able to
complete the HABIT@ analysis for single representative run due to the intensive computation
requirements. Moreover, we were only able to complete the analysis at the scale of the SJNF for
a single indicator species - three-toed woodpecker. Otherwise, all four indicator species were
analyzed for the Columbine District and Hermosa Creek Watershed. As with the landscape
metrics, we computed the 90% coefficient of variation in each species’ Landscape Capability
(LC) index for each landscape. This provided a relative measure of the range of variation in each
species’ habitat that was comparable across species and landscape extents.

Similarly, to assess the effects of landscape context, we compared the range of variability in landscape structure among landscapes at the finest and intermediate scales (i.e., watershed and districts scales). We were unable to assess the relationship between landscape context and wildlife habitat capability as we were only able to complete the HABIT@ analyses for a single extent (i.e., sublandscape) at each scale. Our goal was not to examine every possible context, but rather to examine a range of contexts in order to demonstrate the importance of context (i.e., that no two landscapes are the same). As before, we combined data across simulation runs (i.e., ignored separate runs) for a total of 350 observations (5 runs x 70 snapshots at 10-year intervals) for each district and watershed. We compared the percent similarity in the distribution of each metric among landscapes at each extent. Specifically, we computed the percent similarity as the range of values common to both landscapes divided by the total range of values across both landscapes, multiplied by 100 to convert to a percentage. For example, for a particular metric, if landscape A had a range of 1-10 and landscape B had a range of 2-15, the percent overlap would be 8 (range 2-10 in common) divided by 15 (range 1-15) times 100, or 53%. For each metric, we computed the percent similarity between each pairwise combination of landscapes at each scale and then computed the average percent similarity across comparisons. We interpreted 100% similarity as evidence that the different landscapes experienced identical ranges of variation and, therefore, that landscape context was not important. We interpreted values <100% as evidence that landscapes behaved differently and therefore that landscape context was important. Zero percent similarity indicated that the landscapes had completely different ranges of variation in the corresponding metric.