Objective

The objective was to predict potential catastrophic impacts (both fire and flooding) of buffelgrass (Pennisetum ciliare) invasion into Saguaro National Park and surrounding exurban communities. The specific goals are to compare both pre- and post-fire flood effects based on existing vegetation and projected buffelgrass invaded communities into the next century. Based on those results, areas of greatest concern for asset protection and fuels treatment are to be identified. The model outputs are intended to help inform fire management planning, as well as plan for strategic mitigation of flood runoff following buffelgrass-fueled fires in the park. 

Study Area Location

Orthophotography

The normalized difference vegetation index (NDVI) was generated from the 1 m resolution 4-band (near infrared (NIR), red, green, blue) 2010 National Agricultural Imagery Program (NAIP) orthophotography over Saguaro NP East and the surrounding areas (including Coronado National Forest, Arizona State Trust, BLM, and Private Lands). The NDVI is calculated as:

NDVI = Band 4 - Band 1 / Band 4 + Band 1

where Band 4 is the NIR (NIR, 675-850 nm) and Band 1 is the Red (590-675 nm).

The NDVI data layer is projected in NAD83 UTM Zone 12 (ESPG:26912). 

Impervious Surface Characterization

The Rincons are a metamorphic core complex of granite and schist characterized by prominent rocky escarpments which rise out of the mountainside in periodic layers breaking up the vegetation continuity.The study area is composed of both impervious parent bedrock (gneiss, grano-diorite, granite, schist) as well as permiable soils, regolith, and saprolite rock. In order to define what is pervious versus impervious a multi-tier classification was developed from the available orthophotography and LiDAR derived topography (DEM profile curvature).

Thresholded classifiers were defined where impervious (e.g. exposed bed rock) versus pervious (e.g. soil) surfaces are seen to occur. Several common techniques for isolating vegetation were also used: (1) Normalized Difference Vegetation Index (NDVI); (2) LiDAR intensity; (3) LiDAR return number, and (4) LiDAR canopy height model. To characterize the impervious surfaces I created (5) planimetric, profile, tangential, and cross-sectional curvatures from the 3m bare earth DEM with Slope, Aspect, Curvature in SAGA. 

(1) The NAIP / NDVI dataset was partitioned at 1,200 m amsl where vegetation was observed to become more dense, e.g. woodland and forest, versus more open, e.g. Sonoran upland and desert scrub. Above 1,200 m elevation vegetation is a good indicator of soil presence. Exposed bedrock was easily defined using the NDVI classification at higher elevation. Below 1,200 m where vegetation is sparser and was possibly senecent or leaf-off during the NAIP acquisition the NDVI classification is less precise. Those results were based on the outputs of the 2007 and 2010 NAIP derived NDVI to exclude bare rock and shadowed areas. For elevations above 1,200 meters vegetation NDVI was classified between 0.075 and 0.6. Shadows were generally >0.6 NDVI and barren surfaces (stream channels, bare rock) were <0.075. Potential impervious surfaces were initially classified as <0.075 NDVI above 1200 m and <0.05 below 1200 m. The 1m NDVI image was aggregated to coarser 10m and 30m resolution using the r.resamp.stats command from GRASS in QGIS.

(2) As already discussed above, intensity had issues with its normalization (0-255 digital number) for brightness along different flight lines which hindered its usefulness. It was decided not to include the LiDAR intensity as a classification parameter, though it was useful for qualitative comparisons in the areas where the intensity was properly normalized. Future assessments may be able to use intensity as an additional band for classification. 

(3) The LiDAR point cloud can also be used to characterize vegetation presence/absence. The presence of 2nd through 4th returns in the point cloud suggest the pulse traveled downward through vegetation and was split at least one time. If the pulse strikes a planar surface it is more likely to have only a single 1st return. This is not always the case, as pulses striking the edges of rocky faces may have a second or third return. 

(4) Aspect, slope, planimetric, profile, tangential, and cross-sectional curvature classifications were done with Slope, Aspect, Curvature module with the 2nd degree polynomial method (Zevenbergen and Thorne 1987). With r.param.scale from GRASS default settings were used and the model constrained through the test pixel. The 3m and 10m DEMs were tested to determine which scale was most useful for defining potential bare rock. 

Existing Vegetation Type maps

Vegetation has been variously classified by the LANDFIRE project, as well as by the USGS/NPS, and Arizona FireScape projects (USFS).

The LANDFIRE project has the benefit of also supporting fire behavior model parameters. 

Download data from LANDFIRE 

NPS recently developed a vegetation type model from the 2008 and 2011 LiDAR data. Those data cover all of the NPS lands, but do not include the surrounding USFS, State Trust, and private land which supports areas of potential buffelgrass invasion.  

Burn Severity

Conventionally there are two methods for characterizing the 'burn severity' of how a fire has affected landscapes. Th first is the Normalized Burn Ratio (NBR) derived from remote sensing using the Landsat series of satellites. The second is the an extended approach which includes vegetation recovery and delayed mortality involving field sampling - called the Composite Burn Ratio (Key and Benson 2006).

For the dNBR images burned areas were generally considered to be > 0. Key and Benson (2006) suggest that unburnt areas are generally between -100 and 100 dNBR; however we find that in the upper Sonoran desert vegetation types around Saguaro NP a dNBR value of 50 characterized areas that have burned in the past (e.g. 1999 Box Canyon Fire). This is a critical exception because columnar cactus like Saguaros do not survive any fires. 

Difference Normalized Burn Ratio (dNBR)

One of project-goals is to predict burn severity of future buffelgrass wildfires and the impacts on hydrologic function in both the lower desert and the currently forested higher elevation headwaters.

Existing burn severity maps use the so-called difference normalized burn ratio (dNBR)(Key and Benson 2006) derived from Landsat 5, 7, and 8.

The Monitoring Trends in Burn Severity (MTBS) program website hosts the burn severity data. All large fires after 1990 were downloaded for the Saguaro National Park.

The scaled Normalized Burn Ratio is given from Key and Benson (2006) as: 

NBR = 1000(R_4 - R_7)/(R_4+R_7)

where the index is scaled by 1000 so that a ratio of 0.15 is equal to 150, R_4 is the Landsat spectral band 4 (0.76 - 0.90 x 10^-6 nm) and R_7 is the spectral band 7 (2.08-2.35 x 10^-6 nm). The final calculation of dNBR is derived as:

dNBR = NBRprefire - NBRpostfire

where the pre- and post- fire images are differenced. dNBR changes across vegetation types - a low severity fire in forest vegetation has a much higher dNBR than a high severity fire in a desert vegetation type. 

The classified dNBR severity rating changes for each fire location. The class thresholds for four different Rincon fires dNBRs are shown below along with their median values.

Keane et al. (2010) developed the Keane Burn Severity Index.

All of the MTBS fires that have been mapped for the Rincons were combined into a single dNBR layer. These layers were used to generate a predictive surface model for dNBR across all elevations, slopes and aspects based on a 'northness' function of:  

northness = sin(slope) * cos(aspect)

where slope and aspect are in radians.

There are phenological differences in the dNBR for different fires in the Rincons at similar elevations. For example, the Chiva fire burned in July, but was differenced with pre- and post- imagery from September, the Distillery fire burned in June but used pre- September and post- July imagery.

Fire History Reconstruction

The WFDSS website supports a set of GIS layers broken out by decade that have the history of burned area in the Park. 

A combined fire frequency layer was generated in ArcGIS following a set of instructions found here.

The resulting layer shows for the 20th century where recorded fires were mapped in the Rincons.

• 1937-2005 ignition point map from NPS.
• Fire atlases from NPS and WFDSS.
• NPS prescribed fire perimeters.

Lightning

Raw flash lightning location data was obtained by Saguaro National Park Fire Management for the potential fire season (April to August) over the years 2006 to 2013. These data were aggregated into a single data file. A Heat Map was generated in QGIS with three radii: 500m, 1km, and 2km. The frequency of cloud to ground strikes within each radial area are reported in those surfaces. 

There was a trend of higher density cloud to ground lightning strikes at higher elevations on the mountains. The location of the strikes appear to be somewhat shifted off to the north and north-east of the elevations. This may be an artifact of the predominant wind direction from the southwest driving convective cloud formation on the south aspects and then sheering off toward the northeast.

http://www.unidata.ucar.edu/data/lightning/gln.html

Relative frequency of fire size distribution

The fire recurrence interval (FRI) vs area burned was calculated from the 68 year record as:

Modeling the hydraulic response of burned landscapes based on burn severity

Moody (2011) described an analytical framework for predicting burned watershed response to flooding. The Moody model uses a Hydraulic Functional Connectivity function (Moody et al. 2008) where dNBR is weighted differentially along the flow path and has a variable response based on the pattern of burn severity in the upper vs lower hillslope positions (see Moody 2011's Figure 4 below). Essentially there are two outcomes to a wildfire in a burned watershed: (1) low and mixed-severity fire vary in intensity along the hillslope where the predicted dNBR is related to some combination of elevation and slope position - as described above, or (2) a high severity fire where all vegetation is burned and the surface is given the highest severity dNBR possible. In both cases the prediction of dNBRs should be based on the observed dNBR records from MTBS.