Fire Behavior and Rate-of-Spread Modeling in Python

When a fuel-treatment analyst asks “how fast and how hot would a fire move through this stand under a design-day wind?”, the answer comes from a surface fire-behavior model. This guide implements the Rothermel (1972) surface spread model in Python — the quasi-empirical energy balance that underpins BehavePlus, FlamMap, and most operational fire-behavior tools — and derives rate of spread, Byram fireline intensity, and flame length across a raster grid. It is a workflow within the Fire Risk & Fuel Assessment section, and it produces the hazard side of the picture that complements the likelihood surface built in wildfire risk modeling.

Rothermel’s model predicts the steady-state rate of spread of a surface fire from the fuel bed’s physical properties, its moisture, wind at midflame height, and terrain slope. It is deterministic and fast, which is exactly why it scales to per-pixel evaluation over a landscape.

Prerequisites checklist

Fire-behavior arithmetic is unforgiving about units — the classical coefficients are in imperial units (Btu, ft, lb, min), and mixing them with SI is the single most common mistake.

Python 3.10+ with numpy >= 1.24 and rasterio >= 1.3 for the gridded run
A fuel-model assignment per pixel — a standard model (Anderson 13 or Scott & Burgan 40) or measured fuel-bed parameters from the fuel load mapping workflow
Dead and live fuel moisture by size class, from a moisture model or the Fire Weather Index
Midflame wind speed (not 10 m wind — apply a wind-adjustment factor first)
A slope raster in the same grid as the fuel layer, from Digital Terrain Model Generation
A single, documented unit system — this guide works in imperial internally and converts outputs to metric

Concept background: the Rothermel energy balance

Rothermel frames spread as a ratio of heat received by unburned fuel ahead of the front to the heat needed to raise that fuel to ignition. The rate of spread is:

R = \frac{I _{R} ξ ( 1 + ϕ _{w} + ϕ _{s} )}{ρ _{b} ε Q _{i g}}

where $I_{R}$ is the reaction intensity (heat release rate per unit area), $ξ$ is the propagating flux ratio (the fraction of that heat that reaches adjacent fuel), $ϕ_{w}$ and $ϕ_{s}$ are dimensionless wind and slope factors that add to the no-wind, no-slope case, $ρ_{b}$ is the fuel-bed bulk density, $ε$ is the effective heating number, and $Q_{i g}$ is the heat of preignition. The numerator is heat supply; the denominator is the heat sink.

The wind and slope factors are what make spread directional and terrain-sensitive:

ϕ_{w} = C (U)^{B} (\frac{β}{β _{o p}})^{- E}, ϕ_{s} = 5.275 β^{- 0.3} (tan φ)^{2}

where $U$ is midflame wind speed, $φ$ is slope angle, $β$ is the packing ratio, $β_{o p}$ its optimum, and $C$ , $B$ , $E$ are functions of the surface-area-to-volume ratio $σ$ . Slope enters as $tan^{2} φ$ , so spread accelerates sharply on steep ground.

Two properties of this formulation matter for how you use it. First, spread scales linearly with reaction intensity in the numerator but is divided by the heat sink, so a fuel bed that is both fine (high $σ$ , easy to ignite) and dry (low moisture) spreads fast on both counts — the effects compound. Second, wind and slope enter additively as $(1 + ϕ_{w} + ϕ_{s})$ , which means on a calm day on flat ground the model still predicts a non-zero baseline spread driven purely by radiant and convective heat transfer within the fuel bed. That baseline is small but real, and it is the number you should benchmark against a published no-wind, no-slope table before trusting any wind- or slope-loaded prediction.

Once you have rate of spread, Byram’s (1959) fireline intensity converts it to energy release along the fire front:

I = H w R

with $H$ the fuel heat yield, $w$ the fuel consumed per unit area, and $R$ the rate of spread. Fireline intensity (kW/m) is the physically meaningful hazard number, and it maps to visible flame length through Byram’s relation:

L = 0.0775 I^{0.46}

with $L$ in metres and $I$ in kW/m. Flame length is what crews use to judge suppression tactics — under ~1.2 m, hand tools are effective; over ~2.4 m, direct attack fails.

Fuel, moisture, wind, and slope drive the Rothermel rate of spread, which feeds Byram fireline intensity and flame length — the operational hazard outputs, evaluated per pixel.

Step-by-step Python pipeline

We build the model bottom-up: fuel-bed intermediates, the no-wind reaction intensity, the wind and slope multipliers, then rate of spread, intensity, and flame length — finally vectorized across a raster grid. All internal arithmetic is imperial; conversions happen only at the boundaries.

Step 1 — Fuel-bed intermediates

The packing ratio and its optimum, the net fuel load, and the damping coefficients are derived once from the fuel-model constants and moisture.

import numpy as np

PARTICLE_DENSITY = 32.0     # lb/ft^3, oven-dry wood
TOTAL_MINERAL = 0.0555      # S_T
EFFECTIVE_MINERAL = 0.010   # S_e
HEAT_CONTENT = 8000.0       # Btu/lb, h


def fuel_bed_intermediates(w0, sigma, depth, mf, mx):
    """Packing ratios, damping coefficients, and reaction velocity terms.

    w0    oven-dry fuel load (lb/ft^2)
    sigma surface-area-to-volume ratio (1/ft)
    depth fuel bed depth (ft)
    mf    fuel moisture (fraction)
    mx    moisture of extinction (fraction)
    """
    bulk_density = w0 / depth                       # rho_b
    packing = bulk_density / PARTICLE_DENSITY       # beta
    packing_opt = 3.348 * sigma ** -0.8189          # beta_op
    net_load = w0 / (1.0 + TOTAL_MINERAL)           # w_n

    rm = min(mf / mx, 1.0)
    eta_moisture = 1.0 - 2.59 * rm + 5.11 * rm ** 2 - 3.52 * rm ** 3
    eta_mineral = min(0.174 * EFFECTIVE_MINERAL ** -0.19, 1.0)

    a = 133.0 * sigma ** -0.7913
    gamma_max = sigma ** 1.5 / (495.0 + 0.0594 * sigma ** 1.5)
    gamma = (gamma_max * (packing / packing_opt) ** a
             * np.exp(a * (1.0 - packing / packing_opt)))
    return dict(bulk_density=bulk_density, packing=packing, packing_opt=packing_opt,
                net_load=net_load, eta_moisture=eta_moisture,
                eta_mineral=eta_mineral, gamma=gamma)

Step 2 — Reaction intensity, propagating flux, and heat sink

Reaction intensity is the heat release rate; the propagating flux ratio and heat of preignition close the energy balance.

def reaction_intensity(im):
    """I_R = gamma' * w_n * h * eta_M * eta_s  (Btu/ft^2/min)."""
    return im["gamma"] * im["net_load"] * HEAT_CONTENT * im["eta_moisture"] * im["eta_mineral"]


def propagating_flux_ratio(sigma, packing):
    """xi: fraction of reaction heat reaching adjacent fuel."""
    return np.exp((0.792 + 0.681 * sigma ** 0.5) * (packing + 0.1)) / (192.0 + 0.2595 * sigma)


def heat_sink(im, sigma, mf):
    """rho_b * epsilon * Q_ig  (Btu/ft^3)."""
    epsilon = np.exp(-138.0 / sigma)                # effective heating number
    q_ig = 250.0 + 1116.0 * mf                      # heat of preignition (Btu/lb)
    return im["bulk_density"] * epsilon * q_ig

Step 3 — Wind and slope factors

These dimensionless multipliers add directional and terrain sensitivity. Midflame wind is in ft/min; slope is the tangent of the slope angle.

def wind_factor(sigma, packing, packing_opt, wind_ft_min):
    c = 7.47 * np.exp(-0.133 * sigma ** 0.55)
    b = 0.02526 * sigma ** 0.54
    e = 0.715 * np.exp(-3.59e-4 * sigma)
    return c * wind_ft_min ** b * (packing / packing_opt) ** -e


def slope_factor(packing, slope_tan):
    """phi_s = 5.275 * beta^-0.3 * (tan slope)^2."""
    return 5.275 * packing ** -0.3 * slope_tan ** 2

Step 4 — Rate of spread, intensity, and flame length

Assemble the full Rothermel ratio, then apply Byram for intensity and flame length. Conversions to metric happen here.

FT_MIN_TO_M_MIN = 0.3048
BTU_FT_MIN_TO_KW_M = 0.05767      # Btu/ft/min -> kW/m


def rate_of_spread(w0, sigma, depth, mf, mx, wind_ft_min, slope_tan):
    im = fuel_bed_intermediates(w0, sigma, depth, mf, mx)
    ir = reaction_intensity(im)
    xi = propagating_flux_ratio(sigma, im["packing"])
    sink = heat_sink(im, sigma, mf)
    phi_w = wind_factor(sigma, im["packing"], im["packing_opt"], wind_ft_min)
    phi_s = slope_factor(im["packing"], slope_tan)
    ros_ft_min = ir * xi * (1.0 + phi_w + phi_s) / sink
    return ros_ft_min, ir


def fireline_intensity(ir, ros_ft_min, residence_depth_ft):
    """Byram I = H * w * R, expressed via reaction intensity and residence time."""
    # Flame residence time (Anderson 1969): t_r ~ 384 / sigma; here we use the
    # reaction-intensity form I = I_R * t_r * R for a compact estimate.
    heat_per_area = ir * residence_depth_ft            # Btu/ft^2
    intensity_btu = heat_per_area * ros_ft_min         # Btu/ft/min
    return intensity_btu * BTU_FT_MIN_TO_KW_M          # kW/m


def flame_length_m(intensity_kw_m):
    return 0.0775 * intensity_kw_m ** 0.46             # Byram, metres

Step 5 — Vectorize across a fuel and terrain grid

The scalar functions above are pure numpy arithmetic, so they broadcast directly over raster arrays. Read the fuel-model parameters per pixel (via a lookup keyed on the fuel-model class raster) and the slope raster, then evaluate.

import rasterio

with rasterio.open("fuel_model.tif") as f:
    profile = f.profile
    fuel_class = f.read(1)
with rasterio.open("slope_deg.tif") as s:
    slope_tan = np.tan(np.radians(s.read(1)))

# Map each fuel class to its parameter arrays (broadcast to grid shape).
w0 = np.vectorize(LOAD_LOOKUP.get)(fuel_class).astype("float32")
sigma = np.vectorize(SAV_LOOKUP.get)(fuel_class).astype("float32")
depth = np.vectorize(DEPTH_LOOKUP.get)(fuel_class).astype("float32")
mx = np.vectorize(MX_LOOKUP.get)(fuel_class).astype("float32")

mf = 0.08                                # dead fuel moisture fraction, design day
wind_ft_min = 5.0 * 88.0                 # 5 mph midflame -> ft/min

ros, ir = rate_of_spread(w0, sigma, depth, mf, mx, wind_ft_min, slope_tan)
ros_m_min = ros * FT_MIN_TO_M_MIN
intensity = fireline_intensity(ir, ros, residence_depth_ft=depth)
flame = flame_length_m(intensity)

out = profile.copy(); out.update(dtype="float32", count=1)
with rasterio.open("rate_of_spread_m_min.tif", "w", **out) as dst:
    dst.write(ros_m_min.astype("float32"), 1)

The single-fuel-model, fully worked scalar version — the right place to start before scaling up — is covered in modeling rate of spread with Rothermel in Python.

Configuration: a fuel-model parameter block

Fuel models are just parameter sets. Keep them in a version-controlled structure rather than inline literals; the Scott & Burgan GR2 (low-load, dry-climate grass) values are illustrative.

{
  "fuel_model": "GR2",
  "description": "Low load, dry climate grass",
  "load_1h_lb_ft2": 0.0046,
  "load_live_herb_lb_ft2": 0.0092,
  "sav_1h_1_ft": 2000,
  "sav_live_herb_1_ft": 1800,
  "bed_depth_ft": 1.0,
  "moisture_of_extinction_frac": 0.15,
  "heat_content_btu_lb": 8000
}

Validation and verification

Anchor the implementation against a published benchmark, then guard the physical bounds. BehavePlus results for standard fuel models are the usual reference; a no-wind, no-slope run should reproduce the model’s tabulated baseline rate of spread within a few percent.

# GR2 baseline: no wind, no slope, 6% moisture -> low but non-zero spread.
ros0, ir0 = rate_of_spread(w0=0.0046, sigma=2000.0, depth=1.0,
                           mf=0.06, mx=0.15, wind_ft_min=0.0, slope_tan=0.0)
assert ros0 > 0.0, "no-wind spread must be positive for dry grass"

# Adding wind and slope can only increase spread.
ros_w, _ = rate_of_spread(0.0046, 2000.0, 1.0, 0.06, 0.15, 440.0, 0.30)
assert ros_w > ros0, "wind and slope must raise rate of spread"

# Physical sanity on the grid.
assert np.all(ros_m_min >= 0.0)
assert np.all(flame >= 0.0)

At the moisture of extinction, the damping coefficient drives reaction intensity — and therefore spread — to zero; a fuel wetter than mx should not spread, which is a useful assertion.

For a stronger check, compare a full response surface against BehavePlus output for the same fuel model: sweep midflame wind from 0 to 20 mph and slope from 0 to 60% and confirm your rate-of-spread grid matches the tool within a few percent across the range, not just at a single point. Discrepancies that grow with wind almost always trace to a missing wind limit; discrepancies that grow with slope point to a sign or exponent error in $ϕ_{s}$ . Keep the benchmark as a regression test so a future refactor cannot silently drift the physics.

Failure modes and gotchas

Mixing unit systems. The Rothermel coefficients are calibrated in imperial units. Feeding SI fuel loads or m/s wind into the imperial equations produces numbers that look plausible and are wrong. Convert only at the input and output boundaries.
10 m wind instead of midflame wind. Rothermel needs wind at midflame height. Using the 10 m or 20 ft open wind overpredicts spread badly; apply a wind-adjustment factor for the canopy first.
Moisture above extinction not zeroed. If mf > mx, the moisture damping term can go negative in naive code. Clamp the moisture ratio at 1.0 so reaction intensity floors at zero rather than turning negative.
Unbounded wind factor. At very high winds Rothermel’s $ϕ_{w}$ grows without limit and overpredicts; operational implementations cap effective wind at a reaction-intensity-dependent limit. Apply the wind limit for extreme winds.
Confusing rate of spread with intensity. A fast grass fire can be low-intensity; a slow timber fire can be high-intensity. Flame length follows intensity, not spread — report both.
Live vs dead moisture treated identically. Multi-size-class fuel models weight dead and live fuels differently; collapsing them to one moisture value misrepresents live-dominated fuels.

Performance and scale notes

Broadcast, do not loop. Every equation is elementwise numpy; evaluating the whole raster in one broadcast pass is far faster than a per-pixel Python loop.
Tile large landscapes. Process in windows so a landscape-scale fuel grid never has to fit in RAM, and stream rate of spread, intensity, and flame length to separate output bands.
Precompute fuel-model constants. The packing ratio, optimum packing, and reaction-velocity terms depend only on the fuel model, not on weather — compute them once per class and look them up.
Vectorize the parameter lookup. np.vectorize over a dict is convenient but slow for huge grids; for production, build lookup arrays indexed by fuel-model code and use fancy indexing.
Cache design-day runs. A landscape only needs re-evaluation when weather or fuel changes; persist the outputs and re-run on the scenario, not on every read.

Frequently Asked Questions

What does the Rothermel model actually predict?

The steady-state rate of spread of a surface fire — how fast the flaming front moves through the fuel bed — from the fuel’s physical properties, its moisture, midflame wind, and slope. It is a quasi-empirical energy balance: heat received by unburned fuel ahead of the front divided by the heat needed to bring that fuel to ignition. It models surface fire only; crown fire needs a separate initiation and spread model.

Why do I need midflame wind rather than the 10 m wind?

Rothermel’s wind factor is calibrated to the wind actually pushing the flame, which inside or beneath a canopy is much lower than the open 10 m or 20 ft wind. Using the higher wind overpredicts spread, sometimes by a large factor. Apply a wind-adjustment factor — a function of canopy cover and fuel-bed height — to convert the measured wind to midflame height first.

How does slope change the rate of spread?

Slope enters through the factor phi_s, which scales with the square of the tangent of the slope angle. Because of the tan-squared term, spread accelerates sharply on steep ground — a fire moves much faster uphill than on the flat, which is why steep slopes account for a disproportionate share of entrapments. The model treats upslope spread; downslope and cross-slope require decomposing the wind-slope vector.

What is the difference between rate of spread and fireline intensity?

Rate of spread is how fast the front moves (m/min). Fireline intensity is the rate of energy release per unit length of fire front (kW/m), given by Byram as heat yield times fuel consumed times rate of spread. A fast, light grass fire can be low-intensity; a slow, heavy timber fire can be high-intensity. Flame length — and therefore suppression difficulty — follows intensity, not spread.

How accurate is a Rothermel prediction?

It is a well-validated operational model for surface fire in the fuel types it was built for, and it underpins BehavePlus and FlamMap. But it assumes a uniform, continuous, steady-state fuel bed; real fires with heterogeneous fuels, wind gusts, spotting, or crown involvement can diverge substantially. Treat its output as a physically grounded estimate for planning, calibrated against local observations, not a guarantee.

Fire Risk & Fuel Assessment — the parent section this behavior workflow belongs to
Modeling Rate of Spread with Rothermel in Python — a minimal single-fuel-model implementation
Wildfire Risk Modeling in Python — pairs this hazard model with ignition likelihood
Fuel Load Mapping from LiDAR — supplies the fuel-bed load the model consumes
Computing the Fire Weather Index in Python — provides the fuel-moisture context for a design day

Explore this section

Modeling Rate of Spread with the Rothermel Model in Python (Single Fuel Model)