Short-term fire front spread prediction using inverse modelling and airborne infrared images

Additional keywords: data assimilation, fire behaviour, Rothermel model.

Wildfire phenomena constitute a widespread problem that causes human, environmental and economic losses all over the world every year. One of the key aspects to reduce their impact is to improve and optimise wildfire-fighting capacities by providing the emergency responders with sound information on the upcoming wildfire dynamics. However, the accuracy of the current available operational wildfire models such as FARSITE (Finney 1998) or PHOENIX Rapidfire (Tolhurst et al. 2008) is limited owing to the scarcity of precise data available to initialise them (Finney et al. 2013) and the empirically developed submodels that they contain (Sullivan 2009), which make them unsuitable to be exported to all sorts of different fire events. An alternative approach to wildfire modelling is computational fluid dynamics (CFD)-based models such as FIRETEC (Linn 1997), FOREFIRE (Filippi et al. 2009; Filippi et al. 2014a) or Wildland–Urban Interface Fire Dynamics Simulator (WFDS) (Mell et al. 2007). However, these are restricted to research use, and small- and particular-scale applications owing to the high computational costs and initialising data required (Viegas 2011). A simplified version of a CFD model that couples weather modelling with Rothermel’s fire-spread model (Rothermel 1972) named WRF-SFIRE (Mandel et al. 2014, 2011) was recently reported to have been applied operationally, as a 6-h forecast could be launched with a 30-min run in a multicore processor.

Cruz and Alexander (2013) reviewed 49 fire spread models datasets and found that only 3% of the simulations acceptably predicted the observed rate of spread (RoS) and revealed the need to adjust those models if they were to be used operationally. To deal with this problem, data assimilation has shown great potential to correct a model’s predictions. Mandel et al. (2009) pioneered the use of data assimilation for wildfire modelling, predicting flame temperature and location using an ensemble Kalman filter framework and an atmosphere-coupled wildfire model. Their work showed promising results while raising some concerns about spurious fire corrections and the computing time required. Following their idea, Rochoux et al. (2014a, 2014b, 2015) explored a data-driven wildfire simulator based on parameter and state estimation that assimilates fire front positions and corrects the wildfire forecast by means of a level set model based on Rothermel’s. They explored a parameter and state estimation strategy with stochastically based estimation of the error covariance matrices. The simulator was run using parallel computing and showed great capacity both in synthetically generated data and a small-scale controlled experiment. Lautenberger (2013) also used Rothermel’s model and a data-assimilating approach to calibrate fire propagating parameters to meet observed fire shape. All these authors used an Eulerian level set method to track fire front propagation. Although this approach offers advantages when dealing with crossovers and large perimeters, it also has some drawbacks such as the level set solution dependency on the mesh (Mallet et al. 2009) and the requirement of a flux limiter to filter fire front shapes (Bova et al. 2016).

The alternative approach to propagate the front is the Lagrangian framework, which tracks front markers individually. This approach has been successfully applied on most broadly used fire simulators (FARSITE (Finney 1998), Phoenix RAPIDFIRE (Tolhurst et al. 2008), or PROMETHEUS (Tymstra et al. 2010) as an example). The final front shape equivalence of both approaches was recently demonstrated by Bova et al. (2016) over complex topographical domains.

Some authors also worked to couple Lagrangian-based fire models with data assimilation. Denham et al. (2012) successfully applied genetic algorithms (GA) to find the wind configurations that best resemble observations to launch an improved forecast. A combination of weather and fuel calibration using fire perimeters has also been implemented using FARSITE (Finney 1998) and high-performance computing, showing great improvements and potential for long-term predictions (Wendt et al. 2011; Artés et al. 2014). Combining these ideas, Rios et al. (2014) developed a data-driven algorithm based on Rothermel’s RoS model (Rothermel 1972) and Huygens’ elliptical propagation (Richards 1990; Glasa and Halada 2008) and proved that it provides a short-term highly accurate forecast of wind-driven wildfires when tested with synthetically generated data.

In the present work, we further develop this algorithm to cope with real wildfire data to produce a reliable short-term forecast. We test it by assimilating airborne infrared images captured from two large-scale experiments conducted in South Australia (Cruz et al. 2013) and assess the forecast reliability by means of Jaccard (Jaccard 1901; Filippi et al. 2014c) and Sørensen (Sørensen 1948; Perry et al. 1999) similarity indices, finding great agreement in the available short-term forecast.

The paper is organised as follows: the Methods section presents the complete algorithm rationale and description of real-scale experiments used for validation. The algorithm performance is explored in depth in the Results and discussion section where the forecasting capabilities and system parameters are investigated and discussed. Finally, we provide some concluding remarks and an outline of proposed future work.

Our system combines two main sequential tasks that are envisaged to take place cyclically. The first one involves solving an inverse fire modelling problem by which the parametric space of certain observed fire front behaviour is found. This approach is based on the assumption that there are several variables that remain constant during a certain period of time, and thus we call them invariants. These invariants can be either physical quantities or unknown model parameters that determine the fire front dynamics. With the aid of a proper assimilation framework, these invariants can be resolved using the available data (e.g. fuel characteristics, topography, weather conditions) to reproduce certain fire behaviour observations (i.e. fire perimeters sensed by airborne IR). The second task involves launching a forecast with the forward fire behaviour model, relying on the corrected invariants, which are then entered as input back into the model to deliver a more accurate forecast. The algorithm begins with the input of the observed fire front isochrones, the complementary available data and initialisation of the invariants. Then, an optimisation loop identifies the proper invariants that best suit the fire dynamics observed, and these are finally used to launch the forecast. The algorithm flow is depicted in Fig. 1 and shows the key aspects that will be explained in detail in the following sections.

Fig. 1.  Flow diagram showing the algorithm rationale. Orange rectangular boxes are inputs and green elliptical boxes are outputs. (For colour reproduction of figures, see online version available at http://www.publish.csiro.au/nid/17.htm.) Adapted from Rios et al. (2014).

The first key aspect of the overall system is the forward model. The forward model is the mathematical model that allows one to explain the phenomenon and exhibits a capacity to reproduce the assimilated wildfire dynamics. To be suitable for an operational inverse-modelling scheme, the forward model has to be computationally efficient to ensure the forecast is delivered before the event takes place (i.e. positive lead time). While keeping options open to explore further models, in the present work, we chose a model based on Rothermel (Rothermel 1972) to determine RoS and a Huygens’ expansion-based approach to propagate the front in 2D (Richards 1990, 1993). The Rothermel model (see Eqn 1) is based on a simple energy balance that expresses the fire heat source by means of the reaction intensity (), the propagating flux ratio (), the fuel heat sink by means of the heat of preignition (), the effective heating factor () and the bulk density (). The RoS is corrected by two experimental factors to account for the slope () and wind effects ().

Using Rothermel’s (1972) correlations, the RoS can be ultimately expressed in terms of the fuel depth (), oven-dry fuel loading (), the surface area-to-volume ratio (SAV), the packing ratio (), the moisture of extinction and fuel moisture content (), and the wind and slope factors ():

In the current version of our algorithm, the slope factor is not considered, thus the model applies only to flat terrains.

In order to reduce the number of invariants, we performed a sensitivity study of Rothermel’s model to explore the dependency and influence of the variables on the RoS. One thousand randomly generated fuels models sets (), ranging among physically possible values taken from Scott and Burgan (2005), were used to explore the variables’ correlation with RoS. For each given set, the variables were swept from their minimum to their maximum while calculating the Pearson’s product-moment correlation coefficient to search for linear dependencies. The fuel depth () was found to be the variable that exhibited the greatest linear behaviour as demonstrated by the Pearson’s coefficient histogram for the randomly generated fuel sets shown in Fig. 2.

Fig. 2.  Sensitivity study of Rothermel’s model dependency with fuel depth . The Pearson’s product-moment correlation distribution (a) and the corresponding rate of spread (RoS) vs fuel depth for 1000 model sets (b).

Thus, to determine useful invariants, the Rothermel RoS was factorised as:

where is an invariant that accounts for the overall effect of fuel moisture, moisture of extinction, bulk density, surface area-to-volume ratio and wind factor.

The Huygens’ expansion approach considers every point in the fire front as a virtual ignition point that will expand following an elliptical geometry called a firelet. The curve that envelops all firelets is then the expanded perimeter. The lateral () and the front velocity () of each firelet (Fig. 3) can be derived from geometrical relationships given RoS and the length-to-breadth ratio of the firelet as:

Fig. 3.  Geometrical relationship between a_j, b_j and c_j velocities for the elliptical firelet corresponding to the node p_j of a given perimeter. is the wind–slope direction and the subscript j identifies a node of a given fire front.

where is a front discretisation parameter and is the spreading time.

Note that Eqns 5 and 6 satisfy Eqn 7, as this is the distance between the virtual ignition point on the perimeter and its perpendicular expansion as illustrated in Fig. 3.

To determine the length-to-breadth ratio , we used Anderson’s experimental correlation (Anderson 1983):

where U is the mid-flame wind speed (in m s^–1).

The curve enveloping the firelets is calculated using the partial differential equations derived by Richards (1990, 1993) (Eqns 9 and 10), which are resolved using a predictor–corrector method.

where is the simulation time, and is the slope–wind direction, which corresponds to the direction in which the wind is blowing (i.e. where it is blowing towards) if the terrain is flat. Variables x and y are the coordinates of the nodes that constitute a front at a given time.

In the present implementation, the mid-flame wind speed in Eqn 8 and direction used in Eqns 9 and 10 are also considered invariants to be identified and are named and respectively.

With Eqns 4–10, the final forward model that provides the output of the set of front perimeters for a given simulation time can be defined as a function of the initial observed perimeter , the fuel depth and the invariants vector as:

Owing to the Lagrangian formulation of the forward model in hand, when the front exhibits convex regions the Huygens expansion approach may create some loops or overlaps along the perimeter (Fig. 4). These topological problems are physically meaningless and a loop-clipping filter is required to remove these loops and reclose the remaining front. There are different methods in the literature to tackle this problem. The turning number approach (Barber et al. 2008) can be used to identify the front sections that are internal to the curve and filter them out. However, if small integration steps are used in Eqns 9 and 10, the loops do not interact with the perimeter (see left loop in Fig. 4) but just fold in on themselves (right loop depicted in Fig. 4). In this case, the loops can be filtered out by finding the adjacency matrix corresponding to the nodes p_j (see Fig. 4) and removing all closed paths.

Fig. 4.  Schematic representation of loop formation scenarios and regridding strategy. The loop on the right (red dotted line) corresponds to a double loop that only intersects once with the perimeter (red asterisks). The left loop interacts multiple times with the perimeter. is the new node introduced by applying the regridding algorithm as is larger than the threshold T. (For colour reproduction of figures, see online version available at http://www.publish.csiro.au/nid/17.htm.)

The Lagrangian formulation also yields another problem. As the fire front grows and spreads, the distance between nodes increases and more nodes are needed to avoid formation of sharp unreliable regions (see Fig. 5). This regridding process is conducted by comparing the adjacent angles α_j and α_j+1 (see Fig. 4) of two consecutive segments of the front. If the largest angle exceeds a given threshold T, an extra node is added at the midpoint of the corresponding segment. With this filter, more nodes are automatically generated in the more abrupt and sharp areas, smoothing the forthcoming simulated isochrones.

Fig. 5.  Fire front expansion without loop clipping or regridding filters (a), and with both filters applied (b). The regridding threshold is set to . The marker line (red crosses) is the initial fire perimeter. (For colour reproduction of figures, see online version available at http://www.publish.csiro.au/nid/17.htm.)

To perform a proper assimilation process, we need to define a cost function () to minimise differences between the observed front and the simulated front. We define this cost function as the sum over all simulated isochrones of the averaged Euclidean distance between simulated and observed fronts. These distances are computed by selecting the modelled (m) points () in between simulated nodes () and each point in the corresponding observed (ob) nodes () determined by perpendicular intersection (as depicted in Fig. 6).

Fig. 6.  Cost function geometrical definition between simulated front (dashed line) and observed front (solid line). Red arrows are the Euclidean distances between simulated and observed perimeters. Shaded (light-blue) zone represents the area between observed and simulated fire perimeter, to be used as a convergence controller.

Then, the cost function () can be written in terms of the invariants vector as:

where is the number of nodes in each isochrone, the subscript i is the identifier of the simulated (and assimilated) front at any given time, and is the weight matrix, which can be set to give more weight to the perimeters assimilated later in time, and give less importance to those assimilated a longer time ago. Note that to efficiently apply the optimisation condition () of Eqn 13, we take the square of the Euclidian distance (removing the squared root) as it does not affect the optimisation.

For the examples presented in the present study, the identity matrix is used, which represents a homogeneous weight function distribution.

Once the cost function is defined, the inverse problem to find can be formulated as the following optimisation problem:

for a given initial front with observed coordinates , fuel depth and assimilating period t. The invariants cover a range of values that have a physical meaning, with and as the maximum invariant values for a given scenario. This turns the optimisation problem into a constrained problem and allows particular optimisation techniques. In the present implementation, if any invariant exceeds these ranges during an assimilation loop, its value is reset to the previous estimated value.

The vector-based cost function definition is critical because it will drive the invariant’s convergence to the final optimised value. The actual vector-based method used to quantify the dissimilarity between two curves may entail some converging discrepancy when compared with the actual area between curves (see Fig. 6), as some close curves may give large vectored distances – i.e. large values of – due to the perpendicular definition even though the area difference is kept small (see Fig. 6). This effect can appear when handling large time steps and irregular front shapes (normally caused by highly heterogeneous terrain or fuel distribution). Thus, although we need a point-to-point function to apply our optimising strategy, we evaluate the area between curves at each optimisation loop to control both convergences.

Two main approaches can be used to solve optimisation problems (Onwubolu and Babu 2013): the gradient-free and the gradient-based algorithms. The former involves searching the optimisation domain to find a global minimum, but is computationally expensive if multiple invariants are used, because the cost function needs to be evaluated multiple times for each invariant change and so does the forward model. On the other hand, the gradient-based strategies do not guarantee an absolute minimum; however, if the forward model behaves smoothly and the algorithm is initialised in the vicinity of the solution, they will provide the correct answer with great computational efficiency (Nocedal and Wright 1999).

Which is the more suitable strategy depends on the particularity of the system and its requirements. The problem at hand is a constrained non-linear optimisation that may have to be solved multiple times during an operational event, and thus has to be computationally efficient. Moreover, we want the algorithm to be ready to handle a larger number of invariants. Therefore, we choose a gradient-based approach.

To speed up the optimisation process, we take advantage of a Tangent Linear Model (TLM) to iteratively linearise the forward model on the vicinity of a set of parameters (). The TLM is applied by means of a first-order Taylor expansion around an initial parameter estimation . The vector with the initialisation of the invariant values is chosen by an educated guess based on available data.

The gradient of the linearised cost function is then written as:

which can be rewritten as:

which is a linear system that can be rapidly solved for () by using a QR factorisation algorithm (Nocedal and Wright 1999).

To solve the gradient , we use automatic differentiation (Bischof et al. 2006) because it does not require running the forward model multiple times as more common methods such us central differences (Griewank 2000) do. Adimat software (Bischof et al. 2002) was used to differentiate the code.

The algorithm was tested with data from two large-scale experimental fires performed in South Australia in March 2008. Despite these being within the framework of a scientific experimental burning program, the fires exhibited real wildfire behaviour patterns as they were performed in large plots under extremely severe weather conditions. These experiments were conducted in Ngarkat Conservation Park (35°45′S, 140°51′E), which consists of a characteristic dune and swale system comprising large flat areas (130 m above sea level) of mallee-heath shrublands. Fire behaviour in mallee-heath fuel types is characterised as being discontinuous and highly variable owing to the heterogeneous characteristics of the various fuel layers that comprise mallee-heath fuel complexes (Cruz et al. 2013).

Data used in the present paper come from two experimental burns. These burns were performed in two different sites named hereafter ‘Shrub site’ and ‘Woodland site’ (plots A and AS2E respectively according to Cruz et al. 2013). The burn in the Shrub site was performed in a 9-ha, 9 year-old heath plot on 4 March 2008. This fuel complex was characterised by scattered small-leafed shrubs, organised in clumps, and a discontinuous litter layer partially buried by sand. The main fire-carrying fuel layer was the discontinuous shrub canopy. Ambient conditions were a temperature of 32°C and relative humidity of 25%. Wind speed (10-m open) averaged 15 km h^–1, with gusts up to 30 km h^–1. Wind direction was south–south-westerly. The characteristic dead fuel moisture of the fuel complex was 7.7%. The fire was ignited with a 150-m long line and 2 min after ignition, flame heights were ~2–2.5 m, with flashes up to 4 m. The fire spread vigorously throughout the plot with sustained flames heights of 4–5 m. Eight minutes after ignition, the main flame front hit the northern border of the plot, concluding a 350-m head fire run (Planas et al. 2011).

The burn in the Woodland site was performed on 5 March 2008, in a 25-ha plot of a 22-year-old woodland dominated by mallee eucalypts. The fuel complex in this plot was characterised by a mallee overstorey 2 to 3 m tall and a shrubby understorey that constituted the main fuel layer supporting fire propagation. Air temperature was 36°C and relative humidity 13%. The 10-m open wind speed averaged 17 km h^–1, gusting up to 36 km h^–1. Sampled dead fuel moisture of dead suspended fuels was 7%. The fire was ignited using a 250-m line and initially spread in heath vegetation with flame heights averaging 4 m. As the fire burnt into mallee vegetation, crowning ensued, with flame heights averaging 4.5 m, and peaking at 10 m. A short lull in wind speed 5 min after ignition slowed fire propagation and revealed multiple spot fires along the fire perimeter. As wind speed increased again, the fire made an intense crown fire run in the mallee vegetation, with flame heights between 8 and 10 m and spot fires developing 40–60 m ahead of the fire.

The Woodland site was used to study and characterise aerial suppression efficiency (Plucinski and Pastor 2013). Thus, for the evaluation of the algorithm developed in the present work, only those perimeters before the first aerial drop were used.

Both plots were filmed from a hovering helicopter equipped with an IR camera. This camera operated within the 7.5–13-µm range and stored sequences of IR images (240 × 320-pixel) at an approximate rate of 5 frames s^–1. The helicopter was positioned so that the majority of the plot was in view for the duration of each fire. The images were georeferenced using beacons deployed along the plot as referencing markers and georectified assuming a flat terrain (see Fig. 7). Fire front positions were extracted every 10 s to create isochrones and RoS maps for both plots by applying a methodology for IR analysis described elsewhere (Pérez et al. 2011).

Fig. 7.  Shrub site composite of (a) the visible channel, (b) the raw infrared (IR) image, (c) the orthonormalisation of the image using reference points and (d) the isochrone extractions and representation in a geographic information system (GIS) interface. All the images represent t = 210 s after ignition.

The fire in the Shrub site lasted 370 s, giving 38 perimeters at 10-s frequency, although the initial 130 s are considered to be within the artificial linear ignition and were discarded, leaving 25 perimeters available for the current algorithm validation (left plot in Fig. 8). The fire in the Woodland site lasted 450 s but the first drop discharge took place 240 s after ignition. Thus, taking the perimeters from 60 s on to avoid the effect of the artificial ignition, 19 perimeters at 10-s frequency were left for the purpose of the current study (right plot Fig. 8).

Fig. 8.  Shrub site (a) and Woodland site (b) burnt during the Ngarkat experiments campaign. Isochrones plotted every 10 s are represented on GoogleEarth after georeferencing and orthonormalising. Map data: Google, CNES/Astrium.

The two presented plots are chosen to illustrate the algorithm performance when tested with real data. The wind speed and direction invariants are initialised using approximate values registered during the experiment (simulating the manner the system could be used in a real operation). The invariant initial value (Table 1) is calculated with Eqns 1 and 2. Required parameters are taken from the standard fire behaviour fuel model SH5 (Scott and Burgan 2005) corresponding to ‘Very high load, dry climate shrubs’ and complemented with values of moisture measured during the Ngarkat experiments (Table 1). The moisture content was measured using the oven-dry method. The upper limits for the invariants and are an upper bound needed to set the optimisation problem. Their values were set considering conservative maximum reasonable values for the scenario. In the present examples, those values were never reached.

Table 1.  All values needed to initialise the algorithm and their source

The fuel depth is considered homogeneous. The system is tested using this approximation as in a real-case scenario where an accurate fuel map might not be available. In the following sections, different outputs of the system are separately evaluated.

The first factor to validate the system at hand is the assimilation stage. Only if the optimisation algorithm manages to find the proper invariants that model the observed fronts will it be possible to launch a reliable forecast. The assimilating process is iterative and stops when the percentage difference of the invariants between two consecutives loops (named relative error, ) reaches a threshold of 0.05% as:

where the index k represents each invariant, the superscript i each loop iteration, and is the invariant initialisation value.

As mentioned before, the area between the curves is tracked to check the correctness of the cost function. The absolute areal error is calculated by summing all the enclosing areas between the simulated front, , and the corresponding observed front, :

where is the XOR logical operator, is the absolute areal error and the subscript i corresponds to a given assimilating time.

The match between simulated and observed fronts is evaluated by three similarity indices:

1. SDI′ is the inverted shape deviation index (SDI) (Cui and Perera 2010):

2. Sørensen (Sørensen 1948):

3. Jaccard (Jaccard 1901):

The three similarity indices have a score range of [0–1], 1 being the best score, indicating identical shapes. Their main difference is the denominator identification or, in other words, the normalising factor. Because it is not clear which of the mentioned indexes performs best to quantify simulation correctness (Filippi et al. 2014b), we analyse the three of them.

Fig. 9 shows six assimilated isochrones for the Woodland site. The isochrone captured at 40 s after the ignition is used as the initial perimeter . The subsequent assimilated isochrones span 10 s. The optimisation loop converges after three iterations and the similarity indices are kept over 0.9, showing great resemblance.

Fig. 9.  Assimilating process for the Woodland site: (a) The assimilation is resolved after three optimisation iterations. (b) Detailed areal difference after optimisation convergence for every simulated isochrone. Solid (black) lines are the observed front, dotted (red) lines the simulated front and dashed (green) lines the initial front considered. (c) Similarity indexes. Note that the y-axis is expanded between 0.90 and 1.00. SDI is shape deviation index. (For colour reproduction of figures, see online version available at http://www.publish.csiro.au/nid/17.htm.)

The same procedure is applied to the Shrub site; in this case, six isochrones spanning 10 s are assimilated. The first isochrone corresponds to 140 s after fire ignition. The optimisation loop converges after eight iterations (Fig. 10). The similarity indices are over 0.9 for the three initial fronts and decrease to 0.8 owing to imperfect simulation of the right flank. This flank showed an extremely low RoS that could not be fully replicated by the forward model as only one common lateral RoS () is resolved.

Fig. 10.  Assimilating process for the Shrub site: (a) The assimilation is resolved after eight optimisation iterations. (b) Detailed error after optimisation convergence for every simulated isochrone. The shaded region (blue) represents the areal error for a given simulated and observed isochrone. (c) Similarity indexes. Note that y axis is expanded between 0.50 and 1.00. SDI is shape deviation index. For lines definition, see Fig. 9 caption.

In both cases, the assimilation stage shows great ability to model the front data regardless of the number of assimilated isochrones.

The effectiveness of the algorithm relies on the individual convergence of each invariant. Only if they are all individually resolved will their value be meaningful and useful to run the forecast. For each assimilation, the convergence of the cost function and the invariants’ individual convergence is tracked and displayed in Fig. 11 and Fig. 12.

Fig. 11.  Woodland site convergence: (a) The vector cost function value (dotted blue line) and the areal cost value (solid orange line) during the optimisation loop. (b) The relative convergence for the three invariants. (For colour reproduction of figures, see online version available at http://www.publish.csiro.au/nid/17.htm.)

Fig. 12.  Shrub site convergence: (a) The vector-based (blue dotted line) and the areal absolute error (orange solid line). (b) The invariants' relative error. The invariant fluctuation is due to the nature of angular direction restricted to [0–2π].

The invariant relative error (, Eqn 17) shows how far from the final value the invariant was initialised. For the Woodland site, the invariants converge to values of 0.542 s^–1, 2.651 m s^–1 and 0.569 rad (see Table 1), which represents a change from the initial values of 500%, –80% and 25% (Fig. 11). Applying Eqn 3 and assuming an average fuel depth of 1.4 m, the invariant can be expressed as a RoS with the value of 45.5 m min^–1. This value represents a RoS average for the whole perimeter during the assimilation period and is in agreement with the values reported in Planas et al. (2011) where RoS was found to range between 20 and 60 m min^–1 along the whole perimeter. The wind speed and direction also had reasonable values although more detailed data are not available for comparison.

The vector cost function calculated as Eqn 13 converges to 1.2 m (see Fig. 11), which indicates the average normal distance between simulated and observed fronts. At the same time, the sum of the absolute areal error for all perimeters converges to a value of 1363 m².

Similar results are obtained for the Shrub site (see Fig. 12). In this case, invariants converge to the values 0.447 s^–1, 5 m s^–1 and 5.95 rad, which represent an average RoS of 10.4 m min^–1. The values observed in Planas et al. (2011) ranged from 5 to 40 m min^–1. In this case, the average distance between simulated and observed fronts is ~4 m whereas the areal distance is 2135 m². As mentioned before, both convergence indicators do not follow exactly the same behaviour and the vector-based one finds a minimum before stabilising. This highlights the need for an areal-based cost function instead of a vector-based one.

Once the invariants are resolved, we can launch a forecast using the last assimilated isochrone as the initial perimeter () and propagate it with the forward model while assessing its efficiency. For each case explored, the observed available fronts are split into two groups. The first is used as data assimilation input whereas the second group is reserved for evaluating the correctness of the forecast. In addition to the forecasting error, to evaluate the appropriateness of the system, the lead time must be computed. The lead time is defined as the period of validity of the model once the computational time is subtracted. In order for the model to be operative, the lead time needs to be positive. For the current exploration, we considered the forecast to be valid as long as the SDI′ index was kept over 0.85.

Following the already defined scenarios, when five isochrones are assimilated every 10 s in the Woodland site, the lead time reaches 130 s as the similarity index gradually decreases, reaching the 0.85 threshold (see Fig. 13). The forecast run for the Shrub site, when six fronts are assimilated, yields a larger lead time. Although the similarity indices decrease to 0.8 at 50 s, they stabilise at 0.9 for the rest of the available fronts, as shown in Fig. 14.

Fig. 13.  Forecasting error for Woodland site when six fronts are assimilated: (a) Similarity index for every 10-s spanned front. The lead time is 130 s (until 13th forecasted front). SDI is shape deviation index. For lines definition, see Fig. 9 caption.

Fig. 14.  Forecast fronts and error for shrub Shrub site (six fronts assimilated with 10-s spacing each). Lead time reaches 190 s as all available fronts are predicted within the accepted error, predicting all the available fronts. SDI is shape deviation index. For lines definition, see Fig. 9 caption.

From this point on, the fire behaviour in the right flank changes, accelerating the rate of spread, and the forecasted front is therefore underpredicted. The underlying cause might be due to the fact that Rothermel’s model was initially derived for surface fire spread. As fire propagates through crowns, the model is no longer valid, although a proper assimilation structure extends its validity as fire in the Woodland site was erratically torching and crowning.

It is important to mention that the forecasting step can be improved by adding extra layers of information, particularly those that can change in time and space (i.e. fuels depth map, wind speed and direction, etc.). If those layers are available, the lead time can be improved as the present algorithm cannot resolve time-dependent variables.

If the algorithm is to be used operationally, the consequences of the amount and frequency of assimilated data need to be investigated. The assimilation window (AW) is defined as the number of fronts assimilated before launching a forecast. To explore this, we use the Shrub site to run our assimilation and forecasting system 12 times, changing the number of assimilated fronts. The assimilation window varies from 8 assimilated fronts (i.e. 70 s) up to 19 assimilated fronts (i.e. 180 s) in the last run. The last assimilated front is kept constant for all runs and corresponds to 190 s after ignition. Then, the forecast is launched to generate the 18 fronts corresponding to the 180 s left of observed data. The areal relative error is computed to assess the validity of each forecast.

The results are shown in Fig. 15. As more fronts are assimilated, the SDI′ index tends to stabilise at ~0.87. It is interesting to note that, for the current experiment, there is an optimum AW at ~11–12 assimilated fronts where the lead time is maximised and reaches 180 s. As more than 12 fronts are assimilated, the fronts corresponding to times close to ignition are taken into account. In particular, when 17 fronts are assimilated, the first assimilated isochrone corresponds to 10 s after the ignition. At this stage, the fire is still driven by initial acceleration effects, and thus the forecast slightly overestimates the RoS and the forecasting areal error tends to grow, dramatically decreasing the lead time.

Fig. 15.  Forecasting areal relative error for shrub site and 10 different assimilation windows (from 70 to 160 s, assimilating one front every 10 s). SDI is shape deviation index.

We present a data assimilation framework that enhances a Rothermel-based model to make it capable of forecasting short-term wildfire dynamics even when using it beyond its original applicability scenarios. The assimilation framework previously applied to synthetically generated data is improved to deal with real scenarios and tested in two large-scale shrubland fire experiments conducted in South Australia, yielding similarity index scores over 0.8 and obtaining positive lead times of 130 and 160 s, depending on the scenario. The system performance when changing the AW is explored to conclude that the more fronts are assimilated, the better the forecast validity, up to a point where the initial fire dynamics perturb the forecast. Different forms of the weighting matrix present in the cost function are also explored in depth although improvement is limited.

Even though in the present work the available fire data last a few minutes and extend up to several hundred metres, the algorithm structure is envisaged to be applied at larger time and space scales. At the operational level, the time required to perform image processing tasks such as georeferencing and fire front edge detection must be considered as they will directly affect the lead time. Current work is being done to automate and improve those tasks (Valero et al. 2015).

In the future, the algorithm will be prepared to handle non-flat terrain scenarios by interacting with digital terrain models input beforehand. Time-varying inputs (such as wind speed or moisture content) will be used to enlarge the lead time and take full advantage of the available data. Space-dependent invariants will be explored to handle fuel and fire dynamics heterogeneity. Finally, the algorithm will be implemented recursively as we believe that by running it systematically, the fire spread prediction could be updated and extended to provide a sound operational tool.