Rothermel

Vegetation

Code

// Rothermel (1972) model for wildland fuels

// -------------------------------------------
// Model parameters
// These must be defined below, or included as a user-defined layer
//
// Mx , Moisture content of extinction
// Se , fuel particle effective mineral content, lb silica-free minerals / lb ovendry wood
// S_T , fuel particle total mineral content, lb minerals / lb ovendry wood
// Mf , lb moisture / lb ovendry wood
// rho_p , ovendry particle density, lb/ft^3
// h , fuel particle low heat content, Btu/lb
// sigma , surface area to volume ratio, 1/ft
// delta , fuel depth, ft
// w_o , ovendry fuel loading, lb/ft^2
// elevation , height of terrain above sea level, m
// -------------------------------------------

// Equations
REAL rho_b = w_o / delta; // ovendry bulk density, lb/ft^3
REAL beta = rho_b / rho_p; // Packing ratio
REAL Q_ig = 250 + 1116 * Mf; // Heat of preignition, Btu/lb
REAL epsilon = exp(-138 / sigma); // effective heating number

// Calculating tan(slope) in direction of the normal of the fire front
REAL tan_phi = dot(normal_vector, grad(elevation));

// slope factor, taking max of tan_phi and 0 results in speed increase uphill and no speed change from flat ground to downhill
REAL phi_s = 5.275 * pow(beta, -0.3) * pow(max(tan_phi, 0), 2); 

REAL W_n = w_o / (1 + S_T); // Net fuel loading lb/ft^2
REAL E = 0.715 * exp(-0.000359 * sigma); 
REAL B = 0.02562 * pow(sigma, 0.54);
REAL C = 7.47 * exp(-0.133 * pow(sigma, 0.55));
REAL beta_op = 3.348 * pow(sigma, -0.8189); // Optimum packing ratio
REAL wind_ftpmin = wind * 54.68; // Converting wind from km/h to ft/min
REAL phi_w = C * pow(wind_ftpmin, B) * pow(beta / beta_op, -E); // Wind coefficient

REAL xi = pow(192 + 0.2595*sigma, -1) * exp((0.792 + 0.681 * pow(sigma, 0.5)) * (beta + 0.1)); // Propagating flux ratio

REAL eta_S = 0.174 * pow(Se, -0.19); // Mineral damping coefficient
REAL eta_M = 1 - 2.59 * Mf/Mx + 5.11 * pow(Mf/Mx, 2) - 3.52 * pow(Mf/Mx, 3); // Moisture damping coefficient

REAL A = 1 / (4.77 * pow(sigma, 0.1) - 7.27);
REAL gamma_max = pow(sigma, 1.5) / (495 + 0.0594 * pow(sigma, 1.5)); // Maximum reaction velocity
REAL gamma = gamma_max * pow((beta / beta_op), A) * exp(A*(1 - beta / beta_op)); // Optimum reaction velocity

REAL I_R = gamma * W_n * h * eta_M * eta_S; // Reaction intensity, Btu/ft^2 min

REAL SROS = I_R * xi * (1 + phi_w + phi_s) / (rho_b * epsilon * Q_ig); // Surface rate of spread, ft/min

// Converting spread rate into m/s from ft/min
speed = SROS * 0.00508;