Crown Fire

Modeling passive and active crown fire

Crown fire occurs when fire spreads through the canopy layer of trees. This tutorial covers crown fire modeling in Elmfire.jl.

using Elmfire
using Plots

Crown Fire Basics

Crown fire is classified as: - Passive (torching): Individual trees or groups torch, but fire doesn’t spread continuously through canopy - Active: Fire spreads continuously through the canopy, driven by wind

Canopy Properties

Crown fire requires canopy data:

ncols, nrows = 100, 100
cellsize = 30.0

# Create canopy grid with typical forest values
canopy = CanopyGrid{Float64}(
    fill(0.15, ncols, nrows),   # cbd: Canopy bulk density (kg/m³)
    fill(2.0, ncols, nrows),     # cbh: Canopy base height (m)
    fill(0.6, ncols, nrows),     # cc: Canopy cover (fraction)
    fill(20.0, ncols, nrows)     # ch: Canopy height (m)
)

println("Canopy properties:")
println("  Bulk density: 0.15 kg/m³")
println("  Base height: 2.0 m")
println("  Cover: 60%")
println("  Height: 20 m")

Canopy properties:
  Bulk density: 0.15 kg/m³
  Base height: 2.0 m
  Cover: 60%
  Height: 20 m

Enabling Crown Fire

Configure simulation for crown fire:

config = SimulationConfig{Float64}(
    enable_crown_fire = true,
    enable_spotting = false,
    crown_fire_adj = 1.0,
    critical_canopy_cover = 0.4,
    foliar_moisture = 100.0
)

SimulationConfig{Float64}(true, false, 1.0, 0.4, 100.0, 250.0, 1.0, nothing, false)

Surface Fire vs Crown Fire

Compare simulations with and without crown fire:

fuel_table = create_standard_fuel_table(Float64)
fuel_ids = fill(10, ncols, nrows)  # FBFM10: Timber understory

slope = zeros(Float64, ncols, nrows)
aspect = zeros(Float64, ncols, nrows)

weather = ConstantWeather{Float64}(
    wind_speed_mph = 20.0,
    wind_direction = 270.0,
    M1 = 0.04, M10 = 0.06, M100 = 0.08,
    MLH = 0.50, MLW = 0.80
)

# Surface fire only
state_surface = FireState{Float64}(ncols, nrows, cellsize)
ignite!(state_surface, 50, 50, 0.0)

config_surface = SimulationConfig{Float64}(enable_crown_fire = false)
weather_interp = create_constant_interpolator(weather, ncols, nrows, cellsize)

simulate_full!(state_surface, fuel_ids, fuel_table, weather_interp, slope, aspect,
    0.0, 30.0; config = config_surface)

# Crown fire enabled
state_crown = FireState{Float64}(ncols, nrows, cellsize)
ignite!(state_crown, 50, 50, 0.0)

config_crown = SimulationConfig{Float64}(
    enable_crown_fire = true,
    critical_canopy_cover = 0.4,
    foliar_moisture = 100.0
)

simulate_full!(state_crown, fuel_ids, fuel_table, weather_interp, slope, aspect,
    0.0, 30.0; canopy = canopy, config = config_crown)

# Compare
p1 = heatmap(state_surface.burned',
    title = "Surface Fire Only\n$(round(get_burned_area_acres(state_surface), digits=1)) acres",
    color = :YlOrRd, aspect_ratio = 1)

p2 = heatmap(state_crown.burned',
    title = "Crown Fire Enabled\n$(round(get_burned_area_acres(state_crown), digits=1)) acres",
    color = :YlOrRd, aspect_ratio = 1)

plot(p1, p2, layout = (1, 2), size = (800, 400))

Critical Fireline Intensity

Crown fire initiation depends on fireline intensity reaching the critical threshold:

# Calculate critical intensity for different canopy base heights
cbh_values = 1.0:0.5:10.0  # meters
foliar_moisture = 100.0  # percent

i_crit = [critical_fireline_intensity(cbh, foliar_moisture) for cbh in cbh_values]

plot(cbh_values, i_crit,
    xlabel = "Canopy Base Height (m)",
    ylabel = "Critical Fireline Intensity (kW/m)",
    title = "Crown Fire Initiation Threshold",
    linewidth = 2,
    legend = false
)

Higher canopy base heights require greater fireline intensity to initiate crown fire.

Canopy Bulk Density Effect

Canopy bulk density (CBD) enters the Cruz (2005) model twice, and the two effects pull in opposite directions.

It raises the potential crown spread rate only weakly, as CBD^0.19:

\[ \text{CROSA} = 11.02 \, U_{10}^{0.9} \, \text{CBD}^{0.19} e^{-17 M_1} \]

But it sets the critical rate for a self-sustaining active crown fire much more strongly, as 3/CBD. The crown activity coefficient is the ratio of the two, and the fire goes active once it exceeds 1:

\[ R_0 = \frac{3}{\text{CBD}}, \qquad \text{CAC} = \frac{\text{CROSA}}{R_0} \]

So CBD’s real influence is on whether the fire crowns actively, not on how fast it then runs. The left panel below shows where the two curves cross.

cbds = 10 .^ range(log10(0.01), log10(0.35), length = 120)
winds = [15.0, 25.0, 40.0]
M1 = 0.04

# Cruz (2005) potential active crown spread rate, ft/min
crosa(cbd, ws) = 11.02 * (ws * 1.609 / 0.87)^0.9 * cbd^0.19 * exp(-0.17 * 100 * M1) / 0.3048

# Critical spread rate for active crowning (3/CBD m/min), ft/min
r0(cbd) = (3.0 / cbd) / 0.3048

p1 = plot(xscale = :log10, xlabel = "Canopy Bulk Density (kg/m³)",
    ylabel = "Spread Rate (ft/min)", title = "Potential rate vs. active-crowning threshold",
    ylims = (0, 800), legend = :topright, left_margin = 5Plots.mm, bottom_margin = 6Plots.mm)
for ws in winds
    plot!(p1, cbds, crosa.(cbds, ws), label = "CROSA @ $(Int(ws)) mph", linewidth = 2)
end
plot!(p1, cbds, r0.(cbds), label = "R₀ = 3/CBD", color = :black,
      linestyle = :dot, linewidth = 2)

# Mark where each wind speed crosses into active crown fire
for ws in winds
    lo, hi = 0.005, 0.5
    for _ in 1:60
        mid = (lo + hi) / 2
        crosa(mid, ws) < r0(mid) ? (lo = mid) : (hi = mid)
    end
    scatter!(p1, [(lo + hi) / 2], [r0((lo + hi) / 2)], label = "", color = :black, markersize = 4)
end

# What the model actually spreads at, once the ELMFIRE cap is applied
p2 = plot(xscale = :log10, xlabel = "Canopy Bulk Density (kg/m³)",
    ylabel = "Crown Spread Rate (ft/min)", title = "Rate used by the model",
    ylims = (0, 400), legend = :bottomright, left_margin = 5Plots.mm, bottom_margin = 6Plots.mm)
# The 25 and 40 mph curves sit on top of each other once capped, so offset the
# line widths to keep both visible
for (ws, lw, ls) in zip(winds, [4, 2.5, 1.5], [:solid, :solid, :dash])
    rates = map(cbds) do cbd
        canopy_pt = CanopyProperties{Float64}(cbd = cbd, cbh = 2.0, cc = 0.6, ch = 20.0)
        crown_spread_rate(canopy_pt, 1e9, ws, M1, 20.0).spread_rate
    end
    plot!(p2, cbds, rates, label = "$(Int(ws)) mph", linewidth = lw, linestyle = ls)
end
hline!(p2, [250.0], label = "CROWN_FIRE_SPREAD_RATE_LIMIT",
       color = :red, linestyle = :dash, linewidth = 2)

plot(p1, p2, layout = (1, 2), size = (980, 430))

The right panel is the practical consequence. Below the threshold the fire is passive and its rate is damped by exp(-CAC) — which falls as CBD rises, since the damping tightens faster than the weak CBD^0.19 term grows, so a denser canopy briefly torches more slowly right up until it goes active. Above the threshold the fire is active but immediately runs into ELMFIRE’s CROWN_FIRE_SPREAD_RATE_LIMIT of 250 ft/min. For the CBD range of real conifer stands (roughly 0.05–0.30 kg/m³) every curve is already flat against that cap.

This is why mapping burned area for several CBD values produces four identical pictures: once the fire is actively crowning, the cap has removed CBD from the spread rate entirely. The threshold crossing in the left panel is the effect worth showing.

Foliar Moisture Effect

fm_values = [80, 100, 120, 140]
plots = []

for fm in fm_values
    config_fm = SimulationConfig{Float64}(
        enable_crown_fire = true,
        foliar_moisture = Float64(fm)
    )

    state = FireState{Float64}(ncols, nrows, cellsize)
    ignite!(state, 50, 50, 0.0)

    simulate_full!(state, fuel_ids, fuel_table, weather_interp, slope, aspect,
        0.0, 25.0; canopy = canopy, config = config_fm)

    push!(plots, heatmap(state.burned',
        title = "FM = $fm%",
        color = :YlOrRd, aspect_ratio = 1, colorbar = false
    ))
end

plot(plots..., layout = (2, 2), size = (700, 700),
    plot_title = "Crown Fire vs Foliar Moisture")

Combined Spread Rate

The total spread rate combines surface and crown fire components:

# Get canopy properties
canopy_props = CanopyProperties{Float64}(0.15, 2.0, 0.6, 20.0)

# Range of surface fireline intensities
flin_range = 100:100:5000

surface_ros = 50.0  # ft/min
combined_ros = Float64[]

for flin in flin_range
    cr = crown_spread_rate(canopy_props, Float64(flin), 20.0, 0.06, surface_ros)
    # SpreadResult fields: velocity, vs0, ir, hpua, flin, phiw, phis
    sr = SpreadResult{Float64}(surface_ros, surface_ros, 500.0, 1000.0, Float64(flin), 1.0, 0.0)
    cros = combined_spread_rate(sr, cr)
    push!(combined_ros, cros)
end

plot(flin_range, combined_ros,
    xlabel = "Surface Fireline Intensity (kW/m)",
    ylabel = "Combined Spread Rate (ft/min)",
    title = "Crown Fire Contribution to Spread Rate",
    linewidth = 2,
    legend = false
)
hline!([surface_ros], linestyle = :dash, label = "Surface only")