A high-fidelity ensemble simulation framework for interrogating wildland-fire behaviour and benchmarking machine learning models

Numerical methods

The simulations presented in this work are performed with Swirl-Fire (Wang et al. 2023), which is an open-source large-eddy simulation framework with a fully coupled combustion model to describe the interaction with the atmosphere. In this framework, the gas phase is described by the Favre-filtered conservation equations for mass, momentum, oxygen mass fraction, and potential temperature. The combustion of the solid fuel is modelled using a one-step mixing-limited oxidation reaction, considering the energetic impact of moisture evaporation (Linn et al. 2002), radiative losses, and convective heat transfer between fire and ambience. With this approach, physico-chemical processes involving pyrolysis, vaporisation, and combustion are simplified following an assumption that the energy release is confined to regions where solid fuel is present (Linn 1997). The governing equations are solved with a low-Mach number approach (Wang et al. 2022) to alleviate limitations of the time-step size due to acoustic waves. For further details on the physical models, governing equations, and numerical algorithm, the reader is referred to Wang et al. (2023).

The solver was validated against the Fire-Flux II measurements (Clements et al. 2019), which considered the fire-spread behaviour on a flat surface and prescribed wind profile. To extend this validation to conditions of interest to this study, we performed additional comparisons against the experiments by Morandini and Silvani (2010), considering different slope angles, wind conditions (speed and direction), fuel densities, and moisture content. Results from these simulations and comparisons with experimental data are presented in the Appendix 1 of this study.

To consider complex terrains, we utilise the immersed-boundary (IB) method (Zhang and Zheng 2007) with direct forcing. This method provides a second-order accurate representation of the terrain without the necessity for unstructured meshes or intricate coordinate transformations, thereby retaining the overall efficiency of a structured mesh representation.

The Swirl-Fire simulation framework was implemented in TensorFlow (Abadi et al. 2015) using a just-in-time compilation approach to generate a highly optimised executable that fully exploits the hardware architecture for parallel execution on TPUs. The program is compiled using the Accelerated Linear Algebra (XLA) compiler, which generates a data flow graph with optimised computational and parallel efficiency. Linear weak and strong scalabilities have been demonstrated on canonical flows, showing its effectiveness for large-scale problems (Wang et al. 2022).

All simulations are performed on the fifth generation TPU v5e computing architecture (Jouppi et al. 2020). Each chip in this computing architecture consists of one tensor compute core, with four matrix-multiplication units to simultaneously perform four matrix–matrix multiplications of size 128² on each core, providing a peak throughput of 197 teraflops. The TPU v5e chips are connected through a 2D Torus, with an inter-chip interconnect bandwidth of 1600 Gbps, forming a so-called TPU pod with 128 chips. The TPU v5e architecture inherently facilitates 8-bit arithmetic operations, while high-precision arithmetic operations can be achieved through software emulation. To achieve a balance between computational efficiency and numerical precision, simulations in this study are conducted using single precision.

Supported by this novel hardware architecture, we are able to perform simulations at significantly larger scales and with reduced time and power utilisation compared to those conducted with conventional CPU architectures. In this investigation, simulations on meshes of 1.35 billion mesh points are executed on 128 TPU chips and take approximately 13 h of wall-clock time per simulation, resulting in a total of 1,664 TPU hours to create the FireBench dataset.

Ensemble simulation framework

To enable high-fidelity ensemble simulations, we are leveraging the open-source optimisation platform Vizier (Golovin et al. 2017). Given only a parameter space and a cost function over this space, Vizier efficiently searches for optimal parameters. In addition to providing various optimisation algorithms, Vizier is a system for managing compute-heavy jobs on multiple machines. In the present study, we are exploring the two-dimensional parameter space, consisting of wind speed and slope angle, as uniformly as possible but also taking advantage of the job management capabilities that Vizier provides. To achieve this, we use a cost function that provides no information about the search space to Vizier, thus forcing it to explore the parameter space uniformly.

Vizier batch ensemble simulations are launched as an automated service through a Remote Procedure Call interface (Tay and Ananda 1990). After a simulation request is submitted, Vizier will schedule the simulation runs automatically following configurations specified in the protocol buffer. Given the availability of the computational resources, Vizier distributes runs across multiple machines to execute simulations in parallel. During these runs, the parameter sets being scheduled can be inspected dynamically through a dashboard interface.

In this study, a total of 117 individual ensemble simulations are scheduled with Vizier, and are performed in parallel with 10 simulations running simultaneously. In situ postprocessing scripts are executed by Vizier to provide statistical analysis of the simulations, including 2D inspections of flow-field results and extraction of fire-front location, burn area, and flame length without manual intervention.

Instantaneous flow-field results

To provide a qualitative comparison of the fire-spread behaviour for the different conditions, in Fig. 3 we present instantaneous flow-field results for nine ensembles, corresponding to conditions U₁₀ = {2, 6, 10} m/s and α = {0, 15, 30}° at a representative time instance where the fire front reaches the location of x = 500 m. The results in this figure show volume renderings of the potential temperature, and projections of axial (x–z plane) and vertical (y–z plane) velocity components. Comparing results along columns provides understanding about the effect of wind speed, showing that increasing wind speed results in a transition from plume-dominated to wind-dominated behaviour. In contrast, for increasing slope angle the fire tends to attach closer to the surface, resulting in a reduction of flame height. This observation agrees with the literature, suggesting a transition between plume-dominated and wind-dominated fire modes as a function of wind speed and slope angle (Ju et al. 2019; Liu 2023).

Fig. 3.

3D visualisations of fire propagation for nine ensembles when the mean location of the fire front is x = 500 m: (a) U10S00, (b) U10S15, (c) U10S30, (d) U06S00, (e) U06S15, (f) U06S30, (g) U02S00, (h) U02S15, and (i) U02S30, showing instantaneous volume renderings of the potential temperature, the projection of axial velocity (x–z plane) at y = L_y/2, and vertical velocity (y–z plane) at x = 1100 m.

These results also provide insight about the coupled interaction of the slope and wind speed, clearly showing that these two parameters cannot be treated independently. To examine the coupling further, we extract the instantaneous fire front ξ_f which is defined as the most upstream point with the potential temperature θ = 400 K on the ground level. These results are presented in Fig. 4, showing isochrones of the fire front at discrete time intervals of 20 s. These results provide a direct illustration of the fire-front corrugation due to the coupling between buoyancy and convective heat-exchange processes (Finney et al. 2015) and fire intermittency by the atmospheric turbulent boundary layer (Wang et al. 2023). Specifically, it can be seen that increasing wind speed and slope not only result in a faster fire spread rate, but also in an increase in the amplitude and frequency of the fire-front corrugation.

Fig. 4.

Isochrones of fire front locations ξ_f for nine ensembles: (a) U10S00, (b) U10S15, (c) U10S30, (d) U06S00, (e) U06S15, (f) U06S30, (g) U02S00, (h) U02S15, and (i) U02S30. The fire front is defined as the most upstream point with θ = 400 K. Numbers next to the isochrones are the time of arrival in units of seconds.

Statistical analysis

To examine the results from the ensemble simulation more quantitatively, we proceed with a statistical analysis. For this, we average the flow-field quantities along the homogeneous lateral direction, and compute a mean fire front 〈ξf〉y=Ly−1∫ξfdy and corresponding fluctuation of the fire front ξf′=ξf−〈ξf〉y. The laterally averaged instantaneous ROS 〈R〉_y is then computed as:

Results from this analysis are presented in Fig. 5, corresponding to the nine cases shown in Figs 3 and 4. The mean fire location and ROS are shown by solid lines, with the shaded area and error bars representing the instantaneous fluctuations of the fire front location and ROS.

Fig. 5.

Comparison of the statistical results for (top) mean fire-front location 〈ξ_f〉_y and (bottom) rate-of-spread ROS 〈R〉_y for three different wind speeds (a, d) U₁₀ = 2 m/s, (b, e) U₁₀ = 6 m/s, and (c, f) U₁₀ = 10 m/s and three slope angles (α = {0, 15, 30}°) shown by different colours in each panel.

The results in Fig. 5 show that the fire dynamics can be characterised by three distinct phases. The different phases are more pronounced for fires with larger slope angles. The first phase is associated with the fire transition following the fire-line ignition. In this phase, the fire advances only slowly without significant fire corrugation (see Fig. 4). Following this transition phase, the fire advances with a nearly constant acceleration that is marked by sporadic intermittency. The spatial extent of this quasi-equilibrium phase depends on the wind speed and slope angle and reduces with increasing α and U₁₀. Based on results shown in Fig. 5, we identify this phase as the region where the fire front location is between ξ₀ = 400 m to ξ₁ = 800 m. The dashed lines in Fig. 5d–f show linear fits of the ROS with respect to time,

where t₀ and t₁ are the times when the mean fire front reaches ξ₀ and ξ₁, respectively, 〈R〉_y,0 is the ROS at t₀, a is the linear acceleration of the fire, and the superscript * denotes the quasi-equilibrium period. By integrating Eqn 3, we obtained the mean ROS in the quasi-equilibrium phase as:

The mean ROS and linear acceleration in the quasi-equilibrium period are presented in Fig. 6. From these results, it can be seen that the fire spreads at almost a constant rate at low wind speeds and small slope angles. This is consistent with observations and several empirical models (Sullivan 2009), where wind speed and slope angle are considered as independent parameters, not considering feedback on the fire-spread rate. Comparisons of predicted zero-slope mean ROS with measurements and the phenomenological models by Rothermel and McArthur’s fire danger index (Mk V) (Rothermel 1972; Noble et al. 1980) are presented in Fig. 2b. Our predictions show overall good agreement with observations and the empirical fire-spread models capture the trends, with Rothermel’s model providing a slightly better agreement. All results fall within the range of the experimental observations, while showing a slightly faster trend, which we attribute to the lower bulk-fuel density.

Fig. 6.

Quantitative analysis of ensemble simulations as a function of slope and wind speed, showing (a) mean ROS, 〈R〉_y, and (b) linear acceleration a, computed from Eqn 3 during the quasi-equilibrium phase for ξ ∈ [400,800] m.

Compared to the zero-slope results, Fig. 6b shows that with increasing slope angle the unsteady behaviour becomes more pronounced. This is attributed to the coupling of the fire to the atmosphere by enhanced air entrainment, buoyancy, and heat transfer to the fuel bed upstream of the fire front. The general trend agrees with observations by Cheney and Gould (1997). In particular, in an experiment with a wind speed of 4.5 m/s, the observed average acceleration is close to 0.04 m/s² (Cheney and Gould 1997), which is close to a = 0.02 m/s² from the simulations. Note that the acceleration is underpredicted by the simulation in this case. This may be due to the underrepresentation of the wind gusts that drive the change of wind directions, which has been identified as a key factor for fire acceleration.

As the fire approaches the plateau at the end of the slope, the fire growth increases rapidly due the development of the recirculating flow field on the plateau, which can be seen from the instantaneous streamwise velocity fields in Fig. 3. Note that the fluctuations and corrugations of the fire front increase with both slope angle and wind speed. However, Fig. 6b shows that as the wind speed increases, the dependence in fire growth on the slope angle diminishes, suggesting that the wind speed is the dominating factor for the fire acceleration in the present configurations.

We proceed by quantifying the fireline intensity and fuel consumption rate. These quantities are important properties of line fires, measuring the heat release rate of wildfires with direct implication in fire management (Alexander 1982). For this, we compute the fireline intensity as (Finney et al. 2021):

where H_c is the heat of combustion and mc″ is the fuel consumption per unit area in the flame zone, which is determined as:

where Ω is the simulation domain, Ȧ_burn is the rate of change of the burned area, and ω̇_f is the rate of fuel consumption. With the burn area being determined as the area where the bulk fuel density is below 20% of its initial value, its rate of change at time t_n is computed as Ȧburntn=(Aburn t n + 1−Aburn t n − 1)/2Δt. Fig. 7a shows the heat density Hcmc″, computed from the simulations in units of MJ/m². The mean heat density for the present study is found to be 25.9 MJ/m², with a standard deviation of 3.3 MJ/m² and a maximum difference of 13.1 MJ/m² for all simulations. These results show that the heat density only exhibits a weak dependence on wind speed and slope angle, which is in agrees with observations by Finney et al. (2021).

Fig. 7.

Quantitative analysis of ensemble simulations as a function of slope angle and wind speed, showing (a) Heat density, Hcmc″, and (b) Fireline intensity, I_B.

In contrast to the heat density, the fireline intensity shows a strong dependency on the wind speed and slope angle, as shown in Fig. 7b. Given that the fireline intensity is defined as the product between the heat density and ROS, and the heat density stays relatively constant with respect to the wind and slope angle, the fireline intensity shows the same response with the wind and slope as the ROS. With increasing wind speed and slope angle, the fireline intensity is larger. Similar to the ROS, shown in Fig. 6a, the dependency on the wind and slope factors is nonlinear, which is indicated by the relatively faster growth rate with respect to the slope angle for the cases with a faster ambient wind.

Wind-slope effects on fire-spread rate

To examine the effect of the wind and slope conditions on the fire growth, we compare ROS obtained from the simulations with the Rothermel model following Rothermel’s original definition and Nelson’s modifications. The original definition of the Rothermel model considers a linear combination of the wind and slope factors,

where R₀ is the zero-wind, zero-slope ROS, u = U₁₀, and α is the slope of the fuel bed.

The modification by Nelson (2002) to the model adjusts the effective wind speed geometrically by incorporating the slope angle. The vertical velocity due to buoyancy is computed as

where I_B is computed with the assumption that mc″ is the total available fuel loading. With this, the effective wind speed parallel to the slope, assuming that the wind is aligned with the gradient of the slope, is computed as

The ROS is then computed following Rothermel’s formulation with the slope factor omitted:

Eqns 8–10 are solved using Newton iteration to obtain the ROS.

Fig. 8 shows results over the parameter space corresponding to that of the ensemble simulation (see Fig. 6a). The Rothermel model with Nelson’s modification shows that the ROS has a stronger dependency on the slope compared to the original formulation. Different from the original Rothermel model, in which the slope and wind factors are derived from conditions in the absence of wind and slope, respectively, the modified Rothermel model considers the non-linear coupling of the slope on the wind speed. This coupling effect is shown to be important from the simulation results. Good qualitative agreement is observed between the modified Rothermel model and the simulation results. Fig. 9 shows the mean ROS for the modified wind speed u*. u* of the simulation is computed with Eqns 8 and 9 for direct comparisons with Nelson’s analysis. A positive correlation between these two factors with a relatively small variance is observed, suggesting that the coupling between the buoyancy-induced vertical velocity and the wind speed is important for predicting the ROS on a slope. Note that the mean ROS obtained from the simulations is higher than what is predicted by Nelson’s model for the range of u* considered in this study. This may be caused by the change in wind directions due to turbulence that is induced from the fire–atmosphere interactions, which needs further investigations.

Fig. 8.

Comparison of mean ROS predictions from (a) Rothermel model and (b) Nelson’s modified model.

Fig. 9.

ROS from the simulation data for the effective wind speed computed with Eqn 9.

To further investigate the fire dynamics and effects of buoyancy, we compute the convective Froude number (Sullivan 2007)

where u is computed from Eqn 9 to incorporate the impact of both wind speed and slope. Fig. 10 shows the Froude number parameterised by the wind speed and slope. Literature identified that the Froude number can be used as parameter to quantify the transition between a plume-dominated and a wind-dominated fire (Sullivan 2007; Moinuddin et al. 2018). The critical Froude number for this transition is found to be around 0.5 (Moinuddin et al. 2018; Morvan and Frangieh 2018), which corresponds to the lower-left region in Fig. 10, with slopes smaller than 15° and wind speed less than 8 m/s.

Fig. 10.

Froude number computed from Eqn 11, with the white contour indicating the critical Froude number Fr = 0.5.

The transition of the fire spread from the plume-dominated mode to the wind-dominated mode is demonstrated with the mean potential temperature and streamlines near the fire front in Fig. 11, shown in the slope-oriented coordinate system. The origin of the slope-aligned coordinate ξ is shifted to the location of the fire front before computing the mean across the y-direction. The mean ambient wind speed is removed from the velocity field to improve the illustration of the local interaction of the fire and the flow field. The three cases selected are associated with U₁₀ = 4 m/s, with slope angles being 0°, 10°, and 30° that correspond to Fr < 0.5, ~0.5, and >0.5, respectively. In Fig. 11a we see that the fires are almost vertically oriented, which indicates that the buoyancy is dominating the convection from the upstream wind. The large-scale fluid motion is perturbed by the fire. Additionally, air entrains from upstream of the fire, which reduces the convective heat flux from the flame to the unburned fuel bed. In contrast, in the case with Fr > 0.5, as shown in Fig. 11c, the flame attaches close to the ground, and the streamlines aligned along the vertical direction without large-scale perturbations. As a consequence, the unburned fuel is preheated mainly by the convective heat flux from the flame. When Fr ~ 0.5, as in Fig. 11b, the flame leans closer towards the ground compared to Fr < 0.5 in Fig. 11a, and the flow field is less chaotic with weakened upstream air entrainment. The mode of unburned fuel preheat starts to transition from radiation dominated to wind dominated.

Fig. 11.

Comparison of the potential temperature and flow field at the fire front from (a) U04S00, (b) U04S10, and (c) U04S30.

We find that this plume-dominated regime correlates with the regime where large differences in ROS are found between the simulations and the modified Rothermel model. This is due to the strong coupling between the fire dynamics and the fire-induced turbulence in a plume-dominated fire, which requires more investigations.