Effects of spatial and temporal variation in environmental conditions on simulation of wildfire spread

Environmental Modelling & Software 67 (2015) 118e127

Contents lists available at ScienceDirect

Environmental Modelling & Software journal homepage: www.elsevier.com/locate/envsoft

Effects of spatial and temporal variation in environmental conditions on simulation of wildfire spread J.E. Hilton a, *, C. Miller a, A.L. Sullivan b, C. Rucinski a a b

CSIRO Digital Productivity Flagship, Clayton South, VIC 3169, Australia CSIRO Land and Water Flagship, Black Mountain, ACT 2601, Australia

a r t i c l e i n f o

a b s t r a c t

Article history: Received 13 November 2014 Received in revised form 23 January 2015 Accepted 23 January 2015 Available online

Environmental conditions, such as fuel load and moisture levels, can influence the behaviour of wildfires. These factors are subject to natural small-scale variation which is usually spatially or temporally averaged for modelling fire propagation. The effect of including this variation in propagation models has not previously been fully examined or quantified. We investigate the effects of incorporating three types of variation on the shape and rate of propagation of a fire perimeter: variation in combustion conditions, wind direction and wind speed. We find that increasing the variation of combustion condition decreases the overall rate of propagation. An analytical model, based on the harmonic mean, is presented to explain this behaviour. Variation in wind direction is found to cause the development of rounded flanks due to cumulative chance of outward fluctuations at the sides of the perimeter. Our findings may be used to develop improved models for fire spread prediction. © 2015 Published by Elsevier Ltd.

Keywords: Perimeter propagation Simulation Modelling Fire growth Level set Spark

1. Introduction Accurately predicting the arrival time of a fire front at a specific location is crucial for effective fire suppression, as well as the coordination of safe evacuation operations. Accurate determination of the likely extent and severity of potential wildfires can enable reliable risk assessments to be performed for vulnerable regions. Wildfires are, however, extremely complex systems and modelling their behaviour, even for the most simple scenarios, can be very challenging. The nature of the generally disparate and heterogenous environmental conditions (i.e. vegetation, weather and topography) in which wildfires spread further complicates the task of accurately predicting the behaviour and spread of wildfires. The physical processes governing a wildfire take place over a wide range of temporal and spatial scales (Sullivan, 2009a) from the combustion chemistry at the molecular scale (Sullivan and Ball, 2012), through to atmospheric interactions with the fire on the scale of kilometres (Sun et al., 2009). At spatial scales on the order of centimetres to tens of meters lie local factors governing the behaviour of the fire. These include the distribution of fuel type/ composition and conditions such as moisture content and

* Corresponding author. E-mail address: [email protected] (J.E. Hilton). http://dx.doi.org/10.1016/j.envsoft.2015.01.015 1364-8152/© 2015 Published by Elsevier Ltd.

availability (Anderson and Rothermel, 1965; Catchpole et al., 1989, 1993; Viegas et al., 2013). Atmospheric effects on the order of centimetres to tens of meters also influence fire behaviour, such as surface level wind turbulence (Luke and McArthur, 1978; Albini, 1982; Sullivan et al., 2012). The local topographic slope, and wind interaction with the slope, also plays a large role in fire propagation (Viegas, 2004; Sharples, 2008; Sullivan et al., 2014). The rate of spread of a wildfire is one of the most important parameters for operational purposes, as it can be used to predict arrival times. Other important parameters include the shape and intensity of the flame front as well as the propensity and behaviour of embers lofted into the air stream (fire spotting) which can lead to breakdown of suppression efforts (Ellis, 2011; Sullivan et al., 2012). The dispersal of smoke is also of interest from an environmental perspective (Goodrick et al., 2012). Due to the importance of predicting the arrival time of a fire front, many empirical expressions have been developed to give the one-dimensional rate of spread of the head of a wildfire, the fastest moving part of the perimeter. The Rothermel model (Rothermel, 1972) is one of the earliest examples of such empirical models still in widespread operational use as part of the US BEHAVE system (Andrews, 1986) which continues to be revised and expanded (Scott and Burgan, 2005; Andrews, 2014). In the time since Rothermel's work, a number of more complex empirical models have been introduced, often for specific vegetation types. Recent

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

examples include the Australian Dry Eucalypt Forest Fire Model (Cheney et al., 2012) and the CSIRO Grassland Fire Spread Model (Cheney et al., 1998; Sullivan, 2010). Dynamic spatial fire spread simulation has generally taken two approaches. The first is computational fluid dynamic (CFD) models, which attempt to replicate fire behaviour based on the fundamental combustion and heat transfer processes. The second is perimeter propagation models, which generally apply empirical rate of spread rules, such as the Rothermel expression, to simulate the propagation of the fire perimeter. CFD models are most often based on the compressible NaviereStokes equations with associated auxiliary relationships for factors such as chemical reactions and turbulence. Such models are less reliant upon extensive experimental relations for robustness but are very computationally expensive, can currently still only model a fraction of the complex processes occurring within wildfires, and are unsuited to operational use due to long computation times (Sullivan, 2009a). Perimeter propagation models attempt to model the large-scale propagation of fire across a landscape rather than directly solve the underlying physical relations governing the fire (Sullivan, 2009c). As such, they can be based on simplified versions of more complex physical models, use empirical relationships measured in the field, or be based on a mathematical rule set. The fire perimeter in these models is the interface between regions which are burnt or burning, and unburnt. These models are usually applied on the scale of tens to hundreds of kilometres. These models can be subdivided into front-tracking methods or cellular methods. In both types, the fire perimeter is represented as a two-dimensional interface, giving a considerable saving in computational cost over three-dimensional CFD models. In the front-tracking approach, the fire perimeter is described as a discretised set of line segments that expand according to a given rate of spread rule set. Each point on the fire perimeter is assumed to be a point source for future fire propagation (Knight and Coleman, 1993; Finney, 1994; Richards and Bryce, 1995; Coleman and Sullivan, 1996). The fire shape from these point sources is usually assumed to be elliptical, as ellipses have been found to provide a good approximation to fire perimeters in long-burning wildfires (Anderson et al., 1982; Richards, 1990). The geometry of the ellipse is determined using a one dimensional rate of spread from an empirical model, an empirical relation for the length-tobreadth ratio of the ellipse, and the duration of the simulation time step. The orientation of the ellipse is usually determined from the direction of the wind. The fire perimeter is updated as the maximum extent of all contributing fire ellipses for the time step. Although these models are very fast, limitations arise from the use of only one type of front shape. While elliptical templates provide a close match to many fire propagation scenarios, other template shapes can provide a better match under certain circumstances (Green et al., 1983). For example, different fuels can produce oval or tear-drop shaped fronts (Peet, 1967). Other limitations in such models include the assumption of static conditions at each point on the perimeter for the period of the time step, and the assumption of instantaneous steady-state motion of the fire perimeter from a point source ignition. A large number of topologydependent rules are often required in these models to resolve overlapping, twisting or colliding fronts (Knight and Coleman, 1993; Filippi et al., 2010; Bose et al., 2009). Models using this approach include SiroFire (Coleman and Sullivan, 1995, 1996), Phoenix RapidFire (Tolhurst et al., 2008), Prometheus (Tymstra et al., 2010), Aurora (Johnston et al., 2008) and FARSITE (Finney, 1998, 2004). Cellular methods discretise the domain into an underlying grid over which all input data is prescribed and all calculations are performed. Rule sets based on empirical or physical formula are

119

used to update the state of the grid over time. Such models include static raster implementations (Green et al., 1990), cellular automata models (Encinas et al., 2007; Achtemeier, 2013) and complex irregular stencil-based models (Trunfio et al., 2011). Examples of models using elliptical stencils include FireStation (Lopes et al., 2002), FIREMAP (Vasconcelos and Guertin, 1992), and PYROCART (Perry et al., 1999). A more recent approach, which is used in this study, is perimeter-growth based on the level set method (Sethian, 2001). In this method a local rate of spread can be applied at any point on the fire perimeter. Topological changes, such as breaking and merging of perimeters, are handled without the need for any specialised rule sets to handle colliding interfaces. The method is also not reliant on the application of any pre-defined templates, such as ellipses. Implementations of such methods appear to be in the early stages of development compared to cellular and front-tracking methods. Rehm and McDermott (2009) showed that level set simulation of ignition points on flat homogeneous terrain evolved naturally into an elliptical form, highlighting the potential of the method for realistic fire perimeter simulation. WRF-Fire, a coupled atmosphere-fire model (Mandel et al., 2011a; Coen et al., 2013), combines the Weather Research and Forecasting (WRF) atmospheric model (Skamarock et al., 2008; Kochanski et al., 2013a) with SFIRE (Mandel et al., 2011b), a fire spread sub-model that uses a level set approach employing the empirical rate of spread model of Rothermel (Rothermel, 1972). WRF-Fire utilises a subset of the physical processes for the coupling of fire with the atmosphere from the Coupled AtmosphereWildland Fire Environment Model (CAWFE) (Coen, 2005). WRFFire allows the interaction between the fire and local atmospheric effects caused by the fire to be directly modelled. Application of the model was found to compare well in both fire propagation behaviour and extent with recorded wildfire events (Jordanov et al., 2011; Kochanski et al., 2013b). However, no studies have systematically examined the effect of variation in the resolved wind component on the evolution of the fire perimeter. Spatial and temporal variation in environmental conditions has long been understood to influence the behaviour and spread of wildfires (Frandsen and Andrews, 1979; Cruz and Alexander, 2013). The primary advantage of simulation techniques such as those described above is that they can easily incorporate such variations (primarily fuel and weather), at least on a scale commensurate with the computational limitations of the method employed (e.g. Green et al. (1990); Hargrove et al. (2000); Pimont et al. (2006)). However, little research has been presented on the understanding of the effects of such variations, particularly at smaller scales (i.e. order of metres), on rate of spread outside of simulation results. Fujioka (1985) undertook analytical analysis of the effect of non-uniform fuel attributes on the Rothermel fire spread model considering three averaging techniques; arithmetic mean of spread rates, spread rate based on mean fuel conditions, and harmonic mean of spread rates. Their findings suggested that harmonic mean of spread rates was the most appropriate estimator of spread rate in non-uniform fuel. In this study we use a level set method to examine the effect of small-scale spatial variation in combustion conditions (fuel state and combustibility) and temporal variation in wind (magnitude and direction) on the evolution of a two-dimensional fire perimeter. The flexibility of the level set method approach allows these local small-scale spatial and temporal variations to be directly handled by the model. A brief description of the implementation of the level set method is provided, which is then applied to fire spread simulation using our codebase, called ‘Spark’. The effects of variation in environmental conditions (combustion and wind) on fire perimeter shape and rate of spread are then investigated and discussed.

120

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

2. Formulation The level set method models the motion of an interface G subject to a locally varying normal speed s. Along with wildfires, the method is broadly applicable to a range of other physical systems involving interfaces, including the behaviour of droplets in twophase fluid dynamics systems (Sussman et al., 1994; Sethian and Smereka, 2003), charged bubbles (Hilton and van der Net, 2009), front evolution during combustion (Rhee et al., 1995; Pitsch and Duchamp de Lageneste, 2002), crystalline growth (Sethian and Strain, 1992) and vapour deposition processes (Sethian, 2001). As previously mentioned, the most conceptually straightforward approach for simulating the propagation of an interface is to divide the interface into segments and apply motion to the endpoints of each segment. Such front-tracking methods can quickly result in degenerate behaviour, such as the interface crossing over itself, which requires application of multiple topology-dependent rules to check and correct (Bose et al., 2009). Furthermore, the length of the segments will change as the interface evolves, requiring re-distribution of the segments. The level set method avoids such problems by tracking the distance from the interface to each point on the domain, rather than modelling the evolution of the interface itself. The motion of a one-dimensional closed curve, for example, is tracked as the nearest distance to the curve f on a two-dimensional surface. Positive values of f are outside the curve and negative values lie inside. This is shown schematically in Fig. 1, where the interface G divides the domain into f > 0 and f < 0. The interface is then identified, rather than tracked, as the f ¼ 0 curve over the domain. The method allows interfaces to naturally merge and split depending on the underlying dynamics of the interface without requiring additional rules to account for these behaviours. In this implementation the two-dimensional surface is discretised into a two-dimensional grid with spacing Dx and Dy in the x and y directions, respectively. The distance function f is defined at each point i and j on this grid, fij, where i is the index of the grid point in the x-direction and j is the index of the grid point in the y-direction. The time evolution of scalar distance function is given by:





vf þ s

Vf

¼ 0 vt

(1)

where s can be negative or zero. In a fire perimeter the speed function s is necessarily greater than or equal to zero, as the perimeter cannot shrink. The normal vector to the interface can be straightforwardly determined from:

b¼ n

Vf jVfj

In our implementation the derivatives in Eq. (2) are calculated using central differences. Despite the apparent simplicity of Eq. (1), the equation must be solved using suitable hyperbolic upwind techniques. Here, we use the method given by Sethian (2001) to discretise the gradient term in Eq. (1):



qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi


Vf
x D2 þ D2 x y

  fi;j  fi1;j fi;j  fiþ1;j ; ;0 Dx ¼ max Dx Dx   fi;j  fi;j1 fi;j  fi;jþ1 ; ;0 Dy ¼ max Dy Dy

(3)

(4)

This discretisation is slightly different from the flux limited difference schemes listed by Rehm and McDermott (2009) and the scheme used by Mandel et al. (2011a). As pointed out by Sethian (1998), this scheme is slightly more convenient to work with and is used as the basis for the Fast Marching Method. Interestingly, Mandel et al. (2011a) reported instabilities in the level set method which required an additional diffusion term to stabilise. We found no such instabilities in this implementation below the CFL condition (Courant et al., 1967) and no term of this type was included. Mandel et al. (2011a) also reported instability using a first order Euler method for the time derivative and that a second order method such as Heun's method should be used. For the cases investigated here, we found no instability using the Euler scheme and negligible differences between the Euler scheme, the second order Runge-Kutta scheme used by Mandel et al. (2011a) and the second order Adams-Bashforth scheme. We have, however, noticed instabilities in the Adams-Bashforth scheme for long term large scale test cases using empirical speed coefficients. For robustness and efficiency the Euler scheme was used for the cases reported in this study. At the start of a level set simulation, the signed distance function, f, must be calculated from an initial interface. This requires the distance from the interface to be calculated at every point in the grid, which is a computationally expensive task. The two main methods of carrying out this calculation are the Fast Marching Method, based on a discretised version of Dijkstra's algorithm (Dijkstra, 1959), or solution of Eq. (1) to steady state with a speed of s ¼ 1 everywhere. This is equivalent to the steady state solution of:





vf ¼ 1 

Vf

vt

Fig. 1. The scalar distance function f defined over a Cartesian grid. The interface G is identified as the zero, or level, set of f. Values with f < 0 lie within the interface and values with f > 0 are outside. The method allows the interface to be propagated at a local speed normal to the interface, s. Time integration gives an updated interface (dashed lines) at the next time step.

(2)

(5)

where t is a pseudo-timestep (Sussman et al., 1994). Eq. (5) is solved in the same manner as Eq. (1), with each timestep correcting the distance function outwards from the interface. The derivative in Eq. (5) is implemented using Eq. (3). Although the Fast Marching Method is theoretically faster, with O (N logN) complexity, we find that the solution of Eq. (1) results in fewer sharp gradients and gives a smoother solution near the interface. This initialization step must be re-applied at regular intervals during the simulation to correct for any warping. In practice, it is not required to solve Eq. (5) to steady-state for each reinitialisation as warping is generally confined to the region near the interface. Instead, N time-steps of

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

the time integration of Eq. (5) are used at periodic intervals (Sethian and Smereka, 2003). Reinitialisation conditions for the cases in this study are included in each of the relevant sections.

3. Application to fire spread simulations 3.1. Perimeter dependency on functional form Before moving to a more complex empirical propagation model, it is worth considering a highly simplified representation of perimeter propagation. Major factors governing fire propagation include the local wind and the combustion condition, such as fuel load, structure or moisture level, of the underlying fuel bed. We therefore choose the form of the simplified model to only include these two factors. Realistic fire behaviour, of course, depends on a number of additional factors, such as local topographic or atmospheric conditions, suppression effects, etc. However, this simplification allows any characteristic dependencies on the variation of these factors to be identified and compared to more complex perimeter propagation models, highlighting any common underlying dynamics. There are numerous ways to combine the factors affecting fire propagation into a functional form, and many have been presented in the literature. For example, the Rothermel expression (neglecting topography) has the general form s ¼ b(1 þ w) (Rothermel, 1972). The parameter b represents the local combustion conditions and w is a parameter representing the effect of the wind. The parameter b is determined by the local properties of the fuel such as the mass, structure, flammability, moisture level and heat yield. The form of this expression is common to many empirical fire models (Sullivan, 2009b). Several other common forms of empirical models are given in Table 1. Current empirical equations for fire spread are one-dimensional, based on the speed of the front of the fire. There are currently no two-dimensional empirical spread models, although this is a current focus of research. In many simulation models this is done by employing a template shape with geometry derived from an empirical model (Anderson et al., 1982). In our case the onedimensional functional forms given in Table 1 must be expressed in two dimensions. Here we make the assumption that the outward speed of the fire front is given by the one-dimensional expression using the component of the wind in the normal direction (Clark b $w, where n b is the outward et al., 1996). This is given by w ¼ n normal vector of the fire perimeter and w is the wind vector. As the fire front cannot travel backward on itself we additionally impose a lower bound of zero, giving:

 w¼

b $w n 0

b $w > 0 n b $w  0 n

(6)

121

of 300 m  300 m with zero flux Neumann conditions at the boundaries. These boundary conditions were imposed by setting the derivative of any stencils to zero over the boundary. The wind vector was aligned with the vertical axis of the image. The reinitialisation stage was carried out every 100 time-steps by timeintegrating Eq. (5) for 10 pseudo-timesteps. The fire extent in Fig. 2 has been coloured using alternating dark and light bands so perimeter locations can clearly be identified at certain times. In these examples a perimeter is plotted every 30 s. The simulation shown in Fig. 2a used b ¼ 0.2. In this case, the circular perimeter simply expands uniformly in all directions with a speed of 0.2 m s1. The simulations shown in Fig. 2b and c used b ¼ 0.2 and jwj ¼ 0.3. For consistency with Fig. 2b the non-linear



b b 2 term was imposed as w2 ¼
w
ð n $ wÞ for Fig. 2c. The inclusion of a wind dependency in these additive models stretches the perimeter in the direction of the wind. However, it can clearly be seen that a non-linear dependency on the wind parameter has a large influence on the perimeter shape. The perimeter for the non-linear case (c) evolves from an initially circular shape to an oval and then into a teardrop shape. In contrast, the perimeter for the linear case (b) evolves into an obround shape. The simulation shown in Fig. 2d used b ¼ 0.5 and jwj ¼ 1.0. In this case the front is simply propagated in the wind direction with a speed of 0.5 m s1. No backing spread was observed, but a negligible amount of outward flank spread was found in the simulations, possibly due to numerical truncation error in the calculation of the normal gradient. This spread was around one grid cell to either side of the perimeter by the end of the simulation, showing that only a minimal level of error resulted from using first order time and space numerical schemes in our method. 3.2. Propagation of a fire perimeter with variation We now consider the effect of spatial and temporal variation of the wind and combustion condition. Both of these variables are naturally highly inhomogeneous; for example, small scale turbulent fluctuations in the wind speed and direction as well as variation in the local fuel and moisture conditions. These are usually considered locally homogeneous in propagation models, or reduced to a set of homogenised values at a given scale. The effects of these natural variations were imposed on the model by selecting local wind and combustion conditions from a random distribution rather than simply applying mean values over the entire model. To illustrate the effect of variation the functional form s ¼ b  w is used in this section. As shown in Fig. 2d, this functional form has no backing or outward rate of spread, allowing differences in the evolution of the fire front with and without variation to be easily identified. The initial condition was a rectangular fire perimeter at t=5m

Perimeter shapes resulting from each of the example functional forms given in Table 1 are shown for comparison in Fig. 2. The initial condition for the simulation was a circular fire perimeter 36 m in diameter with parameters of Dx ¼ Dy ¼ 1 m, a time-step of 0.5 s and a total simulation time of 5 min. Simulations were run on a domain

t=4m

w

t=3m t=2m t=1m Initial perimeter

Table 1 Example functional forms for one-dimensional fire perimeter propagation, where s is the front speed, b is a combustion condition and w is a wind parameter. Functional form a b c d

S S S S

¼ ¼ ¼ ¼

b bþw b þ w2 bw

Description Combustion condition only Additive model Non-linear model Multiplicative model

a

b

c

d

Fig. 2. Example of perimeter propagation for the example functional forms given in Table 1, starting from an initial circular perimeter. The wind vector w is aligned with the vertical axis in the image. The perimeter extents every 30 s are coloured by alternating dark and light bands for clarity and example times are given on the righthand side of example (d).

122

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

the base of the simulation domain of width 160 m and 10 m depth. The wind vector was aligned with the y-axis in the simulation, which is the vertical direction in the figure, with a magnitude of 1 m s1. The combustion condition parameter b was set to be 1 over the entire domain. The simulation used parameters of Dx ¼ Dy ¼ 1 m, with a time-step of 0.1 s and was run for a total simulation time of 3 min. The simulation domain was 300 m  300 m with zero flux Neumann boundaries conditions. The level set was reinitialised with N ¼ 10 every 100 time-steps. An example simulation with no variation is shown in Fig. 3. The fire extent simply moves forward in time in the direction of the wind, ending 180 m from the initial location, as expected, with no lateral propagation. As in Fig. 2d, there is no backing or outward rate of spread. It should be noted that the rounding of the corners during the evolution of the perimeter is a geometric, rather than numerical, effect resulting from the requirement for the normal vector to be continuous. This causes regions of high curvature, such as the initial sharp corners, to diffuse into regions of low curvature. To investigate the effect of variation, the value of b was picked randomly from a normal distribution with a mean of 1 and a given standard deviation sF. The values of b were initialised at the start of each simulation for the entire computational grid. The wind vector w was given a magnitude of 1 m s1 with a random angle that was chosen from a normal distribution of width sW. The wind direction was assumed to vary on the timescale of the simulation time-step and a random value was picked whenever Eq. (6) required evaluation. It should be emphasised that the use of a normal form for the wind variation is an over-simplification of the turbulent sub-scale fluctuations found in the atmospheric boundary layer. However, this distribution is used to illustrate the effect of variation in the wind vector, rather than to emulate turbulent fluctuations. The random values were generated using the Mersenne Twister implementation (MT19937) provided in the Cþþ standard library. A matrix of cases showing the perimeter locations for distributions of different widths is given in Fig. 4. The rows consist of a fixed wind vector (top), and variations in the angle of the vector (lower two rows). The variations in the angle were picked from normal distributions with standard deviations of 25 and 50 , respectively. The columns consist of b ¼ 1 with no variation (left) and standard deviations of 0.2, 0.3 and 0.4. These values were chosen to illustrate the effect of variation in this simplified model, rather than to represent the variation of any particular naturally occurring fuel type. The values of b from the distribution were bounded at zero to prevent negative speeds. This is imposed as the perimeter must can only expand for fire growth, so the speed must always be zero or positive. Effects due to variation can clearly be seen for both wind direction and combustion condition. Variation in b created an irregular front shape as different parts of the front propagated forward

200 m 150 m 100 m

y

50 m

x 0

Fig. 4. Effect of variation of the wind vector (rows) and combustion condition parameter b (columns), starting from an initial fire line with a speed function given in Table 1d. The simulations are shown for three minutes duration, where the bands are shaded using the same colour scale as Fig. 3. Variation in b creates an irregular front, and wind variation causes lateral spread in the perimeter. Both types of variation slow the arrival time.

at different rates. More importantly, however, variation was found to slow the propagation of the fire as a whole, as can be seen by comparing the positions of the final dark bands of the top left and top right images in Fig. 4. Variation of the wind direction also slowed the rate of propagation of the fire, as can be seen by comparing the top left and lower left images in Fig. 4, as well as increasing flank propagation rates, giving a more rounded fire perimeter. The reduction in front speed was quite pronounced when wind and fuel amount were both varied. The case with the maximum variations considered, shown in the lower right image, showed a maximum extent of 142 m, a minimum extent of 120 m and a median of 131 m from the initial perimeter after 180 s. The median extent for this case was therefore approximately 73% of the extent of the case with no variation, shown in the upper left image. To investigate the effect of directionality on b, a number of tests were carried out in which b was varied in one direction but remained

t=3m t = 2 m 45 s t = 2 m 30 s t = 2 m 15 s t=2m t = 1 m 45 s t = 1 m 30 s t = 1 m 15 s t=1m t = 45 s t = 30 s t = 15 s Initial perimeter

Fig. 3. Example of perimeter propagation with no variation, starting from an initial fire line aligned with the x axis and the multiplicative speed function given in Table 1d. The perimeter extents every 15 s are coloured by alternating dark and light bands for clarity.

Fig. 5. Effect of combustion condition variation perpendicular to the initial front (upper) and parallel to the initial front (lower). Three trials are shown for both the perpendicular and parallel cases, each with different random variations in b.

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

homogeneous in the other direction. The directions used were along the x-axis, parallel to the initial front, and along the y-axis, perpendicular to the initial front. For the case along the x-axis, one value of b was picked from a distribution of width sF ¼ 0.2 for each x position. Correspondingly for the case along the y-axis, one value of b was picked for each y position using the same distribution width of sF ¼ 0.2. The perimeter positions for three different random normal distributions as a function of time are shown in Fig. 5. It can be seen that variation of b in x created a highly uneven perimeter (Fig. 5, upper), with some parts spreading much farther ahead than neighbouring parts. In contrast, variation of b in y had no effect on the perimeter shape and only a minor effect on the maximum distance travelled by the front. In order to verify the findings of the model, resolution and numerical tests were carried out. Reduction of the time-step gave the same trends observed in Fig. 4. Similarly, varying the spatial discretisation gave the same overall results. A resolution study is shown in Fig. 6 for three different spatial and temporal discretisations for a distribution width of sW ¼ 25 for the wind direction and sF ¼ 0.4 for b. The same random seed and distribution was used for each simulation. The smoothness of the perimeter is dependant on the underlying grid resolution, and the perimeter shape is different due to the different number of random values sampled. However, the median perimeter extent is approximately the same over each of the cases, 135 m for Dx ¼ Dy ¼ 2, 135 m for Dx ¼ Dy ¼ 1 and 134.75 m for Dx ¼ Dy ¼ 0.5. The overall effects of variation on the perimeter extent therefore appear to be independent of the model resolution.

123

2.00 km t = 44 mins

1.75 km 1.50 km 1.25 km 1.00 km 0.75 km 0.50 km 0.25 km

t = 40 mins t = 44 mins t = 40 mins t = 36 mins t = 32 mins t = 28 mins t = 24 mins

t = 36 mins t = 32 mins t = 28 mins t = 24 mins t = 20 mins

t = 20 mins

t = 16 mins

t = 16 mins

t = 12 mins

t = 12 mins t = 8 mins t = 4 mins Initial fire line

t = 8 mins t = 4 mins Initial fire line

0.00 km Fig. 7. Extent of fire using the CSIRO grassland fire spread model, coloured by alternating dark and light bands. Variation in curing and wind are used for the results shown on the left and the results with no variation (unmodified values) are shown on the right. The width of each band is four minutes and the parameters used for the simulations are given in Table 2.

Cm

8 < expð0:108MÞ ¼ 0:684  0:0342M : 0:547  0:0228M

M  12% M > 12%; W  10 km h1 M > 12%; W > 10 km h1

(8)

3.3. Empirical perimeter propagation model We now assess the dependence on wind and fuel values for an operational fire spread model. The speed function in this case was taken directly from a published empirical rate of spread model, where the rate of propagation was calculated from spatially varying values defined on an underlying data surface. We used the CSIRO Grassland Fire Spread Model (Cheney et al., 1998; Cheney and Sullivan, 2008), which is based on the analysis of 121 experimental grassland fires and other field observations (Cheney et al., 1998). The rate of forward spread in m s1, applicable to fire spread in cut or grazed grass, is given by:

( sðWÞ ¼

Cm Cc 1:1=3600 W  5 km h1  . 0:844 3600 W >5 km h1 Cm Cc 1:1þ0:715ðW 5Þ (7)

where Cm is a moisture coefficient, Cc a curing coefficient and W is the wind speed (km h1) at a height of 10 m in the open. The moisture coefficient is given by:

Δt = 0.2 Δx = Δy = 2

Δt = 0.1 Δx = Δy =1

Δt = 0.05 Δx = Δy = 0.5

Fig. 6. Resolution test for a simulations with distribution widths of sW ¼ 25 for the wind direction and sF ¼ 0.4 for b. The irregularity of the perimeter appears higher at low resolutions, but the arrival times are very similar between each of the cases.

where M is the grass fuel moisture content percentage. This can be measured directly or modelled from the temperature and relative humidity of the air as (Sullivan, 2010):

M ¼ 9:58  0:205T þ 0:138H

(9)

where T is air temperature ( C)and H the relative humidity (%). The curing coefficient is given by:

Cc ¼

1:12 1 þ 59:2expð0:124ðc  50ÞÞ

(10)

where c is the percentage of grass curing, 50%  c  100%. The effect of wind direction was incorporated in Eq. (7) as b $wÞ, where n b is the outside normal to the interface, calculated sð n using Eq. (2) and w is the wind vector with magnitude equal to W. In the following simulations the wind speed, wind direction and curing percentage were varied to determine the effect on the propagation of the fire. The simulation time-step was set from the CFL condition (Courant et al., 1967) Dt ¼ ½ min{Dx, Dy}/max{s}, where max{s} was the maximum absolute cell speed in the simulation. A narrow band technique (Sethian and Smereka, 2003) was also used to increase the speed of the level set method by restricting calculations to only the local neighbourhood of the interface. The width of the narrow band was chosen to be 15 cells from the perimeter, and the level set was reinitialised and the band re-calculated when the front came within 5 cells of the edge of the band. The simulation set-up used a domain of size 2 km  0.7 km with Dx ¼ Dy ¼ 2.5 m with zero flux Neumann boundaries conditions. The initial fire perimeter was a rectangular strip 110 m long by 7.5 m deep. This set-up used values for the temperature, wind speed, relative humidity and curing percentage from reported field conditions from an actual event (Dalton, 2012), given in Table 2. These simulation conditions are typical of moderate intensity

124

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

wildfire conditions. As a point of comparison with the twodimensional model, rates of spread calculated from the onedimensional empirical expression Eq. (7) for the mean values used in Table 2 are given in Table 3 (top). The effect of variation was imposed in the same manner as the basic propagation model discussed in the previous section. Here, we investigated the variation in the curing percentage, wind speed and wind direction. The curing percentage was picked from a normal distribution at the start of the simulation and remained fixed throughout. The random values were bounded to the model limits of 50%  c  100%. The wind speed and direction were assumed to be temporally uncorrelated and were picked at each time-step from a random normal distribution of given width. As an example, the extent and positions of the fire perimeter at four minute intervals for the parameters and distribution widths given in Table 2 are shown on the left of Fig. 7. The fire perimeter formed an approximately oval shape during the early stages of propagation before stretching into a sharp tip at the front. The early oval stages qualitatively match experimentally observed grassland fires (Green et al., 1983; Cheney et al., 1993; Cheney and Gould, 1995) lit from point ignitions. As a baseline comparison, the results of the same simulation without variation are shown on the right. It can be seen that including variation results in the broadening of the flanks and a shortening of the overall shape, in the same manner found in the simplified model. The dependency of fire perimeter on the width of the curing distribution is shown in Fig. 8. Changing the curing distribution can be seen to have a large influence on both the rate of propagation and the shape of the perimeter for the same mean value. Rates of spread for the different curing variations, calculated using the maximum extent divided by the total time, are given in Table 3 (bottom). The case with zero variation resulted in the prediction of a highly elongated fire shape. The rate of spread in this case was approximately equal to the value in Table 3 (top) for the mean curing, as expected. For s ¼ 15% the perimeter shape was a rounded oval, which became progressively more pointed in the direction of the wind as the width of the curing distribution was reduced. The maximum predicted extent of the fire at the end of the simulation was found to be around 50% shorter for s ¼ 15% than s ¼ 2%. The dependency of propagation on the variation in wind speed and direction is shown in Figs. 9 and 10, respectively. Variation in wind direction appeared to have no effect on the forward rate of spread, as can be seen by comparing the maximum extents for a directional distribution of zero than s ¼ 15 in Fig. 9. However, the width of the perimeter grew as the width of the directional distribution was increased, in a similar trend to that found in the simplified two-dimensional model considered in the previous section. In contrast, Fig. 10 shows that random variation in wind speed appeared to have no effect in the model on either the extent or width of the fire perimeter. It is interesting to note that experimental fire perimeters have been reported of two types, a ‘broad parabolic shaped headfire where the flanks developed to a width wider than the ignition line’

Table 3 One-dimensional rates of spread for varying curing percentages for the CSIRO Grassland Fire Spread Model given by Eqns. 7e10 (top). Rates of spread computed using the level set model with varying distribution widths (bottom). Curing (%)

Curing distribution width (%)

Empirical rate of spread (km h)1

80 78 75 70 65

e e e e e

2.47 2.12 1.64 1.01 0.59

Curing (%)

Curing distribution width (%)

Computed rate of spread (km h)1

80 80 80 80 80

0 2 5 10 15

2.44 2.45 2.34 1.88 1.37

and ‘a narrow pointed headfire where the width of the fire along the flanks often did not exceed the width of the ignition line' (Cheney et al., 1993). In this study both of these behaviours were qualitatively reproduced by variation of only the fuel distribution (Fig. 8) or variation in wind direction (Fig. 9). This suggests that local variation of these factors may play a significant role in the overall form of the fire perimeter.

Fig. 8. Extent of fire for increasing widths of the curing distribution. The wind speed and direction distribution widths are fixed at s ¼ 3 km h1 and s ¼ 10 , respectively.

Table 2 Base case parameters for the CSIRO Grassland Fire Spread Model given by Eqns. 7e10. The results from a simulation using these parameters is shown in Fig. 7. Parameter

Value

Temperature, T Relative humidity, H Wind mean speed, W Wind standard deviation Wind direction standard deviation Curing mean, c Curing standard deviation

32 25% 27 km h1 3 km h1 5 80% 10%

Fig. 9. Extent of fire for increasing widths of the wind direction distribution. The wind speed and curing distribution widths are fixed at s ¼ 3 km h1 and s ¼ 10%, respectively.

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127 n X 1 s ¼n F i¼0 i

!1



Fig. 10. Extent of fire for increasing widths of the wind speed distribution. The wind direction and curing distribution widths are fixed at s ¼ 10 and s ¼ 10%, respectively.

4. Discussion 4.1. Variation in combustion condition One might expect that using a distribution in combustion condition (such as curing state) would have little effect as the variations would cancel over time, giving propagation behaviour tending towards an average trend. For example, consider a highly simplified one-dimensional spread model where the combustion condition is linearly proportional to the fuel amount. The fire front passes over terrain in which the amount of fuel F has a variance s2 around some mean value. A schematic diagram of this model is shown in Fig. 11, where the fire passes from left to right in a discretised series of steps i ¼ 1...n. All units are non-dimensionalised for simplicity, and the spacing between each step is set to unity. Assuming that the speed of the fire is given by s ¼ F, where we have neglected any constant of proportionality, the speed in each cell i is si ¼ Fi. It may then be assumed that the front will move at a mean P speed s ¼ ni¼0 Fi =n, and that s/F as n / ∞. However, this is not the case and can result in a large overestimation of the front speed. As an example, consider the limiting case where only one cell j > 0 has zero fuel Fj ¼ 0. Assuming only direct cell-to-cell propagation, the front must stop indefinitely within this cell and the arrival time after this cell must therefore tend towards infinity, t / ∞. The mean speed is the total distance travelled, n, divided by P the total time, s ¼ n=t, giving s/0. However, s ¼ ni¼0 Fi =n is clearly non-zero as the only cell with zero fuel is Fj. This apparently inconsistent result arises purely from the assumption that the arithmetic mean can be used. In fact, the arrival time should be the P sum of the time spent within each cell, t ¼ ni¼0 ti , giving Pn s ¼ n= i¼0 ti . The time in each cell is ti ¼ 1/si ¼ 1/Fi, revealing that the harmonic, rather than arithmetic mean should be used:

Fig. 11. A simplified one-dimensional model of fire propagation consisting of n cells, where the ith cell has a random fuel amount F picked from a distribution. A normal distribution is shown on the far right hand side as an example. The speed of the fire, s, within each is assumed to be proportional to F.

125

(11)

where we have used an asterisk in s* to denote the harmonic mean. A similar result was presented by Fujioka (1985) and is, essentially, the zero-lag case derived by Catchpole et al. (1989) using a Markov process, which they regard as a reasonable approximation in realistic situations. The harmonic mean is very sensitive to low values due to the dependency on reciprocals, and recovers the correct behaviour in the degenerate case due to the presence of the single term 1/Fj / ∞. This sensitivity was seen in the simplified two-dimensional model, but is particularly clear in the model using the empirical rate of spread. The curing coefficient Cc used in the model and the reciprocal, 1/Cc is plotted against the curing percentage c in Fig. 12. In a random distribution of c, low values will tend to dominate the harmonic mean through the strong dependence on 1/Cc. Such behaviour clearly gives rise to the significantly reduced front extent in Fig. 8 at high widths of curing distribution. Furthermore, it can be shown that the harmonic mean is always less than or equal to the arithmetic mean (Burk, 1985). Unfortunately, tight bounds on the degree of overestimation is difficult to gauge analytically. Some general bounds on any set of positive non-zero elements {si} are given by Meyer (1984):

s  minfsi g

pffiffiffi 2 ðB  1Þ2  s  s  minfsi g B  1 nðB þ 1Þ

(12)

where the maximum normalised element B in the set {si} is given by B ¼ max{si}/min{si}. The upper bound on the left-hand side of Eq. (12) is difficult to quantify from this expression as it approaches the arithmetic mean as n / ∞. However, any B > 1 leads to the surprising result that using the arithmetic mean in such models will always overestimate the front speed and hence arrival time of the front. It is also interesting that the lower bound does not depend on n and can be re-arranged to give:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  s  s  maxfsi g þ minfsi g  2 maxfsi gminfsi g

(13)

showing that the lowest possible value the harmonic mean can take for any distribution depends only on the maximum and minimum values in the set {si}. If max{si} [ min{si}, the difference between the harmonic and arithmetic mean will be O (max{si}). Use of the arithmetic mean for calculation of speed-dependent terms in a rate

Fig. 12. The curing coefficient Cc given in Eq. (10) and the reciprocal, 1/Cc against the curing percentage, c.

126

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127

of spread expression, such as the simplified one-dimensional model above and the level set models, can therefore result in large error estimations for the rate of fire propagation. 4.2. Variation in wind speed and direction The variation in the wind considered in this study was temporal, in contrast to the spatial variation in combustion condition discussed above. Taking the example for the previous section, we assume that the speed of the fire is dependent on the wind speed W by s ¼ W, where we have non-dimensionalised and neglected any constant of proportionality. The speed in each cell i is si ¼ Wi, and as Wi is independent of position the law of large numbers gives P s/ ni¼0 Wi =n as n / ∞. The front speed will therefore approach the mean wind speed at long times. The approach to the mean value will, of course, depend on the underlying distribution of Wi and bounds on this can be determined using methods such as the law of the iterated logarithm (Bingham, 1986). The lack of dependence on wind speed deviation shown in Fig. 10 is, most likely, explained by this behaviour. The wind vector is re-generated at each grid point at each time step, unlike the curing which remains fixed throughout the simulation. However, it should be noted that the use of a normal distribution for the mean wind speed is a considerable over-simplification of the dynamic properties of the wind around the fire perimeter, which is likely to possess a high degree of local correlation both temporally and spatially (Taylor, 1935). Further work is planned to investigate the distribution and correlation of wind flow in the proximity of the front and to integrate this into the fire propagation model. Variation in local wind direction was found to increase the perimeter width in both the simplified (Fig. 4) and empirical model (Fig. 9). The mechanism for this widening is the random outward fluctuations in the wind direction at the flanks. These outward fluctuations result in a positive outward speed through the dot b $w. In the case considered here, the wind vector is product n aligned with the y-axis and the distribution is assumed normal with standard deviation s. For a front with a normal aligned with b ,w > 0 is equal to 1/2. The amount of the x-axis, the probability of n widening is proportional to the outward component of wind which is, in turn, dependant on the width of the distribution. The average lateral spread of the front, in this restricted case, therefore scales with the standard deviation of the wind direction. A more detailed analysis of flank spread due to wind variation will be addressed in a future study. However, it is evident that the spread is directly related to the variation in wind direction and that this variation will have a significant impact on the flank spread of a fire perimeter. 5. Conclusion We have investigated the effect of variation of factors such as wind and combustion condition on the rate of propagation and perimeter shape predicted by a fire spread simulation method based on the level set technique. The method qualitatively matched observed perimeter shapes for the CSIRO grassland model. The predicted rate of spread was also equal to the empirical rate of spread when no variation was used in the model. A simplified twodimensional model showed the primary characteristics of introducing variation. Firstly, combustion condition variation perpendicular to the net motion of the perimeter creates an irregular fire perimeter and slows the rate of propagation. Secondly, variation of wind speed and direction also slows the rate of propagation and allows the flanks of the fire to expand. The effect of combustion condition can be determined by using a one-dimensional spread rate over a heterogeneous fuel bed. A simplified one-dimensional model showed that it is necessary to use a harmonic, rather than

arithmetic, mean to represent any variables which have a spatial dependency on the speed of the fire front. This finding is of interest from an operational point of view, as using an arithmetic mean over large areas of fuel could greatly overestimate propagation rates. A two-dimensional representation is, however, necessary to truly replicate perimeter behaviour as the irregular shape of the fire perimeter was found to be strongly dependent on the orientation of the fuel distribution. Use of a more complex rate of propagation model showed similar characteristics to the simplified model. Pointed and parabolic fire shapes reported in the experimental literature could be reproduced using the model. Interestingly, the fire shape was found to be only dependent on combustion condition in this study. This suggests that the perimeter behaviour is highly dependent on the homogeneity of the underlying fuel, at least for the assumption of normal distribution used in this study. The directionality of the wind was found to determine the width of the fire, which is crucial for prediction of flank spread. Variations in width of the wind speed distribution, however, may not affect the rate of spread or shape of the fire. Our findings indicate that small-scale variation in combustion condition and wind direction is very important to the rate of spread and perimeter shape, which are the two main outputs required by fire propagation models. This suggests variation must be taken into account in operational fire models to provide more accurate estimates of fire behaviour. Future work will aim to further quantify the effect of variation, as well as to employ ensemble runs to build spatial risk predictions for fire behaviour over heterogeneous landscapes. Acknowledgments We would like to gratefully acknowledge Ryan Fraser and Mahesh Prakash of the CSIRO Digital Productivity Flagship for providing the strategic funding to enable us to carry out this work. References Achtemeier, G.L., 2013. Field validation of a free-agent cellular automata model of fire spread with fire-atmosphere coupling. Int. J. Wildland Fire 22, 148e156. Albini, F.A., 1982. Response of free-burning fires to non-steady wind. Combust. Sci. Technol. 29, 225e241. Anderson, H.E., Rothermel, R.C., 1965. Influence of moisture and wind upon the characteristics of free-burning fires. Symp. Int. Combust. 10, 1009e1019. Anderson, D.H., Catchpole, E.A., de Mestre, N.J., Parkes, T., 1982. Modelling the spread of grass fires. J. Aust. Math. Soc. 23, 451e466. Andrews, P.L., 1986. BEHAVE: Fire Behavior Prediction and Fuel Modeling System e BURN Subsystem, Part 1. USDA Forest Service Gen. Tech. Rep. INT-194. Andrews, P.L., 2014. Current status and future needs of the behaveplus fire modeling system. Int. J. Wildland Fire 23, 21e33. Bingham, N.H., 1986. Variants on the law of the iterated logarithm. Bull. Lond. Math. Soc. 18, 433e467. Bose, C., Bryce, R., Dueck, G., 2009. Untangling the prometheus nightmare. In: 18th World IMACS Congress and MODSIM09 International Congress on Modelling and Simulation, 13e17 July 2009, Cairns, Australia, pp. 240e246. Burk, F., 1985. By all means. Am. Math. Mon. 92, 50. Catchpole, E.A., Hatton, T.J., Catchpole, W.R., 1989. Fire spread through nonhomogeneous fuel modelled as a Markov process. Ecol. Model. 48, 101e112. Catchpole, E.A., Catchpole, W.R., Rothermel, R.C., 1993. Fire behaviour experiments in mixed fuel complexes. Int. J. Wildland Fire 3, 45e57. Cheney, N.P., Gould, J.S., 1995. Fire growth in grassland fuels. Int. J. Wildland Fire 4, 237e247. Cheney, N.P., Sullivan, A.L., 2008. Grassfires: Fuel, Weather and Fire Behaviour. CSIRO Publishing, Collingwood, Australia. Cheney, N.P., Gould, J.S., Catchpole, W.R., 1993. The influence of fuel, weather and fire shape variables on fire-spread in grasslands. Int. J. Wildland Fire 3, 31e44. Cheney, N.P., Gould, J.S., Catchpole, W.R., 1998. Prediction of fire spread in grasslands. Int. J. Wildland Fire 8, 1e13. Cheney, N.P., Gould, J.S., McCaw, W.L., Anderson, W.R., 2012. Predicting fire behaviour in dry eucalypt forest in southern Australia. For. Ecol. Manag. 280, 120e131. Clark, T., Jenkins, M., Coen, J.L., Packham, D., 1996. A coupled atmosphere-fire model: convective feedback on fire-line dynamics. J. Appl. Meteorol. 6, 875e901.

J.E. Hilton et al. / Environmental Modelling & Software 67 (2015) 118e127 Coen, J.L., 2005. Simulation of the Big Elk Fire using coupled atmosphere/fire modeling. Int. J. Wildland Fire 1, 49e59. Coen, J.L., Cameron, M., Michalakes, J., Patton, E.G., Riggan, P.J., Yedinak, K.M., 2013. WRF-fire: coupled weather wildland fire modeling with the weather research and forecasting model. J. Appl. Meteor. Climatol. 52, 16e38. Coleman, J.R., Sullivan, A.L., 1995. SiroFire: the CSIRO bushfire spread simulator. In: Proc. Inst. Foresters Aust. 16th Biennial Conf., Canberra, pp. 309e319. Coleman, J.R., Sullivan, A.L., 1996. A real-time computer application for the prediction of fire spread across the Australian landscape. Simulation 67, 230e240. Courant, R., Friedrichs, K., Lewy, H., 1967. On the partial difference equations of mathematical physics. IBM J. Res. Dev. 11, 215e234. Cruz, M.G., Alexander, M.E., 2013. Uncertainty associated with model predictions of surface and crown fire rates of spread. Environ. Model. Softw. 47, 16e18. Dalton, J., 2012. Operational Review: Westmeadows Grassfire 24 January 2012. Department of Sustainability and Environment, Fire Services Commissioner Victoria. Dijkstra, E.W., 1959. A note on two problems in connexion with graphs. Numer. Math. 1, 269e271. Ellis, P.F.M., 2011. Fuelbed ignition potential and bark morphology explain the notoriety of the eucalypt messmate ‘stringybark’ for intense spotting. Int. J. Wildland Fire 20, 897e907. Encinas, L.H., White, S.H., del Ray, A.M., Sanchez, G.R., 2007. Modelling forest fire spread using hexagonal cellular automata. Appl. Math. Model. 31, 1213e1227. Filippi, J.B., Morandini, F., Balbi, J.H., 2010. Discrete event front-tracking simulation of a physical fire-spread model. Simulation 86, 629e644. Finney, M.A., 1994. Modelling the spread and behaviour of prescribed natural fires. In: Proc. 12th Conf. Fire and Forest Meteorology, pp. 138e143. Finney, M.A., 1998. FARSITE: Fire Area Simulatoremodel Development and Evaluation. USDA Forest Service, RMRS-RP-4. Finney, M.A., 2004. FARSITE: Fire Area Simulator e Model Development and Evaluation. Res. Pap. RMRSRP-4 (revised). USDA Forest Service, Rock Mountain Research Station, Ogden, UT. Frandsen, W.H., Andrews, P.L., 1979. Fire Behaviour in Non-uniform Fuels. USDA Forest Service. INT-232. Fujioka, F.M., 1985. Estimating wildland fire rate of spread in a spatially nonuniform environment. For. Sci. 31, 21e29. Goodrick, S.L., Achtemeier, G.L., Larkin, N.K., Liu, Y., Strand, T.M., 2012. Modelling smoke transport from wildland fires: a review. Int. J. Wildland Fire 22, 83e94. Green, D.G., Gill, A.M., Noble, 1983. I.R., Fire shapes and the adequacy of fire-spread models. Ecol. Model. 20, 33e45. Green, D.G., Tridgell, A., Gill, A.M., 1990. Interactive simulation of bushfires in heterogeneous fuels. Math. Comput. Model. 13, 57e66. Hargrove, W.W., Gardner, R.H., Turner, M.G., Romme, W.H., Despain, D.G., 2000. Simulating fire patterns in heterogeneous landscapes. Ecol. Model. 135, 243e263. Hilton, J.E., van der Net, A., 2009. Dynamics of charged hemispherical soap bubbles. EPL 86, 24003. Johnston, P., Kelso, J., Milne, G.J., 2008. Efficient simulation of wildfire spread on an irregular grid. Int. J. Wildland Fire 17, 614e627. Jordanov, G., Beezley, J.D., Dobrinkova, N., Kochanski, A.K., Mandel, J., Sousedik, B., June 6e10, 2011. Simulation of the 2009 Harmanli fire (Bulgaria). In: 8th International Conference on Large-scale Scientific Computations (Sozopol, Bulgaria). Knight, I., Coleman, J., 1993. A fire perimeter expansion algorithm based on Huygens wavelet propagation. Int. J. Wildland Fire 3, 73e84. Kochanski, A.K., Jenkins, M.A., Mandel, J., Beezley, J.D., Clements, C.B., Krueger, S., 2013a. Evaluation of WRF-SFIRE performance with field observations from the FireFlux experiment. Geosci. Model Dev. 6, 1109e1126. Kochanski, A.K., Jenkins, M.A., Mandel, J., Beezley, J.D., Krueger, S., 2013b. Real time simulation of 2007 Santa Ana fires. For. Ecol. Manag. 294, 136e149. Lopes, A.M.G., Cruz, M.G., Viegas, D.X., 2002. FireStation e an integrated software system for the numerical simulation of fire spread on complex topography. Environ. Model. Softw. 17, 269e285. Luke, R.H., McArthur, A.G., 1978. Bushfires in Australia. Australian Government Publishing Service, Canberra. Mandel, J., Beezley, J.D., Kochanski, A.K., 2011a. Coupled atmosphere-wildland fire modelling with WRF 3.3 and SFIRE 2011. Geosci. Model Dev. 4, 591e610. Mandel, J., Beezley, J.D., Kochanski, A.K., Kondratenko, V.Y., Zhang, L., Anderson, E., Daniels II, J., Silva, C.T., Johnson, C.R., October 2011b. A wildland fire modelling and visualization environment. In: Ninth Symposium on Fire and Forest Meteorology. Palm Springs, CA. Meyer, B., 1984. Some inequalities for elementary mean values. Math. Comput. 42, 193e194. Peet, G.B., 1967. The shape of mild fires in Jarrah forest. Aust. For. 31, 121e127.

127

Perry, G.L.W., Sparrow, A.D., Owens, I.F., 1999. A GIS-supported model for the simulation of the spatial structure of wildland fire, Cass Basin, New Zealand. J. Appl. Ecol. 36, 502e518. Pimont, F., Dupuy, J.L., Scarella, G., Caraglio, Y., Morvan, D., 2006. Effects of small scale heterogeneity of vegetation on radiative transfer in forest fire. For. Ecol. Manag. 234, S88. Pitsch, H., Duchamp de Lageneste, L., 2002. Large-eddy simulation of premixed turbulent combustion using a level-set approach. Proc. Combust. Inst. 29, 2001e2008. Rehm, R.G., McDermott, R.J., 2009. Fire-front Propagation Using the Level Set Method. NIST Technical Note 1611. Rhee, C.W., Talbot, L., Sethian, J.A., 1995. Dynamical study of a premixed V flame. J. Fluid Mech. 300, 87e115. Richards, G.D., 1990. An elliptical growth model of forest fire fronts and its numerical solution. Int. J. Numer. Methods Eng. 30, 1163e1179. Richards, G.D., Bryce, R.W., 1995. A computer algorithm for simulating the spread of wildland fire perimeters for heterogeneous fuel and meteorological conditions. Int. J. Wildland Fire 5, 73e80. Rothermel, R.C., 1972. A Mathematical Model for Predicting Fire Spread in Wildland Fuels. Res. Pap. INT-115. U.S. Department of Agriculture, Forest Service, Intermountain Forest and Range Experiment Station, Ogden, UT. Scott, J., Burgan, R., 2005. Standard Fire Behavior Fuel Models : a Comprehensive Set for Use with Rothermel's Surface Fire Spread Model. USDA Forest Service, Rocky Mountain Research Station. General Technical Report RMRS-GTR-153, 72 pp. Sethian, J.A., 1998. Fast marching methods. SIAM Rev. 41, 199e235. Sethian, J.A., 2001. Evolution, implementation, and application of level set and fast marching methods for advancing fronts. J. Comput. Phys. 169, 503e555. Sethian, J.A., Smereka, P., 2003. Level set methods for fluid interfaces. Annu. Rev. Fluid Mech. 35, 341e372. Sethian, J.A., Strain, J.D., 1992. Crystal growth and dendritic solidification. J. Comput. Phys. 98, 231e253. Sharples, J.J., 2008. Review of formal methodologies for wind-slope correction of wild fire rate of spread. Int. J. Wildland Fire 17, 179e183. Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Duda, M.G., Huang, X.Y., Wang, W., Powers, J.G., 2008. RA Description of the Advanced Research WRF Version 3. NCAR Technical Note TN-475. Sullivan, A.L., 2009a. Wildland surface fire spread modelling, 19902007. 1: physical and quasi-physical models. Int. J. Wildland Fire 18, 349e368. Sullivan, A.L., 2009b. Wildland surface fire spread modelling, 19902007. 2: empirical and quasi-empirical models. Int. J. Wildland Fire 18, 369e386. Sullivan, A.L., 2009c. Wildland surface fire spread modelling, 19902007. 3: simulation and mathematical analogue models. Int. J. Wildland Fire 18, 387e403. Sullivan, A.L., 2010. Grassland fire management in future climate. Adv. Agron. 106, 173e208. Sullivan, A.L., Ball, R., 2012. Thermal decomposition and combustion chemistry of cellulosic biomass. Atmos. Environ. 47, 133e141. Sullivan, A.L., McCaw, W.L., Cruz, M.G., Matthews, S., Ellis, P.F., 2012. Fuel, fire weather and fire behaviour in Australian ecosystems. Flammable Australia: Fire Regimes, Biodiversity and Ecosystems in a Changing World. CSIRO Publishing, Collingwood, Vic, pp. 51e77. Sullivan, A.L., Sharples, J.J., Matthews, S., Plucinski, M.P., 2014. A downslope fire spread correction factor based on landscape-scale fire behaviour. Environ. Model. Softw. 62, 153e163. Sun, R., Krueger, S.K., Jenkins, M.A., Zulauf, M.A., Charney, J.J., 2009. The importance of fireatmosphere coupling and boundary-layer turbulence to wildfire spread. Int. J. Wildland Fire 18, 50e60. Sussman, M., Smereka, P., Osher, S., 1994. A level set approach for computing solutions to incompressible two-phase flow. J. Comput. Phys. 114, 146e159. Taylor, G.I., 1935. Statistical theory of turbulence e Part 1. Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci. 151, 421e444. Tolhurst, K., Shields, B., Chong, D., 2008. Phoenix: development and application of a bushfire risk management tool. Aust. J. Emerg. Manag. 23, 47e54. Trunfio, G.A., D'Ambrosio, D., Rongo, R., Spataro, W., Di Gregorio, S., 2011. A new algorithm for simulating wildfire spread through cellular automata. ACM Trans. Model. Comput. Simul. 22, 1e26. Tymstra, C., Bryce, R.W., Wotton, B.M., Taylor, S.W., Armitage, O.B., 2010. Development and Structure of Prometheus: the Canadian Wildland Fire Growth Simulation Model. Information Report NOR-X-417. Canadian Forest Service, Northern Forestry Centre. Vasconcelos, M.J., Guertin, D.P., 1992. FIREMAP e simulation of fire growth with a geographic information system. Int. J. Wildland Fire 2, 87e96. Viegas, D.X., 2004. Slope and wind effects on fire propagation. Int. J. Wildland Fire 13, 143e156. Viegas, D.X., Soares, J., Almeida, M., 2013. Combustibility of a mixture of live and dead fuel components. Int. J. Wildland Fire 22, 992e1002.