Spread Model

The SpreadModel module provides composable, differentiable components for driving level set fire simulations. Instead of manually constructing a spread rate matrix F, you build a FireSpreadModel from pluggable wind, moisture, and terrain components — each a callable (t, x, y) struct.

Architecture

FireSpreadModel(fuel, wind, moisture, terrain)
       │
       ▼
model(t, x, y) → rate_of_spread(fuel; moisture, wind, slope)
       │
       ▼
simulate!(grid, model) → level set evolution

Each component is a callable struct with signature (t, x, y):

Component	Returns	Example
`wind::AbstractWind`	`(speed, direction)` — speed [km/h], direction [radians]	`UniformWind(speed=8.0)`
`moisture::AbstractMoisture`	`FuelClasses` — moisture fractions	`UniformMoisture(...)`
`terrain::AbstractTerrain`	`(slope, aspect)` — slope [fraction], aspect [radians]	`FlatTerrain()`

Quick Example

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

model = FireSpreadModel(
    SHORT_GRASS,
    UniformWind(speed=8.0),
    UniformMoisture(M),
    FlatTerrain()
)

# Evaluate spread rate at a single point
model(0.0, 100.0, 100.0)

31.119112730486354

# Full simulation
grid = LevelSetGrid(200, 200, dx=30.0)
ignite!(grid, 3000.0, 3000.0, 50.0)
simulate!(grid, model, steps=100, dt=0.5)

fig = Figure()
ax = Axis(fig[1, 1], title="t = $(grid.t) min", aspect=DataAspect(),
    xlabel="x (m)", ylabel="y (m)")
fireplot!(ax, grid)
fig

Directional Spread Models

By default, FireSpreadModel uses cosine blending to distribute the head-fire rate across angles. This produces fires that are wider than typically observed. The elliptical model (Anderson, 1983) produces more realistic elongated fire shapes by computing a length-to-breadth ratio from wind speed.

Pass a directional model as the 5th argument to FireSpreadModel:

CosineBlending() (default) — R(θ) = R_base + (R_head - R_base) · max(0, cos θ)
EllipticalBlending() — R(θ) = R_head · (1 - ε) / (1 - ε · cos θ) where ε is fire eccentricity

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

fig = Figure(size=(700, 300))
for (col, (label, dir)) in enumerate([
    ("Cosine (default)", CosineBlending()),
    ("Elliptical (Anderson)", EllipticalBlending()),
])
    g = LevelSetGrid(200, 200, dx=30.0)
    ignite!(g, 3000.0, 3000.0, 50.0)
    m = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain(), dir)
    simulate!(g, m, steps=100, dt=0.5)
    ax = Axis(fig[1, col], title=label, aspect=DataAspect())
    fireplot!(ax, g)
    hidedecorations!(ax)
end
fig

Length-to-Breadth Formulas

The EllipticalBlending model computes a length-to-breadth ratio (LB) from midflame wind speed using a published formula. Two options are available:

:anderson (default) — Anderson (1983), general-purpose
:green — Green (1983), suitable for grass fuels

EllipticalBlending(formula=:green)  # use Green's formula

EllipticalBlending(:green)

The helper functions length_to_breadth and fire_eccentricity are exported for use in custom spread models:

U = 2.0  # m/s midflame wind speed
LB = length_to_breadth(U)
ε = fire_eccentricity(LB)
(LB=LB, ε=ε)

(LB = 1.5049627492754702, ε = 0.7473164746215192)

Wind

`UniformWind`

Spatially and temporally constant wind field.

wind = UniformWind(speed=10.0, direction=π/4)
wind(0.0, 0.0, 0.0)  # (speed, direction)

(10.0, 0.7853981633974483)

Moisture

`UniformMoisture`

Spatially and temporally constant fuel moisture.

moist = UniformMoisture(FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0))
moist(0.0, 0.0, 0.0)

FuelClasses{Float64}(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

`DynamicMoisture`

Spatially varying fuel moisture that responds to fire-induced drying. As the fire front approaches, radiative heat dries unburned fuel ahead of the front. The 1-hr dead fuel moisture (d1) varies spatially while other size classes remain constant.

The drying model at each unburned cell:

\[\frac{dM}{dt} = -\frac{\text{dry\_rate}}{\phi^2 + 1} + \text{recovery\_rate} \cdot (M_{\text{ambient}} - M)\]

where \(\phi\) is the level set value (approximate distance to the fire front in meters).

grid_d = LevelSetGrid(100, 100, dx=30.0)
ignite!(grid_d, 1500.0, 1500.0, 100.0)

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)
dm = DynamicMoisture(grid_d, M, dry_rate=0.1, recovery_rate=0.001)

DynamicMoisture{Float64, Matrix{Float64}}([0.06 0.06 … 0.06 0.06; 0.06 0.06 … 0.06 0.06; … ; 0.06 0.06 … 0.06 0.06; 0.06 0.06 … 0.06 0.06], FuelClasses{Float64}(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0), 0.06, 0.1, 0.001, 0.03, 30.0, 30.0, 0.0, 0.0)

Parameter guidance:

Parameter	Default	Typical Range	Effect
`dry_rate`	`0.1`	`0.01` – `0.5`	How fast fire dries nearby fuel. Higher values = stronger pre-heating effect, faster fire acceleration.
`recovery_rate`	`0.001`	`0.0001` – `0.01`	How fast fuel rewets toward ambient. Higher values = fuel recovers quickly when fire moves away.
`min_d1`	`0.03`	`0.02` – `0.05`	Minimum 1-hr dead fuel moisture (fraction). Below ~0.03, wildland fuels are essentially oven-dry.

The dry_rate is the most sensitive parameter — it controls the positive feedback loop where fire dries adjacent fuel, increasing spread rate, which dries more fuel. Start with the default (0.1) and adjust based on whether the simulated fire accelerates too quickly or too slowly.

Comparing static vs. dynamic moisture:

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

# Static
grid_s = LevelSetGrid(150, 150, dx=20.0)
ignite!(grid_s, 1500.0, 1500.0, 50.0)
model_s = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain())
simulate!(grid_s, model_s, steps=150, dt=0.5)

# Dynamic
grid_d = LevelSetGrid(150, 150, dx=20.0)
ignite!(grid_d, 1500.0, 1500.0, 50.0)
model_d = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), DynamicMoisture(grid_d, M), FlatTerrain())
simulate!(grid_d, model_d, steps=150, dt=0.5)

fig = Figure(size=(700, 300))
ax1 = Axis(fig[1, 1], title="Static Moisture\n$(count(<(0), grid_s.φ)) cells burned",
    aspect=DataAspect())
fireplot!(ax1, grid_s)
hidedecorations!(ax1)

ax2 = Axis(fig[1, 2], title="Dynamic Moisture\n$(count(<(0), grid_d.φ)) cells burned",
    aspect=DataAspect())
fireplot!(ax2, grid_d)
hidedecorations!(ax2)
fig

Terrain

`FlatTerrain`

Zero slope everywhere.

FlatTerrain()(0.0, 0.0, 0.0)  # (slope, aspect)

(0.0, 0.0)

`UniformSlope`

Spatially constant terrain slope.

slope = UniformSlope(slope=0.3, aspect=0.0)
slope(0.0, 0.0, 0.0)

(0.3, 0.0)

Effect of slope on fire spread:

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

fig = Figure(size=(800, 250))
for (col, s) in enumerate([0.0, 0.3, 0.6])
    g = LevelSetGrid(150, 150, dx=20.0)
    ignite!(g, 1500.0, 1500.0, 50.0)
    m = FireSpreadModel(SHORT_GRASS, UniformWind(speed=5.0), UniformMoisture(M), UniformSlope(slope=s))
    simulate!(g, m, steps=150, dt=0.5)
    ax = Axis(fig[1, col], title="slope = $s", aspect=DataAspect())
    fireplot!(ax, g)
    hidedecorations!(ax)
end
fig

Fuel Breaks and Unburnable Areas

The LevelSetGrid tracks per-cell ignition state via t_ignite. By default all cells are burnable (t_ignite = Inf). Cells marked as unburnable have t_ignite = NaN, which sets their spread rate to zero.

Use set_unburnable! to mark circular regions, or directly set grid.t_ignite[i, j] = NaN for arbitrary shapes. The convenience function burnable(grid) returns a BitMatrix of burnable cells.

This is useful for representing water bodies, roads, fuel breaks, and pre-existing fire scars.

Creating a Fuel Break

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)
model = FireSpreadModel(SHORT_GRASS, UniformWind(speed=0.0), UniformMoisture(M), FlatTerrain())

fig = Figure(size=(700, 300))

# Without fuel break
g1 = LevelSetGrid(200, 200, dx=30.0)
ignite!(g1, 3000.0, 3000.0, 50.0)
simulate!(g1, model, steps=200, dt=0.5)
ax1 = Axis(fig[1, 1], title="No fuel break", aspect=DataAspect())
fireplot!(ax1, g1)
hidedecorations!(ax1)

# With fuel break (unburnable vertical strip)
g2 = LevelSetGrid(200, 200, dx=30.0)
ignite!(g2, 3000.0, 3000.0, 50.0)
xs = xcoords(g2)
for j in eachindex(xs)
    if 3300.0 <= xs[j] <= 3420.0
        g2.t_ignite[:, j] .= NaN
    end
end
simulate!(g2, model, steps=200, dt=0.5)
ax2 = Axis(fig[1, 2], title="With fuel break", aspect=DataAspect())
fireplot!(ax2, g2)
vspan!(ax2, 3300.0, 3420.0, color=(:black, 0.15))
hidedecorations!(ax2)
fig

Marking Pre-Burned Areas

To represent an existing fire scar that should not re-ignite, combine set_unburnable! with setting φ < 0:

grid = LevelSetGrid(200, 200, dx=30.0)

# Pre-existing fire scar: set φ < 0 AND mark unburnable
cx, cy, r = 2500.0, 3000.0, 300.0
xs = xcoords(grid)
ys = ycoords(grid)
for j in eachindex(xs), i in eachindex(ys)
    if hypot(xs[j] - cx, ys[i] - cy) <= r
        grid.φ[i, j] = -1.0
    end
end
set_unburnable!(grid, cx, cy, r)

# New ignition away from the scar
ignite!(grid, 4000.0, 3000.0, 50.0)

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)
model = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain())
simulate!(grid, model, steps=100, dt=0.5)

fig = Figure()
ax = Axis(fig[1, 1], title="Fire scar blocks spread", aspect=DataAspect(),
    xlabel="x (m)", ylabel="y (m)")
fireplot!(ax, grid)
fig

Without set_unburnable!, the fire scar boundary would be treated as an active fire front and expand outward.

Burnout

By default, once a cell ignites it continues to contribute to fire spread indefinitely. In reality, fuel is consumed and the flaming front passes. The burnout keyword in simulate! models this by scaling the spread rate based on how long a cell has been burning.

Three burnout models are available:

Type	Formula	Behavior
`NoBurnout()`	`1.0`	Fire spreads indefinitely (default)
`ExponentialBurnout(τ)`	`exp(-t/τ)`	Smooth exponential decay over residence time `τ`
`LinearBurnout(τ)`	`max(0, 1 - t/τ)`	Linear ramp to zero over residence time `τ`

The residence time τ is computed from the Rothermel fuel model using Anderson’s (1969) formula:

\[t_r = \frac{384}{\sigma_{\text{char}} \cdot 60} \text{ [min]}\]

where \(\sigma_{\text{char}}\) is the characteristic surface-area-to-volume ratio [ft\(^{-1}\)].

using Wildfires.Rothermel: residence_time

# Residence time varies by fuel type
(SHORT_GRASS = residence_time(SHORT_GRASS),
 CHAPARRAL = residence_time(CHAPARRAL),
 TIMBER_UNDERSTORY = residence_time(TIMBER_UNDERSTORY))

(SHORT_GRASS = 0.0018285714285714285, CHAPARRAL = 0.0036797903975260777, TIMBER_UNDERSTORY = 0.003626517214741357)

Using Burnout in Simulation

Pass a burnout model to simulate!:

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

fig = Figure(size=(700, 300))

# Without burnout
g1 = LevelSetGrid(200, 200, dx=30.0)
ignite!(g1, 3000.0, 3000.0, 50.0)
m1 = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain())
simulate!(g1, m1, steps=200, dt=0.5)
ax1 = Axis(fig[1, 1], title="No burnout\n$(count(<(0), g1.φ)) cells burned", aspect=DataAspect())
fireplot!(ax1, g1)
hidedecorations!(ax1)

# With exponential burnout
g2 = LevelSetGrid(200, 200, dx=30.0)
ignite!(g2, 3000.0, 3000.0, 50.0)
t_r = residence_time(SHORT_GRASS)
simulate!(g2, m1, steps=200, dt=0.5, burnout=ExponentialBurnout(t_r))
ax2 = Axis(fig[1, 2], title="Exponential burnout (τ=$(round(t_r, sigdigits=3)) min)\n$(count(<(0), g2.φ)) cells burned", aspect=DataAspect())
fireplot!(ax2, g2)
hidedecorations!(ax2)
fig

Note

For fast-burning fuels like short grass, the residence time is very small (< 0.01 min), so burnout primarily affects the spread rate near the ignition point where the fire has been burning longest. For heavier fuels with longer residence times, the effect is more pronounced.

Custom Components

To create a custom component, define a callable struct that subtypes AbstractWind, AbstractMoisture, or AbstractTerrain.

For example, a wind field that varies in space:

struct GradientWind <: AbstractWind
    base_speed::Float64
    gradient::Float64   # speed increase per meter in x
end

function (w::GradientWind)(t, x, y)
    speed = w.base_speed + w.gradient * x
    direction = 0.0  # wind from the west
    (speed, direction)
end

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)
g = LevelSetGrid(150, 150, dx=20.0)
ignite!(g, 1500.0, 1500.0, 50.0)
m = FireSpreadModel(SHORT_GRASS, GradientWind(5.0, 0.003), UniformMoisture(M), FlatTerrain())
simulate!(g, m, steps=150, dt=0.5)

fig = Figure()
ax = Axis(fig[1, 1], title="Gradient Wind", aspect=DataAspect(), xlabel="x (m)", ylabel="y (m)")
fireplot!(ax, g)
fig

For dynamic components that respond to fire state, also implement update!:

SpreadModel.update!(w::MyDynamicWind, grid::LevelSetGrid, dt) = ...

Simulation Parameters

The simulate! function accepts several keyword arguments that control numerical behavior:

Parameter	Default	Description
`steps`	`100`	Number of time steps to run
`dt`	`nothing`	Time step in minutes. When `nothing`, computed automatically via the CFL condition each step.
`cfl`	`0.5`	Safety factor for the CFL condition (`dt = cfl * min(dx,dy) / max(F)`). Only used when `dt=nothing`.
`reinit_every`	`10`	Reinitialization frequency. Every N steps, the level set is restored to a signed distance function.
`solver`	`Godunov()`	Numerical solver for the level set equation. See Solvers.
`curvature`	`0.0`	Curvature regularization coefficient. When > 0, smooths the fire front via `F_eff = max(F - b·κ, 0)`. See Curvature.
`reinit`	`IterativeReinit()`	Reinitialization method. `IterativeReinit()` (default) or `NewtonReinit()` for machine-precision signed distance. See Reinitialization.

solver: Controls the numerical scheme for advancing the level set. Godunov() (default) uses first-order upwind differencing with forward Euler. Superbee() uses a second-order Superbee flux limiter with RK2 (Heun’s method) for sharper fire fronts. WENO5() uses a 5th-order WENO scheme with SSP-RK3 for the highest accuracy. See the Solvers section for details.

simulate!(grid, model, steps=100)                          # Godunov (default)
simulate!(grid, model, steps=100, solver=Superbee())       # 2nd-order TVD
simulate!(grid, model, steps=100, solver=WENO5())          # 5th-order WENO
simulate!(grid, model, steps=100, curvature=5.0)           # curvature smoothing
simulate!(grid, model, steps=100, reinit=NewtonReinit())   # Newton reinit

cfl: The default 0.5 is conservative — it uses half the maximum stable time step. Values closer to 1.0 are faster but risk instability if the spread rate field changes rapidly. If you see oscillations or divergence, lower this value.

reinit_every: Over many steps, φ drifts from a true signed distance function, degrading the accuracy of the upwind scheme. More frequent reinitialization (e.g., 5) improves accuracy at modest cost; less frequent (e.g., 50) is faster but may introduce artifacts in long simulations. For most cases, the default of 10 works well.

dt: When set to nothing (default), simulate! calls spread_rate_field! first, then computes a CFL-limited time step automatically. If you pass a fixed dt, you are responsible for ensuring it satisfies the CFL condition — violating it causes numerical blowup.

Animation with Trace

The Trace type records snapshots of the level set field during simulation, making it easy to create animations of fire spread. Pass trace=Trace(grid, every) to simulate! to record φ every every steps.

M = FuelClasses(d1=0.06, d10=0.07, d100=0.08, herb=0.0, wood=0.0)

grid = LevelSetGrid(200, 200, dx=30.0)
ignite!(grid, 3000.0, 3000.0, 50.0)
model = FireSpreadModel(SHORT_GRASS, UniformWind(speed=8.0), UniformMoisture(M), FlatTerrain(),
    EllipticalBlending())

trace = Trace(grid, 5)
simulate!(grid, model, steps=150, dt=0.5, trace=trace)

length(trace.stack)  # initial + every 5 steps

Each entry in trace.stack is a (time, φ) tuple. Use firegif to create a GIF:

firegif(joinpath(@__DIR__, "fire_spread.gif"), trace, grid)

"/home/runner/work/Wildfires.jl/Wildfires.jl/docs/fire_spread.gif"

Performance: Float64 vs Float32

The level set grid and all simulation operations support Float32 for reduced memory and potentially faster computation. To create a Float32 simulation, pass Float32 values to all constructors:

Note

The LevelSetGrid constructor determines its element type via promote_type on all keyword arguments (dx, dy, x0, y0). The defaults are Float64, so you must pass Float32 values for all of them — e.g. LevelSetGrid(200, 200, dx=30f0, x0=0f0, y0=0f0).

using BenchmarkTools

function run_simulation(::Type{T}, n=500) where T
    grid = LevelSetGrid(n, n,
        dx=T(10), x0=zero(T), y0=zero(T))
    ignite!(grid, T(2500), T(2500), T(50))
    F = fill(T(10), size(grid))
    for _ in 1:200
        advance!(grid, F, T(0.25))
    end
    reinitialize!(grid)
    grid
end

b64 = @benchmark run_simulation(Float64)
b32 = @benchmark run_simulation(Float32)
nothing

Metric	Float64	Float32	Ratio
Median time	317.84 ms	241.78 ms	1.31x
Memory	397.1 MiB	198.8 MiB	2.00x

References

Rothermel, R.C. (1972). A Mathematical Model for Predicting Fire Spread in Wildland Fuels. Res. Paper INT-115, USDA Forest Service.
Osher, S. & Sethian, J.A. (1988). Fronts propagating with curvature-dependent speed. J. Computational Physics, 79(1), 12-49.