Modelling Bushfire Fuel Hazard Using Biophysical Parameters

Environmental gradients or biophysical parameters such as climate, topography and geology drive landscape-scale vegetation structure, species distribution and productivity. These gradients have the potential to provide detailed, fine-scale spatial prediction of the accumulation of bushfire fuels and hence fire hazard by elucidating patterns in field information in a consistent and repeatable way. Rapid visual assessment of bushfire fuel hazard via ratings provides fire and land management agencies with a measure of the probability of first attack success and general suppression difficulty of bushfires at a location. This study used generalised additive modelling to examine how measures of fuel hazard, recorded for locations in New South Wales, Australia, varied in response to environmental gradients and whether these gradients could be used to predict fuel hazard at a landscape scale. We found that time since last fire, temperature and precipitation were strong predictors of fuel hazard. Our model predictions for fuel hazard outperformed current operational methods; however, both methods tended to overestimate lower fuel hazard and underestimate higher fuel hazard. Biophysical modelling of fuel hazard provides significant advancement for predicting fuel hazard. These models have the capability to be improved and developed as additional fuel hazard data, fire history mapping and remote sensing of environmental variables advance both spatially and temporally.


Introduction
Landscape-scale wildfires are driven by four factors: fuel, its availability to burn, the influence of weather on rate of spread of the fire and sources of ignition [1,2]. Of these, fuel is the single factor that can be manipulated by fire and land management agencies to mitigate the risk of fire to human life, property and the environment. Accurate mapping of fuel characteristics across the landscape is required by land management agencies for strategic planning of where and how to treat fuels, fire behaviour predictions and determining impacts of fire such as ecological effects and smoke emissions [3,4]. Disturbance history is one of the most important factors driving fuel accumulation and structure, and hence current fuel measures used operationally are generalised extrapolations based on vegetation mapping and time since last fire [5][6][7]. Fuel is highly variable, both spatially and temporally, as variations in available fuel are a result of complex and dynamic biological and physical processes that influence species distribution and productivity, including climate, topography and geology. This recognition underlies the necessity for more sophisticated methods to estimate and map fuel characteristics at a landscape scale.
for estimation in a particular setting. These relationships are commonly extended to near-surface, elevated and bark fuels, which are unlikely to conform to the Olsen curve [31]. In contrast, statistical approaches can allow for a variety of functional forms. Generalised additive models (GAMs) represent one such approach. GAMs are extensions of generalised linear models (GLMs) in which the linear predictor includes smooth functions for some or all of the predictor variables. This allows GAMs to flexibly reproduce a wide range of response shapes in a data-driven manner, rather than being constrained to a shape dictated by an a priori model [32]. The ability of GAMS to handle non-linear data can aid in the development of ecological models that better represent the underlying data [33][34][35].
Relationships between biophysical processes and fuel accumulation and decomposition can be used to predict fuel hazard, and the value of this approach is that it provides an ecological context in which to understand and predict fuel dynamics [4]. Biophysical parameters can be topographical (e.g., elevation, slope, aspect), geological (e.g., soil), biological (e.g., productivity, disturbance history) and climatic (e.g., temperature, precipitation, radiation). As outlined previously, the conversion between hazard ratings and ratings to numeric values of fuel loads (t ha −1 ) lacks accuracy and precision. Development of models incorporating environmental gradients that drive landscape structure, productivity and fuel accumulation will allow for explicit, fine-scale predictions of fuel.
To understand the drivers of wildfire fuel hazard across New South Wales forests and woodlands, we used land management agency datasets of visually assessed fuel hazard to address the following aims: (1) determine whether biophysical processes could be used to predict and map landscape fuel hazard; (2) assess the predictive capability of biophysical models and (3) compare biophysical models of fuel hazard to current operational methods.

Data Sources, Collation, Filtering and Manipulation
Fuel in South-Eastern Australian forests, woodlands and shrublands were divided into four strata for the purpose of bushfire risk planning, each based on the position in the vegetation profile, and these are defined as bark, elevated, near-surface and surface fuels. A fuel hazard rating of low, moderate, high, very high or extreme can be assigned to each fuel stratum by assessing it against key attributes such as fuel quantity and continuity [36]. Overall fuel hazard (OFH) is a combination of bark + elevated + combined surface and near-surface fuel hazard, which is used by land and fire management agencies for risk planning [22,36,37]. A similar visual method of fuel hazard assessment for dry eucalypt forests with a litter and shrub understorey has been developed for use in conjunction with the Vesta fire behaviour model [38]. A fuel hazard rating (measured using the overall fuel hazard guide) for each stratum can be converted into Vesta fuel hazard ratings, which are then directly related to predicted fire behaviour [15]. We examined both fuel hazard and near-surface cover (as a proxy for hazard) from three large datasets for New South Wales: (1) NSW Office of Environment and Heritage, National Parks and Wildlife Service, OFH assessment database (n = 3399); (2) NSW Office of Environment and Heritage, National Parks and Wildlife Service, Enhanced Bushfire Management Plan dataset (n = 98) and (3) Centre for Environmental Risk Management of Bushfires, University of Wollongong dataset (n = 766) ( Figure 1 and Table S1).
A total of 4263 sites were included for fuel hazard measures across forested vegetation classes of dry sclerophyll forest (n = 3201), wet sclerophyll forest (n = 720), grassy woodland (n = 243), rainforest (n = 25) and semi-arid woodlands (n = 74). Analysis was restricted to sites with a known fire history and a valid date of fuel assessment. Sites spanned a broad range of elevations (0-1650 m asl.), mean annual precipitation levels (307-1883 mm) and mean annual temperatures (6.7-19.7 • C). All fuel hazard ratings were assessed using the Overall Fuel Hazard Guide [22] and translated into numerical Vesta fuel hazard ratings [36] for statistical analysis. This guide was developed for dry sclerophyll forests, hence the bias towards these vegetation types in the datasets. Additionally, as these hazard ratings were collected for Forests 2020, 11, 925 4 of 13 operational purposes, there was an inherent bias towards locations assumed to pose a potential risk due to proximity to property or age of fuel. The distribution of data varied between fuel strata ( Figure 2). A total of 4263 sites were included for fuel hazard measures across forested vegetation classes of dry sclerophyll forest (n = 3201), wet sclerophyll forest (n = 720), grassy woodland (n = 243), rainforest (n = 25) and semi-arid woodlands (n = 74). Analysis was restricted to sites with a known fire history and a valid date of fuel assessment. Sites spanned a broad range of elevations (0-1650 m asl.), mean annual precipitation levels (307-1883 mm) and mean annual temperatures (6.7-19.7 °C). All fuel hazard ratings were assessed using the Overall Fuel Hazard Guide [22] and translated into numerical Vesta fuel hazard ratings [36] for statistical analysis. This guide was developed for dry sclerophyll forests, hence the bias towards these vegetation types in the datasets. Additionally, as these hazard ratings were collected for operational purposes, there was an inherent bias towards locations assumed to pose a potential risk due to proximity to property or age of fuel. The distribution of data varied between fuel strata ( Figure 2).
Biophysical parameters were sourced from spatially referenced data (Tables S2 and S3) and compiled using ArcMap ver. 10.4.1 (Esri, Redlands, CA, USA). Data layers were selected for their potential influence on forest distribution and productivity and hence fuel accumulation.

Statistical Analysis
All data processing and statistical analyses were undertaken using the R statistical programming language version 3.3.1 [39]. Collinearity between predictors was assessed using Pearson's correlation coefficient and a subset of predictors was discarded from further analysis if pairwise correlations  Biophysical parameters were sourced from spatially referenced data (Tables S2 and S3) and compiled using ArcMap ver. 10.4.1 (Esri, Redlands, CA, USA). Data layers were selected for their potential influence on forest distribution and productivity and hence fuel accumulation.

Statistical Analysis
All data processing and statistical analyses were undertaken using the R statistical programming language version 3.3.1 [39]. Collinearity between predictors was assessed using Pearson's correlation coefficient and a subset of predictors was discarded from further analysis if pairwise correlations exceeded 0.7 (given in Table S3).
Overall fuel hazard ratings were converted from hazard ratings (low, moderate, high, very high or extreme) to integers of 1 to 5 and treated as pseudo-continuous, whilst all other surface, near-surface, elevated and bark fuel hazard ratings were converted to numeric equivalents according to the Vesta fire behaviour model [15,36].
Generalised additive models (GAM, mgcv package) [40] were used to identify which combination of biophysical parameters best predicted each of the fuel measures examined and allow for the possibility of a non-linear relationship between parameter and fuel measure (Table S2). Setting each fuel measure in turn as the response variable, models were fitted for all combinations of predictor variables using the dredge function in the MuMIN package [41], which returns a list of models representing every possible combination of predictor variable and ranked in order of model performance based on the Akaike information criterion (AIC; the best model being that with the lowest AIC value). Models with ∆AIC < 4 points of the top model (i.e., with the smallest AIC) were considered as being strongly supported [42].
Continuous variables were included as smooth terms using thin-plate regression splines, while categorical variables were included as random effects. All response variables were treated as pseudo-continuous and models were fitted with a Gaussian distribution and identity link function. Interactions between predictors were not considered in this study.
The performance of the best-ranked biophysical model was then compared to a reference model using only variables currently used by operationally. Fire and land management agencies use a fuel type map overlaid with time since last fire to calculate the current fuel level based on the fuel accumulation parameters [6]. Generalised additive models were fitted using only fuel group (as a random effect) and time since last fire (as a smooth term) as variables to provide a reference model.
The predictive accuracy of both the selected biophysical model and the reference model was estimated using leave-one-out cross-validation to calculate mean absolute error (MAE).
For each continuous predictor included in the selected model, its marginal contribution to the predicted fuel measure estimate was graphed to show the shape and direction of the effect.

Results
Across all of the fuel hazard ratings, GAMs based on biophysical variables had substantially better predictive performance than the corresponding reference models representing current operational practice (Table 1). For overall fuel hazard, the best-ranked biophysical model explained 43% of the deviance, compared to 27.8% for the reference model. For bark fuel hazard, the best-ranked biophysical model explained 36.3% of the deviance, compared to 17.1% for the reference. For elevated fuel hazard, the best-ranked biophysical model explained 31.2% of the deviance, compared to 13.8% for the reference. For near-surface cover, the best-ranked biophysical model explained 31.2% of the deviance, compared to 15.4% for the reference. For surface fuel hazard, the best-ranked biophysical model explained 48.5% of the deviance, compared to 28.5% for the reference (Table 1).

Model Prediction Error
The performance of the selected models for all fuel hazard ratings was best at the moderate (2) and high (3) hazard ratings, as indicated by the prediction error overlapping the dashed line ( Figure 3). There was a general tendency to underestimate very high (4) and extreme (5) and to overestimate low (1) hazard ratings, indicated by greater deviation from the dashed line ( Figure 3). Overall fuel hazard has a mean absolute error of 0.83 over the entire dataset (using leave-one-out cross-validation), compared to 0.96 in the reference model (Table 1). All other fuel hazard rating best-ranked models had lower mean absolute error rates than the corresponding reference models, with values ranging from 0.58 to 0.76. Prediction error ( Figure 3) indicated that, across the entire model, the number of underand over-predictions were similar and that, at overall fuel hazard ratings < 3, a greater proportion were over-predicted whilst, at rating levels > 4, most ratings were under-predicted.

General Trends in Biophysical Parameters
Time since fire was a key predictor in all models ( Table 2). All fuel hazard measures increased immediately following fire and then plateaued; the rate of increase varied with fuel strata. Overall fuel hazard increased for the first~15 years following fire ( Figure S1a), bark fuel hazard increased for the first~8 years ( Figure S2a), elevated fuel hazard increased for the first~15 years ( Figure S3a), near-surface cover increased for the first~7 years ( Figure S4a) and surface fuel hazard increased for the first~12 years following fire ( Figure S5a).  (Table 1). All other fuel hazard rating bestranked models had lower mean absolute error rates than the corresponding reference models, with values ranging from 0.58 to 0.76. Prediction error (Figure 3) indicated that, across the entire model, the number of under-and over-predictions were similar and that, at overall fuel hazard ratings < 3, a greater proportion were over-predicted whilst, at rating levels > 4, most ratings were underpredicted.

General Trends in Biophysical Parameters
Time since fire was a key predictor in all models ( Table 2). All fuel hazard measures increased immediately following fire and then plateaued; the rate of increase varied with fuel strata. Overall fuel hazard increased for the first ~15 years following fire ( Figure S1a), bark fuel hazard increased for the first ~8 years ( Figure S2a), elevated fuel hazard increased for the first ~15 years ( Figure S3a), nearsurface cover increased for the first ~7 years ( Figure S4a) and surface fuel hazard increased for the first ~12 years following fire ( Figure S5a). Table 2. Full model terms for the best-ranked biophysical model for fuel hazard. p-values of terms included in best-ranked models (GAM) are indicated by * p < 0.05, ** p < 0.01, *** p < 0.00. Model terms not included in the best-ranked model are indicated by a dash (−).

Model Terms Overall Bark Elevated Near-Surface Cover Surface
Smooth terms Time since last fire *** *** *** *** *** Mean annual temp.   Table 2. Full model terms for the best-ranked biophysical model for fuel hazard. p-values of terms included in best-ranked models (GAM) are indicated by * p < 0.05, ** p < 0.01, *** p < 0.00. Model terms not included in the best-ranked model are indicated by a dash (−).

Model Terms Overall Bark Elevated Near-Surface Cover Surface
Smooth terms Time since last fire *** *** *** *** *** Mean annual temp. Fire type * 0.11 *** ** *** Soil order ** *** *** *** *** Region *** *** *** *** *** Climatic variables were included in all best-ranked models ( Table 2). Mean precipitation during the warmest quarter and mean annual precipitation (Figures S1b-S5b) showed a similar trend, with the greatest fuel hazard in moderate precipitation ranges and reduced hazard at lower and higher precipitation levels. Mean annual temperature and mean temperature of the warmest quarter were included in the best-ranked models (Figures S1c-S5c). Overall and elevated fuel hazard showed a weak negative trend with mean annual temperature, and near-surface cover showed a negative trend in mean temperature of the warmest quarter. All other fuel measures showed variable trends with mean temperature of the warmest quarter.
Other significant predictor variables included elevation, solar radiation, woody foliage cover, net primary productivity and CTI (Compound topographic index) (Figures S1c-i-S5c-i), although relationships were varied and weak. Random factors included in all best-ranked models ( Table 2) included fire type (of previous fire wildfire or prescribed fire), soil order and region.

Discussion
We developed biophysical models of fuel hazard across forests and woodlands of New South Wales. Our results show that the predictions of these models substantially outperformed those of reference models that considered only time since fire and fuel type, as per current operational approaches. Time since last fire, precipitation, temperature and elevation were all strong predictors of fuel hazard measures, consistent with the notion that climate and more specifically moisture availability is the key environmental determinate of the structure and distribution of Australian forests and woodlands [29,43]. Studies examining fuel accumulation are often limited in either spatial or temporal scope [44,45]. Biophysical models provide a substantial improvement in predictive capability for fuel mapping; with the increased quality and availability of spatial data, there needs to be a consistent and reproducible quantitative approach to fuel hazard mapping.

Time since Last Fire
Our analysis showed that time since last fire was a strong predictor in all selected biophysical models. The relationship between fuel hazard and time since last fire followed a similar form to that traditionally used to predict the accumulation of litter and fuel over time [30,42]. A study by McColl-Gausden, Bennett [13] in Victoria, Australia found that relationships with time since fire varied with strata, surface and bark strata showing a similar plateau, whilst, in contrast to our study, near-surface and elevated fuel hazard increased to~10-20 years of time since fire then gradually decreased until 100 years since fire. The advantage of using generalised additive models is that the form of the relationship between each predictor variable and the fuel hazard measure being modelled is data-driven rather than based on a priori assumptions of functional form, therefore allowing for examination of these complex relationships.
Variations in these fuel hazard components after extended time since last fire could be attributed to shrub senescence and die-back in the absence of fire [46][47][48], leading to a reduction in near-surface cover and elevated fuels. Our study showed little evidence of this, with all strata following a relationship comparable to the Olson curve currently used for growth and fuel modelling [30]. This is in contrast to with findings from Gould, McCaw [14], which suggested that, when surface fuel and understorey layers have reached an effective steady state, the hazard rating of fibrous barked trees will continue to increase. This has important implications in operational terms as bark fuel hazard is not considered when predicting fire rate of spread or flame height using fire behaviour models; however, it is considered when predicting spotting distance [15], which has been associated with house loss [49]. However, bark hazard is incorporated into the overall fuel hazard rating used for determining chances of extended first attack success [36]. Understanding the underlying factors that drive an increase in overall fuel hazard is of the utmost importance if hazard reduction treatments are to be appropriate and effective.

Climatic Variables
Landscape variation in fuel is a function of the local water and energy budgets available for the production and desiccation of plant biomass and therefore a function of climate [50]. Traditional fuel accumulation curves use time since last fire along with constants for initial fuel load, litter fall and rates of decomposition. Our models take this a step further to explore which climatic and landscape variables influence fuel accumulation, without applying a priori constants. Eucalypt open forests predominately occupy areas with relatively high moisture availability, whilst less productive woodlands are common in drier environments; therefore, fuel quantity can be linked with moisture availability [2]. In our study, mean precipitation during the warmest quarter was present in all selected models, with the exception Forests 2020, 11, 925 9 of 13 of near-surface cover, in which mean annual precipitation was a stronger predictor. These models showed a similar trend between fuel hazard and precipitation during the warmest quarter, with the greatest hazard in areas with around 300 mm precipitation and reduced hazard at lower precipitation levels for overall, bark and surface fuel hazard models, suggesting that areas with dry summer periods had reduced fuel loads. The trend between near-surface cover and mean annual precipitation was similar, with reduced near-surface cover in areas with less than 500 mm precipitation, which agrees with previous studies that suggest that higher moisture availability implies higher growth and hence higher fuel loads [51]. Our results also showed a reduction in fuel hazard for some fuel measures at higher precipitation levels, suggesting that, above certain thresholds, moisture may also limit fuel accumulation.
Temperature was present as a variable in all selected models in our study. This is consistent with Gordon et al. [34], who found that mean annual temperature was the strongest predictor of above-ground carbon (biomass), with a strong negative and near-linear relationship with total above-ground carbon, tree and litter carbon. Similarly, Thomas et al. [52] found negative relationships between mean annual temperature and steady state litter load across eucalypt forests of New South Wales. Our study suggested weak negative trends in mean annual temperature for overall and elevated fuel hazard and negative trends in mean temperature of the warmest quarter and near-surface cover. All other fuel measures showed variable trends with mean temperature of the warmest quarter. The negative relationships of temperature with overall and elevated fuel hazard and near-surface cover agree with previous work by Mathews et al. [53] in dry eucalypt forests of New South Wales, which reported that warmer and drier climates produce low amounts of fine fuel. Interestingly, bark fuel hazard showed the opposite trend of increasing hazard in areas with hotter and drier climates. This is possibly a reflection of eucalypt species distribution and bark type.
Elevation was also included as a predictor variable in all best-ranked models, with little agreement in trends. In a similar analysis, Rollins, Keane [10] found that elevation was commonly included as an important variable in most fuel models and suggested that this is likely due to the high accuracy of mapped elevation relative to other biophysical gradients. Further research on the biophysical parameters driving bark hazard are required. Models of the hazard associated with each bark type (ribbon/candle, stringy/fibrous barks and other bark types) would be desirable but would likely require additional data collection in order to ensure adequate coverage of spatial and biophysical variation.

Other Parameters and Limitations
Although biophysical models performed better than current operational models, the deviance explained was still low (between 48.5 and 31.2% deviance explained). All models performed best in the mid-ranges (2)(3)(4) and worst at the upper (5) and lower (1) levels of hazard. There are several limitations with the datasets that may explain the deviation.
Twelve National Parks and Wildlife Service regions were examined across the three datasets. For this analysis, as region was a categorical variable, we included region as a random effect. We suggest that park region influenced fuel hazard measures for a number of reasons. The distribution of data among national park regions was not equal. Some regions contained very little data, i.e., for overall fuel hazard, several regions contained less than 50 sites, equating to approximately 1% of the sites examined, whilst other regions contained greater than 20% of the sites. Models would be significantly improved by a more stratified sampling scheme, with effort focused on sampling regions which currently have insufficient data.
It is well established that considerable observer variation occurs in the visual estimation of vegetation cover [54][55][56] and fuel hazard rating [27,28]. Visual estimations are considerably quicker and more cost effective than direct sampling methods; however, this may come at a cost to accuracy and consistency. Variability in observations can be due to differing individual perceptions and difficulties with identifying boundaries between fuel strata layers [28]. We agree with previous studies that propose that ongoing training and calibration may minimise this uncertainty in visual observation methods [28,54,55].
Another potential source of variability in estimating hazard ratings or cover scores is the timing of the assessment (i.e., season). In eucalypt dominated forests, maximum litter fall generally occurs in summer and minimum in winter [57][58][59] and hence this may influence the observed surface fuel hazard. For example, Gould et al. [21] observed seasonal litter fall patterns in jarrah forest and noted that this pattern was consistent across a range of forest types and hence influenced the litter bed and resultant hazard rating. Despite these limitations, visual assessment of fuel hazard is routinely used by fire and land management agencies to assess fuels, assess risk and populate fire behaviour models.
Soil type was also included as a random effect in all the selected models. McColl-Gausden, Bennett [13] showed that soil variables explained the greatest variation in fuel hazard in almost all fuel hazard strata examined in Victorian forests, with soil bulk density as one of the strongest predictors. Soil fertility is known to influence the distribution of species [60,61] and interactions between moisture availability and soil texture/depth are also likely to play an important role [62]. Soil texture, rainfall and nutrient status of the soil have also been shown to influence the species composition of fire-prone communities [63].
The type of last fire (prescribed or wildfire) was also included in all selected models as a random effect. This suggests that wildfire and prescribed burn not only consume different amounts of biomass, and therefore the starting point of fuel accumulation after fire is likely different between the fire types, but also potentially influence regeneration patterns. Vivien et al. [64] reported that higher and lower severity sites were characterised by differences in the retention of the overstorey, the percentage cover of litter and shrubs and the geology. The development of future models could be greatly improved by examination of both soil variables and fire type; however, this would require additional data collection regarding spatial variation and accurate mapping of fire severity classes, soil texture and nutrient status.

Conclusions
There is a need for comprehensive, consistent and comparable spatial fuels data for fire and land management agencies. Modelling landscape-scale fuel requires significant resources for adequate data collection and collation, comprehensive geographical information system (GIS) processing and complex statistical analysis [4]. The integration of biophysical gradients with fuel data improved the prediction of landscape-scale fuel hazard when compared to models using time since last fire and fuel group as the only variables, as per current operational practice. This analysis of the individual fuel hazard components allowed insight into the biophysical variables driving fuel hazard. Our analysis of more than 4000 observations of fuel found that time since last fire, climatic parameters, region, soil and fire type drive fuel hazard across New South Wales. Biophysical models provide a framework that can be improved over time with increased data collection of fuel hazard measures and advancement of spatial data. This modelling could help to improve preparedness and operational decisions such as risk planning and fire behaviour prediction.
Supplementary Materials: The following are available online at http://www.mdpi.com/1999-4907/11/9/925/s1, Figure S1. Plots of GAM smooth terms for overall fuel hazard rating of best-ranked model; Figure S2. Plots of GAM smooth terms for bark fuel hazard rating of best-ranked model; Figure S3. Plots of GAM smooth terms for elevated fuel hazard rating of best-ranked model; Figure S4. Plots of GAM smooth terms for near-surface fuel cover of best-ranked model; Figure S5. Plots of GAM smooth terms for surface fuel hazard of best-ranked model; Table S1. Fuel measures, units and categories; Table S2. Candidate predictor variables considered for modelling; Table S3. Predictor variables considered for analysis; however, these were discarded due to high correlation with other selected variables.