Fuels, fire weather, and Rx fire behavior in North Dakota

Finalizing figures and results

Methods

Regression analysis

Figures generated by the script below and included in this document are available on GitHub, as well.

For analysis, data are filtered to \(\textsf{Rate of spread < 40 m/s}\) and \(\textsf{Maximum temperature > 40}^\circ \textsf{C}\).

Imputation of missing values

We used the multiple imputation method in the mice package to fill in missing values. The method creates 5 datasets, each with different but reasonable values for the missing data based on patterns in existing data, and performs the regression analysis on each. These results are pooled into a single set of results.

Model specification

Mixed-effect regression models fit with lme4::glmer and constructed as follows:

\[y \approx RH + WS + FM + FL + (1|location/year/plot),~ \text{family = gamma, link = log}\]

where

\(RH\) = Relative humidity,

\(WS\) = Average wind speed,

\(FM\) = Fuel moisture (from clipped samples),

\(FL\) = Fuel load (t/ha from clipped samples)

Weather variables \(RH\) and \(WS\) were taken from NDAWN stations for the hour in which the fire behavior observation occurred.

All variables were scaled to a common range prior to analysis, and units for each predictor variable are not reported.

Response variables included

• Maximum temperature (C) 15 cm above soil surface (mean of three thermocouples)
• Maximum temperature (C) at soil surface (single thermocouple)
• Rate of spread (m/s) through three points of a 1 m equilateral triangle

Thermocouple data recorded with the open-source FeatherFlame system.

Multivariate analysis

Principal Components Analysis (PCA) on fire behavior responses fit with vegan::rda. Post-hoc group (location) and gradient (fire weather) analysis with vegan::envfit.

An additional composite variable, vapor pressure deficit (VPD), was calculated as

\(\text{VPD} = e - e_s\)

where

\(e = 6.11 \cdot 10^{\frac{7.5 \cdot \text{Dewpoint}}{237.3 + \text{Dewpoint}}}\)

and

\(e_s = 6.11 \cdot 10^{\frac{7.5 \cdot \text{Air temperature}}{237.3 + \text{Air temperature}}}\)

Results

Fuel, weather, and fire behavior across locations

Variables associated with fire behavior

Rate of spread

Summary statistics for Rate of spread.

term	t	P
Relative humidity	-2.67	0.02
Wind speed	5.09	0
Fuel moisture	-2.78	0.03
Fuel load	-0.28	0.78

Aboveground temperatures

Summary statistics for flame temperature 15 cm above soil surface.

term	t	P
Relative humidity	-1.45	0.15
Wind speed	0.28	0.78
Fuel moisture	-1.95	0.08
Fuel load	1.26	0.27

Soil surface temperature

Summary statistics for soil surface temperature.

term	t	P
Relative humidity	-2.3	0.02
Wind speed	0.13	0.9
Fuel moisture	-0.19	0.86
Fuel load	0.46	0.66

Multivariate analysis

Eigenvalues and proportion explained by composite axes (principal components).

	PC1	PC2	PC3
Eigenvalue	1.574	1.002	0.4243
Proportion Explained	0.5246	0.334	0.1414
Cumulative Proportion	0.5246	0.8586	1

## Centroids:
##                PC1     PC2     PC3
## locationCG  0.0497  0.0106 -0.0348
## locationH  -0.1148 -0.0244  0.0805
## 
## Goodness of fit:
##              r2 Pr(>r)  
## location 0.0207  0.065 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Blocks:  strata 
## Permutation: free
## Number of permutations: 199

No difference in fire behavior patterns between Hettinger and Central Grasslands; this is expected.

## 
## ***VECTORS
## 
##                   PC1      PC2      PC3     r2 Pr(>r)   
## MaxWindSpeed -0.23784  0.87424 -0.42323 0.0388  0.298   
## AirTemp       0.64969 -0.75916  0.03972 0.0433  0.223   
## dpC           0.90056 -0.14982 -0.40809 0.0757  0.059 . 
## RH            0.60329  0.47220 -0.64270 0.1112  0.003 **
## VPD           0.30924 -0.89685  0.31627 0.0519  0.079 . 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## Blocks:  strata 
## Permutation: free
## Number of permutations: 999

dpC, RH, and VPD had statistically-significant associations with variation in fire behavior, but even these relationships were weak.

Principal Components Analysis of fire behavior data (response variables in red) for prescribed burns on rangeland at Hettinger (H) and Central Grasslands (CG) Research Extension Centers. No difference between locations and only the vapor pressure deficit (VPD) and relative humidty (RH) gradients had significant associations with variation in fire behavior.

It is interesting that the gradient in temperature at the soil surface (soil °C) is orthogonal to both canopy temperature and rate of spread (canopy °C and ROS, respectively). This confirms trends observed in the regression analysis that the conditions most conducive to fast, high-intensity fires are not the same that result in surface heating.

Managers might be able to use high windspeeds to mitigate soil heating when fuels are dry, or focus on conditions with low RH to ensure litter burns down to the soil surface.

Script

# Packages 
  pacman::p_load(tidyverse, readr, mice, broom.mixed, vegan,
                 lubridate, gridExtra, wesanderson)
set.seed(42)
# Data
  fp = url('https://raw.githubusercontent.com/devanmcg/SpatialFireBehavior/main/data/AnalysisDataKY.csv') 
  AnalysisData <- read_csv(fp) 
# Extra scripts
  source("https://raw.githubusercontent.com/devanmcg/IntroRangeR/master/11_IntroMultivariate/ordinationsGGplot.R")
  source('https://raw.githubusercontent.com/cran/mice/master/R/mipo.R')
###
### N O T   R U N
###
 AllData <-  
    read_csv(paste0(FilePath, "/data/fromMZ/CompiledData2.csv")) %>%
    filter(location != "OAK") %>%
      mutate(date = as.Date(date, format = "%m/%d/%Y"),
               L = str_remove(location, "REC"), 
               B = str_sub(block, 1,3), 
               Ps = str_replace(pasture, "[.]", ""), 
               Ps = str_sub(Ps, 1,2), 
               patch = str_replace(patch, "[.]", ""),
               y = format(date, "%y")) %>%
        unite("FireCode", c(L,B,Ps,patch,y), sep=".") %>%
    mutate(time = str_remove(MaxTempTime, "[.]+[0-9]"))%>%
    unite(timestamp, c(date, time), sep = " ") %>%
    mutate(timestamp = as.POSIXct(timestamp, format = "%Y-%m-%d %H:%M:%S")) %>%
        select(FireCode, timestamp, plot, array, TC, MaxC, 
               AirTemp, RH, dpC, WindSpeed, 
               LAI, FMC, KgHa) 
# Isolate soil surface temperature (TC 4)
  SoilTemp <-
    filter(AllData, TC == 4) %>%
      select(FireCode, plot, array, MaxC) %>%
        rename(SoilC = MaxC)
# Summarize array-level data 
  DataMeans <- 
    AllData %>%
      filter(TC %in% c('1', '2', '3')) %>% 
      select(-timestamp) %>%
      pivot_longer(cols = c(MaxC:KgHa), 
                   names_to = "var",
                   values_to = "value") %>%
      group_by(FireCode, plot, array, var) %>%
        summarize(Mean = mean(value) ) %>%
      ungroup() %>%
      pivot_wider(names_from = var, 
                  values_from = Mean) 
# Calculate Vapor Pressure Deficit
  DataMeans <- 
    DataMeans %>%
      mutate(e  = 6.11*(10^((7.5*dpC)/(237.3+dpC))), 
             es = 6.11*(10^((7.5*AirTemp)/(237.3+AirTemp))), 
             VPD = es - e) %>%
      select(-e, -es)
# Calculate ROS 
  D = 1   # Distance between thermocouples (m)
  ROS <- 
    AllData %>%
      filter(TC %in% c('1', '2', '3')) %>%
      mutate(timestamp = format(timestamp, "%H:%M:%S"), 
             ArrivalTime = seconds(hms(timestamp)) ) %>%
    select(FireCode, plot, array, ArrivalTime) %>%
    group_by(FireCode, plot, array) %>%
    arrange(ArrivalTime, .by_group = TRUE) %>% 
    mutate(position = order(order(ArrivalTime, decreasing=FALSE)), 
           position = recode(position, "1"="a", "2"="b", "3"="c"), 
           ArrivalTime = as.numeric(ArrivalTime) /60 ) %>%
    spread(position, ArrivalTime)  %>%
    ungroup %>% 
    mutate( theta_rad = atan((2*c - b - a) / (sqrt(3)*(b - a))), 
            ros = case_when(
              a == b ~ (sqrt(3) / 2) / (c - a) , 
              a != b ~  (D*cos(theta_rad) / (b - a) ) 
            )) %>%
    select(-a, -b, -c, -theta_rad)
# Combine, filter into final tibble for analysis 
  AnalysisData  <- 
    full_join(DataMeans, ROS) %>%
              left_join(SoilTemp) %>%
      filter( ros <= 40, 
              MaxC >= 40) %>%
      rename(FuelMoisture = FMC, 
             SoilMaxC = SoilC) %>%
        mutate(FuelMoisture = ifelse(FuelMoisture >= 0, 
                                      FuelMoisture, NA), 
               FuelMoisture = FuelMoisture * 100) %>%
      separate(FireCode, into = c("location", "block", "pasture", 
                           "patch", "year"), 
               remove = F)
###
###
###
# Calculate imputed datasets on scaled data
  imp_sc <- AnalysisData %>% 
          select(-LAI, -JD) %>%
          mutate_at(vars(AirTemp:tHa), ~as.numeric(scale(., center=F))) %>%
                    mutate(across(location:array, as.factor)) %>%
                     mice(m=5, seed = 23109, print=F)
# Reduce mids object to tibble for multivariate analysis    
    imp_raw <- 
      complete(imp_sc, 'long') %>%
      as_tibble() %>%
      unite("TreeID", c(location, block, pasture,
                        year, plot, array), sep = ".") %>%
      select(-patch, -.id, -.imp, -FireCode) %>%
      pivot_longer(names_to = "response", 
                   values_to = "values", 
                   -TreeID) %>%
      group_by(TreeID, response) %>%
      summarize(value = median(values)) %>%
      ungroup() %>%
      pivot_wider(names_from = response, 
                  values_from = value) %>%
      separate(TreeID, c("location", "block", "pasture",
                         "year", "plot", "array")) %>% 
      mutate(across(location:array, as.factor))
# Mixed-effect regression models 
#
# Rate of spread
#   
  # Fit model 
    ros_imp <- with(imp_sc, suppressMessages(
                lme4::glmer(ros ~ RH + WindSpeed + 
                              FuelMoisture + tHa + 
                              (1|location/block/year/plot), 
                            family=Gamma(link = "log"), 
                            control=lme4::glmerControl(optimizer="bobyqa", 
                                                       optCtrl=list(maxfun=100000)) )) )
    ros_imp_pooled <- pool(ros_imp)$pooled  %>%
                          as_tibble()
    # Get terms
    ros_terms <- 
      full_join(
      summary(pool(ros_imp)) %>% 
        as_tibble() %>%
        rownames_to_column("row"), 
  
      confint.mipo(pool(ros_imp)) %>%
        as_tibble() %>%
        rownames_to_column("row") ) 
    

#
# Maximum canopy temperature  
#
# Fit model 
canopy_imp <- with(imp_sc, suppressMessages(
              lme4::glmer(MaxC ~ RH + WindSpeed +  
                            FuelMoisture + tHa + 
                            (1|location/block/year/plot), 
                          family=Gamma(link = "log"), 
                          control=lme4::glmerControl(optimizer="bobyqa", 
                                                     optCtrl=list(maxfun=100000)) )) )

  canopy_imp_pooled <- pool(canopy_imp)$pooled  %>%
                      as_tibble()
  
  canopy_terms <- 
    full_join(
      summary(pool(canopy_imp)) %>% 
        as_tibble() %>%
        rownames_to_column("row"), 
      
      confint.mipo(pool(canopy_imp)) %>%
        as_tibble() %>%
        rownames_to_column("row") ) 

#
# Maximum soil surface temperature  
#
# Fit model 
  soil_imp <- with(imp_sc, suppressMessages(
                  lme4::glmer(log(SoilMaxC+1) ~ RH + WindSpeed +  
                      FuelMoisture + tHa + 
                      (1|location/block/year/plot), 
                    family=Gamma(link = "log"), 
                    control=lme4::glmerControl(optimizer="bobyqa", 
                                               optCtrl=list(maxfun=100000)) )) )

soil_imp_pooled <- pool(soil_imp)$pooled  %>%
  as_tibble()

soil_terms <- 
  full_join(
    summary(pool(soil_imp)) %>% 
      as_tibble() %>%
      rownames_to_column("row"), 
    
    confint.mipo(pool(soil_imp)) %>%
      as_tibble() %>%
      rownames_to_column("row") ) 


response_CIs <-   
  bind_rows(
    mutate(soil_terms, response = 'Temp. at surface'), 
    mutate(canopy_terms, response = 'Temp. above surface'),
    mutate(ros_terms, response = "Rate of spread") 
  ) %>%
  as_tibble() %>%
  rename(lwr = `2.5 %`, upr = `97.5 %`) %>%
  filter(term != '(Intercept)') %>%
  mutate(term = recode(term, tHa = 'Fuel load', 
                             WindSpeed = 'Wind speed', 
                             RH = 'Relative humidity', 
                             FuelMoisture = 'Fuel moisture')) 
response_CIs %>%
  filter(response == "Rate of spread") %>%
  rename(t = statistic, 
         P = p.value) %>%
    select(term, t, P) %>%
  mutate(across(t:P, ~round(., 2))) %>%
  pander::pander("Summary statistics for Rate of spread.") 
response_CIs %>%
  filter(response == "Temp. above surface") %>%
  rename(t = statistic, 
         P = p.value) %>%
    select(term, t, P) %>%
  mutate(across(t:P, ~round(., 2))) %>%
  pander::pander("Summary statistics for flame temperature 15 cm above soil surface.") 
response_CIs %>%
  filter(response == "Temp. at surface") %>%
  rename(t = statistic, 
         P = p.value) %>%
    select(term, t, P) %>%
  mutate(across(t:P, ~round(., 2))) %>%
  pander::pander("Summary statistics for soil surface temperature.") 
# Fire behavior PCA 
  fb_d <- 
    imp_raw %>%
    select(MaxC, ros, SoilMaxC) 
  
  fb_pca <- rda(fb_d ~ 1, 'euc', scale = T)
# Testing for differences between locations 
pca_fc <-  envfit( fb_pca ~ location, imp_raw, 
                   choices = c(1:3), 
                   strata = imp_raw$year, 
                   199)$factors
pca_fc
# Testing fire weather against PCA
pca_ev <- envfit(fb_pca ~ MaxWindSpeed+AirTemp+dpC+RH+VPD, 
                 data = imp_raw, 
                 choices = c(1:3), 
                 strata = imp_raw$location)
pca_ev