What is the overall ease of understanding? o What is the overall ease of understanding?

or truncation?

the legend?

Possible design requirements include the accommodation of colour vision deficiencies and other technical considerations such as the adequate appearance on low-quality displays, colour printers or standard photocopiers. Possible design requirements include the accommodation of colour vision deficiencies and other technical considerations such as the adequate appearance on low-quality displays, colour printers or standard photocopiers.

#### *2.4. Step 4: Designing the Categorization and Visualization Scheme 2.4. Step 4: Designing the Categorization and Visualization Scheme*

There are four design options for categorizing and colouring the values in the scale (Figure 2): There are four design options for categorizing and colouring the values in the scale (Figure 2):


#### 2.4.1. Gradient Continuity

The alternatives for gradient continuity are either to use the original values or categorize them (Figure 2, parts 1a and 1b). For a continuous gradient, the percentage of colour saturation is proportional to the raw datum magnitude. Continuous gradients are therefore precise and accurate but harder to interpret using the legend compared to categorized gradients.

#### 2.4.2. Hue Selection

Colours have strong psychological associations that affect the inferred meaning of information and its speed and ease of interpretation. Blue and green are associated with relaxation, calm and hope [21], and red is associated with danger [20] (Figure 2, part 2f). Historically, variations of a blue–green–yellow–orange–red sequence (which alludes to water, growing vegetation, dried vegetation and flame) have been used to represent escalating fire danger (see examples in Table 2). Accommodating colour vision deficiencies reduces the colour combination choices, particularly most of those in the traditional fire danger sequence [31]. Tools are available to assist with the testing of colour palettes for accessibility [31,32]. The design task is to choose colours appropriate for the implications of the data magnitudes, particularly the degree of attention or alarm.

Regarding the gradient design, there are a few distinct alternatives [33] (Figure 2, parts 2a–2e). A single sequential gradient is for data ranging over a meaning of zero or neutral to bad or good. A divergent sequential gradient is for data ranging over a meaning of good through neutral to bad. In Figure 2, parts 2b and 2c grade from a calming blue through white to an alarming red. The left and right ends each have a single hue. Part 2b is neutral in the middle, whereas part 2c has compressed and expanded ends. Part 2d, multiple sequential, is analogous to part 2a except that part 2d has multiple hues, which are in a rainbow spectrum in the example. Compared to a single hue, multiple hues provide more contrast over the range, making it easier to match the legend and signal levels of attention or alarm. Part 2e, multiple divergent, is analogous to part 2b except that part 2e has multiple hues for each side. The R software [34] package "inlmisc" [35,36] is useful for constructing continuous or categorized gradients.

A key concern is how these many design alternatives support or oppose the design goals. Only the single sequential and single divergent gradients (Figure 2, parts 2a and 2b) have the accuracy of continuous gradient continuity (Figure 2, part 1a), but they have a difficult interpretability. The remaining gradient alternatives require matching the legend to identify the magnitudes, but this can become quick and easy to interpret with familiarity and an effective use of colour psychology.

#### 2.4.3. Range Completeness and Scale Linearity

These are described together because an incomplete range is an extreme form of nonlinearity. The alternatives for range completeness (Figure 2, part 3) are whether to show the full range of the data or truncate the top or bottom and group the truncated data. All resolution is lost beyond the truncation points. The alternatives for scale linearity (Figure 2, part 4) are for two components, independently:

	- For categorized gradient continuity, a nonlinear scale is used to vary the resolution over the range
	- For categorized gradient continuity, a nonlinear colour gradient is used to communicate the varying meaning or importance of the information over the range

For a continuous gradient, the same result can be achieved from nonlinearity in either of the above two components.

Nonlinear-systematic (part 4b) methods use a smooth function such as log or power to transform the output, while nonlinear-irregular methods use a non-smooth progression such as Jenks [37]. Nonlinear colouring requires the referral to the legend to understand the magnitudes. This is a trade-off between the goals of drawing attention to where it is needed and improving the speed and ease of understanding.

#### 2.4.4. The Design Process

There are copious settings and combinations of alternatives for the four design options. Getting to a result is an iterative process, with analysts and subject matter experts trying alternatives and making trade-offs, hopefully avoiding the anchoring to tradition or early trials. We cannot recommend a path through the four design options other than saying it is iterative and concurrent. We do, however, recommend a starting point or baseline, which is the extreme of displaying all the data completely and accurately and ignoring the goals of quick and easy understanding. The baseline has continuous gradient continuity, an achromatic colour gradient of white through greys to black and a linear scale with no truncation (Table 3).

**Table 3.** The baseline case alternatives for the design options. This is a starting point that presents complete and accurate information, while ignoring the goals of a quick and easy understanding of the information.


If categories are used, determining the number and their boundaries are fundamental design decisions [38]. The number of categories is a trade-off of accuracy (requiring more) and speed and ease of understanding (requiring fewer). Ideally, individual categories have no operationally significant physical differences, while adjacent categories do. In practice, all considerations require compromise. An example of determining categories for FWI System outputs based on physical differences is given by [9].

#### *2.5. Step 5: Evaluating and Revising the Scheme*

Once the design process is done and implemented, an essential further step is the ongoing work with decision-makers to evaluate the outputs and revise the design as necessary.

#### **3. Case Study: Designing the Scheme for Ontario's FOP Models**

We now describe the application of the above method to categorizing and visualizing FOP data.

#### *3.1. Case Study—Step 1: Understanding and Scoping the Data*

The data are outputs from process and statistical FOP models for Ontario, so we begin with their description. The lightning- and the human-caused FOP models have been used operationally since the mid-2000s and 2015, respectively. The daily lightning-caused fire occurrence is modelled as two separate processes [3]: the probability of a lightning strike will lead to a holdover ignition and the probability that an existing ignition "arrives" (is reported). The ignition model is mainly driven by the forest floor organic layer's moisture content, which determines the sustainability of smouldering and the survivability of the ignition. Additional factors are other moisture indicators, ecoregional modifiers and lightning strike polarity. The "arrival" model, which is conditional on a holdover ignition being present, is influenced by the surface litter moisture, organic layer moisture, wind

speed and ecoregional differences. Ontario's human-caused fire prediction system uses a set of logistic generalized additive models to model inherently nonlinear relationships with key drivers of human-caused fire occurrence, including seasonal and spatial patterns, fuel moisture and the characteristics of human land use [1]. The models are stratified regionally and by cause categories to account for different seasonal patterns in fire occurrence.

Both the lightning- and human-caused FOP models produce outputs for each of the 2574 cells in a grid that spans the province's approximately 91.9 million ha wildland fire management area (Figure 3). Most of the cells are about 20 km × 20 km or 40,000 ha, with some fractional cells at the boundaries. The units of the FOP output data are interpreted as the expected number of fires per cell, regardless of the cell size. The model calculations and outputs have a 64-bit floating point precision, but the data are transferred to the mapping software via a text file holding up to 12 significant digits. *Fire* **2021**, *4*, x FOR PEER REVIEW 10 of 20

**Figure 3.** The approximately 91.9 million ha extent of Ontario's fire management area and the 20 km × 20 km resolution of the fire occurrence prediction (FOP) grid for Ontario. **Figure 3.** The approximately 91.9 million ha extent of Ontario's fire management area and the 20 km × 20 km resolution of the fire occurrence prediction (FOP) grid for Ontario.

Regarding the range and frequency distribution, we analysed historical FOP data for each cell for each day from May 15 to August 31; 2016–2018 for human-caused and 1992– 2006 for lightning-caused fires. For this analysis, the start and end dates in each fire season were chosen to avoid the variability in spring and fall snow-free conditions, when the models are not making predictions for the entire province. The lower limit of the range is zero; there is no theoretical upper limit. Table 4 presents summary statistics for the data, with zeros excluded to characterize the important data more clearly. Most of the distribution statistics of the human- and lightning-caused data differ by about an order of magni-Regarding the range and frequency distribution, we analysed historical FOP data for each cell for each day from 15 May to 31 August; 2016–2018 for human-caused and 1992–2006 for lightning-caused fires. For this analysis, the start and end dates in each fire season were chosen to avoid the variability in spring and fall snow-free conditions, when the models are not making predictions for the entire province. The lower limit of the range is zero; there is no theoretical upper limit. Table 4 presents summary statistics for the data, with zeros excluded to characterize the important data more clearly. Most of the

of fires/cell) generated from May 15–August 31 for 2016–2018 (human-caused) and 1992–2006

Human 607,383 0.00003 0.00013 0.00107 0.00048 1.37168 Lightning 1,026,860 0.00090 0.00260 0.00955 0.00750 3.89129

Figure 4 shows the empirical probability distributions of the non-zero data. Those data were mostly clustered close to zero in both models, so we log-transformed the data for illustration (untransformed data are required for operational use). The magnitudes of

**Quartile Median Mean Third** 

**Quartile** 

**Maximu m** 

**First** 

(lightning-caused), with zeros removed.

**Number of Observations** 

tude.

distribution statistics of the human- and lightning-caused data differ by about an order of magnitude.

**Table 4.** Summary statistics of historical fire occurrence prediction model data (expected number of fires/cell) generated from 15 May–31 August for 2016–2018 (human-caused) and 1992–2006 (lightning-caused), with zeros removed.


*Fire* **2021**, *4*, x FOR PEER REVIEW 11 of 20

Figure 4 shows the empirical probability distributions of the non-zero data. Those data were mostly clustered close to zero in both models, so we log-transformed the data for illustration (untransformed data are required for operational use). The magnitudes of data between the lower tails of the two distributions differ by orders of magnitude. This presents a challenge for categorizing and colouring the data on a common scale for mapping. data between the lower tails of the two distributions differ by orders of magnitude. This presents a challenge for categorizing and colouring the data on a common scale for mapping.

**Figure 4.** Empirical probability distributions of the non-zero, human- and lightning-caused fire occurrence predictions (FOPs) for Ontario, May 15–August 31, for 2016–2018 (human-caused) and 1992–2006 (lightning-caused). FOPs for the two causes have distinctly different ranges and central tendencies because of the spatio-temporal processes of ignition. The data are transformed by log 10 for visualization here, but untransformed data are required for operational use. **Figure 4.** Empirical probability distributions of the non-zero, human- and lightning-caused fire occurrence predictions (FOPs) for Ontario, 15 May–31 August, for 2016–2018 (human-caused) and 1992–2006 (lightning-caused). FOPs for the two causes have distinctly different ranges and central tendencies because of the spatio-temporal processes of ignition. The data are transformed by log 10 for visualization here, but untransformed data are required for operational use.

#### *3.2. Case Study—Step 2: Understanding the Decision-Making Uses of the Information 3.2. Case Study—Step 2: Understanding the Decision-Making Uses of the Information*

We used a variety of methods to understand the FOP information needed and how it is used for daily decision-making: reviewing documentation, observing operational decision-making and (for some) working part-time in operational functions where FOP information is used. Most importantly, we held a series of engagements with fire management agency personnel. For example, we hosted a workshop in 2017 attended by agency personnel including regional and provincial Fire Intelligence Officers, external researchers and students. The workshop's topics included the purposes and methods of subjective FOPs by experts. We used a variety of methods to understand the FOP information needed and how it is used for daily decision-making: reviewing documentation, observing operational decision-making and (for some) working part-time in operational functions where FOP information is used. Most importantly, we held a series of engagements with fire management agency personnel. For example, we hosted a workshop in 2017 attended by agency personnel including regional and provincial Fire Intelligence Officers, external researchers and students. The workshop's topics included the purposes and methods of subjective FOPs by experts.

To understand the agency's use of FOP information, it is necessary to outline the agency's hierarchical structure and the responsibilities of each level. Ontario's fire management area (Figure 3) has two main parts, the Northeast and Northwest Regions, each To understand the agency's use of FOP information, it is necessary to outline the agency's hierarchical structure and the responsibilities of each level. Ontario's fire management area (Figure 3) has two main parts, the Northeast and Northwest Regions, each of

of which is divided into six or seven Fire Response Sectors. There is also a Provincial level. Each level has a sole or shared responsibility for various decisions; most are made with

sible for strategic fire response decisions and management, and the Sectors are primarily responsible for tactical fire response decisions and operations. The Province is primarily responsible for adjusting the near-term (1- to 21-day) capacity according to the demand via temporary commercial hiring and inter-provincial and international resource sharing. Several decisions are directly dependent on the potential number of fires anticipated

Prevention: for example, escalated and targeted messaging, temporary fire bans;

in various parts of the province. These decisions are made at the stated levels:

made by the Province, Regions and Sectors

which is divided into six or seven Fire Response Sectors. There is also a Provincial level. Each level has a sole or shared responsibility for various decisions; most are made with consultation or coordination between adjacent levels. The Regions are primarily responsible for strategic fire response decisions and management, and the Sectors are primarily responsible for tactical fire response decisions and operations. The Province is primarily responsible for adjusting the near-term (1- to 21-day) capacity according to the demand via temporary commercial hiring and inter-provincial and international resource sharing.

Several decisions are directly dependent on the potential number of fires anticipated in various parts of the province. These decisions are made at the stated levels:


The various FOP-dependent decisions have diverse needs in terms of the spatial extent and resolution of FOP information. For example, detection route designers can use relatively fine resolution information on the order of kilometres, while the Province needs only aspatial, numeric FOPs by region for resource-sharing decision-making.

#### *3.3. Case Study—Step 3: Specifying the Design Goals and Requirements*

The primary and conflicting goals are of course to display complete and accurate information and have the information quickly and easily understood. In twice-daily briefings, decision-makers have limited time (minutes) to view, interpret and absorb each of many information items regarding, for example, weather values, FWI and FBP System outputs, FOP, active fires, logistics and personnel. Completeness and accuracy are important because the information supports the many decisions described above, which are made under uncertainty and have potentially significant consequences.

Regarding specific information requirements:


Decision-makers expressed strong preferences for the number of categories, ranging from three to many categories, and they desired a familiar colour sequence (blue–green– yellow–red). There was also a strong preference for integers for all numbers related to FOP. Finally, we wished to accommodate colour vision deficiencies.

#### *3.4. Case Study—Step 4: Designing the Categorization and Visualization Scheme*

Design has been described as a messy process with a tidy outcome. We do not detail our circuitous journey but show and describe some key alternatives, stages and considerations. Figure 5a illustrates the baseline alternative (Table 3) applied to the humancaused FOP and actual fire arrivals for a selected day. The range of the colour gradient is 0–3 fires/cell, the upper limit of which is between the maximum human- and lightningcaused FOP (Table 4). The map area looks mostly white, with three small, pale grey patches; there is little useful information, especially considering the six actual fires that day, which is a low-to-moderately busy day for this cause. Adding more hues alone would make no meaningful difference because the data are clustered very near zero. Categorizing at this stage would make it worse. Truncating the upper limit to a low value somewhat close to

zero would add resolution and colour to the human-caused FOP here, but such truncation would lose all resolution of the rarer but critically important high lightning-caused FOP. *Fire* **2021**, *4*, x FOR PEER REVIEW 13 of 20

**Figure 5.** The human-caused fire occurrence prediction of June 10, 2018, mapped using the categorization and colouring schemes indicated in the legends. The black dots are the human-caused fires reported later that day. The scales range from zero to 3 fires/cell: (**a**) the baseline scheme (Table 3), which shows little useful information; (**b**) 4 + 1 categories; (**c**) 10 + 1 categories; (**d**) 20 + 1 categories. The additional category in (**b**), (**c**) and (**d**) are for a "no forecast model" or true zero. Using more categories and hues greatly increases the information portrayed but makes matching the colour to the legend more difficult. The spatial pattern in (**d**) corresponds with roads and settlements. **Figure 5.** The human-caused fire occurrence prediction of 10 June 2018, mapped using the categorization and colouring schemes indicated in the legends. The black dots are the human-caused fires reported later that day. The scales range from zero to 3 fires/cell: (**a**) the baseline scheme (Table 3), which shows little useful information; (**b**) 4 + 1 categories; (**c**) 10 + 1 categories; (**d**) 20 + 1 categories. The additional category in (**b**–**d**) are for a "no forecast model" or true zero. Using more categories and hues greatly increases the information portrayed but makes matching the colour to the legend more difficult. The spatial pattern in (**d**) corresponds with roads and settlements.

Our original solution was to use separate scales for the two fire occurrence causes (Figure 6). The lightning-caused FOP scale was linear. For the human-caused FOP scale, we led subject matter experts through a scenario for an area of Ontario that has a relatively high occurrence. When the FFMC (a strong indicator of sustainable ignition) was 90, that area was considered to have an elevated concern suitable for a High classification. We used the corresponding FOP magnitude (0.4 fires/cell) as the upper limit for the then 10 category scale. We determined Moderate similarly and interpolated with equal linear steps. In Figure 6, this scale is shown by the blue line, except that categories 1 to 10 in the original have been mapped to a 0 to 20 scale for comparability. Scaling the original totalfires map required a creative logic. The maps by individual cause were acceptable to decision-makers, but the inconsistency between the causes was ultimately unacceptable (and motivated the present work). Our original solution was to use separate scales for the two fire occurrence causes (Figure 6). The lightning-caused FOP scale was linear. For the human-caused FOP scale, we led subject matter experts through a scenario for an area of Ontario that has a relatively high occurrence. When the FFMC (a strong indicator of sustainable ignition) was 90, that area was considered to have an elevated concern suitable for a High classification. We used the corresponding FOP magnitude (0.4 fires/cell) as the upper limit for the then 10-category scale. We determined Moderate similarly and interpolated with equal linear steps. In Figure 6, this scale is shown by the blue line, except that categories 1 to 10 in the original have been mapped to a 0 to 20 scale for comparability. Scaling the original total-fires map required a creative logic. The maps by individual cause were acceptable to decision-makers, but the inconsistency between the causes was ultimately unacceptable (and motivated the present work).

**Figure 6.** The original and unified fire occurrence prediction (FOP) classification scales. Because of order-of-magnitude differences in FOP, a separate scale was originally used for each cause: linear for lightning (black) and subjectively determined, irregular, nonlinear for human (blue). The unified scale (red) is a systematic, nonlinear one generated by a power function; a cube root in this example. **Figure 6.** The original and unified fire occurrence prediction (FOP) classification scales. Because of order-of-magnitude differences in FOP, a separate scale was originally used for each cause: linear for lightning (black) and subjectively determined, irregular, nonlinear for human (blue). The unified scale (red) is a systematic, nonlinear one generated by a power function; a cube root in this example.

Several alternatives were considered for a unified scale that accommodated the conflicting needs of fine resolution at the low end and a high upper limit. Discussions with decision-makers indicated that concern increases relatively quickly as the likelihood of fire rises from zero. Providing a high resolution at the low end while retaining a high upper limit would require a great many colours (for example, like those of precipitation radar maps). That would be unfamiliar and confusing and would not correspond with psychological colour associations nor accommodate colour vision deficiencies. Piecewise linear and irregular scales were explored, but their abrupt changes made them difficult to interpret. We wanted a smooth, systematic progression of category boundaries and tested logarithmic and power functions. Those functions can be made to match fairly closely, but the power function had a more suitable shape at the low end. A power function takes the general form () = , with the shape controlled by parameters and . Our desired behaviour for the scale to increase quickly for low FOP but then increase progressively more slowly is provided when > 0 and 0 < < 1, since this family of power functions is monotonically increasing and concave down. The FOP scale is , and the category scale is (). We developed a parametric scaling tool with three inputs to generate and plot boundaries: shape parameter, 1 ; the upper limit of the FOP scale, ; and the number of categories, . The boundary for the top of category is Several alternatives were considered for a unified scale that accommodated the conflicting needs of fine resolution at the low end and a high upper limit. Discussions with decision-makers indicated that concern increases relatively quickly as the likelihood of fire rises from zero. Providing a high resolution at the low end while retaining a high upper limit would require a great many colours (for example, like those of precipitation radar maps). That would be unfamiliar and confusing and would not correspond with psychological colour associations nor accommodate colour vision deficiencies. Piecewise linear and irregular scales were explored, but their abrupt changes made them difficult to interpret. We wanted a smooth, systematic progression of category boundaries and tested logarithmic and power functions. Those functions can be made to match fairly closely, but the power function had a more suitable shape at the low end. A power function takes the general form *f*(*x*) = *ax<sup>b</sup>* , with the shape controlled by parameters *a* and *b*. Our desired behaviour for the scale to increase quickly for low FOP but then increase progressively more slowly is provided when *a* > 0 and 0 < *b* < 1, since this family of power functions is monotonically increasing and concave down. The FOP scale is *x*, and the category scale is *f*(*x*). We developed a parametric scaling tool with three inputs to generate and plot boundaries: shape parameter, <sup>1</sup>/*b*; the upper limit of the FOP scale, *FOPMax*; and the number of categories, *NumCat*. The boundary for the top of category *Cat* is

$$\text{Boundary\\_Fog\\_FOPMax}\_{\text{Cat}} \overset{\text{FOPMax}}{=} \text{FOPMax} \overset{\text{cart}}{\text{(\*)}} \overset{\text{bart}}{\text{(\*)}} \overset{\text{bart}}{\text{(\*)}} \overset{\text{bart}}{\text{(\*)}} \overset{\text{bart}}{\text{(\*)}} \overset{\text{u}}{\text{Cdf}} \overset{\text{u}}{=} \text{I}\_{\text{'}}^{\downarrow \text{b}} \overset{\text{u}}{\text{Cdf}} \overset{\text{u}}{=} \text{I}\_{\text{'}}^{\downarrow \text{u}} \overset{\text{u}}{\text{Cdf}} \overset{\text{cst}}{\text{5}} \dots \overset{\text{u}}{\text{ NumCat.}} \tag{1}$$

added if required. For convenience, we parameterized 1 , for which we tested values in the range of 1.5 to 4.5; for example, 2 yields a square root shape. works out to be . The red curve in Figure 6 illustrates the boundaries for 1 = 3, = 3 and = 20; for example, the boundary for = 11 (the top of category 11) is ~0.5. The tool lists the boundaries and graphs them as in Figure 6. Any FOP > *FOPMax* stays in the highest category. A category for true zero can be added if required. For convenience, we parameterized 1/*b*, for which we tested values in the range of 1.5 to 4.5; for example, 2 yields a square root shape. *a* works out to be *NumCat FOPMax<sup>b</sup>* . The red curve in Figure 6 illustrates the boundaries for <sup>1</sup>/*<sup>b</sup>* = 3, *FOPMax* = 3 and *NumCat* = 20; for example, the boundary for *Cat* = 11 (the top of category 11) is ~0.5. The tool lists the boundaries and graphs them as in Figure 6.

While working directly with subject matter experts, the tool facilitated the joint testing of alternatives for the number of categories, truncation and nonlinearity design options. These alternatives plus colouring needed to be adjusted simultaneously when trading off the goals because their effects interact. The FOP outputs for a set of representative days were mapped using R [35] for candidate sets of boundaries and colouring. Table 5 lists ways in which the alternatives for the number of categories, the amount of truncation While working directly with subject matter experts, the tool facilitated the joint testing of alternatives for the number of categories, truncation and nonlinearity design options. These alternatives plus colouring needed to be adjusted simultaneously when trading off the goals because their effects interact. The FOP outputs for a set of representative days were mapped using R [35] for candidate sets of boundaries and colouring. Table 5 lists

ways in which the alternatives for the number of categories, the amount of truncation and nonlinearity and number of hues generally interact in affecting several attributes related to completeness and accuracy of information or speed and ease of understanding. There are exceptions for some combinations and edge conditions. Every alternative for the design options improves some attributes and worsens others. The trade-off behaviour is more straightforward to work with than it may seem because the tool shows most of the trade-offs immediately. The difficulty lies in subjectively assessing the results and compromising on the attributes and goals.

**Table 5.** Tabulation of how alternatives for the number of categories, amount of truncation and nonlinearity and number of hues generally interact in affecting several attributes related to completeness and accuracy of information or speed and ease of understanding. There are exceptions for some combinations and boundary conditions. Every alternative improves some attributes (blue-grey shading) and worsens others (orange-tan shading).


#### *3.5. Case Study—Step 4: Results of the Design Process*

We state our current category and colouring design and give the rationale for the trade-offs made. The pressure to have a small number of categories and colours (~4) could not accommodate the need for fine resolution at the low end; we used 20 categories (plus a true zero if needed). Figure 5b–d show the same FOP data as Figure 5a but with 4, 10 and 20 categories, respectively. Only the largest number of categories reveals the meaningful network pattern of lines and nodes that corresponds roughly with roads and settlements.

In addition, a highly nonlinear scale was needed to show the patterns in the data. We used a nonlinear-systematic scale with boundaries obtained using Equation (1) with parameters <sup>1</sup>/*<sup>b</sup>* = 3, *FOPMax* = 3 and *NumCat* = 20. The boundaries are given in Table 6, stated in units of fires/cell and cells/fire. The final shape of the nonlinear scaling corresponds to parameter settings of *a* = 20·3 <sup>−</sup> <sup>1</sup>/<sup>3</sup> and *b* = <sup>1</sup>/3.

$$f(\mathbf{x}) = \begin{cases} \ 20 \cdot 3^{-1/3} \text{ x}^{1/3} \\ \ 20 \end{cases}, \quad 0 \le \mathbf{x} < 3 \tag{2}$$

which is illustrated by the red curve in Figure 6.


**Table 6.** Categories and colours used for mapping fire occurrence prediction model outputs. Note that the map legend has a highly simplified integer scale and broad adjective categories.

> Regarding truncation, we considered the rarity and operational importance of extreme FOP magnitudes. The 20th category ends at 3 fires/cell according to Equation (1), but that category is used for all higher FOP magnitudes. The maximum FOP in Table 4 is ≈3.9 fires/cell.

() = > 0 0 < < 1 () 1 ⁄ = ∙ ( ) 1 ⁄ , = 1, 2, 3, … , . Regarding colouring, we used a nonlinear, divergent scale with one hue for the low end and multiple hues for the high end (Table 6). The hues transition from light blue to yellow through orange to red, which mostly follows the traditional blue-to-red progression. Avoiding green in that sequence accommodates some types of colour vision deficiency [31]. Even though there are 20 categories, having four main colours is easy to interpret and consistent with other CFFDRS outputs (for example, Table 1). The gradient within each main colour is difficult to match with the legend, but nonetheless reveals meaningful spatial patterns in the maps. Figure 7a–c show the final categorization and colouring scheme in example daily FOP maps of human- and lightning-caused fires and total fires, respectively. We intentionally show maps as formatted for operational use. They are intended for display on large monitors, but these reduced versions still show the colouring and categorization results adequately. Larger versions are provided in the Supplementary Material.

= 20 = 11

  ⁄ = 3 = 3

1 

1 ⁄

**Figure 7.** Examples of fire occurrence prediction maps that use the categorization and colouring in Table 6: (**a**) human-caused; (**b**) lightning-caused; (**c**) total. The legend uses units of fires/cell and cells/fire to make the magnitudes integer. Larger versions are in the supplementary material. **Figure 7.** Examples of fire occurrence prediction maps that use the categorization and colouring in Table 6: (**a**) human-caused; (**b**) lightning-caused; (**c**) total. The legend uses units of fires/cell and cells/fire to make the magnitudes integer. Larger versions are in the Supplementary Material.

Note that a sequential scale is logical for FOP because any non-zero FOP is "bad" in this context. But the calming colours are assigned to very low magnitudes, and this provides a slightly greater distinction between the remaining colours. A much greater distinction could easily be achieved by adding more colours—for example, magenta– purple–black for the highest categories, which would also draw more attention to the critical extremes. This, however, would not accommodate colour vision deficiencies.

Broad adjective categories were added to emulate the familiar four- or five-category pattern and simplify interpretation. Those boundaries fit the general association between Ontario's fire arrival density and fire situation severity. The three-significant-digit category boundaries were replaced in the legend by integers for selected category midpoints or boundaries (Table 6). Note that the units change from fires/cell to cells/fire to show integer magnitudes, which are far more meaningful than fractions.

#### *3.6. Case Study—Step 5: Evaluating and Revising the Scheme*

Several significant revisions were done in arriving at the current design in Figure 7. The first lightning-caused FOP maps had a linear scale with four equal-interval categories using traditional blue–green–yellow–red. When human-caused FOP was added, the scale was changed to 10 categories using a smooth green–yellow–orange–red transitional gradient and later to a special blue–yellow–orange–red to accommodate some colour vision deficiencies. As stated above, the inconsistent subjective scales were replaced completely in early 2020 per our method. The maps originally showed the expected number of fires/cell but were changed to show the density of the expected number of fires/unit area because of fractional cells. Note that the map legends intentionally omit this complication; density seems to be the automatic, intuitive interpretation. Additional revisions are planned, and further evaluation is always ongoing.

#### **4. Discussion**

While implementing the FOP models for operational evaluation, we arrived at and applied this method, which may seem well structured and straightforward now, but it was far from that during the design work. The method emerged as a by-product, having evolved during iterative deliberations. As illustrated in the introduction and identified in [9], the categorization of the FWI System outputs could benefit from this approach. We have since applied the method to categorize and colour outputs from other models for operational display [39–41], examples of which are in the Supplementary Material.

The many considerations discussed in this paper need to be addressed to ensure that the models are interpreted and used appropriately. We emphasize that any schemes including those built into software applications need this careful consideration. Incorporating the outputs from scientific models into operational decision-making is not straightforward. The key for a successful application is that researchers and practitioners work together closely throughout the process, from problem identification through to implementation and evaluation [1]. Through this collaborative approach, outcomes tend to have a higher acceptance and usefulness.

An additional factor not addressed by our method is that categories need to be meaningful for more than just map design, because there is a tendency to extend the use categories to guidelines and standard operating procedures and vice versa [9]. The classifications and their boundaries can also be used as unwritten mental shortcuts or heuristics in the place of a more deliberative consideration of complex information. Well-designed, science-based classifications can be consistent with less time-constrained situational analyses, while poor classifications may lead to suboptimal decision-making.

Design considerations for presenting quantitative data for decision support go beyond the categorization and colouring of numerical scales. Also important are options for spatial resolution and the use of simulated three-dimensional displays, examples of which are given in the Supplementary Material [42].

We presented a method of designing classification and colouring schemes that consider many complex factors, interactions and trade-offs. These satisfied the ultimate goals of showing decision-makers complete and accurate quantitative information that was understood quickly and easily.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/fire4030050/s1: other applications of the model (Figures S1 and S2), classification of numbers, additional considerations for the spatial display of data (Figures S3–S5) and larger versions of Figure 7a–c (Figures S6–S8).

**Author Contributions:** Conceptualization: D.B., J.E., C.C.H., C.B.M., D.G.W. and M.W. (Mike Wotton); methodology: D.B., C.B.M. and D.G.W.; software: J.E., A.S. and D.G.W.; formal analysis: D.B., C.C.H., C.B.M. and M.W. (Melanie Wheatley), D.G.W. and M.W. (Mike Wotton); investigation: D.B., J.E., C.C.H., C.B.M., A.S. and M.W. (Melanie Wheatley), D.G.W. and M.W. (Mike Wotton); writing—original draft preparation: D.B., C.C.H., C.B.M., M.W. (Melanie Wheatley) and D.G.W.; writing—review and editing: D.B., J.E., C.C.H., C.B.M., A.S., M.W. (Mike Wotton), D.G.W. and M.W. (Melanie Wheatley); visualization: D.B., J.E., C.C.H., C.B.M., M.W. (Melanie Wheatley), D.G.W., M.W. (Mike Wotton) and A.S.; project administration: D.B., C.B.M. and D.G.W. All authors have read and agreed to the published version of the manuscript.

**Funding:** We acknowledge the support of the Natural Sciences and Engineering Research Council of Canada (NSERC) [RGPIN-2015–04221].

**Acknowledgments:** Thank you to the Ministry of Northern Development, Mines, Natural Resources and Forestry, Aviation Forest Fire and Emergency Services staff for important discussions and management for support. Thank you to Barry Graham, Dan Johnston, Dan Leonard and David Martell for insightful advice, Darren McLarty for the review, Jim Caputo, Natasha Jurko and Benito Russo for technical support and the Fire Intelligence Officer cadre for their regular engagement and feedback. A special thank you goes to Darryl Pajunen for supporting the model's development. We thank the anonymous reviewers for their careful reading and helpful comments that improved the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Ujjwal K. C. 1,**<sup>∗</sup> **, Jagannath Aryal 2 , James Hilton <sup>3</sup> and Saurabh Garg 1**


**Abstract:** Rapid estimates of the risk from potential wildfires are necessary for operational management and mitigation efforts. Computational models can provide risk metrics, but are typically deterministic and may neglect uncertainties inherent in factors driving the fire. Modeling these uncertainties can more accurately predict risks associated with a particular wildfire, but requires a large number of simulations with a corresponding increase in required computational time. Surrogate models provide a means to rapidly estimate the outcome of a particular model based on implicit uncertainties within the model and are very computationally efficient. In this paper, we detail the development of a surrogate model for the growth of a wildfire based on initial meteorological conditions: temperature, relative humidity, and wind speed. Multiple simulated fires under different conditions are used to develop the surrogate model based on the relationship between the area burnt by the fire and each meteorological variable. The results from nine bio-regions in Tasmania show that the surrogate model can closely represent the change in the size of a wildfire over time. The model could be used for a rapid initial estimate of likely fire risk for operational wildfire management.

**Keywords:** fire spread models; surrogate modeling; sensitivity analysis; global sensitivity analysis

#### **1. Introduction**

Wildfires are a major threat in fire-prone areas such as Australia, Mediterranean regions in Europe, and the United States. These can destroy homes and infrastructure causing millions of dollars of damage [1,2] as well as threaten lives and ecologically sensitive areas. Consequently, fire suppression and mitigation costs have significantly increased over the years [3]. For disasters such as wildfires, where the fire conditions rapidly change with time, any models that can give realistic fire predictions in as little time as possible can be critical for effective preparedness, early warning, fire suppression, and evacuation efforts. Computational wildfire simulations using different fire propagation models are an effective means of predicting the wildfire behaviors and risk [4]. However, such models are typically deterministic and may neglect uncertainties inherent in factors driving the fire [5]. Modeling these uncertainties can more accurately predict risks associated with a particular wildfire, but require a large number of simulations with a corresponding increase in computation time. Under operational constraints, this computation time may be longer than the time window available to prepare and respond to the wildfire. Consequently, computationally cheaper models that can rapidly estimate risks could be very useful for operational wildfire management.

Wildfire behavior models typically fall into one of three categories [6]: those based on physical processes [7–10], empirical models [11–13], or simulation and mathematical models [14–16]. Mathematical concepts like Ordinary Differential Equations (ODEs), polynomial chaos, Gaussian process, and spatial process-based models like Cellular Automata, and Level Set methods have been used to model the spread of wildfires [17–22]. Currently,

**Citation:** K.C., U.; Aryal, J.; Hilton, J.; Garg, S. A Surrogate Model for Rapidly Assessing the Size of a Wildfire over Time. *Fire* **2021**, *4*, 20. https://doi.org/10.3390/fire4020020

Academic Editor: James R. Meldrum

Received: 30 March 2021 Accepted: 22 April 2021 Published: 23 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

operational fire behavior tools such as Spark [14], FARSITE [15], and Phoenix [23] are used to predict the wildfire behavior. These tools predict fire behavior by simulating the complex relationships between different contributing factors in a one-to-one deterministic manner without incorporating uncertainty in the prediction results [24]. In reality, there are uncertainties for each of the contributing factors which consequently introduce uncertainty in the resulting model output.

Quantifying such uncertainties in operational fire models requires identifying, classifying, and assessing their sources and influence on the model output [25]. Kaschek et al. [26] highlighted the importance of parameter value estimation for the experimental data in understanding the model dynamics for any permissible value of the parameters and subsequent uncertainty quantification. However, such uncertainties have not been well quantified in current operational fire models. As such, Sensitivity Analysis (SA) is a currently active area of research, in which the spatial and temporal variability in the model outputs is quantified as a function of the variation in the input parameters for an optimized model [27]. Note that uncertainty and sensitivity analysis are different as uncertainty analysis assesses the uncertainty in the model output due to uncertainties in input parameters while sensitivity analysis assesses the contribution of input parameters on the uncertainty of the model output [28]. Sensitivity analyses of input parameters for different fire models have been carried out in several studies [5,29–33], where the model outputs have been considered for a fixed period. Sensitivity analyses may also help quantify such uncertainties and determine the influence of each input parameter on fire progression but such studies for operational fire models, to our knowledge, have not been well-explored.

On the other hand, the model performance and the effectiveness of risk management achieved with such models is determined to a large degree by how well such uncertainties are understood and communicated [34]. Therefore, calculating more accurate risk metrics quantifying all these uncertainties can be a complex task [24,29]. Risk calculation can require a large number of simulations to be run under different scenarios, which can take several hours to days to complete [35]. Such time-consuming analyses are currently impractical for operational management. Consequently, developing computationally efficient methods, such as surrogate models, can be an alternative to such models for operational wildfire management.

Surrogate models, also known as metamodels, are the models of the outcomes that mimic the behavior of the simulation model as closely as possible while being computationally inexpensive. Surrogate models have been extensively used to replace, or as an alternative to, computationally expensive models, as they represent input–output relationships derived from complex models. These have been used in natural hazard models including flooding [36,37] and storms [38,39]) as well as engineering design problems [40–42] to rapidly estimate the associated risks and concerns. The surrogate model introduced in this study can facilitate rapid estimation of wildfire risks based on the data from high-fidelity operational fire simulation models.

Several studies have shown that not all the parameters in wildfire models have significant influence on the spread rate, and consequently fire behaviors can be explained by prioritizing highly significant parameters [43,44]. Sharples et al. [44] demonstrated that the spread of wildfires could be expressed in a "universal" fire spread model based on only temperature, relative humidity, and wind speed input parameters. Taking this as a starting point, we use these input parameters to construct a computationally efficient surrogate model. Next, we present a one-at-a-time (OAT) sensitivity analysis of input parameters for the variability of fire dynamics in the wildfire prediction tool Spark, whereby the influence of a parameter on the model output is assessed one at a time instead of multiple parameters simultaneously. Based on the results obtained from the analyses, the nature of the contribution of the parameters is determined and used to define the functional form for the surrogate model. Finally, we use a data fitting technique to determine the values of unknown factors in the functional form to define a complete surrogate model to explain fire behaviors based on initial weather conditions. We test the efficacy of the approach by applying it to the entire Tasmanian region and its nine different bio-regions. The specific contributions of this study are as follows:


The rest of the paper is organized as follows. Section 2 explains the general methods for the mathematical formulation of the surrogate model and time-based sensitivity analysis. Section 3 discusses the investigative results obtained in our experiments. Section 4 concludes the paper and discusses possible future work.

#### **2. Methodology**

#### *2.1. Study Area*

We chose Tasmania for our study due to the frequent occurrences of wildfires in the region, high-quality land data sets as maintained by Tasmania Fire Service (TFS) and State Emergency Service (SES), and a well-studied and systematic grid configuration for possible fire start locations in the region. Tasmania has a total area of 68,401 sq. km. The grid configuration for fire simulation, as maintained by TFS, places each possible fire start location at an interval of 1 sq. km irrespective of land classification. Any start locations falling on water bodies are shifted to the nearest land locations. There is a total of 68,048 possible fire start locations in the grid configuration for the entire Tasmanian region. We followed Interim Biogeographic Regionalization for Australia Version 7 (IBRA 7) [45] for nine different bio-regions (see Figure 1) to run fire simulations.

**Figure 1.** Nine bio-regions of Tasmania. Sourced from the work in [46].

#### *2.2. Fire Simulations—Spark*

Spark [14] is a computational wildfire modeling framework that predicts the spread and potential impact of wildfires based on different behavior fire models. The framework simulates the progression of the fire over time in different fire conditions and fuels, resulting in different rates-of-spread (RoS). The fire simulations in Spark require many input data sets,

including the information on land classification, fuel type, topography, and meteorological data. The simulations can be run for any number of distinct fire perimeters and can model factors such as firebreaks, spot fires, and coalescence between different parts of the fire over time. The simulations use several different empirical fire models for fuels found in Tasmania. Different fuel types found in Tasmania are mapped from the TasVeg [47] vegetation types to several empirical fire spread models. These include models for eucalypt forest [48,49], buttongrass moorland [50], heathland [51], and grasslands [52]. A detailed description of these models is included in the Appendix A.

The models require meteorological data inputs. For this study, only temperature, relative humidity, and wind speed were used [44]. Other input parameters in the fire simulation such as fuel load and land topography were used as per the configuration and records maintained by TFS and Tasmanian Government [53]. The fuel load was based on a Tasmania fire history data set allowing the time since the fire to be calculated and populated using Olson fuel load curves [54–56]. Our current study covered only the westerly winds as they are the most common winds in Tasmania [57,58]. The configuration and input files can be made available from the authors upon request. Based on the historical records and as considered in [5], the ranges for the data inputs are fixed as follows:


The fire simulations were started at locations in each of the nine bio-regions ensuring minimal obstructions (lakes/water bodies) irrespective of land classification locations, as shown in Figure 2. The simulations were run for a total of five hours and the area of each simulation was recorded every half-hour. The simulation data used in this study are a subset of the data available in [59].

**Figure 2.** Fire start locations for different bio-regions in Tasmania. The locations are chosen ensuring minimal obstructions due to water bodies irrespective of land classification.

#### *2.3. Sensitivity Analysis*

To quantify the interactions of meteorological inputs on fire dynamics, we focus on the parametric sensitivity analysis of fire simulations over time. For SA, we follow an OAT approach where the influence of a parameter is assessed by running the fire simulations for different values of the parameter within the permissible range while keeping the other parameters constant as explained in the work [60]. The influence of a parameter on the fire area is assessed by comparing the fire values obtained for the maximum and minimum values of the parameter (with minimum and maximum values of other parameters, respectively) against the fire areas obtained for the maximum ( *all-high*) and minimum values (*all-low*) of all the parameters. The variability of the output of the fire simulations can be analyzed at different time steps by considering the model outputs in all these time steps to determine the scale of influence of that input parameter on the fire area. However, for our analysis to determine the cumulative influence of each parameter on the fire size, we use the total fire size obtained at the end of the total simulation time. The quantification of the variability in the fire area is carried out for each bio-region and the entirety of Tasmania as well. Similar analyses can be done to understand the influence of each parameter on the fire spread during a particular duration of time by assessing the values of the fire area.

#### *2.4. Surrogate Modeling*

For surrogate modeling of fire simulations, we first establish a mathematical model based on how fires spread over time when they ignite under some initial environmental conditions. Then, unknown values in the mathematical model are determined by fitting the simulation data by considering the results obtained from the sensitivity analysis. The mathematical relationship and simulation data fitting are discussed in detail as follows.

#### 2.4.1. Mathematical Foundation

For a fire starting at a particular location, the shape of the fire is ideally assumed to be elliptical [61]. Consequently, the total area burnt by the fire, represented by *A*, is directly proportional to the square of the time (*t*) for which the fire is burning [61]. We focus our mathematical foundation on the same principle but replace the squared power of the time with a variable *n* to account for any dependency on fire geometry (ideally, the value of *n* in our mathematical foundation must be less than or equal to 2). As such, we mathematically express *A* as

$$A = \mathfrak{a}t^n\tag{1}$$

where *α* is a fitting factor that depends on the initial values of temperature (*T*0), relative humidity (*R*0), and wind speed (*W*0). This fitting factor does not account for factors such as daily temperature changes or wind changes during the progression of the fire. However, such factors are accounted for by analyzing the fire sizes for different sets (original and changed) of initial conditions and assessing the difference. The value of *α* is determined by defining a functional form for the three input parameters based on the results of sensitivity analysis. The values of *n* and *α* are determined by fitting simulation data.

#### 2.4.2. Simulation Data Fitting

The unknown factors in Equation (1) can be determined by fitting a curve through experimental analyses. We use nonlinear least squares method [62] for data fitting. For experimental analyses, the fire simulations are run for five hours over a set of points in all nine different bio-regions in Tasmania. The simulation of fire behavior in Spark is graphically represented in Figure 3, which shows that the propagation of the fire over time can be different for different bio-regions. The steps for surrogate modeling are carried out initially for the entire Tasmanian region and then to specific bio-regions one at a time.

#### (**a**) Southern Ranges (**b**) West

**Figure 3.** Illustration of simulation in Fire dynamics in Spark (*T* = 40, *R* = 90, *W* = 60). The wildfire propagates differently based on the fuel load and land topography and consequently, the total area burnt by the fire varies based on the bio-region. Different colors represent the area burnt by the fire over time in half hourly steps. The purple color represents the area burnt by the fire in the first half hour, while the red color represents the area burnt from 4.5–5 h mark.

#### *2.5. Surrogate Model Validation*

For the validation of the dynamic model constructed in this study, we use the hold-out cross-validation method as explained in [63]. This method splits the simulation data into 70% training set and 30% test set for cross-validation. Then, to quantify the accuracy of the dynamic model, Pearson's correlation coefficient [64] is calculated between the model predicted values and the actual values obtained from the simulations. The value of the correlation coefficient closer to one (1) represents a more accurate model.

#### **3. Results and Discussions**

#### *3.1. Sensitivity Analysis*

We analyzed the effect of each input parameter on the fire area size at different instants of time as given by the fire simulations in Spark run over different bio-regions in Tasmania. A plot of the variation in fire area for different bio-regions caused by variability in temperature (◦C) is shown in Figure 4a. The fire area generally grows bigger over time for the maximum value of the temperature. In the analysis, the fire burns a greater area at *T* = 40 ◦C and minimum values of other parameters. For *T* = 40 ◦C, the increase in the area burnt by the fire is steady until 2 h, but after that, the increase is exponential to reach a maximum value of 4121.190 ha at the five-hour mark. For a lower value of T, the increase in the fire area over time is steady and close to a linear relationship. For the entire Tasmania, the variability in burnt fire area caused by the variation of temperature is as high as 4115.340 ha. The respective variations in the fire area for different bio-regions are listed in Table 1.

The variation of the fire area with the variability of relative humidity (%) is shown in Figure 4b. It is clear from the plots that relative humidity has a significant influence on the variation of the fire area. The maximum values of the area burnt by the fire are achieved with the value of relative humidity at 10% and maximum values of other parameters for all the locations. For this particular combination of parameter values, the increase in the fire area size with time is brisk. There is a sharp increase in the rate of increase in the fire area after 2 h in general for most of the bio-regions. The maximum area burnt by fire over 5 h is 28,682.6 ha. The variation in the fire area caused by the variability of relative humidity for the entire Tasmanian region is as high as 28,680.08 ha, which is significantly higher and signifies the greater impact of the parameter. The respective variations in the fire area for different bio-regions are listed in Table 1.

The variation in the fire area caused by the variability in the wind speed (kmh−<sup>1</sup> ) is shown in Figure 4c. The area burnt by the fire is the greatest at *W* = 60 kmh−<sup>1</sup> and minimum values of other parameters for all the locations. The fire areas for other cases are

relatively low. The increase in the fire area with time is steady until 2 h after which the fire area starts to increase at a higher rate. The area burnt by the fire over 5 h reaches 11,543.2 ha. The variation in the fire area caused by the variability of wind speed is approximately 11,539.96 ha for the entirety of Tasmania, which indicates the significant influence of the wind on the fire area. The respective variations in the fire area for different bio-regions are listed in Table 1.

**Figure 4.** Variation of fire area (in *ha*) with the variability of different input parameters, shown as ∆*A* on the plots. *All high* is the plot for fire area obtained with maximum values of all parameters, while *all low* is with minimum values. The variation of fire area caused by relative humidity is the highest while that of the temperature is the lowest.

In the analyses, we derived the upper and lower bounds for the values of the fire area for each input parameter considering the values of other parameters as well. The minimum value of an input parameter did not necessarily signify the lower bound for the fire area given by the parameter, and vice versa. For the entire Tasmanian region, the variability in fire area brought about by extreme values of temperature for extreme values of other parameters is ~4115.340 ha in five hours (see Figure 4a). The variability in fire area in the same duration caused by relative humidity is ~28,680.08 ha (see Figure 4a), while the variability caused by wind speed stands at 11,539.96 ha (see Figure 4a). As such, the extent of variability in fire area caused by relative humidity and wind speed is approximately 7 and 3 times greater than the same caused by temperature, respectively. Thus, it can be concluded that relative humidity has the highest effect on the area burnt by the fire at different instants of time for the entire Tasmanian region given by the fire simulations by Spark while the temperature has the least effect. For Northern Slopes, the variations in fire area caused by temperature and wind speed are comparable, which points out the impact caused by fuel load and land topography. Similar analyses can be done to understand the influence of each parameter on the fire spread during a particular duration of time. For example, during the second hour of a wildfire (from 1 to 2-h mark), the variation caused in the fire area by the variability in the values of temperature, relative humidity, and wind speed in the entire Tasmanian region are 739.26 ha, 5378.58 ha, and 2095.47 ha, respectively. Such information about the influence of the input parameters on the model output with consideration of the time can lead to a better quantification of the feedback of meteorological inputs on the fire dynamics.

We used an OAT method to measure the influence of each input parameter on the output of the fire simulations. To understand the variation in fire area in different bioregions, we carried out the analyses for all the regions independently. In our analyses, we derived the upper and lower bounds for the values of the fire area for each input parameter considering the values of other parameters as well. For all the bio-regions, the variation in the fire area caused by the variability of relative humidity is the highest, while the variation caused by that of the temperature is the lowest. The variations in fire area caused by the variability of relative humidity and wind speed are as high as about 11 and 5 times more when compared to the same caused by the variability of temperature, respectively. For Northern Slopes, the variations in fire area caused by temperature and wind speed are comparable, which points out the impact caused by fuel load and land topography. The quantitative measure of the variability in fire area brought by each parameter during the OAT method is conclusive. Our findings established relative humidity as the parameter with the greatest impact on the spread of the fire over time and temperature as the one with the least impact. This finding can be attributed to the fire models in Spark that are highly sensitive to fuel moisture content, which in turn is highly sensitive to relative humidity compared to temperature (for example, the Dry Eucalypt model equation for fuel moisture content, *MC*, has a greater dependency on relative humidity than temperature). However, the extent of the influence of each parameter on the variability of the fire area size is not necessarily the same for all the bio-regions even though they have similar relative comparisons.


**Table 1.** Variability in fire area caused by variability in input parameters.

#### *3.2. Surrogate Modeling*

Following the methods defined for surrogate modeling (Equation (1)), we first incorporated the meteorological influence and adapted a functional form for temperature, relative humidity, and wind speed to determine the value of *α*. The functional form was developed by further analyzing the results of the sensitivity analysis in which the growth of fire area was investigated at different time instants (every half-hour mark) against the highest values of the input parameters (Figure 5). As can be seen in the figure, the growth

of the fire area is close to linear over time for temperature and relative humidity, whereas an exponential fit is more suitable for wind speed. As such, the functional form for input parameter to define the factor in our surrogate model (Equation (1)) is

$$\alpha = k\_1 T\_0 + k\_2 R\_0 + e^{k\_3} W\_0 \tag{2}$$

After fitting the simulation data obtained from all the Tasmanian regions, the surrogate model is defined as follows:

$$A(t) = \left(16.72T\_0 - 26.80R\_0 + e^{2.03}W\_0\right) \text{ } t^{1.76} \tag{3}$$

Note that the functional form used for the surrogate model is slightly different from the formulation in the work of Sharples et al. [44] as the surrogate model is based on the findings from the uncertainty quantification.

**Figure 5.** Rate of increase of fire area with time for Southern Ranges bio-region (The rate of increase in the fire area is in hectares per minute (ha min −1 ). The contributions in the rate of change of fire area made by temperature (blue line) and relative humidity (red) are steady over time while the same made by the wind speed (green) increases exponentially over time).

The simulation and predicted data are shown in Figure 6. The value of the correlation coefficient is 0.897. Interestingly, the value of power to which the time is raised in the model is close to 2 (1.76), indicating fire growth scales approximately as the expected square of the time. Despite the high value of the correlation coefficient, the model over fits the fire area for some extreme values of the input parameters. Although not ideal, an overfitted result produces a more conservative estimate of larger fire sizes, rather than a possibly dangerous underestimate. The overfitting in the model is possibly due to the significant impact of the fuel loads, which can be different for different bio-regions, and the interaction between the input parameters, which are not considered in the study. The fitting parameters for each individual region in Tasmania are given in Table 2 and plotted in Figure 7.

**Figure 6.** Simulation Data vs. model predicted data for entire Tasmania (The area burnt by the fire is in hectares (*ha*). The red line represents the ideal line where the scatter plots should align for the most accurate surrogate prediction data. The calculated value of correlation coefficient is 0.897, which indicates that the surrogate model closely predicted the actual simulation data.)


**Table 2.** Fitting parameters in the surrogate model for nine bio-regions in Tasmania.

**Figure 7.** Bio-region-specific surrogate model predicted data compared against simulation data. The fire area is in hectares (*ha*). The red line represents the ideal line where the scatter plots should align for the most accurate surrogate prediction data. The surrogate model constructed for the Southern Ranges bio-region has the highest value of correlation coefficient of 0.996 while that of Northern Slopes is the lowest with a value of 0.910.

Our surrogate model based on fire simulations has reasonable success with promising findings. The mathematical equation defined in the surrogate model in terms of the initial values of different meteorological data gave satisfactory information about the progression of the fire over time as the value of the calculated correlation coefficient for the entire Tasmanian region is about 0.897. This finding indicates the fact that constructing a surrogate model to represent all the bio-regions with a diverse range of fuel loads and land topography may not be so efficient. Consequently, we constructed a single surrogate model for each bio-region and calculated the correlation coefficient between the simulation data and the model predicted data, as shown in Figure 7. Moreover, as listed in Table 2, the values of unknown parameters also represent the respective influence of the parameters in the forest area burnt by the fire. The findings of the surrogate models are consistent with the results obtained from our sensitivity analysis carried out by considering different instants of time. The negative values of *k*<sup>2</sup> indicate that the fire grows rapidly under lower values of relative humidity when compared to the higher values, while the higher values of *k*<sup>2</sup> reflect the strong influence of relative humidity on the fire area. The negative value of *k*<sup>2</sup> as determined in our study is quite consistent with the functional form defined by Shaples et al. in [44] for spread index, where the relative humidity in the denominator indicates a higher spread rate (consequently higher fire area) at lower values when compared to

the same at higher values. Further, the surrogate model constructed to represent the dynamics of fire in terms of initial environmental conditions closely represents the actual behavior of the fire over time as the values of the correlation coefficient are well over 0.987. The surrogate model incorrectly identifies small fires, possibly due to the presence of a lake or other water bodies near the fire start location. However, the model shows a good match for larger fires in most cases. This is useful for operational fire management, as recent studies have shown the opposite trend for fire simulations [65], where fire predictions using forecast weather in fire simulation tools under-predicted the fire size in extreme weather conditions. Considering other location-specific information in the model can improve the closeness of the model estimation with the actual values. For the **Southern Ranges** bio-region, the model has a correlation coefficient of 0.996, which signifies the fact that the model closely matches the fire dynamics from the simulations in the region. Such a good prediction by the model appears to be due to the region being relatively free of obstructions (lakes/built-up areas) where the fire would otherwise be stopped. This finding indicates that surrogate models, like the one built here in the study, can be helpful for quickly estimating fires under extreme fire weather conditions for large open areas. The lowest value of the correlation coefficient for the model is 0.910 in the **Northern Slopes** region, which is likely due to the presence of wetlands in the region and due to nonlinear interactions between the input parameters. The nonlinear interactions between the input parameters are not considered for this study. This particular aspect of the surrogate model will be strengthened in the future to improve the closeness of the data predicted by the model with the actual simulation data.

The constructed surrogate model can represent the change in the size of the fire with time-based on the initial values of meteorological inputs (temperature, relative humidity, and wind speed). This could be used to rapidly estimate the wildfire risks, although it should be emphasized that the model can only give an estimate of the size of the fire rather than the shape, growth, or possible impact on certain areas. The functional form to estimate the fire size obtained in this study can be an alternative form to the one explained in [44]. Additionally, our attempt to investigate the path of fire growth has shown that in practical scenarios, the fire area does not follow an ideal *t* 2 relationship as the values of *n* as obtained in our experimental analyses for all nine bio-regions are less than 2.

#### **4. Conclusions and Future Works**

In this paper, we analyzed the contribution of different parameters on the area (geographical extent) burnt by the fire over time in a conventional one-at-a-time (OAT) fashion by considering nine different bio-regions in Tasmania. Relative Humidity was found to be the parameter with the greatest impact on the spread of fire over time, while Temperature was found to have the least effect. Moreover, our analyses were time-dependent, allowing new ways to quantify the feedback of meteorological inputs on the fire dynamics, and vice versa. Information about the sensitivity of parameters in wildfire models could facilitate new operational approaches for risk management by excluding less sensitive parameters from the models allowing risk to be calculated in a shorter time-frame. We used simulations to construct a computationally-efficient surrogate model for the area of a fire over time in terms of initial starting conditions. The surrogate model considered the nature and the extent of the influence of each meteorological input on the fire dynamics. The region-specific surrogate models constructed in this study represented the dynamic fire behaviors except for some extreme cases. Such surrogate models may be helpful in operational wildfire management during emergencies to quickly make informed decisions for response activities by providing a rapid estimate of potential fire size.

This study is one of the preliminary works to explore sensitivity analyses at different instants of time and how results of such analyses can be integrated into natural hazard modeling systems. In the future, we will take the nonlinear interactions between the input parameters into account to strengthen the surrogate model. The study can be extended for other regions as well with sufficient fire simulation data in hand.

**Author Contributions:** Conceptualization, U.K. and J.A.; Formal analysis, U.K.; Investigation, U.K.; Methodology, U.K. and J.A.; Software, J.H.; Supervision, J.A., J.H. and S.G.; Validation, U.K. and J.H.; Visualization, U.K.; Writing—original draft, U.K. and J.H.; Writing—review and editing, U.K., J.A., J.H., and S.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank our academic colleagues who helped in improving the quality of the manuscript at various stages of the work.

**Conflicts of Interest:** The authors declare no conflicts of interest.

#### **Appendix A**

*Appendix A.1. McArthur Grassland Fire Danger Meters*

McArthur [48] described his research into grassland fire behavior in terms of the Grassland Fire Danger Index (GFDI), based on which the rate of spread was calculated. The equations for GFDI and rate of spread (R) are listed as follows:

$$\text{GFDI} = 2 \exp(-23.6 + 5.01lc(\mathbb{C}) + 0.0281T - 0.226\sqrt{RH} + 0.633\sqrt{\mathcal{U}\_{10}})$$

*C* is the degree of curing (%), *T* is the air temperature (*oC*), *RH* is the relative humidity (%), and *U*10 is the wind speed (km/h) measured at a height of 10-m in the open.

$$R = 0.13GFDI$$

*R* is the headfire rate of spread (km/h).

#### *Appendix A.2. Dry Eucalypt Model (Cheney et al.)*

This wildfire model is used for predicting the spread of fire behaviors in dry eucalypt forest. The model was developed from a sequence of experiments called "project Vesta" [49] carried out in south-western Australia. The mathematical equations for the model are listed as follows:

$$\mathcal{R} = \begin{cases} 30 \text{ } \phi M\_f & \text{, } \mathcal{U}\_{10} \le 5 \text{ km/h} \\ \left[ 30 + 1.531 (\mathcal{U}\_{10} - 5)^{0.858} FHS\_{\text{s}}^{0.93} (FHS\_{\text{ns}}H\_{\text{ns}})^{0.637} B\_1 \right] \phi M\_f & \text{, } \mathcal{U}\_{10} > 5 \text{ km/h} \end{cases}$$

*Note that the term B*<sup>1</sup> *is the model correction for bias and is taken as 1.03 for this. FHSns is the near-surface fuel hazard score, FHSs is the surface fuel hazard score, and Hns is near-surface height and are derived from fuel age using the tall shrubs regression equations as explained in [66]*.

*R* is the rate of spread (m/h), *U*<sup>10</sup> is the average 10-m open wind speed (km/h), and *φM<sup>f</sup>* is the fuel moisture function.

$$\phi M\_f = 18.35 \,\,\,\text{MC}^{-1.495}\_{-},$$

*MC* is the moisture content (%) and taken from [67] and [68].

$$\text{MC} = 2.76 + 0.124 \text{ } RH - 0.0187 \text{ } T$$

*T* is the air temperature (*oC*) and *RH* is the relative humidity (%).

#### *Appendix A.3. CSIRO Grassland Model*

Based on the experimental results obtained from the burning project in the Northern Territory of Australia, Cheney et al. [52] described a relationship for the rate of fire spread on grasslands. This model was developed to settle down the confusion created after the calculation of different Grassland Fire Danger Index (GFDI) values under the same conditions thereby questioning the actual effect of fuel load on fire spread rate. The mathematical equations in the model are explained as follows.

$$R\_{c\mu} = \frac{1}{3.6} \begin{cases} (0.054 + 0.209 \text{ }\text{U}\_{10}) & \phi M \phi \text{C} \\ (1.1 + 0.715 \text{ }(\text{U}\_{10} - 5) ^{0.844}) & \phi M \phi \text{C} \end{cases}, \text{\$\mathcal{U}\_{10} > 5\$ km/h\$}$$

*Rcu* is the cut/grazed rate of spreed in m/s, *U*10 is the 10-m open wind speed (km/h), *φM* is the fuel moisture coefficient, and *φC* is the curing coefficient.

$$\phi M = \begin{cases} \exp(-0.108 \text{ MC}) & \text{, } \text{MC} < 12\% \\ 0.684 - 0.0342 \text{ MC} & \text{, } \text{MC} \ge 12\% \end{cases}, \text{\textdegree C} \ge 12\% \text{, } \text{\textdegree U}\_{10} < 10 \text{ km/h} $$

*MC* is the dead fuel moisture content (% oven-dry weight basis).

$$\text{MC} = 9.58 - 0.205 \,\,\, T + 0.138 \,\, \,\, RH$$

*T* is the temperature (◦C and *RH* is the relative humidity (%).

$$\phi \mathbb{C} = \begin{cases} \frac{1.036}{1 + 103.99 \exp(-0.0996(\mathbb{C} - 20))} & \text{, } \mathbb{C} > 20\%\\ 0 & \text{, } \mathbb{C} \le 20 \end{cases}$$

*C* is the degree of grass curing (%).

*Note that the original equation for φC proposed by Cheney et al. was modified by Cruz et al. as the curing level can be as low as 20% and the damping effects of the fuel in grassland is less.*

#### *Appendix A.4. Marsden–Smedley and Catchpole Buttongrass Model*

Marsden–Smedly and Catchpole [50] described the fire behaviors for buttongrass moorlands in Tasmania for fire danger rating and fire behavior prediction. The mathematical equations for the models are explained as follows:

$$R = 0.678 \text{U}\_2^{1.312} \exp(-0.0243 \text{MC}) (1 - \exp(-0.116 \text{AGE})) $$

whee *U*<sup>2</sup> is the wind speed (km/h) measured at a 2-m height, *MC* is the dead fuel moisture content (%), and *AGE* is the time since the last fire (years).

$$\text{MC} = \exp(1.66 + 0.0214RH - 0.0292T\_{\text{dew}})$$

*Tdew* is the dew point temperature (*oC*) and *RH* is the relative humidity (%).

#### *Appendix A.5. Anderson et al. Heathland Model*

Anderson et al. [51] used a dataset that covered a wider range of heathland and shrubland species to develop a model for healthland as an extension to the work [69]. The equations for the model are given as follows:

$$R = \begin{cases} [R\_\vartheta + 0.2(5.67(5WF)^0.91 - R\_\vartheta) \mathcal{U}\_{10}] H^{0.22} \exp(-0.076 \text{MC}) & \mathcal{U}\_{10} < 5 \\ 5.67(WF \mathcal{U}\_{10}^0.91) H^{0.22} \exp(-0.076 \text{MC}) & \mathcal{U}\_{10} \ge 5 \end{cases}$$

where *R* is the rate of spread (m/min) with vegetation height and without live fuel moisture content, *U*<sup>10</sup> is the 10-m open wind speed (km/h), *H* is the average vegetation height (m), *WF* is the wind adjustment factor, and *MC* is the dead fuel moisture content (%).

$$\text{MC} = 4.37 + 0.161RH - 0.1(T - 25) - \Delta 0.027RH$$

∆ = 1 for sunny days from 12:00 to 17:00 from October to March and 0 elsewhere. *T* is the ambient air temperature (◦C and *RH* is the relative humidity (%).

#### **References**


*Article*
