# An introduction to firebehavioR

## 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:

devtools::install_github("EcoFire/firebehavioR")

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

library("firebehavioR")

## 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.

data("coForest")
print(coForest)
#>      site status trees_perha basalArea_m2ha qmd_cm height_m sfl_kgm2
#> 1    Heil    pre         418           19.1   24.1      9.8     1.04
#> 2    Heil   post         295           14.5   25.0     10.1     0.76
#> 3     Unc    pre         504           30.0   27.5     22.1     1.20
#> 4     Unc   post         314           17.6   16.7     25.5     1.06
#> 5    Pike    pre         934           24.8   18.4     15.8     1.33
#> 6    Pike   post         352            9.5   18.5     16.6     1.33
#> 7    RedF    pre         516           14.1   18.7     14.2     1.16
#> 8    RedF   post         181            9.0   19.0     15.0     0.61
#> 9  Kaibab    pre         487           25.6   25.9     20.9     0.36
#> 10 Kaibab   post         203           19.7   35.2     25.3     0.29
#> 11  MessG    pre         685           22.7   20.5     14.9     0.41
#> 12  MessG   post         269           10.6   22.4     16.6     0.62
#> 13   Zuni    pre         327           19.3   27.4     14.3     0.24
#> 14   Zuni   post          63            7.5   38.9     18.2     0.24
#>    cbd_kgm3 cbh_m cfl_kgm2
#> 1      0.14  3.63    0.837
#> 2      0.10  3.69    0.628
#> 3      0.15  4.35    2.653
#> 4      0.08  4.14    1.674
#> 5      0.15  2.90    2.096
#> 6      0.06  3.60    0.830
#> 7      0.08  2.60    0.918
#> 8      0.05  3.40    0.610
#> 9      0.13  4.90    2.235
#> 10     0.10  6.70    1.811
#> 11     0.14  3.50    1.619
#> 12     0.06  3.60    0.810
#> 13     0.09  3.83    0.914
#> 14     0.03  4.60    0.396

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.

ex1 = cfis(fsg = coForest$cbh_m, u10 = 20, effm = 6, sfc = coForest$sfl_kgm2 * 10,
cbd = coForest$cbd_kgm3, id = 1) print(ex1) #> type pCrown cROS sepDist #> 1 active 99.30 40.54 381.88 #> 2 active 89.63 38.03 358.23 #> 3 active 98.84 41.08 386.92 #> 4 passive 99.00 13.79 129.91 #> 5 active 99.58 41.08 386.92 #> 6 passive 99.32 17.31 163.03 #> 7 passive 99.66 13.79 129.91 #> 8 passive 91.40 19.13 180.17 #> 9 active 78.55 39.97 376.54 #> 10 active 50.50 38.03 358.23 #> 11 active 90.82 40.54 381.88 #> 12 passive 90.21 17.31 163.03 #> 13 active 88.67 37.28 351.13 #> 14 passive 81.92 22.36 210.59 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. ex2 = aggregate(x = ex1$cROS, by = list(treatmentStatus = coForest$status), FUN = mean) print(ex2) #> treatmentStatus x #> 1 post 23.70857 #> 2 pre 36.32571 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. ex3 = coForest[rep(seq_len(nrow(coForest)), 11), ] ex3$u10 = sort(rep(10:20, 14))
#>   site status trees_perha basalArea_m2ha qmd_cm height_m sfl_kgm2 cbd_kgm3
#> 1 Heil    pre         418           19.1   24.1      9.8     1.04     0.14
#> 2 Heil   post         295           14.5   25.0     10.1     0.76     0.10
#> 3  Unc    pre         504           30.0   27.5     22.1     1.20     0.15
#> 4  Unc   post         314           17.6   16.7     25.5     1.06     0.08
#> 5 Pike    pre         934           24.8   18.4     15.8     1.33     0.15
#> 6 Pike   post         352            9.5   18.5     16.6     1.33     0.06
#>   cbh_m cfl_kgm2 u10
#> 1  3.63    0.837  10
#> 2  3.69    0.628  10
#> 3  4.35    2.653  10
#> 4  4.14    1.674  10
#> 5  2.90    2.096  10
#> 6  3.60    0.830  10

ex3$type = cfis(ex3$cbh_m, ex3$u10, 6, ex3$sfl_kgm2 * 10, ex3$cbd_kgm3, 1)$type
ex4$id = 0:10 library("ggplot2") ggplot(ex4, aes(sepDist, id)) + geom_ribbon(aes(ymin = pmin(sepDist, id), ymax = 0), fill = "red", col = "red", alpha = 0.5) + geom_line(size = 1.5) + labs(x = "Critical seperation distance (m)", y = "Ignition delay (min)") + annotate("text", x = 450, y = 3, label = "Spot fire ignition zone") + theme_classic() 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. data(fuelModels, fuelMoisture) fuelModels['A10',] #> fuelModelType loadLitter load1hr load10hr load100hr loadLiveHerb #> A10 S 0 6.73 4.48 11.23 0 #> loadLiveWoody savLitter sav1hr sav10hr sav100hr savLiveHerb #> A10 4.48 0 6561.68 358 98 0 #> savLiveWoody fuelBedDepth mxDead heat description #> A10 4921.26 30.48 25 18622 Timber (litter and understory) #> source #> A10 Anderson (1982) exampSurfFuel = fuelModels['A10',] fuelMoisture['D1L1',] #> fmLitter fm1hr fm10hr fm100hr fmLiveHerb fmLiveWoody #> D1L1 3 3 4 5 30 60 #> description #> D1L1 Very dry dead FM, fully cured herb exampFuelMoisture = fuelMoisture['D1L1',] exampCrownFuel = data.frame( CBD = coForest$cbd_kgm3,
FMC = 100,
CBH = coForest$cbh_m, CFL = coForest$cfl_kgm2
)

exampEnviro = data.frame(
slope = 10,
windspeed = 40,
direction = 0,
waf = 0.2
)
exampSurfFuel = sapply(exampSurfFuel, rep, 14)
exampFuelMoisture = sapply(exampFuelMoisture, rep, 14)
exampEnviro = sapply(exampEnviro,rep,14)
ex5 = rothermel(exampSurfFuel, exampFuelMoisture, exampCrownFuel, exampEnviro)

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

head(ex5$fireBehavior) #> Type of Fire Crown Fraction Burned [%] Rate of Spread [m/min] #> 1 surface 8.62 4.64 #> 2 surface 4.95 4.64 #> 3 conditional 0.00 4.64 #> 4 conditional 0.00 4.64 #> 5 passive 23.16 4.64 #> 6 surface 3.62 4.64 #> Heat per Unit Area [kJ/m2] Fireline Intensity [kW/m] Flame Length [m] #> 1 18733.15 1446.50 2.31 #> 2 17967.30 1387.37 2.22 #> 3 17388.16 1342.65 2.13 #> 4 17388.16 1342.65 2.13 #> 5 26435.24 2041.23 2.98 #> 6 17948.02 1385.88 2.20 #> Direction of max spread [deg] Scorch height [m] Torching Index [m/min] #> 1 0 15.66 35.76 #> 2 0 15.16 36.47 #> 3 0 14.77 44.26 #> 4 0 14.77 41.78 #> 5 0 20.42 27.12 #> 6 0 15.15 35.41 #> Crowning Index [km/hr] Surfacing Index [km/hr] #> 1 24.70 24.70 #> 2 31.81 31.81 #> 3 23.43 23.43 #> 4 37.51 37.51 #> 5 23.43 23.43 #> 6 46.26 46.26 #> Effective Midflame Wind [km/hr] Flame Residence Time [min] #> 1 8.79 0.22 #> 2 8.50 0.22 #> 3 8.11 0.22 #> 4 8.11 0.22 #> 5 9.95 0.22 #> 6 8.39 0.22 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. head(ex5$critInit)
#>   Fireline Intensity [kW/m] Flame length (m) Surface ROS [m/min]
#> 1                   1164.96             1.99                4.02
#> 2                   1193.96             2.02                4.12
#> 3                   1528.21             2.26                5.28
#> 4                   1418.90             2.18                4.90
#> 5                    831.85             1.71                2.87
#> 6                   1150.54             1.98                3.97
#>   Canopy base height [m]
#> 1                   3.99
#> 2                   3.99
#> 3                   3.99
#> 4                   3.99
#> 5                   3.99
#> 6                   3.99
head(ex5$critActive) #> Canopy bulk density [kg/m3] ROS, crown (R'active) [m/min] #> 1 0.07 21.43 #> 2 0.07 30.00 #> 3 0.07 20.00 #> 4 0.07 37.50 #> 5 0.07 20.00 #> 6 0.07 50.00 head(ex5$critCess)
#>   Canopy base height [m] O'cessation [km/hr]
#> 1                  17.06                5.43
#> 2                  17.06                5.60
#> 3                  17.06                7.40
#> 4                  17.06                6.84
#> 5                  17.06                3.30
#> 6                  17.06                5.35

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.

exampEnviro = data.frame(
slope = 30,
windspeed = 60,
direction = 0,
waf = 0.2
)

exampCrownFuel = data.frame(
CBD = truncnorm::rtruncnorm(1000, a=0, coForest$cbd_kgm3[1], sd = coForest$cbd_kgm3[1]*.2),
FMC = rep(100,1000),
CBH = truncnorm::rtruncnorm(1000, a=0, coForest$cbh_m[1], sd = coForest$cbh_m[1]*.2),
CFL = truncnorm::rtruncnorm(1000, a=0, coForest$cfl_kgm2[1], sd = coForest$cfl_kgm2[1]*.2)
)

ex6 = rep(NA,200)
for (i in 1:200){
ex6[i] = rothermel(exampSurfFuel[1,], exampFuelMoisture[1,], exampCrownFuel[i,], exampEnviro, cfbForm = "sr")$fireBehavior[,3] } ggplot(data.frame(ex6), aes(ex6)) + geom_density(fill='blue',alpha = 0.5,adjust = 1/5) + xlab("Rate of spread (m/s)") + theme_classic() + coord_cartesian(expand=FALSE) In another example, we explore how rate of spread for a grass fire varies with wind direction and topographic slope. fuelModels['A5',] #> fuelModelType loadLitter load1hr load10hr load100hr loadLiveHerb #> A5 S 0 2.25 1.12 0 0 #> loadLiveWoody savLitter sav1hr sav10hr sav100hr savLiveHerb #> A5 4.48 0 6561.68 358 98 0 #> savLiveWoody fuelBedDepth mxDead heat description source #> A5 4921.26 61 20 18622 Brush Anderson (1982) exampSurfFuel = fuelModels['A5',] fuelMoisture['D1L1',] #> fmLitter fm1hr fm10hr fm100hr fmLiveHerb fmLiveWoody #> D1L1 3 3 4 5 30 60 #> description #> D1L1 Very dry dead FM, fully cured herb exampFuelMoisture = fuelMoisture['D1L1',] exampCrownFuel = data.frame( CBD = Inf, FMC = 100, CBH = Inf, CFL = Inf ) exampEnviro = data.frame( slope = 40, windspeed = 60, direction = 0, waf = 0.3 ) ex7 = expand.grid(winddir=seq(-180,180,10),slope = seq(0,30,5)) ex7$ros = NA
for (i in 1:259){
exampEnviro$direction = ex7$winddir[i]
exampEnviro$slope = ex7$slope[i]
ex7$ros[i] = rothermel(exampSurfFuel, exampFuelMoisture, exampCrownFuel, exampEnviro)$fireBehavior[,3]
}
ggplot(ex7,aes(winddir,ros,color=slope,group=slope)) +geom_path() + theme_classic() + labs(x='Wind direction (deg)',y='Rate of spread (m/s)')

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().

ba = c(10, 15)
ht = c(12, 20)
tph = c(100, 300)
type = c("df", "pp", "mc","lp")
cbind(type,canFuel(ba, ht, tph, type))
#>   type    cfl    cbd    cbh
#> 1   df 0.2862 0.0334  5.327
#> 2   pp 0.5354 0.1143  8.729
#> 3   mc 0.2371 0.0234  5.733
#> 4   lp 0.4231 0.0578 11.430

### 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)).

data(rrRAWS)
#>               dateTime temp_c rh windSpeed_kmh precip_mm
#> 1 04/01/2017 00:35 MDT   -5.0 98          0.00         0
#> 2 04/01/2017 01:35 MDT   -4.4 98          1.63         0
#> 3 04/01/2017 02:35 MDT   -5.6 96          3.23         0
#> 4 04/01/2017 03:35 MDT   -5.6 96          6.44         0
#> 5 04/01/2017 04:35 MDT   -6.1 97          4.83         0
#> 6 04/01/2017 05:35 MDT   -6.1 96          8.03         0

ff = rbind(
data.frame(ffm = ffm("simard",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="simard"), data.frame(ffm = ffm("wagner",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="wagner"),
data.frame(ffm = ffm("anderson",rrRAWS$rh, rrRAWS$temp_c)$fm1hr,method="anderson") ) ff$dateTime = rep(rrRAWS$dateTime,3) ggplot(ff, aes(ffm, color = method, fill = method)) + geom_density(alpha = 0.5) + xlab("Fine fuel moisture (%)") + theme_classic() 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}$ rrRAWS_ma = subset(rrRAWS, rh <= 70 & temp_c >= 10 & rh > 42.5 - 1.25 * temp_c & rh < 94.5 - 1.35 * temp_c) ff.ma = data.frame(ffm = ffm("mcarthur", rrRAWS_ma$rh, rrRAWS_ma$temp_c)$fm1hr, method = "mcarthur",
dateTime = rrRAWS_ma$dateTime) ff = rbind(data.frame(ffm = ffm("simard", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "simard"),
data.frame(ffm = ffm("wagner", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "wagner"), data.frame(ffm = ffm("anderson", rrRAWS$rh, rrRAWS$temp_c)$fm1hr, method = "anderson"))
ff$dateTime = rep(rrRAWS$dateTime, 3)
ff = rbind(ff, ff.ma)
ff$dateTime <- strptime(ff$dateTime, "%m/%d/%Y %H:%M")
ff$dateTime <- as.POSIXct(ff$dateTime)
ffm.plot = ggplot(ff, aes(x = dateTime, y = ffm, color = method)) + geom_smooth(span = 0.1,
method = "loess", se = F) + theme_classic() + labs(x = "Time", y = "Fuel moisture (%)")
print(ffm.plot)

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.

rh = rrRAWS$rh temp = rrRAWS$temp_c
month = as.numeric(format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%m")) hour = as.numeric(format(strptime(rrRAWS$dateTime, "%m/%d/%Y %H:%M"), "%H"))
ff.fbo = data.frame(ffm = ffm(method = "fbo", rh, temp, month, hour, asp = "N", slp = 10,
bla = "l", shade = "n")$fm1hr, method = "FBO", dateTime = strptime(rrRAWS$dateTime,
"%m/%d/%Y %H:%M"))
ff = rbind(ff, ff.fbo)
ffm.plot + scale_colour_grey() + geom_smooth(data = ff.fbo, aes(x = dateTime, y = ffm),
span = 0.1, method = "loess", se = F, color = "blue") + guides(color = FALSE) +
labs(x = "Time", y = "Fuel moisture (%)")

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:

waf(forestHt = 10, cr = 40, cc = 40)
#> [1] 0.209721

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$$.

waf(forestHt = 10, sheltered = "y")
#> [1] 0.1531586

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

waf(fuelDepth = 1)
#> [1] 0.2166844

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

waf(fuelDepth = 1, forestHt = 10, cr = 10, cc = 10)
#> [1] 0.2166844

### 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.

fireChart("Valley Fire",1000,50) +xlim(0,12000)

## 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) ###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)

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 %