An introduction to firebehavioR

Justin Ziegler

; Forest & Rangeland Stewardship, Colorado State University, Fort Collins, CO, USA

2019-02-28

1.0 About firebehavioR

This vignette is a brief tutorial to firebehavioR. This package is a compilation of different mathematical models of fire behavior or fire hazard indices. The former make explicit, measurable predictions of fire behavior while the later describe the relative hazard posed by current weather and/or fuel conditions. In both cases, these models express the potential for fire behavior posed by fuels, weather and/or topography (the three sides of the wildland fire behavior triangle). This is not a package for statistical methods in fire research. Most of these models have been incorporated, in some form, within other software and often in a more user-friendly graphical interface. My rationale for developing this package was twofold: 1) to take advantage of the workflow efficiency in R; for example, users can wrangle data, directly perform complex statistical analyses on outputs, and produce visualizations in one environment; and 2) users can view the source code; the transparency of fire modelling benefits researchers and educators. However, this tutorial will not focus on either of these as needs can vary greatly across potential users. Each model has its own sets of limitations and assumptions for which the user is responsible for understanding.

To provide benefit to the most users this tutorial assumes that the user has only a novice level understanding of the R programming language as well as a recent version of R installed (>=3.4). Users are encouraged to play around with the code and firebehavioR’s functions on their own to find how this package can benefit their needs. This tutorial also assumes a basic understanding of variables used in fire modelling. Appendix A contains a listing of all the variables used in this package; it is critical that users understand the measurement units used in this package.

1.1 Citing firebehavioR

Maintenance of firebehavioR is done on my own time and I do it to benefit the fire research community. That said, please cite firebehavioR appropriately if you use it. Not only does it give me credit for developing this package, but it also signals to me that this package is useful and deserves regular maintenance and revision.

To cite this package, there is an upcoming technical note that I am preparing for a scientific journal. When this is published, that will be the preferred reference. In the interim you can cite the development version:

Ziegler, J. P. 2018. firebehavioR. https://github.com/EcoFire/firebehavioR.

1.2 Installing firebehavioR

The firebehavioR package is not yet on the CRAN repository.

The development version is available at https://github.com/EcoFire/firebehavioR. You can download this from with R with the following:

Once firebehavioR is installed, you can load it just as any other package:

2.0 Fire models

2.1 Crown Fire Initiation & Spread model

Crown Fire Initiation & Spread (CFIS) predicts crown fire behavior. The arguments include the fuel stratum gap, 10-m open wind speed, fine fuel moisture, the load of surface fuel consumed, the canopy bulk density and spot ignition delay. To demonstrate this function, we will use data(coForest). This data frame contains summary information on canopy and surface fuels taken from different forest stands in the Southern Rocky Mountains, USA. As with all objects, more information can found by entering in ?coForest.

In this example, we will enter in the measured fuel stratum gap and canopy bulk density while assuming that the surface fuel consumption is 100% of the measured surface fuel load. We also assume the 10-m open wind speed is 20 km/hr and the effective fine fuel moisture is 6%. Note that sfl_kgm2 are surface fuel loads in units of kg/m2 and must be converted to units of Mg/ha.

The output returns the type of fire (surface, passive crown or active crown), probability of crown fire, crown fire rate of spread, and separation distance (described later).

coForest contains information on seven forest stands sampled before and after silvicultural treatment. Let"s explore how treatments affect average crown fire rate of spread ceteris paribus.

The results state that mean crown fire rate of spread is lower post-treatment than pre-treatment.

In this next example, we will explore the interaction between effect of treatment and open wind speed on fire type. To do this we need to replicate our 7 forest stands, each with an increasing level of open wind speed.

The results show that as wind speeds increase, pre-treatment stands move from a mix of active and surface fire types to a mainly active fire type. In contrast, post-treatment stands move from a mostly surface fire type to a mostly passive fire as wind speeds increase. A simple interpretation on treatment effectiveness would state that treatments reduce the propensity for the most severe class of crown fires (active), regardless of wind speed, but do not preclude passive crown fire.

2.1. Spot fire calculation

The id argument for cfis() does not effect fire type, probability of crown fire occurrence or crown fire rate of spread. This argument is used to predict the critical separation distance for a spot fire.
Spot fires result from the generation of embers from a source fire, transportation aloft, and ignition at a target location ahead of the fire. During this time, the source fire is advancing according to its rate of spread. Therefore, the time delay between generation and ignition determines how far an ember must travel in order to avoid landing in the advancing source fire. In this example, we will look at the relationship between ignition delay and critical travel distance.

Requisite conditions for successful spot fire ignition given the projected crown fire rate of spread.

Requisite conditions for successful spot fire ignition given the projected crown fire rate of spread.

Given the predicted crown fire rate of spread, the above figure demonstrates the set of ignition delays and separation distances under which an ember could ignite a spot fire in advance of a source fire.

2.2 Rothermel

The Rothermel modelling system, invoked with rothermel(), predicts fire behavior using linked models of surface fire spread, crown fire initiation, and crown fire spread. rothermel() takes four arguments input as data frames. In each of these data frames, the column order is important.

  1. Surface fuel characteristics
    • Fuel model type, static or dynamic
    • Litter load
    • 1-hr load
    • 10-hr load
    • 100-hr load
    • Herbaceous load
    • Shrub load
    • 1-hr SAV
    • 10-hr SAV
    • 100-hr SAV
    • Herbaceous SAV
    • Shrub SAV
    • Fuel bed depth
    • Moisture of extinction
    • Heat content
  2. Surface fuel moisture, for:
    • Litter
    • 1-hr
    • 10-hr
    • 100-hr
    • Herbaceous
    • Shrub
  3. Crown fuel characteristics
    • Canopy bulk density
    • Foliar moisture content
    • Canopy base height
    • Canopy fuel load
  4. Environmental characteristics
    • Topographic slope
    • 10-m Open windspeed
    • Wind direction, from uphill
    • Wind adjustment factor

Additional arguments include a crown fire rate of spread multiplication factor, the method for determining crown fraction burned, and whether or not a foliar moisture effect is calculated.

This example will use a set of data objects to populate the input arguments. data(fuelModels) is a data frame with rows of stylized surface fuel models. We will use Fuel Model 10, represented timber litter. The next argument is a data frame of fuel moisture scenarios, data(fuelMoisture), of which we will choose D1L1. Again, we will use data(coForest) for some canopy characteristics. Next, we will assign hypothetical values to environmental variables and leave rosMult, cfbForm and folMoist as defaults. Note that input data frames are repeats since rothermel() predicts fire behavior for each entered row of input argument data frames and coForest contains 14 different forest stands.

rothermel() outputs a list of data frames. The first describes the finalized outputs.

Intermediate outputs used for calculating surface and crown fire behavior are also in the output.

head(ex5$detailSurface)
#>   Potential ROS [m/min] No wind, no slope ROS [m/min] Slope factor [0-1]
#> 1                  4.64                          0.43               0.18
#> 2                  4.64                          0.43               0.18
#> 3                  4.64                          0.43               0.18
#> 4                  4.64                          0.43               0.18
#> 5                  4.64                          0.43               0.18
#> 6                  4.64                          0.43               0.18
#>   Wind factor [0-1] Characteristic dead fuel moisture [%]
#> 1              9.68                                  3.08
#> 2              9.68                                  3.08
#> 3              9.68                                  3.08
#> 4              9.68                                  3.08
#> 5              9.68                                  3.08
#> 6              9.68                                  3.08
#>   Characteristic live fuel moisture [%] Characteristic SAV [m2/m3]
#> 1                                 10.33                    5789.31
#> 2                                 10.33                    5789.31
#> 3                                 10.33                    5789.31
#> 4                                 10.33                    5789.31
#> 5                                 10.33                    5789.31
#> 6                                 10.33                    5789.31
#>   Bulk density [kg/m3] Packing ratio [-] Relative packing ratio [-]
#> 1                 8.83              0.02                       2.35
#> 2                 8.83              0.02                       2.35
#> 3                 8.83              0.02                       2.35
#> 4                 8.83              0.02                       2.35
#> 5                 8.83              0.02                       2.35
#> 6                 8.83              0.02                       2.35
#>   Reaction intensity [kW/m2] Heat source [kW/m2] Heat sink [kJ/m3]
#> 1                    1330.81              697.75           9036.08
#> 2                    1330.81              697.75           9036.08
#> 3                    1330.81              697.75           9036.08
#> 4                    1330.81              697.75           9036.08
#> 5                    1330.81              697.75           9036.08
#> 6                    1330.81              697.75           9036.08
head(ex5$detailCrown)
#>   Potential ROS [m/min] No wind, no slope ROS [m/min] Slope factor [0-1]
#> 1                 40.97                           1.5               0.18
#> 2                 40.97                           1.5               0.18
#> 3                 40.97                           1.5               0.18
#> 4                 40.97                           1.5               0.18
#> 5                 40.97                           1.5               0.18
#> 6                 40.97                           1.5               0.18
#>   Wind factor [0-1] Characteristic dead fuel moisture [%]
#> 1             26.08                                  3.08
#> 2             26.08                                  3.08
#> 3             26.08                                  3.08
#> 4             26.08                                  3.08
#> 5             26.08                                  3.08
#> 6             26.08                                  3.08
#>   Characteristic live fuel moisture [%] Characteristic SAV [m2/m3]
#> 1                                  5.17                    5788.49
#> 2                                  5.17                    5788.49
#> 3                                  5.17                    5788.49
#> 4                                  5.17                    5788.49
#> 5                                  5.17                    5788.49
#> 6                                  5.17                    5788.49
#>   Bulk density [kg/m3] Packing ratio [-] Relative packing ratio [-]
#> 1                 8.84              0.02                       2.35
#> 2                 8.84              0.02                       2.35
#> 3                 8.84              0.02                       2.35
#> 4                 8.84              0.02                       2.35
#> 5                 8.84              0.02                       2.35
#> 6                 8.84              0.02                       2.35
#>   Reaction intensity [kW/m2] Heat source [kW/m2] Heat sink [kJ/m3]
#> 1                    1399.37             1843.02           9050.41
#> 2                    1399.37             1843.02           9050.41
#> 3                    1399.37             1843.02           9050.41
#> 4                    1399.37             1843.02           9050.41
#> 5                    1399.37             1843.02           9050.41
#> 6                    1399.37             1843.02           9050.41

Some users may also find critical values useful. This includes critical values for the initiation of crown fire, the transition from passive to active crown fires, and the cessation of active crown fire.

Users can extend the functionality of the Rothermel fire modelling system. Since crown fuels measurements come from sampling, we could ask what the consequences of uncertainty are on potential fire rate of spread.

Effect of uncetain canopy fuels on predicted rate of spread.

Effect of uncetain canopy fuels on predicted rate of spread.

In another example, we explore how rate of spread for a grass fire varies with wind direction and topographic slope.

Distribution of hourly fuel moistures predicted by three methods using `ffm()`.

Distribution of hourly fuel moistures predicted by three methods using ffm().

We see here that the influence of wind direction increases with slope. Clocking from 0 to 90 degrees from the upslope direction, wind direction countervails the effect of slope.

3.0 Helper Functions

The fire models require an array of inputs describing the fire behavior triangle. Some of these inputs are not observable using direct measurement techniques and these must be indirectly estimated. While the user ultimately decides the estimation method, firebehavioR contains some useful functions.

3.1 Estimating canopy characteristics

canFuel() estimates base height, fuel load, and bulk density of the canopy using often measured forest inventory variables. This example will estimate the canopy characteristics for stands of different forest cover types and otherwise identical canopy structures. The required inputs are basal area, canopy height, trees per hectare and forest cover type. These outputs can be used as inputs for either cfis() or rothermel().

3.2 Estimating fine fuel moisture

There are many ways to use meteorological observations to indirectly estimate surface fuel moisture. These all estimate 1-hr time-lag class woody surface fuels (0-0.635 mm particle diameter; i.e. fallen twigs). This functions assumes that larger diameter surface fuel components can be simply estimated by adding additional moisture content.

In the first example, we will estimate 1-hr fuel moisture content using three methods, "simard", "wagner", and "anderson", in ffm. The input data we will use contains hourly weather observations from a Remote Automated Weather Station during one calendar year"s fire season (data(rrRAWS)).

Distribution of hourly fuel moistures predicted by three methods using `ffm()`.

Distribution of hourly fuel moistures predicted by three methods using ffm().

The resulting plot shows the distribution of hourly, 1-hr fine fuel moisture content by method. Methods yield slightly differing distributions, particularly with regard to centrality and right-tail skewness.

Now we will add on a fourth method and view the results as a smoothed time series. When method="mcarthur", it is advisable to estimate moisture content only under a set of weather conditions: \[ \begin{eqnarray} & rh \leqq70\\ & temp \geqq10\\ & 42.5-1.25temp < rh < 94.5 - 1.35temp \end{eqnarray} \]

Hourly fuel moistures predicted by four methods using `ffm()`

Hourly fuel moistures predicted by four methods using ffm()

While the prior methods require just temperature and relative humidity, the Fire Behavior Officer“s Table (fbo) method requires also month, hour, topographic aspect, percent slope, and shading (whether the fine fuels of interest are shaded by cloud cover or forest canopy). Further, there is a modification to state whether weather observations were taken higher than, lower than, or within 305 m +/- of the fuel”s elevation. Here we will assume that the aspect is north, the slope is 10%, the fuel is near the same elevation as the RAWS station and the fuels are unshaded.

Hourly fuel moistures predicted by `method = 'fbo'`, in blue, relative to other methods, in grayscale.

Hourly fuel moistures predicted by method = 'fbo', in blue, relative to other methods, in grayscale.

The figure here illustrates that the fbo method, in blue, follows similar trends as the simard method, except that peaks and troughs are relatively dampened.

3.3 Estimating wind adjustment factor

The rothermel() requires the user to enter in a wind adjustment factor (WAF; ratio of mid-flame wind speed to 20-ft open wind speed). waf() contains separate equations for unsheltered (i.e. no forest canopy) fuelbeds and sheltered fuelbeds. If the fuelbed is sheltered, the forest canopy height, crown ratio (average tree crown length [tree height minus crown base height] divided by tree height), and canopy cover are needed:

The crown ratio and canopy cover are used internally to derive a quantity known as crown fill, \(CrownFill = CrownRatio * CanopyCover/300\). If crown ratio and canopy cover are not known, but the sheltered WAF equation is desired, the waf() assumes \(CrownFill = 0.1\).

If the fuelbed is unsheltered, the only required argument is the surface fuel depth:

It is worth noting that if crown fill is less than 0.05, the unsheltered equation will be used selected For example:

3.4 Visualizing fire behavior ourputs

To make quick assessments of fire behavior, some fire behavior analysts refer to a fire characteristics chart. These charts plot heat per unit area against rate of spread. To aid interpretation, these charts also contain predicted flame lengths and images that refer to the recommended method of direct fire suppression. To expedite fire behavior assessments with fire characteristics charts, fireChart() instructs ggplot with a pre-populated set of arguments and uses the following arguments from the user: the name of each observation, observations" heat per unit area and observations" rate of spread.

Fire characteristics chart

Fire characteristics chart

4.0 Fire weather indices

As opposed to the fire behavior models described Section 2.0, fire weather indices describe the relative fire hazard posed by environmental conditions; most indices ignore the topography leg of the fire behavior triangle and, at most, minimally account for the fuels. This package divides indices into two categories, “static indices” which can use arguments with a minimum length of 1, and “cumulative indices”" that require continuous data with multiple observations. To elaborate, if one used two days of weather observations to derive a fire hazard index, a static index“s calculation would not require the prior day”s index value whereas a cumulative index would.

###Static indices In this package, instantaneous indices are found in fireIndex(). This function uses vectors of air temperature, wind speed, relative humidity, grass fuel load (optional) and grass curing level (optional) to produce the Angstrom Index, Hot, Dry, Windy Index, the Fuel Moisture Index, the Fosberg Fire Weather Index, the MacArthur Grassland Mark 4 Index, the MacArthur Grassland Mark 5 Index, and the Chandler Burning Index. We will use the data(rrRAWS) to illustrate the function, calculating the daily indices at 14:30 MDT. Then we will normalize the scales to compare trends across indices.

rrRAWS.daily =   rrRAWS[format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H:%M")=="14:35",]
indices = fireIndex(temp=rrRAWS.daily$temp_c, u= rrRAWS.daily$windSpeed_kmh, rh = rrRAWS.daily$rh)
normalize=function(x,low=0,high=1){low+(x-min(x,na.rm=T))*(high-low)/(max(x,na.rm=T)-min(x,na.rm=T))}
indices = data.frame(sapply(indices,normalize),  Date = strptime(rrRAWS.daily$dateTime, "%m/%d/%Y %H:%M"))
indices = setNames(reshape2::melt(indices,id="Date"),c("Date","Index","Value"))
ggplot(indices,aes(Date,Value,group=Index,color=Index))+geom_smooth(span = .1,  method = "loess", se = F) +theme_classic()+coord_cartesian(expand=F)
Continuous fire index values applied to `data(rrRAWS)`.

Continuous fire index values applied to data(rrRAWS).

###Cumulative indices Cumulative fire weather indices explicitly use the prior day"s hazard in order to account for seasonal trends. The indices covered by the fireIndexKBDI() include the Keetch-Bryam Drought Index (KBDI), the Drought Factor, MacArthur Forest Mark 5 Index, KBDI-modified Fosberg Fire Index, KBDI-modified Fuel Moisture Index, the original and modified Nesterov Indices, and the Zdenko Index. This example will also use data(rrRAWS), which has a mean annual precipitation of 610 mm. Again, we will normalize indices for comparison.

data(rrRAWS)
daily.precip = rrRAWS
daily.precip$Date = strptime(daily.precip$dateTime, "%m/%d/%Y")
daily.precip = setNames(aggregate(daily.precip$precip_mm, by = list(as.character(daily.precip$Date)), 
    FUN = sum), c("Date", "DailyPrecip"))
rrRAWS.daily = rrRAWS[format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H:%M") == 
    "14:35", ]
rrRAWS.daily$DailyPrecip = daily.precip$DailyPrecip
indices = fireIndexKBDI(temp = rrRAWS.daily$temp_c, precip = rrRAWS.daily$DailyPrecip, 
    map = 610, rh = rrRAWS.daily$rh, u = rrRAWS.daily$windSpeed_kmh)
indices = data.frame(sapply(indices,normalize), Date = strptime(rrRAWS.daily$dateTime, 
    "%m/%d/%Y %H:%M"))
indices = setNames(reshape2::melt(indices, id = "Date"), c("Date", "Index", "Value"))
ggplot(indices, aes(Date, Value, group = Index, color = Index)) + geom_smooth(span = 0.1, 
    method = "loess", se = F) + theme_classic() + coord_cartesian(expand = F)
Static fire index values applied to `data(rrRAWS)`.

Static fire index values applied to data(rrRAWS).

There are three time periods during which most indices agree on high fire hazard, later April, late June, and early-to-mid September.Because the cumulative indices track accumulating droughtiness and incorporate precipitation events, these indices account for the late Spring showers and Summer monsoon rains typical of this region.

The estimation methods in fireIndexKBDI() vary by their requisite arguments. Below is a chart that maps needed arguments to each Index method. fireIndexKBDI() will output results only for those methods that have been supplied arguments.
Method Inputs
1 KBDI temp, precip, map
2 Drought factor temp, precip, map
3 Forest Mark 5 temp, precip, map, u, rh
4 Fosberg-KBDI temp, precip, map, u, rh
5 Fuel moisture-KBDI temp, precip, map, u, rh
6 Nesterov temp, precip, rh
7 Nesterov-Modified temp, precip, rh
8 Zdenko temp, precip, rh

##Conclusion This vignette covers the current state of firebehavioR. This package is under active development and may change in subsequent years. Some changes may include incorporation of meteorological data to predict dynamic fuel moisture, such as adsorption, desorption and wetting events, prediction of live fuel moisture, incorporation of the National Fire Danger Rating System, and linkages with other packages to enhance analyses.

Appendix A: Input quantities

The following quantities are used as input variables for this package“s functions. Use”?" to see documentation of individual functions.
Variable Units
1 Fuel stratum gap m
2 10-m Open wind speed km/hr
3 Fuel moisture content %
4 Canopy bulk density kg/m3
5 Seperation distance m
6 Basal area m2/ha
7 Average stand tree heights m
8 Trees per hectare trees/ha
9 Relative humidity %
10 Temperature degrees C
11 Month Month of the year (1-12)
12 Hour Hour of day (1-24)
13 Topographic aspect N,S,W, or E
14 Topographic slope %
15 Heat per unit area kJ/m2
16 Fire rate of spread m/min
17 surface fuel load Mg/ha
18 Surface area to volume m2/m3
19 Fuelbed depth cm
20 Moisture of extinction %
21 Heat content J/g
22 Canopy fuel load kg/m2
23 Wind direction 0-360
24 Wind adjustment factor Ratio of 20-ft open wind speed to midflame wind speed
25 Crown ratio %
26 Canopy cover %