Mathematical Models
This section describes the mathematical models implemented in the wildfire level-set solver.
Fire Spread Models
Multiple fire spread model families are implemented. Select the active
model with fire_spread_model:
rothermel— Rothermel (1972) empirical surface fire spread (default)balbi— Balbi (2009) physics-based radiation-driven spreadcheney_gould— Cheney & Gould (1995/1998) Australian grassland empirical modelfbp_o1a/fbp_o1b— Canadian FBP grass fuel types (matted / standing)fbp_s1/fbp_s2/fbp_s3— Canadian FBP slash fuel typeslautenberger— Lautenberger (2013) semi-physical Eulerian modelcruz_crown— Cruz et al. active crown fire rate of spread
Rothermel (1972) Fire Spread Model
The Rothermel fire spread model is the foundation for computing the rate of spread (ROS) of wildfire. The model computes the no-wind, no-slope rate of spread and then applies wind and slope correction factors.
Basic Rate of Spread
The no-wind, no-slope rate of spread \(R_0\) (ft/min) is given by:
where:
\(I_R\) is the reaction intensity (BTU/ft²/min)
\(\xi\) is the propagating flux ratio
\(\rho_b\) is the oven-dry bulk density (lb/ft³)
\(\epsilon_h\) is the effective heating number
\(Q_{ig}\) is the heat of preignition (BTU/lb)
Reaction Intensity
The reaction intensity \(I_R\) is computed as:
where:
\(\Gamma'\) is the optimum reaction velocity (1/min)
\(w_n\) is the net fuel loading (lb/ft²)
\(h\) is the heat content of fuel (BTU/lb)
\(\eta_M\) is the moisture damping coefficient
\(\eta_s\) is the mineral damping coefficient
The optimum reaction velocity is:
where:
\(\sigma\) is the surface-area-to-volume ratio (ft²/ft³)
\(\beta\) is the packing ratio (dimensionless)
\(\beta_{op}\) is the optimum packing ratio
Packing Ratio
The packing ratio is defined as:
where:
\(\rho_b = w_0 / \delta\) is the oven-dry bulk density (lb/ft³)
\(\rho_p\) is the oven-dry particle density (lb/ft³)
\(w_0\) is the oven-dry fuel loading (lb/ft²)
\(\delta\) is the fuel bed depth (ft)
Moisture Damping
Single-class path (aggregate)
When no per-class loads are specified (the default), the moisture damping coefficient \(\eta_M\) is computed from a single composite moisture content \(M_f\) and extinction moisture \(M_x\):
where \(r_m = \min(M_f / M_x, 1)\), \(M_f\) is the fuel moisture content, and \(M_x\) is the moisture of extinction.
Multi-class path (per size class)
When per-class fuel loads are provided (via a fuel model database entry or
explicit rothermel.w_d1, etc. inputs), the full Rothermel (1972)
multi-class formulation is used. Each size class has its own oven-dry fuel
load \(w_i\), surface-area-to-volume ratio \(\sigma_i\), and moisture
content \(M_i\).
The fixed SAV values for coarser dead-fuel classes are: \(\sigma_{d10} = 109\ \text{ft}^{-1}\) and \(\sigma_{d100} = 30\ \text{ft}^{-1}\).
The per-class weighting factor is:
Dead-fuel composite moisture ratio (Rothermel 1972, Eq. 86):
Live-fuel moisture of extinction (Rothermel 1972, Eq. 88 approximation):
Live-fuel composite moisture ratio:
Both dead and live moisture damping coefficients use the cubic polynomial (Eq. 29).
Reaction intensity sums dead and live contributions:
where \(w_{n,\text{dead}} = \sum_{i \in \text{dead}} w_i (1-S_T)\) and \(w_{n,\text{live}} = \sum_{i \in \text{live}} w_i (1-S_T)\).
The characteristic SAV \(\sigma_c\) (load-weighted mean across all classes) is used in place of the single aggregate \(\sigma\) for \(\Gamma'\), \(\xi\), \(\epsilon_h\), \(Q_{ig}\), and the wind-factor coefficients \(C\), \(B\), \(E\).
Mineral Damping
The mineral damping coefficient \(\eta_s\) is:
where \(S_e\) is the effective mineral content.
Propagating Flux Ratio
The propagating flux ratio \(\xi\) is:
Heat of Preignition
The heat of preignition \(Q_{ig}\) (BTU/lb) depends on fuel moisture:
Wind Factor
The wind factor \(\phi_w\) multiplies the base rate of spread to account for wind effects:
where:
\(U\) is the midflame wind speed (ft/min)
\(C = 7.47 \exp(-0.133 \sigma^{0.55})\)
\(B = 0.02526 \sigma^{0.54}\)
\(E = 0.715 \exp(-3.59 \times 10^{-4} \sigma)\)
The rate of spread with wind is:
Slope Factor
The slope factor \(\phi_s\) accounts for terrain slope:
where \(\theta\) is the slope angle. The slope can be decomposed into x and y components:
Combined Wind and Slope
When both wind and slope effects are present, the total rate of spread is:
Balbi (2009) Physical Fire Spread Model
The Balbi model provides an alternative rate-of-spread calculation rooted in
radiation physics rather than empirical curve fits. It is activated by setting
fire_spread_model = balbi and automatically replaces the Rothermel model in
both the level-set advection path and the pre-computed \(R\) field used by
FARSITE. The model is fully compatible with all six Viegas wind-terrain options
and FARSITE landscape file inputs.
Physical Basis
The rate of spread \(R\) [m/s] is derived from the balance between radiant heat flux from a tilted flame and the energy required to ignite the unburned fuel ahead:
where:
\(A_{\text{coeff}} = \chi\,\sigma_m\,\delta_m / (2\,\tau_0\,B^*)\) [m/s] — fuel radiative amplitude
\(\chi = r_{00}\,\sigma_m / (1 + r_{00}\,\sigma_m)\) — forward radiation fraction
\(B^* = (C_{pf}\,(T_i - T_a) + M_f\,\Lambda) / h\) — dimensionless ignition energy
\(\tan\alpha = U / v_b + \tan\theta\) — flame tilt from wind and slope
\(v_b = \sqrt{g\,\delta_m\,(T_f - T_a) / T_a}\) — buoyancy velocity scale [m/s]
Inputs from the LCP / Fuel Database
The Balbi model shares fuel parameters with the Rothermel model:
Parameter |
Source |
Description |
|---|---|---|
\(\sigma\) [ft⁻¹→m⁻¹] |
Fuel database / LCP |
Surface-area-to-volume ratio |
\(\delta\) [ft→m] |
Fuel database / LCP |
Fuel bed depth |
\(w_0\) [lb/ft²→kg/m²] |
Fuel database / LCP |
Oven-dry fuel load |
\(\rho_p\) [lb/ft³→kg/m³] |
Fuel database / LCP |
Particle density |
\(h\) [BTU/lb→J/kg] |
Fuel database / LCP |
Heat of combustion |
\(M_f\) [−] |
|
Fuel moisture fraction |
slope, aspect |
LCP / terrain file |
Terrain slope components |
Wind \(U\) [m/s] |
|
Wind speed |
Extra Inputs via Parmparse (balbi.*)
Parameter |
Default |
Description |
|---|---|---|
|
300.0 K |
Ambient temperature |
|
1000.0 K |
Mean flame temperature |
|
600.0 K |
Ignition temperature |
|
2.26e6 J/kg |
Latent heat of water vaporisation |
|
1800.0 J/(kg·K) |
Specific heat of dry fuel |
|
2.5e-4 m |
Radiation length scale |
|
75591.0 s/m |
Residence-time coefficient |
Auto-Generated Balbi Fuel Parameter Table
When fire_spread_model = balbi, the solver automatically prints a table of
pre-computed Balbi parameters for every fuel model in the active database
(FBFM13 or FBFM40) at startup. Columns are: fuel code, short name,
\(\sigma_m\), \(\delta_m\), \(\chi\), \(B^*\), \(v_b\),
\(A_{\text{coeff}}\), and the predicted ROS at 5 m/s wind on flat ground.
Reference
Balbi, J.-H., Rossi, J.-L., Marcelli, T., and Santoni, P.-A. (2009). “A physical model for wildland fires.” Combustion and Flame, 156(12), 2217–2230. https://doi.org/10.1016/j.combustflame.2009.07.010
Cheney & Gould (1995 / 1998) Grassland Fire Spread Model
The Cheney–Gould model is a purely empirical rate-of-spread (ROS) formula calibrated against a large number of experimental grassland fires conducted in Australia. It is activated by setting fire_spread_model = cheney_gould and generally outperforms the Rothermel model in open grassland fuels.
Head-Fire Rate of Spread
The head-fire rate of spread \(R\) (km/h) is a piecewise-linear function of the 10-m open wind speed \(U_{10}\) (km/h), a moisture correction factor \(f_{MC}\), and a curing factor \(CF\):
The result is clamped to \(R \ge 0\) — in the high-wind regime the intercept \(-0.020\) makes the base ROS slightly negative for extremely low wind speeds just above 5 km/h, which is physically meaningless — and then converted to m/s for the simulation:
The 10-m wind speed \(U_{10}\) is obtained by converting the simulation wind field (m/s) via \(U_{10} = |\mathbf{u}| \times 3.6\).
Moisture Correction Factor
where \(MC\) is the dead fine fuel moisture content [%]. At \(MC = 0\) (bone-dry fuel) \(f_{MC} = 1\); at \(MC \approx 22\%\) the ROS is halved relative to dry conditions.
Curing Factor
\(CF \in [0,\,1]\) represents the degree of curing of the grass: \(CF = 1\) is fully cured (dry standing grass), while \(CF = 0\) is completely green (no spread). The factor is clamped to the \([0,1]\) interval at runtime.
Note on Terrain Slope
Terrain slope is not accounted for in the original Cheney–Gould empirical formulation. The model was calibrated on flat or gently sloping Australian grasslands and the slope correction is intentionally omitted in this implementation. For slope effects, use the Rothermel or Balbi models instead.
Canadian Forest Fire Behaviour Prediction (FBP) System
The Canadian FBP system provides empirical rate-of-spread equations calibrated for specific fuel types. The solver implements five fuel types: O1a (matted grass), O1b (standing grass), S1 (Jack/Lodgepole Pine slash), S2 (White Spruce/Balsam slash), and S3 (Coastal Cedar/Hemlock/Douglas Fir slash).
Select with fire_spread_model = fbp_o1a (or fbp_o1b, fbp_s1, fbp_s2, fbp_s3).
Grass fuel types (O1a / O1b)
The fine-fuel moisture factor follows Van Wagner (1985):
The wind factor and Initial Spread Index (ISI) are:
where \(U_{10}\) is the 10-m wind speed (km/h) and \(M_f\) is the fine fuel moisture content (%). The curing factor (CF) accounts for the fraction of dead fuel:
Rate of spread [m/s]:
Coefficients per fuel type:
O1a (matted grass): \(C_c = 190/60,\; k_c = 0.031\)
O1b (standing grass): \(C_c = 250/60,\; k_c = 0.035\)
Slash fuel types (S1 / S2 / S3)
Same ISI computation; \(\mathrm{CF} = 1\). Rate of spread [m/s]:
Coefficients:
S1: \(C_s = 75/60,\; k_s = 0.110\)
S2: \(C_s = 200/60,\; k_s = 0.062\)
S3: \(C_s = 320/60,\; k_s = 0.010\)
Parameters: fbp.fuel_type, fbp.curing (curing [%], O1a/O1b only), fbp.moisture (dead fine fuel moisture [%]).
Reference: Forestry Canada Fire Danger Group (1992). Development and Structure of the Canadian Forest Fire Behavior Prediction System. Information Report ST-X-3.
Lautenberger (2013) Physics-Based Fire Spread Model
Select with fire_spread_model = lautenberger.
The model derives rate of spread from a semi-physical fuel-energy balance calibrated against US wildfire data (Lautenberger 2013):
where \(\sigma_m\) is the fuel SAV ratio [m-1] (converted from ft-1), \(M_f\) is fuel moisture (fraction), \(A_L\) is a pre-factor [m²/s], and \(B_L\) is the moisture sensitivity.
Wind and slope corrections are applied linearly:
Parameter |
Default |
Description |
|---|---|---|
|
|
Pre-factor \(A_L\) [m²/s] |
|
|
Moisture sensitivity \(B_L\) [-] |
|
|
Wind correction \(C_L\) [(m/s)-1] |
|
|
Slope correction \(D_L\) [-] |
Reference: Lautenberger, C. (2013). Wildland fire modeling with an Eulerian level set method and automated calibration. Fire Safety Journal, 62, 289–298.
ROS Floor (FARSITE Stall Threshold)
A FARSITE-style stall threshold is applied across all spread models: cells where
the computed ROS is below WildfireConst::ROS_MIN_MS (1 × 10⁻⁴ m/s) are
forced to zero. This prevents numerical drift in extremely low-moisture conditions
and matches FARSITE’s internal practice of not propagating the fire front when
spread is negligible.
Parameters: None (constant in constants.H).
Non-Burnable Cell Masking
Fuel model codes that correspond to non-burnable landscape types are explicitly zeroed in the ROS kernel. This prevents fire from numerically creeping across water bodies, bare rock, or urban areas due to floating-point noise in zero-fuel cells.
Affected codes (FBFM40 / FBFM13): 91 (Urban/Developed), 92 (Snow/Ice), 93 (Agriculture), 98 (Open Water), 99 (Bare Ground).
When rothermel.landscape_file is provided and the per-cell fuel code maps to
one of these values, the kernel returns R = 0 for that cell without entering
the wind-slope-fuel computation.
Parameters: None. Automatic when rothermel.landscape_file is used.
Propagation Methods
Three propagation methods are available and selected with propagation_method:
levelset (WENO5-Z level-set advection), farsite (FARSITE Huygens-wavelet
expansion), and mtt (Minimum Travel Time fast-march).
Level-Set Method
The level-set function \(\phi(x,y,t)\) represents the fire front as the zero level set \(\phi = 0\). The evolution equation is:
where \(V\) is the local rate of spread computed from the fire models above.
Numerical Scheme
The level-set equation is discretized using:
Upwind finite differences for spatial derivatives
Explicit time stepping (e.g., forward Euler or RK2)
Reinitialization to maintain signed distance property
The gradient \(|\nabla\phi|\) is computed using:
Terrain-Corrected Gradient
When a terrain file is present, the gradient is computed on the terrain surface rather than on the flat horizontal plane. The horizontal grid spacing \(\Delta x\) corresponds to a surface arc length of \(\Delta s_x = \Delta x \sqrt{1 + (\partial z/\partial x)^2}\), and similarly in \(y\). Dividing finite differences by these arc-length spacings gives the surface gradient components:
where \(s_x = \partial z / \partial x\) and
\(s_y = \partial z / \partial y\) are the terrain slope components stored
in the two-component slopes MultiFab.
The terrain-corrected gradient magnitude is then:
This correction is applied inside godunov_norm_grad_phi
(src/numerical_schemes.H) using effective spacings
\(\Delta x_\text{eff} = \Delta x \sqrt{1 + s_x^2}\) and
\(\Delta y_\text{eff} = \Delta y \sqrt{1 + s_y^2}\) in both the WENO3
and the first-order fallback stencils. Reinitialization uses the default
flat-terrain path (slope = 0) because the signed-distance property is defined
in the horizontal plane.
For steep terrain (e.g., slope magnitude 1.2, equivalent to ~50°) the effective spacing is 57 % larger than the flat value, so the surface gradient magnitude is correspondingly smaller. Without this correction the gradient would be over-estimated on steep flanks, causing the fire to propagate too slowly upslope.
FARSITE Elliptical Expansion Model (Richards 1990)
The FARSITE model uses an elliptical expansion pattern to represent the anisotropic spread of fire under the influence of wind and terrain.
Richards’ Coefficients
The elliptical fire spread is characterized by three coefficients:
\(a\) - head fire coefficient (maximum spread, downwind)
\(b\) - flank fire coefficient (crosswind spread)
\(c\) - backing fire coefficient (minimum spread, upwind)
These coefficients relate the directional rate of spread to the base rate:
where \(\theta\) is the angle between the fire normal and the wind direction.
Anderson Length-to-Width Ratio
The Anderson (1983) model relates the ellipse shape to wind speed:
where \(L/W\) is the length-to-width ratio and \(U\) is the wind speed (mph). This ratio is then converted to the Richards’ coefficients for the elliptical spread calculation.
Alternative Fire Shape Models
Four fire shape models are available for the FARSITE propagation path. All share the same three user coefficients \(a\) (head), \(b\) (flank), and \(c\) (backing), and the base rate of spread \(R_{head}\) from the active fire spread model (Rothermel, Balbi, etc.). Select the shape with:
farsite.fire_shape_model = richards # default
# or: catchpole_demestre | wilson | lemniscate
Richards (1990) – Linear Blend (default)
The default model is a simple linear combination of head, flank, and backing spread distances:
Enabled by farsite.fire_shape_model = richards (default).
Catchpole & de Mestre (1986) – True Double-Ellipse
Each Huygens wavelet is composed of two half-ellipses joined at the crosswind axis. The head half-ellipse has semi-axes \(a\) and \(b\); the backing half has semi-axes \(c\) and \(b\).
The polar spread distance from the fire-front origin is:
This model is geometrically exact for the double-ellipse shape and gives:
\(r(0) = a\,R_{head}\,\Delta t\) (head fire rate)
\(r(\pi/2) = b\,R_{head}\,\Delta t\) (flank rate)
\(r(\pi) = c\,R_{head}\,\Delta t\) (backing rate)
Enabled by farsite.fire_shape_model = catchpole_demestre.
Reference: Catchpole, E.A. & de Mestre, N.J. (1986). Australian Forestry, 49(2), 102–111.
Wilson (1988) – Single Ellipse from Rear Focus
Wilson’s model places the fire origin at the rear focus of a single ellipse. The ellipse geometry is derived from the head (\(a\)) and backing (\(c\)) coefficients only; the flank spread follows naturally.
Defining:
the ellipse parameters are:
The conic-section polar equation from the rear focus gives:
At the principal angles:
\(r(0) = R_H = a\,R_{head}\,\Delta t\) (head fire)
\(r(\pi) = R_B = c\,R_{head}\,\Delta t\) (backing fire)
\(r(\pi/2) = \ell = \dfrac{2R_H R_B}{R_H + R_B}\) (harmonic mean — flank)
The \(b\) coefficient is not used by this model; the flank spread is determined solely by \(a\) and \(c\).
Enabled by farsite.fire_shape_model = wilson.
Reference: Wilson, A.A.G. (1988). Canadian Journal of Forest Research, 18(6), 682–687.
Alexander et al. – Lemniscate (Limaçon)
The lemniscate model represents the fire wavelet as a Limaçon (cardioid family):
Principal values:
\(r(0) = \left(b + \tfrac{a-c}{2}\right) R_{head}\,\Delta t\)
\(r(\pi/2) = b\,R_{head}\,\Delta t\) (flank rate equals \(b\) exactly)
\(r(\pi) = \left(b - \tfrac{a-c}{2}\right) R_{head}\,\Delta t\)
When the natural choice \(b = (a + c)/2\) is used, the formula reduces to:
giving \(r(0) = a\), \(r(\pi) = c\), and a flank rate that is the arithmetic mean of head and backing. Values below zero are clamped to zero.
Enabled by farsite.fire_shape_model = lemniscate.
Reference: Alexander, M.E., Stocks, B.J., & Lawson, B.D. (1991). Information Report NOR-X-310, Canadian Forest Service.
Minimum Travel Time (MTT) Propagation
The Minimum Travel Time method (propagation_method = mtt) pre-computes a scalar arrival-time field \(T(i,j,k)\) — the time at which fire first reaches cell \((i,j,k)\) — using a Dijkstra fast-marching sweep over the precomputed ROS field. The level-set field is then updated analytically at each simulation time \(t\):
giving:
\(\phi < 0\) — burned (arrival before current time)
\(\phi = 0\) — fire front (arriving now)
\(\phi > 0\) — unburned (fire not yet arrived)
Arrival Time Computation
Starting from all ignition cells (where initial \(\phi < 0\), \(T = 0\)), the arrival time propagates to neighbouring cells via:
where \(d_s\) is the cell-centre-to-cell-centre distance and \(R_{\text{interface}}\) is the harmonic mean of the two adjacent ROS values:
The harmonic mean is chosen to penalise paths through low-ROS cells and prevents numerical blow-up when one cell has very small ROS.
A standard min-heap priority queue (Dijkstra) guarantees that the globally minimum arrival time is always processed next, making the result exact for any non-negative, spatially-varying ROS field.
Advantages and Limitations
Advantages relative to the level-set and FARSITE methods:
The arrival-time field is monotonically increasing; no reinitialization is needed.
The phi update is trivially parallel (one subtraction per cell per step).
The fast-march is run only once, regardless of the number of time steps.
Limitations:
The ROS field is frozen at the start; time-varying wind, moisture, or diurnal models are not captured after the initial sweep.
The method is isotropic (spread depends only on local ROS magnitude, not direction). Wind direction effects are included through the directional ROS model but fire does not curve around obstacles post-computation.
Select MTT with:
propagation_method = mtt
It is compatible with all fire spread models (rothermel, balbi,
cheney_gould, cruz_crown).
Reference: Finney, M.A. (2002). “Fire growth using minimum travel time methods.” Canadian Journal of Forest Research, 32(8), 1420–1424. https://doi.org/10.1139/x02-068
Crown Models
Van Wagner (1977) Crown Fire Initiation
Crown fire initiation occurs when sufficient heat is available to raise the canopy fuel to ignition temperature.
Crown Fire Initiation Criterion
The critical surface fire intensity \(I_0\) (kW/m) required for crown fire initiation is:
where:
\(CBH\) is the canopy base height (m)
\(M_c\) is the canopy moisture content (%)
\(h\) is the crown fuel consumption depth (m)
If the actual surface fire intensity \(I > I_0\), then crown fire initiation occurs.
Active Crown Fire Criterion
For active crown fire spread, the critical wind speed \(U_{active}\) is:
where \(CBD\) is the canopy bulk density (kg/m³).
Rothermel (1991) Active Crown Fire ROS
When crown.use_rothermel1991_crown = 1, active-crown-fire cells use a simple
multiplicative relationship derived by Rothermel (1991):
This is applied only in cells where crown_activity == 2 (active crown fire),
and replaces (or is blended with) the surface ROS.
Parameters: crown.use_rothermel1991_crown = 1.
Reference: Rothermel, R.C. (1991). Predicting Behavior and Size of Crown Fires in the Northern Rocky Mountains. USDA Forest Service Research Paper INT-438.
Van Wagner (1977) Passive-Crown Blending
When crown.use_passive_blend = 1 the transition from surface to active crown
fire is smoothed via Van Wagner’s (1977) passive crown-fire factor:
The blended ROS is:
where \(I_B\) is the Byram fireline intensity (kW/m) and \(I_o\) is the Van Wagner crown fire initiation threshold (kW/m):
Parameters: crown.use_passive_blend = 1, crown.CBH, crown.FMC.
Reference: Van Wagner, C.E. (1977). Conditions for the start and spread of crown fire. Canadian Journal of Forest Research, 7(1), 23–34.
Scott & Reinhardt (2001) Crown Fire Surface Area (CFSA)
The Scott & Reinhardt (2001) Crown Fire Surface Area (CFSA) model computes the effective burning surface area per unit ground area during active crown fire. This accounts for the 3-D structure of the canopy layer.
Model equation:
where:
\(\mathrm{CBD}\) is the canopy bulk density [kg/m³]
\(\mathrm{CH}\) is the canopy height [m]
\(\mathrm{CBH}\) is the canopy base height [m]
\(k_{sa}\) is the surface area coefficient [m²/kg] (typical: 4-8 m²/kg)
\(\mathrm{CFSA}_{\max}\) is the maximum CFSA [dimensionless] (typical: 2-4)
The surface area coefficient \(k_{sa}\) relates fuel loading to effective burning surface area, accounting for foliage clumping, branch geometry, and crown architecture.
Applications: Crown fire heat release calculations, emissions modeling, radiation heat flux from active crown fires, plume rise and smoke production.
Parameters: crown.cfsa_k_sa (default: 6.0), crown.cfsa_max (default: 3.0).
Reference: Scott, J.H. & Reinhardt, E.D. (2001). Assessing Crown Fire Potential by Linking Models of Surface and Crown Fire Behavior. USDA Forest Service Research Paper RMRS-RP-29.
Scott & Reinhardt (2001) Bisection-Based TI and CI
The full Scott & Reinhardt (2001) Torching Index (TI) and Crowning Index (CI) are
defined as the minimum 20-ft (6.1 m) open-wind speed (km/h) at which passive torching or
active crowning can occur. These are output as the plotfile fields
torching_index_kmh and crowning_index_kmh.
Enable with scott_reinhardt_full.enable = 1.
Torching Index is found by bisection on the wind speed \(U\) that satisfies:
where \(I_B(U) = H_l \, w_n \, R(U) / 60\) is the Byram fireline intensity (kW/m) computed from the Rothermel surface ROS \(R(U)\) [m/min], net fuel load \(w_n\) [kg/m²], and heat of combustion \(H_l\) [kJ/kg].
Crowning Index is found by bisection on the wind speed \(U\) that satisfies:
Reference: Scott, J.H. & Reinhardt, E.D. (2001). Assessing Crown Fire Potential by Linking Models of Surface and Crown Fire Behavior. USDA Forest Service Research Paper RMRS-RP-29.
Spotting Models
Firebrand Spotting
Firebrand spotting generates new ignition points ahead of the main fire front. Two independent spotting models are available; they can be enabled separately or together.
Probability-Based Spotting (spotting.*)
A stochastic model driven by wind speed, fire intensity, and fuel moisture. For each fire-front cell a spotting probability is computed and a Bernoulli draw decides whether a firebrand is generated. If so, the landing distance is sampled from a log-normal or exponential distribution and the spot is placed downwind with optional lateral dispersion.
Albini (1983) Spotting with 2-D Trajectory (albini_spotting.*)
A physics-based spotting model that couples Albini’s thermal-plume lofting formula with a 2-D horizontal particle trajectory integrated through the simulation wind field.
Stage 1 – Lofting height (Albini 1983)
Byram’s fire line intensity is computed from the Rothermel rate-of-spread:
The maximum height reached by a firebrand above the fire front is:
where \(I_B\) is in kW/m. Cells with \(I_B < I_{B,\min}\) do not generate firebrands.
Stage 2 – 2-D horizontal trajectory
Each lofted firebrand descends at a constant terminal velocity \(v_t\) (m/s). The flight time is:
During the flight the horizontal position is integrated with forward Euler over \(n_{\text{traj}}\) sub-steps using bilinear interpolation of the cell-centred velocity field \((u, v)\):
where \(\Delta t = t_f / n_{\text{traj}}\).
The landing position is the final \((x, y)\) after the full trajectory.
A circular ignition zone of radius albini_spotting.spot_radius is then
imposed by setting \(\phi\) to a negative signed-distance value in that
neighbourhood.
Diagnostic output fields
Each plot file contains four Albini-specific scalar fields:
albini_Hz– lofting height \(H_z\) at fire-front source cells.albini_count– number of firebrands launched from each source cell at the last spotting step.albini_dist– maximum landing distance from each source cell.albini_active– flag (1) at cells that received a spot ignition.
Vorticity-Enhanced Spotting (optional)
When both the Albini spotting model (albini_spotting.enable = 1) and the
Weise & Biging fire whirl model (weise_biging.enable = 1) are active,
vorticity-enhanced spotting can be enabled via weise_biging.enhance_spotting = 1.
This couples the fire whirl characteristics to the Albini spotting distance,
capturing the enhanced lofting and transport caused by fire whirls.
The enhancement factor multiplies the base Albini spotting distance:
where:
\(d_{\text{base}}\) is the baseline Albini spotting distance
\(\alpha\) is the empirical enhancement coefficient (
weise_biging.alpha, default: 1.0)\(v_\theta\) is the maximum tangential velocity at the whirl core edge [m/s]
\(\Omega\) is the angular velocity [rad/s]
\(r_w\) is the whirl core radius [m]
\(U\) is the ambient wind speed [m/s]
The enhancement factor is clamped to the range [1.0, 5.0] to avoid extreme values. When fire whirl characteristics are negligible (low intensity or wind), the enhancement factor defaults to 1.0 (no enhancement).
Configuration: To enable, set all three flags:
weise_biging.enable = 1
weise_biging.enhance_spotting = 1
albini_spotting.enable = 1
weise_biging.alpha = 1.0 # Optional: default is 1.0
Flux-Based Ember Cascade Model (ember_cascade.*)
The flux-based ember cascade model complements the single-brand Albini trajectory approach by treating dense firebrand showers as a continuous landing-flux field rather than tracking individual particles. It is appropriate when fire intensity is high enough to generate large, simultaneous ember cascades from the convective plume — a regime where the Monte-Carlo single-brand treatment becomes computationally or physically inadequate (Sardoy et al. 2007).
Stage 1 – Plume lofting height (Albini 1983)
Same as the trajectory model:
where \(I_B\) [kW/m] is the Byram fireline intensity. Cells with \(I_B < I_{B,\min}\) produce no embers.
Stage 2 – Ember source flux
The volumetric firebrand production rate per unit fire-line length is parameterised as a power-law of intensity:
where k_flux, I_B_ref, and flux_exp (\(\alpha\)) are user
parameters.
Stage 3 – Wind-advected Gaussian landing kernel
Flight time and mean transport displacement:
Turbulent diffusion from the plume broadens the ember cloud:
The landing-flux density at cell \((x, y)\) from source \(s\) is:
where \(\Delta x_s = x - (x_s + \delta x_s)\). Only cells within \(n_{\sigma,\text{cut}} \times \sigma_s\) of the mean landing centre receive contributions (Gaussian truncation for efficiency). The total landing flux is the sum over all active source cells:
Stage 4 – Stochastic Poisson ignition
During each check interval of duration \(\Delta t_{\text{check}}\), each cell is treated as an independent Poisson receiver:
A uniform random draw on \([0,1]\) is compared against \(P_{\text{ign}}\);
cells that pass receive a spot-fire ignition of radius spot_radius.
N_min_density [embers/m²/s] controls the ignition sensitivity.
3-D wind transport
When ember_cascade.use_3d_wind = 1 (or require_3d_wind = 1), the
height-averaged horizontal wind from a
massconsistent_amr plotfile
is used for \((u, v)\) in the transport displacement. Setting
require_3d_wind = 1 causes the simulation to abort with a diagnostic message
if no valid 3-D wind file is provided.
Diagnostic output fields
ember_cascade_flux– accumulated ember landing-flux density \(F_{\text{land}}\) [embers/m²/s] at each cell.ember_cascade_ignition– flag (1.0) at cells that received a spot-fire ignition during the most recent check interval.
Input parameters (prefix ember_cascade.):
Parameter |
Default |
Description |
|---|---|---|
|
0 |
1 to activate flux-based ember cascade model |
|
10.0 |
Minimum Byram fireline intensity [kW/m] to emit embers |
|
1.0 |
Firebrand terminal descent velocity [m/s] |
|
1.0 |
Ember source flux coefficient |
|
100.0 |
Reference intensity [kW/m] for flux scaling |
|
1.0 |
Intensity exponent \(\alpha\) |
|
50.0 |
Minimum Gaussian spread radius [m] |
|
0.1 |
Spread growth per metre of plume height [m/m] |
|
4.0 |
Gaussian kernel truncation (multiples of \(\sigma\)) |
|
1.0×10⁻³ |
Ignition threshold [embers/m²/s] |
|
5.0 |
New ignition zone radius [m] |
|
5 |
Run cascade check every N timesteps |
|
0 |
RNG seed (0 = system time) |
|
0 |
1 to use height-averaged wind from plt file |
|
0 |
1 to abort if no valid 3-D plt wind data provided |
|
(empty) |
Path to massconsistent_amr plt directory |
References: Albini, F.A. (1983). Potential spotting distance from wind-driven surface fires. USDA Forest Service Research Paper INT-309. Sardoy, N., Consalvi, J.L., Porterie, B. & Fernandez-Pello, A.C. (2007). Modeling transport and combustion of firebrands from burning trees. Combustion and Flame 150(3):151–169.
Spotting Suppression Inside Retardant Drop Zones
Inside active retardant drop zones the spotting probability is scaled by the same suppression factor as ROS:
This ensures that cells covered by retardant do not generate long-range spot fires, consistent with FARSITE’s suppression model.
Parameters: Existing retardant_file (no new parameters).
Model Enhancements
These modules augment the primary fire spread models. They are activated independently via ParmParse flags and are compatible with any of the fire spread models unless otherwise noted.
Barrier Polygons / Firebreaks
The barrier polygon feature (barrier_files) allows one or more CSV files of barrier vertex coordinates to be read at startup. At every time step, any grid cell whose centre falls nearest to a barrier vertex and that is currently burning (\(\phi < 0\)) is extinguished by setting:
This models firebreaks, fuel breaks, roads, and other physical barriers without requiring the AMReX Embedded Boundary (EB) framework.
Nearest-Cell Mapping
For each vertex \((X_v, Y_v)\) in a barrier CSV file, the nearest cell index is found by:
with the result clamped to the domain bounds. The set of barrier cells is de-duplicated once at startup and stored on the GPU device for efficiency.
Usage:
# Two CSV firebreak files (space-separated list)
barrier_files = road_break.csv ridge_break.csv
CSV format:
# X Y (one vertex per line; lines starting with # are ignored)
330100.0 3775300.0
330200.0 3775300.0
The barrier is applied after every propagation step regardless of the
propagation method (levelset, farsite, or mtt).
Andrews (2018) Wind Adjustments for Rothermel
Andrews (2018) documents two critical wind adjustments that improve the physical accuracy of the Rothermel surface fire spread model when wind input is from NWP or WRF (20-ft / 10-m height). Both are optional flags that can be enabled independently.
Multi-Layer Canopy Wind Profile
The multi-layer canopy wind profile model (Massman 1997, Inoue 1963) provides a physically realistic representation of wind speed variation with height through forest canopies. This improves accuracy for surface vs. crown fire transitions and layered fuel structures.
Within-canopy exponential profile (z < h):
where:
\(u_h\) is the wind speed at canopy top h [m/s]
\(\alpha\) is the attenuation coefficient [dimensionless]
\(z\) is height above ground [m]
\(h\) is canopy height [m]
Above-canopy logarithmic profile (z > h):
where \(d = 0.7h\) is the displacement height and \(z_0 = 0.1h\) is the roughness length.
The attenuation coefficient \(\alpha\) is related to leaf area index (LAI) via \(\alpha \approx 2.5 \times \mathrm{LAI}^{0.5}\) (Massman 1997).
Parameters: canopy_wind.enable, canopy_wind.n_layers (default: 5),
canopy_wind.alpha (default: 2.5), canopy_wind.LAI,
canopy_wind.d_ratio (default: 0.7), canopy_wind.z0_ratio (default: 0.1).
References: Massman, W.J. (1997). An analytical one-dimensional model of momentum transfer by vegetation of arbitrary structure. Boundary-Layer Meteorology, 83(3), 407-421. Inoue, E. (1963). On the turbulent structure of airflow within crop canopies. Journal of the Meteorological Society of Japan, 41, 317-326.
Wind Adjustment Factor (WAF)
The Rothermel model requires wind speed at midflame height, but field measurements and WRF/NWP output are typically at 20 ft (6.1 m) above open terrain. The WAF converts 20-ft open wind to midflame height.
Two formulas are available, selected via rothermel.waf_formula:
Andrews logarithmic formula (waf_formula = "andrews", default)
Derived from Albini & Baughman (1979). Used for all fuel types; the sheltered (closed-canopy) variant is applied automatically when canopy data are available.
Unsheltered (open / shrub fuel beds):
where \(h\) is the fuel bed depth [ft]. Typical values: 0.42 for FM4 (\(h = 2\) ft), 0.35 for FM9 (\(h = 3\) ft).
Sheltered (closed-canopy, FARSITE-style):
When canopy cover fraction \(f_c \geq 0.5\) and canopy height \(h_c > h\):
BehavePlus linear formula (waf_formula = "behaviorplus")
A simpler linear approximation used in BehavePlus for open and shrub fuel models:
where \(h_\text{in} = 12 \times h_\text{ft}\) is the fuel bed depth in inches. Typical values for common fuels:
Fuel |
h [ft] |
h [in] |
WAF (linear) |
|---|---|---|---|
FM1 |
1.0 |
12 |
0.408 |
FM2/FM3 |
2.5 |
30 |
0.480 |
FM4 |
6.0 |
72 |
0.648 |
FM6 |
2.5 |
30 |
0.480 |
FM9 |
0.2 |
2.4 |
0.370 |
For closed-canopy (sheltered) cells the BehavePlus path uses an exponential Beer–Lambert-style canopy attenuation:
where \(\alpha_c\) is the canopy attenuation coefficient
(rothermel.waf_canopy_alpha, default 1.5) and \(f_c\) is the canopy
cover fraction. At \(f_c = 1\) and \(\alpha_c = 1.5\) the canopy
reduces the midflame wind to about 22 % of the above-canopy value.
Enable via rothermel.use_waf = 1. When a landscape file is active and
provides canopy cover and canopy height data, WAF is computed per cell using
each fuel model’s depth and the local canopy structure.
Maximum Effective Wind Speed (MEWS)
Rothermel’s empirical wind factor \(\phi_w\) grows without bound, producing unrealistically large ROS at high wind speeds outside the model’s calibration range. The MEWS cap limits the effective wind factor by inverting Rothermel’s wind-factor formula at a maximum wind factor value equal to 90% of the reaction intensity (matching units: \(I_R\) in BTU/ft²/min, \(\phi_w\) dimensionless via the original Rothermel calibration):
Wind speed used in the ROS formula is capped at \(U_{\max}\) before
evaluating the wind factor. Enable via rothermel.use_wind_limit = 1.
Reference
Andrews, P.L. (2018). The Rothermel Surface Fire Spread Model and Associated Developments: A Comprehensive Explanation. Gen. Tech. Rep. RMRS-GTR-371. USDA Forest Service. https://doi.org/10.2737/RMRS-GTR-371
Albini, F.A. & Baughman, R.G. (1979). Estimating Windspeeds for Predicting Wildland Fire Behavior. USDA For. Serv. Res. Pap. INT-221.
Massman, W.J. (1987). A comparative study of some mathematical models of the mean wind structure and aerodynamic drag of plant canopies. Boundary-Layer Meteorology, 40(1–2), 179–197.
Heat Flux MultiFab and Fire-Induced Wind
A spatially-varying (or uniform) heat flux field \(Q\) [W/m²] represents the fire heat release rate at each cell. It drives two WindNinja-style wind corrections that are applied after any terrain-based modification.
Upward Convective Velocity
The fire plume drives a vertical velocity proportional to the buoyancy of the hot gas column above the fuel bed:
where:
\(g = 9.81\) m/s² – gravitational acceleration
\(Q\) – heat flux [W/m²] at the cell
\(H_{\text{ref}}\) – reference height [m] (
heat_flux.ref_height, default 10 m)\(\rho_{\text{air}}\) – air density [kg/m³] (
heat_flux.rho_air, default 1.2)\(C_p\) – specific heat of air [J/(kg·K)] (
heat_flux.Cp_air, default 1005)\(T_a\) – ambient temperature [K] (
heat_flux.T_a, default 300)\(k_{\uparrow}\) – upward velocity coefficient (
heat_flux.k_upward, default 1.0)
In a 3-D build, \(w_{\uparrow}\) is added directly to the vertical wind component. In a 2-D build the upward motion is projected onto the horizontal plane as an inflow opposite to the ambient wind direction.
Induced Horizontal Inflow
The fire-plume entrainment also creates a horizontal inflow toward the fire perimeter:
directed anti-parallel to the level-set gradient \(\nabla\phi\) (i.e., toward decreasing \(\phi\), which is the interior of the fire). This term is only applied outside the fire perimeter (\(\phi \ge 0\)).
Balbi Buoyancy Velocity Augmentation
For the Balbi spread model, the fire heat flux additionally augments the fuel-derived buoyancy velocity:
The combined \(v_{b,\text{eff}}\) is then used in the Balbi tilt-angle calculation. A larger buoyancy velocity makes the flame more vertical (smaller tilt angle), reducing radiant heat transfer to unburned fuel and thereby reducing ROS. This is physically consistent with observations of strongly buoyant plumes above high-intensity fires.
Input Parameters
Parameter |
Default |
Description |
|---|---|---|
|
0.0 |
Uniform heat flux [W/m²]. Set > 0 to activate. |
|
“” |
ASCII file (X Y Q columns) for spatially-varying heat flux (2D only). |
|
1.2 |
Air density [kg/m³] |
|
1005.0 |
Specific heat of air [J/(kg·K)] |
|
300.0 |
Ambient temperature [K] |
|
10.0 |
Reference height for buoyancy velocity [m] |
|
1.0 |
Upward velocity coefficient |
|
0.5 |
Induced horizontal inflow coefficient |
|
0 |
1 to enable upward velocity term |
|
0 |
1 to enable induced horizontal inflow term |
FMC Seasonal / Phenological Schedule
The foliar moisture content (FMC) used by the Van Wagner (1977) crown fire initiation threshold is updated each timestep from a seasonal schedule. Two modes:
CSV file (
fmc_schedule.file): two-column CSV of (day-of-year, fmc_pct).Built-in parametric curve (
fmc_schedule.use_farsite_curve = 1): piecewise linear phenology with configurable spring green-up and autumn curing breakpoints.
Default breakpoints match Southern California chaparral (Rothermel 1983):
Phase |
DOY range |
FMC |
|---|---|---|
Winter dormant |
1–89 |
85 % |
Green-up |
90–149 |
rises |
Summer peak |
150–239 |
140 % |
Curing |
240–300 |
falls |
Post-cure |
301–365 |
85 % |
Parameter |
Default |
Description |
|---|---|---|
|
0 |
1 to enable FMC seasonal schedule |
|
“” |
Path to two-column CSV (doy, fmc_pct) |
|
0 |
1 to use built-in parametric curve |
|
1 |
Day-of-year at simulation t = 0 |
|
90 |
DOY when green-up begins |
|
150 |
DOY when FMC reaches maximum |
|
240 |
DOY when curing begins |
|
300 |
DOY when curing ends |
|
85.0 |
Dormant / cured FMC [%] |
|
140.0 |
Peak green FMC [%] |
Precipitation-Driven Dead Fuel Moisture Wetting
When diurnal_moisture.enable = 1 and a rain rate is specified, dead fuel
moisture is driven by an exponential wetting / drying model:
Wetting (rain > threshold): \(M_\text{new} = M_\text{sat}(1 - e^{-\Delta t / \tau_\text{wet}}) + M_\text{old}\,e^{-\Delta t / \tau_\text{wet}}\)
Drying (no rain): \(M_\text{new} = M_\text{EMC} + (M_\text{old} - M_\text{EMC})\,e^{-\Delta t / \tau_\text{dry}}\)
Per size-class time constants (Rothermel 1983):
Class |
τ_wet |
τ_dry |
|---|---|---|
1-hr |
1 hr |
4 hr |
10-hr |
2 hr |
24 hr |
100-hr |
4 hr |
10 days |
Parameter |
Default |
Description |
|---|---|---|
|
0.0 |
Constant rain rate [mm/hr] |
|
“” |
CSV of (time_s, rain_mm_hr) for time-varying rain |
|
0.25 |
Minimum rain rate to trigger wetting [mm/hr] |
|
1.20 |
Saturation moisture content [fraction] |
Fine Fuel Moisture Time-Lag Differential Equations
The time-lag differential equation model (Nelson 2000, Viney 1991) provides physically-based fuel moisture dynamics with continuous evolution driven by ambient conditions.
Governing equation:
where:
\(M\) is fuel moisture content [fraction]
\(M_e\) is equilibrium moisture content (EMC) [fraction]
\(\tau_{\text{eff}}\) is the effective time-lag constant [hours]
\(P(t)\) is the precipitation wetting term [1/hour]
Equilibrium moisture content with hysteresis (Nelson 2000):
Adsorption (wetting):
Desorption (drying):
where \(H = \mathrm{RH}/100\) is relative humidity as a fraction.
Temperature correction:
where \(T\) is temperature [°C], \(T_{\text{ref}} = 20°\text{C}\), and \(k = 0.015\) °C⁻¹. Higher temperature accelerates drying.
Time-lag constants:
Class |
τ [hours] |
|---|---|
1-hr |
1 |
10-hr |
10 |
100-hr |
100 |
1000-hr |
1000 |
Parameters: fuel_moisture_de.enable, fuel_moisture_de.method = "timelag_de",
fuel_moisture_de.T_ambient, fuel_moisture_de.RH,
fuel_moisture_de.tau_1hr (default: 1.0), fuel_moisture_de.tau_10hr (default: 10.0),
fuel_moisture_de.use_hysteresis (default: 1), fuel_moisture_de.temp_correction (default: 1).
References: Nelson, R.M. (2000). Prediction of diurnal change in 10-h fuel stick moisture content. Canadian Journal of Forest Research, 30, 1071-1087. Viney, N.R. (1991). A review of fine fuel moisture modelling. International Journal of Wildland Fire, 1(4), 215-234.
Per-Cell Live Canopy Moisture from FMS File
A FARSITE fuel moisture scenario (.fms) file can be supplied via fms_file to
provide per-fuel-model dead and live moisture values. These are stored in a spatially-
varying MultiFab and override the global rothermel.M_d1, …, rothermel.M_lw on
a per-cell basis.
Requires: A landscape file (rothermel.landscape_file) to supply per-cell fuel
model codes.
Input parameter: fms_file = path/to/fire.fms
Per-Fuel-Model Burnout (Residence) Time
When a landscape file is active, the residence time for each fuel code is derived from the Rothermel (1983) formula:
where \(\alpha\) = 3600 s·ft⁻¹/(lb/ft³), \(\rho_p\) = 32 lb/ft³, and
\(\sigma_i\) = 1739 ft⁻¹. When no landscape file is active the global
farsite.tau_residence scalar is used (backward compatible).
Parameters: Automatic when rothermel.landscape_file is present.
Live Fuel Moisture Conditioning
During the pre-simulation conditioning period (conditioning.n_days > 0),
live fuel moistures (M_lh, M_lw) are ramped linearly from their initial
values toward equilibrium dead-fuel targets over the conditioning steps,
matching FARSITE’s behaviour:
M_lh_target = M_d100 × 1.5 (live herbaceous ≈ 150% of 100-hr dead)
M_lw_target = M_d100 × 2.0 (live woody ≈ 200% of 100-hr dead)
Parameters: Uses conditioning.n_days (existing parameter).
Burn-Period Controls (Diurnal Active-Spread Window)
When burn_period.enable = 1, the computed rate-of-spread field is zeroed
outside the specified local clock window [start_hour, end_hour),
preventing any level-set, FARSITE, or MTT advance during inactive hours.
All other processes (moisture evolution, spotting diagnostics, ecology metrics)
continue normally.
Midnight-crossing windows (e.g. start_hour = 22, end_hour = 6)
are handled correctly.
Parameter |
Description |
|---|---|
|
1 to enable burn-period gating (default: 0) |
|
Local hour (decimal) when fire becomes active (default: 10.0) |
|
Local hour (decimal) when fire becomes inactive (default: 20.0) |
|
Local clock hour at simulation t=0 (default: 0.0) |
FARSITE Topographic Horizon Shading
solar_radiation.H computes per-cell terrain shading from the local surface
normal (incidence-angle test) and canopy shading from cover fraction. Neither
method detects shadows cast by terrain tens or hundreds of cells away.
FARSITE addresses this with a full topographic horizon scan: for each cell,
march 8 compass directions to find the maximum elevation angle to any visible
terrain point, then shade the cell whenever the solar elevation is below that
angle. Enable via solar_radiation.use_topographic_horizon = 1.
Precomputation (one-time)
For each locally-owned cell \((i, j)\) and each of the 8 compass directions \(d\), the function ray-marches outward:
where \(\Delta_d\) is the physical step distance in direction \(d\) (cell spacing for axis-aligned directions, \(\sqrt{\Delta x^2 + \Delta y^2}\) for diagonals). The result is stored in an 8-component MultiFab.
Per-timestep shading check (GPU)
The GPU kernel interpolates the horizon angle for the current solar azimuth between the two bounding compass directions:
A cell is marked as fully shaded (\(f = 1\)) when \(\theta_s < \theta_\text{hz}\), where \(\theta_s\) is the solar elevation angle.
The total shade fraction applied to fuel temperature is:
f = max(f_local_terrain, f_canopy, f_topographic_horizon, cloud_cover)
Parameter |
Default |
Description |
|---|---|---|
|
|
Set to |
|
|
Maximum ray-march distance [m]. |
Warning
The topographic horizon scan is expensive (one-time O(N³) CPU sweep plus
MPI global elevation gather). Use horizon_scan_max_dist_m to bound cost
on large domains. Disabled by default.
Diagnostic Models
The following models can be run alongside any firespread model. They do not alter fire-front propagation but write additional diagnostic fields to every plotfile.
Byram Fire Behavior Diagnostics
Fireline Intensity
Byram’s (1959) fire line intensity \(I_B\) quantifies the rate of heat release per unit length of fire front:
where \(H\) is the low heat of combustion, \(w_a\) is the available (net) fuel load per unit area, and \(R\) is the rate of spread. The available fuel load accounts for the mineral content correction:
(where \(w_0\) is in lb/ft² and 4.8824 is the lb/ft² to kg/m² conversion).
Flame Length
Byram’s (1959) empirical relationship between fireline intensity and flame length:
These fields (fireline_intensity and flame_length) are computed at every time
step and written to each plotfile.
Fire Intensity Classification
The fireline intensity \(I_B\) is classified into six discrete classes following the FARSITE / NFFL convention:
Class |
Intensity (kW/m) |
Flame Length (m) |
Description |
|---|---|---|---|
I |
< 10 |
< 0.5 |
Creeping / smouldering |
II |
10-175 |
0.5-1.5 |
Low intensity surface fire |
III |
175-500 |
1.5-2.5 |
Moderate intensity |
IV |
500-2000 |
2.5-5 |
High intensity; spotting begins |
V |
2000-10000 |
5-11 |
Very high; active crown fire |
VI |
> 10000 |
> 11 |
Extreme; firestorm |
The classification is automatically computed and written to every plotfile as fire_intensity_class.
To selectively output only this field:
plot_vars = fire_intensity_class fireline_intensity
Implementation: src/fire_intensity_class.H provides GPU-compatible classification functions.
Byram Convective Number
Byram’s (1959) convective number \(N_c\) classifies fire behavior as wind-driven vs. plume-dominated:
where:
\(g = 9.81\) m/s² — gravitational acceleration
\(I_B\) [kW/m] — Byram fireline intensity
\(\tau_{\text{res}}\) [s] — flame residence time
\(\rho_{\text{air}} = 1.2\) kg/m³ — air density
\(C_p = 1005\) J/(kg·K) — specific heat of air
\(T_a = 300\) K — ambient temperature
\(U\) [m/s] — wind speed magnitude
When \(N_c > 1\), fire is plume-dominated (strong vertical convection, increased spotting). When \(N_c < 1\), fire is wind-driven (horizontal spread dominates).
Output variable: convective_number
Flame Tilt Angle
Flame tilt angle from vertical due to wind and buoyancy forces:
where \(v_b = \sqrt{(g \times I_B) / (\rho_{\text{air}} \times C_p \times T_a)}\) is the buoyancy-driven vertical velocity.
Output in degrees from vertical (0° = vertical flame, 90° = horizontal).
Output variable: flame_tilt
Packing Ratio
Rothermel (1972) fuel bed packing ratio and relative packing ratio:
where \(w_0\) is fuel loading [lb/ft²], \(\delta\) is fuel bed depth [ft], \(\rho_p\) is particle density [lb/ft³], and \(\sigma\) is the surface-area-to-volume ratio [ft⁻¹]. Values of \(\beta/\beta_{\text{opt}}\) near 1.0 indicate near-optimal fuel bed compaction for fire spread.
Output variables: packing_ratio, relative_packing
Flame Front Depth
Spatial thickness of the actively burning zone:
where \(R\) [m/s] is the rate of spread and \(\tau_{\text{res}}\) [s] is the flame residence time. This diagnostic is useful for heat transfer calculations, backing fire behavior, and fire line construction planning.
Output variable: flame_depth
McArthur Forest Fire Danger Index
The McArthur Mark 5 Forest Fire Danger Index (FFDI) is the Australian standard for forest fire danger assessment (McArthur 1967, Noble et al. 1980):
where:
\(DF\) — Drought Factor (0–10, default: 5.0 for moderate drought)
\(H\) [%] — Relative humidity (clamped to 1–100%)
\(T\) [°C] — Air temperature
\(V\) [km/h] — Wind speed
FFDI interpretation:
0–5: Low–Moderate fire danger
5–12: High fire danger
12–24: Very High fire danger
24–50: Severe fire danger
50–100: Extreme fire danger
100+: Catastrophic (Code Red)
Output variable: mcarthur_ffdi
References: Byram (1959), McArthur (1967), Rothermel (1972), Noble et al. (1980)
Advanced Physics Features
The following advanced physics models enhance the realism of fire behavior simulations.
Radiation-Driven Preheating Distance
The distance ahead of the fire front where fuel is preheated by radiant energy from flames:
where \(L_f\) is the Byram flame length [m], \(\theta_{\text{tilt}}\) is the flame tilt angle from vertical [radians], and \(F_v\) is the radiative view factor (typically 0.3-0.5).
The flame tilt angle is computed from wind and buoyancy:
where \(U_{\text{eff}}\) is the effective wind speed [m/s] and \(v_{\text{buoy}} = \sqrt{g \times L_f}\) is the buoyant rise velocity scale.
This affects ignition timing and spread rate in non-uniform fuels.
Reference: Butler et al. (2004), Weber (1991)
Implementation: src/radiation_preheating.H
Fuel Particle Temperature Evolution
Tracks the surface temperature of fuel particles using a lumped-capacitance energy balance:
where:
\(T\) is fuel particle temperature [K]
\(Q_{\text{rad}}\) is radiative heat flux from flames [W/m²]
\(Q_{\text{conv}}\) is convective heat flux from hot gases [W/m²]
\(Q_{\text{loss}} = \varepsilon \sigma (T^4 - T_a^4)\) is re-radiation loss [W/m²]
\(C_{\text{thermal}} = \rho_p c_p V_{\text{char}}\) is thermal capacitance [J/K/m²]
Ignition occurs when \(T \geq T_{\text{ignition}}\) (typically ~600 K for wood).
This determines actual ignition timing (not instantaneous) and affects spotting ignition probability.
Reference: Dupuy & Maréchal (2011), Pickett et al. (2010), Drysdale (2011)
Implementation: src/fuel_temperature.H
Fire Line Intensity Rate of Change
The temporal derivative of fireline intensity provides an indicator of fire acceleration/deceleration:
Interpretation:
\(dI/dt > 0\): Fire is intensifying (accelerating)
\(dI/dt = 0\): Fire is steady-state
\(dI/dt < 0\): Fire is weakening (decelerating)
Large positive values (e.g., > 50 kW/m/s) indicate dangerous “blow-up” conditions.
Useful for suppression resource allocation and identifying rapid fire growth.
Implementation: src/intensity_rate_of_change.H
Intensity-Dependent Flame Residence Time
Enhances the constant residence time model with intensity dependence:
where \(\tau_{\text{base}}\) is the baseline residence time [s] at reference intensity \(I_{\text{base}}\) [kW/m], \(I_{\text{byram}}\) is the current fireline intensity [kW/m], and \(\alpha\) is the intensity exponent (typically 0.3-0.5).
Higher intensity fires burn more vigorously and complete combustion faster.
Reference: Byram (1959), Rothermel (1972)
Implementation: src/intensity_residence_time.H
Wind-Fuel Interaction Feedback
Adjusts effective wind speed based on fuel canopy sheltering:
where \(U_{\text{ref}}\) is the reference wind speed [m/s], \(a\) is the extinction coefficient (typically 0.5-2.0), LAI is the leaf area index [dimensionless], \(z\) is the height of interest [m], and \(h\) is canopy height [m].
LAI is estimated from fuel bed properties:
where \(k_{\text{LAI}}\) is an empirical coefficient (typically 3-6), \(\delta\) is fuel bed depth [m], and \(\beta\) is packing ratio [dimensionless].
Dense fuels shelter the surface from wind, reducing effective wind at fuel bed level.
Reference: Wilson (1985), Albini & Baughman (1979)
Implementation: src/wind_fuel_interaction.H
Fuel Moisture of Extinction Gradient
Makes moisture of extinction (\(M_x\)) depend on fuel surface-area-to-volume ratio rather than using a constant 0.30:
where \(\sigma_{\text{fuel}}\) is the fuel SAV [ft⁻¹] and \(\sigma_{\text{ref}}\) is the reference SAV (typically 1739 ft⁻¹ for 1-hr fuels).
Typical values:
Grass (σ ~ 2000 ft⁻¹): \(M_x \approx 0.25\)
1-hr fuels (σ ~ 1739): \(M_x \approx 0.30\)
10-hr fuels (σ ~ 109): \(M_x \approx 0.40\)
More accurate fire/no-fire transitions with variable \(M_x\).
Reference: Rothermel (1972), Anderson (1970), Albini (1976)
Implementation: src/moisture_of_extinction.H
Flame Intermittency Factor
Accounts for pulsating/intermittent flames rather than steady burning:
where \(\gamma \in [0, 1]\) is the intermittency factor, \(I_{\text{byram}}\) is Byram fireline intensity [kW/m], and \(k\) is an empirical coefficient (typically 0.00325 [m/kW]).
The effective radiative heat flux is:
For low intensity fires: \(\gamma \to 0\) (highly intermittent)
For high intensity fires: \(\gamma \to 1\) (continuous flame)
Real flames flicker and pulse, affecting heat transfer efficiency and spotting (firebrands released in pulses).
Reference: Finney et al. (2015), Frankman et al. (2013)
Implementation: src/flame_intermittency.H
Plume Entrainment Momentum Feedback
Fire plumes create horizontal inflow winds that strengthen ROS:
where \(k_{\text{momentum}}\) is the momentum coupling coefficient (typically 0.1-0.3), \(w_{\text{plume}}\) is the upward plume velocity [m/s], \(I\) is fireline intensity [kW/m], and \(I_{\text{ref}}\) is reference intensity [kW/m] (typically 500-1000).
The upward plume velocity is estimated from buoyancy:
The total effective wind is:
Important for fire whirl formation and rapid fire growth.
Reference: Finney et al. (2015), Linn et al. (2002)
Implementation: src/plume_momentum_feedback.H
Spatially Varying Fuel Bed Bulk Density
Allows spatially varying fuel loading from landscape files:
where \(w_{0,\text{base}}\) is the base fuel loading [lb/ft²] and \(m(x, y)\) is a spatial multiplier field [dimensionless].
Real landscapes have patchy fuel loading due to terrain, moisture, past fires, etc. Higher local loading leads to higher intensity and potentially higher ROS.
The multiplier field can be loaded from:
LCP auxiliary data
Separate ASCII file with X Y multiplier columns
Generated randomly for ensemble runs
Reference: Finney (1998), Andrews (2018)
Implementation: src/fuel_loading_variation.H
Critical Heat Flux for Ignition
Replaces binary ignition (φ < 0) with threshold-based ignition from incident heat flux:
where \(q_{\text{base}}\) is the base critical heat flux for dry fuel [kW/m²] (typically 10-15), \(k_{\text{moisture}}\) is the moisture sensitivity coefficient (typically 2-4), and \(M_{\text{fuel}}\) is fuel moisture content [fraction].
Ignition criterion:
Typical values:
Dry fuel (M=0.05): \(q_{\text{crit}} \approx 11\) kW/m²
Moderate (M=0.15): \(q_{\text{crit}} \approx 16\) kW/m²
Wet fuel (M=0.30): \(q_{\text{crit}} \approx 25\) kW/m²
More physically realistic than geometric ignition; wet fuels require more heat to ignite.
Reference: Drysdale (2011), Quintiere (2006), Balbi et al. (2009), Lautenberger (2013)
Implementation: src/critical_heat_flux.H
Viegas (2004) Eruptive Fire Diagnostics
The Viegas (2004) model is an optional parallel diagnostic that characterises
eruptive (blow-up) fire behaviour on steep terrain. Unlike the primary spread
models it does not alter the fire front propagation; instead it writes five
diagnostic fields to every plotfile so that potential under-prediction by
Rothermel’s quadratic slope factor can be identified. Enable with
viegas.enable = 1.
Physical Basis
On steep slopes the Rothermel quadratic slope factor \(\phi_s\) significantly under-predicts the rate of spread compared to field and laboratory observations of eruptive fires (Viegas 2004). Viegas introduces an exponential slope enhancement factor that captures the runaway acceleration observed at critical slope angles.
Exponential Slope Enhancement Factor
where:
\(a_V\) is the Viegas slope coefficient (default 1.83, calibrated to laboratory fuels; dimensionless)
\(\varphi\) is the terrain slope angle; \(\tan\varphi = \sqrt{(\partial z/\partial x)^2 + (\partial z/\partial y)^2}\)
Viegas Rate of Spread
The Viegas ROS uses the Rothermel no-wind, no-slope rate \(R_0\) and the Rothermel wind factor \(\phi_w\) as a baseline and replaces only the slope factor:
Eruptive Regime Flag
A cell is flagged as being in the eruptive regime when the terrain slope exceeds a critical value \(\tan\varphi_c\) (Viegas 2004, Section 3):
The default critical slope is \(\tan\varphi_c = 0.4\) (approximately 22°). Flagged cells indicate locations where the Rothermel model is expected to under-predict the rate of spread.
ROS Excess Ratio
The signed excess ratio quantifies the relative difference between the two models:
Positive values indicate that Viegas predicts a higher spread rate than Rothermel; the larger \(\varepsilon\), the greater the potential under-prediction by Rothermel on steep terrain.
Flame-Tilt Angle (Hazard Assessment)
The combined wind–slope flame-tilt angle \(\alpha\) is computed for hazard assessment purposes only; it is not fed back into the rate-of-spread calculation:
where the buoyancy velocity scale \(v_b\) (m/s) is:
and:
\(U\) – wind speed magnitude (m/s)
\(g = 9.81\) m/s² – gravitational acceleration
\(\delta_m\) – fuel bed depth (m), converted from the Rothermel database value in ft
\(T_f\) – mean flame temperature (K); default 1000 K
\(T_a\) – ambient temperature (K); default 300 K
Merging Strategy
The following table summarises which couplings between the Viegas and Rothermel models are permitted in this implementation:
Coupling |
Status |
Rationale |
|---|---|---|
Rothermel \(R_0\) used as Viegas baseline |
✅ compatible |
Shared no-wind, no-slope ROS |
Same slope/wind/fuel inputs as Rothermel |
✅ compatible |
Consistent environmental state per cell |
Viegas flame-tilt for hazard assessment only |
✅ compatible |
Read-only diagnostic output |
Viegas induced wind fed into Rothermel |
❌ excluded |
Would invalidate Rothermel calibration |
Viegas eruptive acceleration fed into Rothermel ROS |
❌ excluded |
Incompatible ROS formulations |
Viegas critical slope replacing Rothermel \(\phi_s\) |
❌ excluded |
Different slope-factor definitions |
Viegas dynamic ROS inside Rothermel’s static formula |
❌ excluded |
Conceptually incompatible |
Diagnostic Output Fields
Five fields are written to every plotfile when viegas.enable = 1:
Field |
Units |
Description |
|---|---|---|
|
m/s |
\(R_0 (1+\phi_w)\exp(a_V \tan\varphi)\) |
|
– |
1 where \(\tan\varphi > \tan\varphi_c\) |
|
– |
\((R_V - R_{\text{Rothermel}})/R_{\text{Rothermel}}\) |
|
rad |
\(\arctan(U/v_b + \tan\varphi)\) (hazard only) |
|
– |
\(\exp(a_V \tan\varphi)\) |
Input Parameters
Parameter |
Default |
Description |
|---|---|---|
|
0 |
Enable (1) or disable (0) the Viegas diagnostics |
|
1.83 |
Exponential slope coefficient (dimensionless) |
|
0.4 |
Critical slope \(\tan\varphi_c\) for eruptive-regime flag (≈ 22°) |
|
300.0 K |
Ambient temperature for buoyancy velocity |
|
1000.0 K |
Mean flame temperature for buoyancy velocity |
Reference
Viegas, D.X. (2004). “Slope and wind effects on fire propagation.” International Journal of Wildland Fire, 13(2), 143–156. https://doi.org/10.1071/WF03046
Viegas-Balbi Coupled Diagnostic
When fire_spread_model = balbi, the Viegas diagnostic uses the Balbi
amplitude coefficient \(A\) and buoyancy velocity \(v_b\) instead of
the Rothermel \(R_0\) and wind factor \(\phi_w\).
The Viegas-Balbi rate of spread is:
where:
\(\alpha_w\) is the wind-only flame tilt angle: \(\tan\alpha_w = U / v_b\) (no terrain slope in the tilt angle — the slope enters only through \(\Phi_{s,V}\))
\(A\) is the Balbi amplitude coefficient [m/s] pre-computed from fuel properties (radiation term, moisture correction, and fuel geometry)
\(\Phi_{s,V} = \exp(a_V \tan\varphi)\) is the same Viegas exponential slope enhancement factor used in the Rothermel baseline
This formulation is physically consistent: Balbi’s radiation-driven ROS is modulated by the same exponential slope term that Viegas calibrated from laboratory and field eruptive fire experiments. The ROS excess ratio becomes
Weise & Biging (1996) Fire Whirl Model
The Weise & Biging (1996) fire whirl model is an optional diagnostic sub-model
that computes fire whirl characteristics from existing fire behavior outputs. It
does not modify the primary fire spread model; instead it runs alongside any
spread model to predict whirl formation when enabled via weise_biging.enable = 1.
Physical Basis
Fire whirls are columnar vortices generated by fire-induced updrafts interacting with ambient wind and terrain. The model characterises whirl geometry and kinematics from Byram fireline intensity and flame length.
Flame Tilt Angle
The Weise & Biging (1996) empirical flame tilt angle \(\theta\):
where \(Fr = U^2 / (g\,L_f)\) is the modified Froude number, \(U\) is wind speed [m/s], \(g = 9.81\) m/s², \(L_f\) is Byram flame length [m], and \(\beta\) is the terrain slope angle.
Vertical Flame Height
Fire Whirl Geometry (Rankine-Vortex Scaling)
Configuration via Parmparse (weise_biging.*)
Parameter |
Default |
Description |
|---|---|---|
|
0 |
Enable (1) or disable (0) the fire whirl model |
|
0.1 |
Whirl core radius-to-height ratio (dimensionless) |
|
1.0 kW/m |
Minimum Byram fireline intensity threshold |
|
0 |
Enable (1) vorticity-enhanced Albini spotting |
|
1.0 |
Spotting enhancement coefficient (dimensionless) |
Diagnostic Output Fields
Six fields are written to every plotfile when enabled:
Field |
Units |
Description |
|---|---|---|
|
m |
Vertical component of tilted flame |
|
rad |
Flame tilt angle from vertical |
|
m |
Fire whirl column height |
|
m |
Fire whirl core radius |
|
rad/s |
Whirl angular velocity |
|
m/s |
Maximum tangential velocity at core edge |
Reference
Weise, D.R. and Biging, G.S. (1996). “Effects of wind velocity and slope on flame properties.” Canadian Journal of Forest Research, 26(10), 1849–1858. https://doi.org/10.1139/x26-210
Bulk Fuel Consumption
The bulk fuel consumption fraction \(f_c\) represents the fraction of available fuel consumed behind the fire front:
where \(\tau_{res}\) is the residence time. The consumption affects the fire intensity and heat release rate.
The consumption fraction is bounded by configurable minimum and maximum values
(farsite.f_consumed_min and farsite.f_consumed_max):
Post-Processing Diagnostics
In addition to fireline intensity and flame length, the solver always writes four fire-ecology fields and three emissions fields to every plotfile. The equations below describe how each field is computed.
Scorch Height (Van Wagner 1973)
Van Wagner’s (1973) empirical relationship between fireline intensity and the height to which foliage is killed by convective heat:
where \(I_B\) [kW/m] is the Byram fireline intensity. Cells where \(I_B = 0\) receive
\(H_s = 0\). Output field: scorch_height.
Probability of Ignition (Anderson 1970)
The probability that a firebrand landing on fine fuel will produce sustained ignition (Anderson 1970, Rothermel 1983):
where:
\(T_{\text{fuel,F}} = T_{a,\text{F}} + \Delta T_{\text{solar}}\) — fine-fuel temperature in °F (ambient temperature converted from °C plus a solar-heating increment; default \(+25\) °F)
\(M_{\%}\) — 1-hr dead fuel moisture content [%]
Output field: prob_ignition.
Tree Mortality (Ryan & Reinhardt 1988)
A logistic mortality model based on the fraction of the live crown that is scorched.
Crown Scorch Ratio (CSR):
where \(H_{\text{CBH}}\) is the canopy base height [m] (crown.CBH) and
\(H_{\text{tree}}\) is the tree height [m] (fire_ecology.tree_height).
Mortality probability (only when \(H_s > H_{\text{CBH}}\)):
The logistic curve gives ≈50 % mortality when 60 % of the live crown is scorched.
Output field: tree_mortality.
Crown Activity Class (Van Wagner 1977)
The crown activity class categorises each cell as surface, passive crown, or active crown fire:
Class |
Label |
Criterion |
|---|---|---|
0 |
Surface fire |
\(I_B < I_o\) |
1 |
Passive crown |
\(I_B \ge I_o\) and \(R < R'_{SA}\) |
2 |
Active crown |
\(I_B \ge I_o\) and \(R \ge R'_{SA}\) |
Crown initiation intensity (Van Wagner 1977):
Critical active-crown ROS (Van Wagner 1977):
where \(C_{\text{BD}}\) [kg/m³] is the canopy bulk density (crown.CBD) and
\(F_{\text{MC}}\) [%] is the foliar moisture content (crown.FMC).
Output field: crown_activity.
Fire Emissions (Seiler & Crutzen 1980)
Per-cell cumulative emissions for three species, following the WRF-Fire convention:
where:
\(W_{\text{fuel}}\) — oven-dry fuel load [kg/m²] (converted from lb/ft²)
\(f_c\) — fuel consumption fraction from the bulk model, or
emissions.default_consumed_frac(default 0.9) for burned cells without the bulk model active\(\text{EF}\) — emission factors [kg species per kg dry fuel consumed]:
CO₂: 1.570 kg/kg (Seiler & Crutzen 1980)
CO: 0.102 kg/kg
PM2.5: 0.0162 kg/kg (Andreae & Merlet 2001)
Output fields: co2_emissions, co_emissions, pm25_emissions.
Smoke Plume Rise with Pasquill-Gifford Stability
Integrates Pasquill-Gifford atmospheric stability classification into the Briggs plume rise model, allowing FARSITE-compatible adjustment of smoke plume rise based on atmospheric conditions.
The base Briggs (1965/1969) formula computes:
With Pasquill-Gifford correction:
This accounts for:
Unstable conditions (A-C): Enhanced vertical mixing → greater plume rise
Stable conditions (E-F): Suppressed vertical motion → reduced plume rise
Stability Classes
Class |
Conditions |
Correction Factor |
Effect on Plume Rise |
|---|---|---|---|
A |
Extremely unstable (strong solar radiation, light winds) |
1.20 |
+20% enhanced rise |
B |
Moderately unstable |
1.10 |
+10% enhanced rise |
C |
Slightly unstable |
1.05 |
+5% enhanced rise |
D |
Neutral (overcast or moderate wind) |
1.00 |
No correction |
E |
Slightly stable (clear night, light winds) |
0.85 |
-15% suppressed rise |
F |
Moderately stable (clear night, very light winds) |
0.70 |
-30% suppressed rise |
Usage:
smoke_plume.enable = 1
smoke_plume.use_stability_correction = 1
smoke_plume.stability_class = A # Options: A, B, C, D, E, or F
Default (without correction) is equivalent to:
smoke_plume.use_stability_correction = 0 # or omit
smoke_plume.stability_class = D
Implementation: src/smoke_plume_rise.H
References: Briggs, G.A. (1965/1969). Plume rise models; Pasquill, F. (1961); Gifford, F.A. (1961); Turner, D.B. (1970). Workbook of atmospheric dispersion estimates.
Always-Present Output Fields
Every plotfile produced by the solver contains the following fields regardless of which models are enabled. Optional-model fields (Albini, Viegas, Weise & Biging, stochastic spotting) are described in their respective model sections.
Core simulation fields
Field |
Units |
Description |
|---|---|---|
|
m |
Level-set signed-distance function. Negative inside the fire perimeter, positive outside. The zero contour \(\phi = 0\) is the active fire front. |
|
m/s |
Wind velocity component in the x-direction at cell centres. |
|
m/s |
Wind velocity component in the y-direction at cell centres. |
Terrain and fuel fields
These fields carry the static landscape inputs and are constant throughout the simulation. They are written once per plotfile for easy co-location with fire behaviour outputs.
Field |
Units |
Description |
|---|---|---|
|
m |
Ground elevation above mean sea level, read from the terrain file or landscape file (LCP). Zero when neither is provided. |
|
– |
Terrain slope magnitude \(\sqrt{(\partial z/\partial x)^2 + (\partial z/\partial y)^2}\). Dimensionless (rise/run). |
|
rad |
Terrain aspect (downslope azimuth from north, clockwise). Used internally for anisotropic slope corrections. |
|
– |
Integer fuel model index for each cell. Values correspond to FBFM13 (1–13) or FBFM40 (91–204) identifiers, or the uniform model code when no landscape file is provided. |
Crown fire field
Field |
Units |
Description |
|---|---|---|
|
– |
Fraction of the total rate of spread contributed by crown fire
(\(R_{\text{crown}} / R_{\text{total}}\) when crown fire is initiated).
Zero for pure surface fire, one for pure crown fire. Computed only when
the crown fire initiation model is active ( |
Bulk fuel consumption field
Field |
Units |
Description |
|---|---|---|
|
– |
Fraction of oven-dry fuel consumed behind the fire front (\(f_c \in [0,1]\), see Bulk Fuel Consumption). Zero in unburned cells and in cells ahead of the front. |
Stochastic spotting fields
These four fields are written when the probability-based spotting model is active
(spotting.enable = 1; see Firebrand Spotting for the full model description).
Field |
Units |
Description |
|---|---|---|
|
– |
Probability that a firebrand generated at this cell produces a secondary ignition. Computed from wind speed, fuel moisture, and fire intensity (Bernoulli spot-probability model). |
|
– |
Number of spot ignitions generated from this source cell at the most recent spotting step. |
|
m |
Maximum spotting distance from this source cell at the most recent spotting step (landing distance of the most distant active firebrand). |
|
– |
Flag (1) at cells that received at least one spot ignition during the most recent spotting step; 0 elsewhere. |
Spatial Fuel Moisture Fields
These five fields are written when the Nelson (2000) diurnal EMC or FARSITE
.fms spatial scenario file is active; otherwise they contain the global
Rothermel moisture values uniformly.
Field |
Units |
Description |
|---|---|---|
|
– |
1-hr dead fuel moisture [fraction] |
|
– |
10-hr dead fuel moisture [fraction] |
|
– |
100-hr dead fuel moisture [fraction] |
|
– |
Live herbaceous fuel moisture [fraction] |
|
– |
Live woody fuel moisture [fraction] |
References
Van Wagner, C.E. (1973). “Height of crown scorch in forest fires.” Canadian Journal of Forest Research, 3(3), 373–378.
Anderson, H.E. (1970). “Forest fuel ignitibility.” Fire Technology, 6(4), 312–319.
Ryan, K.C. & Reinhardt, E.D. (1988). “Predicting postfire mortality of seven western conifers.” Canadian Journal of Forest Research, 18(10), 1291–1297.
Seiler, W. & Crutzen, P.J. (1980). “Estimates of gross and net fluxes of carbon between the biosphere and the atmosphere from biomass burning.” Climatic Change, 2(3), 207–247.
Andreae, M.O. & Merlet, P. (2001). “Emission of trace gases and aerosols from biomass burning.” Global Biogeochemical Cycles, 15(4), 955–966.
Wind Models
Several wind modelling approaches are available, from a simple constant uniform wind to spatially-varying fields and terrain-following feedback models. The selection is made through a combination of input parameters as described below.
Uniform Wind
The simplest wind specification sets constant velocity components over the entire domain for the whole simulation:
u_x = 2.0 # x-wind component [m/s]
u_y = 0.5 # y-wind component [m/s]
The default is u_x = 0.25, u_y = 0.0. This is the
fallback when no wind file or terrain-following solver is used.
Wind File — Steady (Spatially-Varying)
A spatially-varying but time-constant wind field can be supplied as a CSV file:
velocity_file = wind_data.csv
The file must contain one sample point per line with columns
X Y U V (in metres and m/s respectively). The solver interpolates
the wind onto the computational grid using inverse-distance weighting (IDW)
at initialisation. When velocity_file is set, it overrides
u_x / u_y / u_z.
CSV format:
# X [m] Y [m] U [m/s] V [m/s]
0 0 2.0 0.5
1000 0 2.1 0.4
0 1000 1.9 0.6
1000 1000 2.0 0.5
Wind File — Unsteady (Time-Dependent)
When use_time_dependent_wind = 1, the solver loads a sequence of wind
snapshot files and performs temporal linear interpolation between them. This
is only available in 2-D builds.
velocity_file = wind_t0.csv # base snapshot (t = 0)
use_time_dependent_wind = 1
wind_time_spacing = 3600.0 # seconds between snapshots
Snapshot files follow the naming convention:
<velocity_file>— time index 0<base>_1.csv— time index 1<base>_2.csv— time index 2
At each simulation step the solver:
Determines which two snapshots bracket the current time.
Loads both files (IDW spatial interpolation).
Performs temporal linear interpolation.
For example, with wind_time_spacing = 3600.0 s:
At t = 0 s : uses snapshot 0
At t = 1800 s : 50 % blend of snapshots 0 and 1
At t = 3600 s : uses snapshot 1
Turbulent Wind Perturbation (turb_wind.*)
Adds stochastic velocity perturbations to the ambient wind field at every timestep. The unperturbed base wind is preserved in a separate MultiFab; the perturbed field is always computed as base + perturbation so the perturbation never drifts away from the intended background.
Two models are available.
Ornstein-Uhlenbeck Process (turb_wind.model = ou_process)
The OU process produces temporally correlated wind gusts whose auto-correlation decays exponentially in time (Uhlenbeck & Ornstein 1930).
Discrete-time OU update (exact Euler–Maruyama)
where:
\(\theta\) [s⁻¹] is the reversion rate (
turb_wind.theta); decorrelation time = \(1/\theta\).\(\sigma\) [m/s] is the stationary standard deviation of each cell’s perturbation (
turb_wind.sigma).\(\eta_u, \eta_v\) are unit-variance noise terms (white when
L_c = 0, spatially correlated whenL_c > 0; see below).
The effective wind used for all ROS computations is
Domain-uniform mode (L_c = 0)
A single pair \((u', v')\) is drawn and added to every cell. This corresponds to a gust event coherent over the entire domain — appropriate when the domain is small relative to the turbulent integral scale.
Spatially correlated mode (L_c > 0)
Per-cell OU states \(u'(x,y)\) and \(v'(x,y)\) are maintained in a MultiFab. The noise field \(\eta(x,y)\) is constructed as:
Draw independent \(\mathcal{N}(0,1)\) white noise \(\xi(x,y)\) at every cell.
Apply a separable 1-D Gaussian convolution in \(x\) then \(y\) with kernel standard deviation
\[\sigma_k = \frac{L_c}{\Delta x} \quad\text{[cells]}\]using Neumann (nearest-cell) clamping at domain and subdomain boundaries.
Rescale the smoothed field to unit RMS so the per-cell stationary standard deviation remains exactly \(\sigma\).
The resulting noise field has the approximate 2-D autocovariance
i.e.a Gaussian envelope with e-folding length \(L_c\) (turb_wind.L_c).
Random Fourier Feature Spectral Noise (turb_wind.model = spectral_noise)
The spectral noise model combines a physically correct spatial power spectrum with OU temporal evolution, following the Random Fourier Feature (RFF) methodology (Kraichnan 1970; Rahimi & Recht 2007).
Initialisation – wavenumber and phase sampling
At startup \(N\) wavenumber pairs and phases are drawn once and remain fixed for the entire simulation:
where \(L_c\) is the user-specified spatial correlation length
(turb_wind.L_c [m]) and \(N\) is the number of modes
(turb_wind.N_modes). Sampling \(k_x, k_y\) independently from
\(\mathcal{N}(0, 1/L_c^2)\) corresponds to a 2-D isotropic Gaussian
power spectral density \(S(k) \propto \exp(-k^2 L_c^2/2)\), whose
Fourier transform is the Gaussian autocorrelation function.
Temporal update – OU amplitude evolution (CPU, every timestep)
Each mode carries a pair of scalar OU amplitude coefficients \((A^u_n, A^v_n)\) with unit stationary variance:
where \(\xi^u_n, \xi^v_n \sim \mathcal{N}(0,1)\) are independent each step.
Field reconstruction (GPU, every timestep)
The perturbation field is reconstructed on the GPU as a weighted cosine superposition:
The scale factor \(\sigma\sqrt{2/N}\) ensures that the stationary variance of each cell’s perturbation equals \(\sigma^2\) [m²/s²]:
Properties
Temporal decorrelation time = \(1/\theta\) [s] (same as
ou_process).Approximate 2-D spatial autocovariance (by the Bochner theorem):
\[C(r) \approx \sigma^2 \exp\!\left(-\frac{r^2}{2 L_c^2}\right).\]Energy is distributed across all resolved wavenumbers according to the Gaussian PSD — unlike the grid-smoothing approach, no artificial truncation of the spectrum occurs.
GPU-accelerated: the inner-mode loop in the reconstruction kernel runs on the device; CPU workload per step is only \(O(N)\) for the OU updates.
References
Kraichnan, R.H. (1970). Diffusion by a random velocity field. Phys. Fluids 13(1):22–31.
Rahimi, A. & Recht, B. (2007). Random features for large-scale kernel machines. NIPS 2007.
Bounded Direction Random Walk (turb_wind.model = direction_walk)
At each timestep a Gaussian angular increment is drawn and accumulated:
The entire base wind field is then rotated by the cumulative angle \(\vartheta\):
Wind speed is preserved exactly; only direction fluctuates. Parameters:
turb_wind.sigma_theta [rad/step] and turb_wind.theta_max [rad].
References
Uhlenbeck, G.E. & Ornstein, L.S. (1930). On the theory of the Brownian motion. Phys. Rev. 36(5):823–841.
Finney, M.A. et al. (2011). Role of wind, fuel moisture, and terrain in controlling fire movement. Ecosphere 2(3):art17.
Cruz, M.G. et al. (2019). Uncertainty in wildfire behaviour research. Current Forestry Reports 5:155–172.
Wind-Terrain Feedback Models (Options 1–6)
The wind_terrain.model ParmParse key selects how terrain-induced or
fire-feedback winds modify the effective wind seen by the fire spread model.
Option 1 (``none``) is the default and preserves all existing behaviour.
Options 2–6 progressively couple Viegas or terrain-channel wind physics into
the actual fire propagation. Option 7 (windninja_ridge_canyon) is described
in the following section.
Options Summary
Option |
|
Effect on fire spread model |
|---|---|---|
1 |
|
No modification — default behaviour |
2 |
|
\(R_{\text{final}} = \max(R_\text{Rothermel}, R_\text{Viegas})\) in eruptive cells (\(\tan\varphi > \tan\varphi_c\)) |
3 |
|
Add buoyancy-induced upslope wind \(U_\text{ind} = v_b\tan\varphi\) in eruptive cells only |
4 |
|
Scale ambient wind by \((1 + k_\text{canyon}\tan\varphi)\) everywhere (Rothermel 1983) |
5 |
|
Add upslope buoyancy wind \(U_\text{ind} = v_b\tan\varphi\) at every cell (no threshold) |
6 |
|
Scale ambient wind by \(\exp(k_\text{pimont}\tan\varphi)\) everywhere (Pimont et al. 2009) |
9 |
|
Viegas’ (2004) canyon eruptive acceleration model by dynamically coupling front orientation to a local terrain-slope induced draft velocity field |
Selecting viegas_ros, viegas_wind, viegas_neto, or canyon_eruptive automatically enables
the Viegas (2004) diagnostic model (viegas.enable = 1). Terrain slopes must be
provided (via rothermel.terrain_file or rothermel.landscape_file) for any
slope-dependent option to have an effect.
Physical Equations
Option 2 — Viegas ROS override (requires terrain slopes):
Options 3 and 5 — Buoyancy-induced upslope wind:
where \(\hat{\mathbf{n}}_s = \nabla z / |\nabla z|\) is the unit upslope
vector, \(\delta_m\) is fuel bed depth [m], and \(T_a\), \(T_f\)
are taken from viegas.T_a / viegas.T_f. Option 3 applies
\(\Delta\mathbf{U}\) only where \(\tan\varphi > \tan\varphi_c\);
Option 5 applies it everywhere.
Option 4 — Canyon wind (Rothermel 1983):
Option 6 — Pimont exponential correction (Pimont et al. 2009):
Option 9 — Canyon eruptive acceleration model (Viegas 2004):
Couples the fire front orientation vector \(\hat{\mathbf{n}} = \nabla\phi / |\nabla\phi|\) to the terrain-slope induced draft velocity field:
where \(\hat{\mathbf{u}}_\text{uphill} = \nabla z / |\nabla z|\) is the uphill unit vector, and \(v_b\) is the buoyancy velocity scale.
Configuration via Parmparse (wind_terrain.*)
Parameter |
Default |
Description |
|---|---|---|
|
|
Model: |
|
1.0 |
Terrain channeling coefficient for |
|
0.5 |
Exponential slope coefficient for |
Viegas parameters (viegas.a_V, viegas.tan_phi_c, viegas.T_a,
viegas.T_f) are shared between the Viegas diagnostic model and the
buoyancy-wind options (2, 3, 5); see
Viegas (2004) Eruptive Fire Diagnostics above.
References
Rothermel, R.C. (1983). How to Predict the Spread and Intensity of Forest and Range Fires. USDA Forest Service General Technical Report INT-143.
Viegas, D.X. & Neto, L.P.S. (1994). “Wind tunnel study of the convective air flow of slope fires.” Annual Report, Project COMOESTAS.
Pimont, F., Dupuy, J.-L., Linn, R.R. & Dupont, S. (2009). “Validation of FIRETEC wind-flows over a canopy and fuel-break.” International Journal of Wildland Fire, 18(7), 775–790. https://doi.org/10.1071/WF07130
WindNinja Ridge/Canyon Terrain Speed-up (Option 7)
The WindNinja empirical ridge/canyon model (Option 7,
wind_terrain.model = windninja_ridge_canyon) accounts for terrain-driven
wind acceleration observed in ridge and canyon topography (Forthofer 2007).
Wind-Slope Alignment
The alignment between the ambient wind and the upslope direction is:
where \(\hat{\mathbf{n}}_s = \nabla z / |\nabla z|\) is the unit upslope vector and \(\mathbf{U}\) is the horizontal wind velocity. The alignment satisfies \(a \in [-1,\,1]\):
\(a > 0\): wind has an upslope component → ridge acceleration
\(a < 0\): wind has a downslope component → canyon channeling
\(a = 0\): wind is perpendicular to slope, no modification
Ridge Speed-up
When the wind climbs a slope, terrain convergence amplifies the wind speed:
Canyon Channeling
When the wind flows down a slope (drainage flow channeling):
In both cases the effective wind is \(\mathbf{U}_{\text{eff}} = f \,\mathbf{U}\) with \(f \ge 1\) (the model never suppresses wind speed).
Parameters:
\(k_{\text{ridge}}\) – ridge speed-up coefficient (default 1.0)
\(k_{\text{canyon,WN}}\) – canyon channeling coefficient (default 0.5)
Reference: Forthofer, J.M. (2007). Modeling Wind in Complex Terrain for Use in Fire Spread Prediction. Colorado State University MS thesis.
FARSITE Wind Stream Simulation (Option 8)
The FARSITE wind stream simulation model (Option 8,
wind_terrain.model = farsite_wind) implements FARSITE’s internal simplified
empirical wind flow model around terrain. Unlike full CFD approaches, this model
uses terrain curvature analysis to identify ridge and valley features, then
applies empirical wind speed adjustments and direction deflections based on
field observations and wind tunnel experiments (Finney 1998).
Terrain Curvature Analysis
The model computes terrain curvature (Laplacian of elevation) to classify topographic features:
where:
\(\kappa > 0\): ridge (convex terrain)
\(\kappa < 0\): valley (concave terrain)
\(|\kappa| < \kappa_{\text{min}}\): flat terrain, no modification
Wind Speed Modifications
Ridge acceleration (positive curvature, wind upslope):
where \(a = (\mathbf{U} \cdot \hat{\mathbf{n}}_s) / |\mathbf{U}|\) is the wind-slope alignment. Maximum acceleration occurs when wind blows directly upslope (\(a = 1\)).
Lee-side sheltering (negative curvature, wind downslope):
Applies in protected areas downwind of ridges when wind has a downslope component. Floor at \(f_{\text{shelter}} \ge 0.2\) prevents zero wind.
Valley channeling (negative curvature, cross-wind or upslope):
Valley funneling is strongest when wind flows perpendicular to the valley axis (\(a \approx 0\)).
Wind Direction Deflection
Wind direction is deflected to follow terrain contours:
where \(\alpha_s\) is the terrain aspect (direction of steepest ascent), \(\theta_{\text{wind}}\) is the wind direction, and \(\varphi\) is the slope angle. Maximum deflection occurs when wind crosses terrain contours at 90°. Deflection is limited to ±45° per cell to preserve numerical stability.
Effective Wind
The final effective wind incorporates both speed and direction modifications:
Parameters
Parameter |
Default |
Description |
|---|---|---|
|
|
Set to |
|
1.5 |
Ridge speed-up coefficient |
|
0.6 |
Lee-side sheltering coefficient |
|
0.8 |
Valley channeling coefficient |
|
0.3 |
Wind direction deflection coefficient |
|
0.0001 |
Minimum curvature threshold [m⁻¹] |
Example Configuration
# FARSITE wind stream simulation on complex terrain
wind_terrain.model = farsite_wind
wind_terrain.k_ridge_farsite = 1.5
wind_terrain.k_shelter = 0.6
wind_terrain.k_valley = 0.8
wind_terrain.k_deflection = 0.3
Reference: Finney, M.A. (1998). FARSITE: Fire Area Simulator—Model Development and Evaluation. USDA Forest Service Research Paper RMRS-RP-4.
Compact Wind Direction Schedule
A compact three-column CSV schedule (wind_dir_schedule_file) updates the
uniform wind field at every timestep without requiring a full spatial wind grid
per snapshot:
wind_dir_schedule_file = wind_schedule.csv
CSV format: time_s, speed_ms, dir_deg where dir_deg is the direction
from which the wind blows (270° = westerly wind blowing eastward). The
schedule is linearly interpolated with circular direction averaging to avoid
0°/360° wrap-around artefacts.
When this file is set it overrides u_x / u_y and the
use_time_dependent_wind file-grid path. Compatible with turbulent wind
perturbation (perturbations are added on top of the scheduled base wind).
Input parameter: wind_dir_schedule_file = path/to/wind_schedule.csv
Multiple Weather Stations with Spatial IDW Interpolation
A new multi_wtr_file input enables loading multiple FARSITE .wtr weather
station files and producing spatially-varying wind, temperature, and relative
humidity via inverse-distance-weighting (IDW). At each timestep, wind U/V
components, temperature, and relative humidity are interpolated to every grid
cell, creating spatially-varying fuel moisture content when used with the
diurnal moisture model.
Station list CSV format (multi_wtr_file):
# station_id, x_m, y_m, wtr_file
1, 330000.0, 3775000.0, station1.wtr
2, 335000.0, 3775000.0, station2.wtr
3, 332500.0, 3780000.0, station3.wtr
At each timestep, the following fields are IDW-interpolated to every grid cell:
where \(d_i\) is the distance from the cell to station \(i\), \(p\) is the IDW power exponent (default 2.0), and \(V_i\) is the value at station \(i\) (wind component, temperature, or RH).
Interpolated Fields:
Wind U and V components (Cartesian interpolation)
Temperature \(T\) [°C] (scalar interpolation)
Relative Humidity RH [%] (scalar interpolation)
When diurnal_moisture.enable = 1 is set, the spatially-interpolated T and
RH create per-cell equilibrium moisture content (EMC) via the Nelson (2000)
model, resulting in realistic spatial moisture gradients.
Parameter |
Description |
|---|---|
|
Path to station list CSV (default: |
|
IDW exponent \(p\) (default: 2.0) |
Test: regtest/moisture/multi_wtr_spatial/
Spot Fire Ignition Delay
Models the delay between firebrand landing and actual ignition based on fuel moisture content, following FARSITE’s approach to more realistically represent fire spread timing.
The delay formula is:
where:
\(\tau_{\text{base}}\) = base delay time at zero moisture [s] (default: 120 s)
\(M_x\) = extinction moisture content [fraction]
\(M_f\) = current fuel moisture content [fraction]
At extinction moisture (\(M_f \to M_x\)), delay approaches infinity (no ignition). At very dry conditions (\(M_f \to 0\)), delay approaches \(\tau_{\text{base}}\).
Usage:
spotting.enable = 1
spotting.enable_delay = 1
spotting.tau_base = 120.0 # base delay in seconds
Implementation: src/spot_ignition_delay.H
FARSITE Temporal Fire Acceleration Model
Implements the full FARSITE temporal acceleration model (McAlpine & Wakimoto 1991 / VanWagner’s equation) with per-cell temporal tracking and wind-onset time-lag capability. Models the delayed approach to quasi-steady-state spread rate in small fires or after wind changes.
Two models are available:
Size-Based Model (Catchpole et al. 1992)
where \(r_{\text{fire}} = \sqrt{A_{\text{burned}} / \pi}\). Global scaling based on current fire size. Fire achieves ~63% of QSS at \(r_{\text{fire}} = L_{\text{acc}}\), ~95% at \(r_{\text{fire}} = 3 \times L_{\text{acc}}\).
FARSITE Temporal Model (McAlpine & Wakimoto 1991)
where:
\(A = A_{\text{point}}\) (0.115 1/min) for small fires (perim < 402.3 m)
\(A = A_{\text{line}}\) (0.300 1/min) for large fires (perim >= 402.3 m)
\(t\) = elapsed time since entering current acceleration state [s]
Per-cell temporal tracking with automatic switching between point and line acceleration constants.
Wind-Onset Time-Lag Extension
When \(R_E\) changes (wind change):
Optional wind-lag model for realistic wind response. Fire ROS ramps up exponentially after wind changes, controlled by time constant tau_wind (default: 180 s = 3 min).
Usage - Size-based model (backward compatible):
acceleration.enable = 1
acceleration.L_acc = 50.0
Usage - FARSITE temporal model:
acceleration.enable = 1
acceleration.use_temporal = 1
acceleration.A_point = 0.115
acceleration.A_line = 0.300
acceleration.perim_limit = 402.3
Usage - Enable wind-onset time-lag:
acceleration.enable = 1
acceleration.use_temporal = 1
acceleration.enable_wind_lag = 1
acceleration.tau_wind = 180.0
Parameter |
Default |
Description |
|---|---|---|
|
0 |
Enable acceleration model (1=yes, 0=no) |
|
0 |
0=size-based, 1=FARSITE temporal |
|
50.0 m |
Length scale for size-based model |
|
0.115 1/min |
Point ignition acceleration constant |
|
0.300 1/min |
Line ignition acceleration constant |
|
402.3 m |
Perimeter threshold (20 chains) |
|
0 |
Enable wind-onset time-lag (1=yes, 0=no) |
|
180.0 s |
Wind response time constant (3 min) |
Tests: regtest/surface_spread/fire_acceleration/
Implementation: src/fire_acceleration.H
References: McAlpine, R.S. & Wakimoto, R.H. (1991). The acceleration of fire from point source to equilibrium spread. Forest Science, 37(5), 1314–1337.
Cell Size Effects Correction
Implements empirical correction factors for fire spread rates based on grid cell resolution. Accounts for the physical reality that finer computational grids better resolve fire-front geometry, leading to higher effective spread rates when compared to coarser grids.
Physical Basis
Fire spread models like Rothermel assume a representative fuel bed at a given scale. When implementing these models on computational grids:
Coarser grids: Numerical diffusion and averaging effects cause under-prediction of fire spread
Finer grids: Better resolution of fire-front topology leads to higher effective spread rates
FARSITE addresses this through empirical cell size correction factors.
Correction Formula
The rate of spread is adjusted as:
where:
\(R_{\text{base}}\) = rate of spread without correction (m/s)
\(dx_{\text{ref}}\) = reference cell size (default: 30 m)
\(dx_{\text{min}}\) = actual minimum cell size in domain = min(dx, dy, dz)
\(\text{exponent}\) = power-law exponent (default: 0.1)
Examples
With dx_ref = 30.0 and exponent = 0.1:
Cell Size (m) |
dx/dx_ref |
Correction Factor |
Effect |
|---|---|---|---|
10 |
3.0 |
(3.0)^0.1 = 1.116 |
+11.6% faster |
15 |
2.0 |
(2.0)^0.1 = 1.072 |
+7.2% faster |
20 |
1.5 |
(1.5)^0.1 = 1.041 |
+4.1% faster |
30 |
1.0 |
(1.0)^0.1 = 1.000 |
No change |
60 |
0.5 |
(0.5)^0.1 = 0.933 |
-6.7% slower |
100 |
0.3 |
(0.3)^0.1 = 0.878 |
-12.2% slower |
Usage
Cell size correction is designed for FARSITE elliptical propagation and applies only when
propagation_method = farsite:
propagation_method = farsite
cellsize.enable = 1
cellsize.dx_ref = 30.0
cellsize.correction_exponent = 0.1
Cell size correction is not applied to level-set or MTT propagation, even if
cellsize.enable = 1. The level-set method naturally incorporates resolution effects
through its CFL-based timestep constraint.
Parameter |
Default |
Description |
|---|---|---|
|
0 |
0: Disable, 1: Enable empirical correction |
|
30.0 |
Reference cell size for correction [m] |
|
0.1 |
Power-law exponent (range: [0.0, 1.0]) |
Tests: regtest/surface_spread/cell_size_correction_*/
Implementation: src/cell_size_correction.H
References: Finney, M.A. (1998). FARSITE: Fire Area Simulator. USDA Forest Service RMRS-RP-4.
Advanced Physics Features
This section describes 10 advanced physics features that have been implemented as standalone header-only modules with GPU-compatible kernels. They can be used independently or combined to enhance fire behavior modeling.
Radiation-Driven Preheating Distance
File: src/radiation_preheating.H
Purpose: Calculate the distance ahead of the fire front where fuel is preheated by radiant energy from flames.
Key Functions:
compute_preheating_distance()- Compute preheating distance field from flame length and windcompute_preheating_distance_scalar()- Variant using scalar wind speed
Physics:
where:
\(L_f\) is the flame length (m)
\(\theta_{\text{tilt}}\) is the flame tilt angle (rad)
\(F_v\) is the view factor (dimensionless)
Value: Affects ignition timing and spread rate, important for understanding fire behavior in non-uniform fuels.
Fuel Particle Temperature Evolution
File: src/fuel_temperature.H
Purpose: Track surface temperature of fuel particles ahead of and behind the fire front.
Key Functions:
update_fuel_temperature()- Time integration of energy balance equationcheck_thermal_ignition()- Check for ignition based on temperature thresholdinit_fuel_temperature()- Initialize temperature fieldcompute_thermal_capacitance_field()- Compute thermal properties
Physics:
where:
\(Q_{\text{rad}}\) is the radiative heat transfer (kW/m²)
\(Q_{\text{conv}}\) is the convective heat transfer (kW/m²)
\(Q_{\text{loss}}\) is the radiative losses (kW/m²)
\(C_{\text{thermal}}\) is the thermal capacitance (kJ/m³/K)
Value: Determines actual ignition timing (not instantaneous), affects spotting ignition probability.
Fire Line Intensity Rate of Change
File: src/intensity_rate_of_change.H
Purpose: Compute temporal derivative of fireline intensity for blow-up detection.
Key Functions:
compute_intensity_rate_of_change()- Compute dI/dt fieldupdate_previous_intensity()- Store previous intensity for next time stepclassify_fire_behavior()- Classify fire as intensifying/steady/weakening
Physics:
where:
\(I\) is the fireline intensity (kW/m)
\(\Delta t\) is the timestep (s)
Value: Simple diagnostic for identifying dangerous “blow-up” conditions, useful for suppression resource allocation.
Flame Residence Time from Fire Intensity
File: src/intensity_residence_time.H
Purpose: Improve the burnout_time model with intensity-dependent residence time.
Key Functions:
compute_intensity_residence_time()- Compute intensity-dependent residencecompute_intensity_residence_field()- Field versioncompute_intensity_burnout_time()- Enhanced burnout time calculationclassify_combustion_rate()- Classify combustion as smoldering/flaming/intense
Physics:
where:
\(\tau_{\text{base}}\) is the base residence time (s)
\(I_{\text{byram}}\) is the Byram fireline intensity (kW/m)
\(\alpha\) is the intensity exponent (dimensionless)
Value: Higher intensity fires burn faster; refines existing burnout model.
Wind-Fuel Interaction Feedback
File: src/wind_fuel_interaction.H
Purpose: Adjust effective wind speed based on fuel structure (sheltering, penetration depth).
Key Functions:
estimate_LAI()- Estimate leaf area index from fuel propertiescompute_wind_reduction_factor()- Exponential attenuation through canopycompute_sheltered_wind_speed()- Combined functioncompute_sheltered_wind_field()- Field versionapply_wind_fuel_reduction()- In-place modification of wind
Physics:
where:
\(U_{\text{eff}}\) is the effective wind speed (m/s)
\(U_{\text{ref}}\) is the reference wind speed (m/s)
\(a\) is the extinction coefficient
\(\text{LAI}\) is the leaf area index
\(z\) is the height (m)
\(h\) is the canopy height (m)
Value: Dense fuels shelter surface from wind; current WAF models are height-based only, this adds fuel-structure feedback.
Fuel Moisture of Extinction Gradient
File: src/moisture_of_extinction.H
Purpose: Make moisture of extinction (\(M_x\)) depend on fuel properties, not constant 0.30.
Key Functions:
compute_moisture_of_extinction()- SAV-dependent M_xcompute_M_x_rothermel()- Original Rothermel formulation with mineral contentcompute_M_x_field()- Field versionclassify_fuel_ignitability()- Classify fuels by M_x
Physics:
where:
\(\sigma_{\text{fuel}}\) is the fuel surface-area-to-volume ratio (m²/m³)
\(\sigma_{\text{ref}}\) is the reference surface-area-to-volume ratio (m²/m³)
Value: M_x varies by fuel type (grass vs. timber), more accurate fire/no-fire transitions.
Flame Intermittency Factor
File: src/flame_intermittency.H
Purpose: Account for pulsating/intermittent flames rather than steady burning.
Key Functions:
compute_flame_intermittency()- Intensity-dependent intermittency factorcompute_flame_intermittency_field()- Field versionapply_flame_intermittency()- Apply to radiative heat fluxclassify_flame_stability()- Classify as flickering/pulsating/continuous
Physics:
where:
\(\gamma\) is the intermittency factor (0-1)
\(k\) is the intermittency coefficient
\(Q_{\text{rad,eff}}\) is the effective radiative heat flux (kW/m²)
Value: Real flames flicker and pulse; affects heat transfer efficiency and spotting.
Plume Entrainment Feedback on Surface Wind
File: src/plume_momentum_feedback.H
Purpose: Enhance heat_flux_model with horizontal momentum effects.
Key Functions:
compute_plume_velocity()- Upward plume velocity from buoyancycompute_plume_wind_increment()- Horizontal wind increment from plumecompute_plume_wind_field()- Field versionapply_plume_momentum_feedback()- In-place wind modificationapply_directional_plume_feedback()- Vectorial 2D feedback
Physics:
where:
\(k_{\text{momentum}}\) is the momentum coupling coefficient
\(w_{\text{plume}}\) is the plume updraft velocity (m/s)
\(I\) is the fireline intensity (kW/m)
Value: Fire creates inflow winds; can strengthen ROS by feeding fire with fresh air, important for fire whirl formation.
Fuel Bed Bulk Density Spatial Variation
File: src/fuel_loading_variation.H
Purpose: Allow spatially varying fuel loading from landscape files.
Key Functions:
apply_fuel_loading_multiplier()- Apply spatial multiplier to fuel loadinginit_uniform_multiplier()- Initialize with constant valueload_fuel_multiplier_from_file()- Load from ASCII file with IDW interpolationgenerate_random_fuel_variation()- Create random variation for ensembles
Physics:
where:
\(w_{0,\text{base}}\) is the base fuel loading (kg/m²)
\(m(x, y)\) is the spatial multiplier (dimensionless)
Value: Current models use uniform fuel properties; real landscapes have patchy fuel loading affecting local ROS and intensity.
Critical Heat Flux for Ignition
File: src/critical_heat_flux.H
Purpose: Replace binary ignition (φ < 0) with threshold-based ignition from incident heat flux.
Key Functions:
compute_critical_heat_flux()- Moisture-dependent critical fluxcompute_critical_heat_flux_field()- Field versioncheck_heat_flux_ignition()- Binary ignition checkcompute_ignition_probability()- Probabilistic ignition near thresholdcheck_moisture_dependent_ignition()- Combined function
Physics:
where:
\(q_{\text{crit}}\) is the critical heat flux (kW/m²)
\(q_{\text{base}}\) is the base critical flux (kW/m²)
\(k_{\text{moisture}}\) is the moisture sensitivity coefficient
\(M_{\text{fuel}}\) is the fuel moisture content (%)
Value: More physically realistic than geometric ignition; wet fuels need more heat. Used in Balbi and Lautenberger models.
Design Principles
All advanced physics features follow these design principles:
Header-only: Implemented as header-only libraries for easy integration
GPU-compatible: All kernels use AMReX GPU macros (
AMREX_GPU_DEVICE,ParallelFor)Standalone: Each feature can be used independently
No core modifications: Features don’t require changes to main solver loop
Well-documented: Extensive comments with physics explanations and references
Integration Example
Features can be integrated into the main solver by:
Including the appropriate header file
Allocating MultiFab fields as needed
Calling the compute functions in the time-stepping loop
Adding fields to plotfile output
Example integration pattern:
#include "radiation_preheating.H"
// Allocate field
MultiFab preheating_dist_mf(ba, dm, 1, 0);
// Compute in time loop
compute_preheating_distance(preheating_dist_mf, flame_length_mf,
vel_x_mf, vel_y_mf, view_factor);
// Add to plotfile
WriteSingleLevelPlotfile("plt", preheating_dist_mf, ...);
References
Each feature is supported by scientific literature:
Butler et al. (2004), Weber (1991) - Preheating
Dupuy & Maréchal (2011), Pickett et al. (2010) - Fuel temperature
Byram (1959), Rothermel (1972) - Intensity dynamics
Byram (1959), Rothermel (1972) - Residence time
Wilson (1985), Albini & Baughman (1979) - Wind-fuel interaction
Rothermel (1972), Anderson (1970) - Moisture of extinction
Finney et al. (2015), Frankman et al. (2013) - Flame intermittency
Finney et al. (2015), Linn et al. (2002) - Plume entrainment
Finney (1998), Andrews (2018) - Spatial variation
Drysdale (2011), Quintiere (2006), Balbi (2009) - Critical heat flux
Canadian Forest Fire Weather Index (FWI) System
File: src/canadian_fwi.H
Purpose: Full implementation of the Canadian FWI System for fire danger rating and fuel moisture tracking.
Key Functions:
compute_ffmc()- Fine Fuel Moisture Code (1-2 hour timelag)compute_dmc()- Duff Moisture Code (10-20 day timelag)compute_dc()- Drought Code (seasonal timelag, 50+ days)compute_isi()- Initial Spread Index (combines FFMC and wind)compute_bui()- Buildup Index (combines DMC and DC)compute_fwi()- Fire Weather Index (overall fire danger rating)
Physics: Van Wagner (1987) equations for moisture code evolution and fire danger indices. Provides standardized fire danger rating used operationally in Canada and many other countries.
Value: Operationally-tested indices for fire danger assessment, fuel moisture estimation, and fire weather monitoring.
Two-Fuel Model Blending
File: src/two_fuel_blending.H
Purpose: Blend rate of spread from two different fuel models for heterogeneous landscapes.
Methods: Linear, harmonic mean, maximum, and Finney (1998) FARSITE-style weighted quadratic blending of ROS values.
Value: Handles fuel transitions (grass/timber), heterogeneous fuel distributions, and multi-model scenarios.
Duff Moisture and Smoldering Combustion
File: src/duff_moisture_smoldering.H
Purpose: Track deep organic layer moisture and smoldering combustion following fire passage.
Physics: Exponential decay smoldering model with moisture-dependent time constants. Tracks long-duration heat release from deep organic layers.
Value: Captures post-frontal heat release, spotting from smoldering logs, and smoke production. Similar to FARSITE/FLAMMAP duff moisture tracking (DMC).