Dry Eucalypt Forest (Vesta Mk 2)

Vegetation

Fuel inputs

• temp
• rel_hum
• WAF
• Ws
• Hu
• drought_fac

• Temperature (degrees C)
• Relative humidity (%)
• Wind Adjustment Factor
• Surface and near surface fuel load (t/ha)
• Understorey height (m)
• McArthur drought factor [0 – 10]

Code

// Dry Eucalypt model - Vesta Mk 2

// -------------------------------------------
// Model parameters
// 0. Wind vector at 10m height (km/h), 'wind_vector' (input)
// 1. Temperature (degrees Celsius), 'temp' (input)
// 2. Relative humidity (%), 'rel_hum' (input)
// 3. Wind Adjustment Factor, 'WAF'
// 4. Surface & Near surface fuel load (t/ha), 'Ws'
// 5. Understorey height (m), 'Hu'
// 6. Drought factor, 'drought_fac'
// -------------------------------------------

// Temporary inputs - To be removed when real data is provided
REAL WAF = 1/0.4; // Update with real data when available
REAL drought_fac = 10; // Update with real data when available
REAL Ws = 1.1 * 10; // Update with real data when available
REAL Hu = 0.5; // Update with real data when available

//// Calculating U10 windspeed from wind vector
REAL wind_speed_10 = length(wind_vector);

// Applying wind reduction factor to get wind at 2m height
REAL wind_speed_2 = wind_speed_10 / WAF;

//// Calculating backing and flanking coefficients compared to head fire ROS ////
REAL wdot = dot(normalize(wind_vector),advect_normal_vector);

// Calculate length-to-breadth ratio (LBR).
REAL LBR = 1.0;
if (wind_speed_10 < 5){
    LBR = 1.0;
} else if (wind_speed_10 < 25){
    LBR = 0.9286 * exp(0.0505 * wind_speed_10); 
} else { 
    LBR = 0.1143 * wind_speed_10 + 0.4143; 
} 

// Determine coefficient for flank rank of spread, Rf = cf * Rh, using elliptical LBR equations 
REAL cc = sqrt(1.0-pow(LBR, -2.0)); 
REAL cb = (1.0-cc)/(1.0+cc); REAL a_LBR = 0.5*(cb+1.0); 
REAL cf = a_LBR/LBR; 

// Determine shape parameters for normal flow 
REAL f = 0.5*(1.0+cb); 
REAL g = 0.5*(1.0-cb); 
REAL h = cf; 

// Now calculate a speed coefficient using normal flow formula
REAL speed_fraction = (g*wdot+sqrt(h*h+(f*f-h*h)*wdot*wdot)); 

//// Estimating fuel moisture content using Gould et al. (2007) and Matthews et al. (2010) //// 
// Initialising the fuel moisture variable 
REAL MC; 

// Calculating fuel moisture between 12:00 and 16:59 (valid for sunny days from October to March) 
if (hour > 11 && hour < 17){
    MC = 2.76 + (0.124*rel_hum) - (0.0187*temp);}

// Calculating fuel moisture for other daylight hours (from 9:00 to 11:59 and 17:00 to 19:59 in this example)
else if ((hour < 12 && hour > 8) || (hour > 16 & hour < 20)){
    MC = 3.6 + (0.169*rel_hum) - (0.045*temp);}

// Calculating fuel moisture for night time hours (from 20:00 to 8:59) in this example
else{
    MC =  3.08 + (0.198*rel_hum) - (0.0483*temp);}  

//// Calculate moisture coefficients from Cruz et al. (2021) ////
// Using the polynomial fit
REAL moisture_coeff = 0.9082 + 0.1206*MC - 0.03106*MC*MC + 0.001853*MC*MC*MC - 0.00003467*MC*MC*MC*MC;

if (MC < 4.1) {moisture_coeff = 1;} 
else if (MC > 24) {moisture_coeff = 0;}

// Calculating fuel availability factor and using it to adjust moisture coefficient
REAL fuel_availability = 1.008 / (1 + 104.9 * exp(-0.9306*drought_fac));

moisture_coeff = moisture_coeff * fuel_availability;

//// Calculate slope effect - Needs to be done earlier for transitional behaviour
REAL slope_in_normal_dir = degrees(atan(dot(normal_vector,grad(elevation))));
slope_in_normal_dir = min(max(slope_in_normal_dir,-20),20);
REAL slope_coeff = pow(2.0, 0.1*fabs(slope_in_normal_dir));
REAL slope_factor;

if (slope_in_normal_dir >= 0)
   slope_factor = slope_coeff;
else
   slope_factor = slope_coeff/(2*slope_coeff-1.0);

//// Calculating head fire RoS for all phases
REAL ROS1_head = (0.03 + 0.05024*pow(wind_speed_2 - 1, 0.92628) * pow(Ws * 0.1, 0.79928)) * moisture_coeff * slope_factor;
if (wind_speed_2 < 2){ROS1_head = 0.03 * moisture_coeff * slope_factor;}

REAL ROS2_head = 0.19591 * pow(wind_speed_2, 0.8257) * pow(Ws * 0.1, 0.4672) * pow(Hu, 0.495) * moisture_coeff * slope_factor;

REAL ROS3_head = 0.05235 * pow(wind_speed_10, 1.19128) * moisture_coeff * slope_factor;

//// Calculating the probability of transitioning to phase 2
REAL G2 = -23.9315 + 1.7033*wind_speed_2 + 12.08822*moisture_coeff + 0.9524*Ws;
REAL P2;
if (Ws < 1){
    P2 = 0;}
else {
    P2 = 1 / (1 + exp(-G2));}

//// Calculating the probability of transitioning to phase 3
REAL G3 = -32.3074 + 0.2951*wind_speed_10 + 26.8723*moisture_coeff;
REAL P3;

if (ROS2_head < (0.3 / 3.6)){
    P3 = 0;}
else {
    P3 = 1 / (1 + exp(-G3));}

//// Combining phases probabilistically to obtain head fire RoS
REAL head_speed;
if (P2 < 0.5){
    head_speed = ROS1_head*(1-P2) + ROS2_head*P2;}
else {
    head_speed = (ROS1_head*(1-P2) + ROS2_head*P2)*(1-P3) + ROS3_head*P3;}

// Converting from km/h to m/s
head_speed = head_speed / 3.6;

//// Adjust for calculated speed coefficient for fire flanks
speed = head_speed * speed_fraction;