Modeling Rate of Spread with Rothermel in Python

This page gives a minimal, correct Python implementation of the Rothermel (1972) surface rate-of-spread equation for a single fuel model, with the wind and slope factors, returning rate of spread in metres per minute. It is the focused, copy-pasteable core of the broader fire behavior and rate-of-spread modeling guide, which sits in the Fire Risk & Fuel Assessment section. Start here when you want one function you can trust before scaling it across a raster; the parent guide covers the gridded, multi-fuel version.

The whole model reduces to one ratio — heat supplied to unburned fuel over heat needed to ignite it — scaled up by wind and slope. Everything below is the arithmetic to compute each term for one homogeneous fuel bed.

Getting a single fuel model exactly right is the foundation for everything downstream: the same function, broadcast over a raster of fuel-model codes and slopes, becomes the landscape rate-of-spread surface described in the parent guide, and its output feeds Byram fireline intensity and flame length. So it pays to implement it once, test it against a known benchmark, and only then vectorize. The version here deliberately handles one dead-fuel size class; that keeps every term visible, which is what you want when you are learning the model or debugging a discrepancy with BehavePlus.

When to use a single-fuel implementation

A one-fuel-model function is the right tool for calibration, teaching, and unit tests; reach for a package or the gridded version when you have spatial or multi-fuel-class data.

Minimal reproducible example

The function works entirely in imperial units internally — the units the Rothermel coefficients were fit in — and converts the result to metric at the end. It implements the standard equations for a single dead-fuel class (e.g. Scott & Burgan GR1/GR2 grass, using the fine dead-fuel load).

import numpy as np

PARTICLE_DENSITY = 32.0      # lb/ft^3
TOTAL_MINERAL = 0.0555       # S_T
EFFECTIVE_MINERAL = 0.010    # S_e
HEAT_CONTENT = 8000.0        # Btu/lb (h)


def rothermel_ros(w0, sigma, depth, mf, mx, wind_mid_mph, slope_frac):
    """Rothermel surface rate of spread for one fuel model.

    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, e.g. 0.08)
    mx            moisture of extinction       (fraction, e.g. 0.15)
    wind_mid_mph  midflame wind speed          (mph)
    slope_frac    slope as rise/run            (tan of slope angle)

    Returns rate of spread in metres per minute.
    """
    # --- fuel-bed geometry ---
    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

    # --- damping coefficients ---
    rm = min(mf / mx, 1.0)                           # capped moisture ratio
    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)

    # --- reaction intensity I_R (Btu/ft^2/min) ---
    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)))
    reaction_intensity = gamma * net_load * HEAT_CONTENT * eta_moisture * eta_mineral

    # --- propagating flux ratio xi ---
    xi = (np.exp((0.792 + 0.681 * sigma ** 0.5) * (packing + 0.1))
          / (192.0 + 0.2595 * sigma))

    # --- wind factor phi_w (wind in ft/min) ---
    wind_ft_min = wind_mid_mph * 88.0
    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)
    phi_w = c * wind_ft_min ** b * (packing / packing_opt) ** -e

    # --- slope factor phi_s ---
    phi_s = 5.275 * packing ** -0.3 * slope_frac ** 2

    # --- heat sink ---
    epsilon = np.exp(-138.0 / sigma)                # effective heating number
    q_ig = 250.0 + 1116.0 * mf                      # heat of preignition (Btu/lb)
    heat_sink = bulk_density * epsilon * q_ig

    ros_ft_min = reaction_intensity * xi * (1.0 + phi_w + phi_s) / heat_sink
    return ros_ft_min * 0.3048                       # ft/min -> m/min


if __name__ == "__main__":
    # Scott & Burgan GR2 (low-load dry grass), 8% moisture, 5 mph midflame, 20% slope.
    ros = rothermel_ros(w0=0.0046, sigma=2000.0, depth=1.0,
                        mf=0.08, mx=0.15, wind_mid_mph=5.0, slope_frac=0.20)
    print(f"Rate of spread: {ros:.2f} m/min")

How the terms combine

Reading the function top to bottom traces the physics. The fuel-bed geometry block converts the fuel model into a packing ratio $\beta$ (how tightly the fuel is stacked) and compares it to the optimum packing ${\beta }_{op}$ — fuel that is too sparse or too dense burns slower than fuel at the optimum. The damping block penalizes reaction intensity for moisture and mineral content: as fuel moisture $M_f$ approaches the moisture of extinction $M_x$, the cubic damping polynomial drives the moisture coefficient toward zero. The reaction-intensity block multiplies the optimum reaction velocity $\gamma$, the net (mineral-free) load, the heat content, and the two damping coefficients into $I_R$, the heat release per unit area.

The propagating flux ratio $\xi$ is the fraction of that heat that actually reaches the unburned fuel ahead of the front — it depends on the fuel-bed geometry, not the weather. The wind and slope factors are pure multipliers on the no-wind, no-slope spread: ${\varphi }_w$ grows with midflame wind, and ${\varphi }_s$ with the square of the slope. Finally the heat sink — bulk density times effective heating number times heat of preignition — is the denominator, the energy the fire must spend to bring the next parcel of fuel to ignition. Divide supply by sink, scale by $(1+{\varphi }_w+{\varphi }_s)$, and convert feet-per-minute to metres-per-minute. Every line of the function maps to one of those physical quantities, which is exactly why a single-fuel implementation is the best place to build intuition.

Parameter reference

Every argument below changes the result; the fuel-bed geometry terms come from the fuel model, the moisture and wind terms from the day’s conditions.

Expected output and verification

For the GR2 example — dry grass, 8% moisture, 5 mph midflame wind, 20% slope — the function returns a rate of spread of a few metres per minute, consistent with BehavePlus for that fuel model. The most useful checks are monotonicity and the extinction boundary:

base = rothermel_ros(0.0046, 2000.0, 1.0, 0.08, 0.15, 0.0, 0.0)
windy = rothermel_ros(0.0046, 2000.0, 1.0, 0.08, 0.15, 10.0, 0.0)
steep = rothermel_ros(0.0046, 2000.0, 1.0, 0.08, 0.15, 0.0, 0.50)

assert base > 0.0, "dry grass must spread with no wind"
assert windy > base, "more wind must increase spread"
assert steep > base, "steeper slope must increase spread"

# At/above the moisture of extinction, spread collapses to ~zero.
dead = rothermel_ros(0.0046, 2000.0, 1.0, 0.15, 0.15, 5.0, 0.0)
assert dead < base, "fuel at extinction moisture barely spreads"

When mf equals mx, the capped moisture ratio drives the damping coefficient — and thus reaction intensity — to its floor, so rate of spread should nearly vanish. That behavior is the cleanest single test that your damping term is coded correctly.

Common pitfalls

• Unit inconsistency. The coefficients are imperial. Pass fuel load in lb/ft², SAV in 1/ft, depth in ft, and wind in mph (converted to ft/min inside). Feeding SI values gives a plausible-looking but wrong number with no error — convert at the boundary only.
• Uncapped moisture damping. If mf > mx and you do not cap mf / mx at 1.0, the damping polynomial goes negative and reaction intensity turns negative. Always min(mf / mx, 1.0).
• 10 m wind in place of midflame wind. Passing the open 10 m or 20 ft wind instead of the midflame wind overpredicts spread badly. Apply a wind-adjustment factor first.
• No wind limit at extreme wind. Rothermel’s ${\varphi }_w$ grows without bound; at very high winds it overpredicts. For extreme-wind scenarios cap the effective wind at the reaction-intensity-based limit rather than trusting the raw factor.

Related

• Fire Behavior and Rate-of-Spread Modeling — the parent guide with the gridded, multi-fuel pipeline and intensity/flame-length steps
• Fire Risk & Fuel Assessment — the section this task belongs to
• Wildfire Risk Modeling in Python — where this hazard estimate joins ignition likelihood
• Computing the Fire Weather Index in Python — supplies the fuel-moisture inputs for a design day

Frequently Asked Questions