# Towards predictive data-driven simulations of wildfire spread – Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation

M.C. Rochoux, S Ricci, D Lucor, B Cuenot, A Trouvé

► To cite this version:

M.C. Rochoux, S Ricci, D Lucor, B Cuenot, A Trouvé. Towards predictive data-driven simulations of wildfire spread – Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation. Natural Hazards and Earth System Sciences, 2014, 14 (11), pp.2951-2973. 10.5194/nhess-14-2951-2014. hal-01332351

HAL Id: hal-01332351

https://hal.sorbonne-universite.fr/hal-01332351

Submitted on 15 Jun 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Towards predictive data-driven simulations of wildfire spread – Part I: Reduced-cost Ensemble Kalman Filter based on a Polynomial Chaos surrogate model for parameter estimation

M. C. Rochoux<sup>1,2,3,4</sup>, S. Ricci<sup>1,2</sup>, D. Lucor<sup>5</sup>, B. Cuenot<sup>1</sup>, and A. Trouvé<sup>6</sup>

*Correspondence to:* M. C. Rochoux (melanie.rochoux@graduates.centraliens.net)

Received: 19 March 2014 – Published in Nat. Hazards Earth Syst. Sci. Discuss.: 9 May 2014

Revised: 2 September 2014 – Accepted: 23 September 2014 – Published: 10 November 2014

**Abstract.** This paper is the first part in a series of two articles and presents a data-driven wildfire simulator for forecasting wildfire spread scenarios, at a reduced computational cost that is consistent with operational systems. The prototype simulator features the following components: an Eulerian front propagation solver FIREFLY that adopts a regional-scale modeling viewpoint, treats wildfires as surface propagating fronts, and uses a description of the local rate of fire spread (ROS) as a function of environmental conditions based on Rothermel's model; a series of airborne-like observations of the fire front positions; and a data assimilation (DA) algorithm based on an ensemble Kalman filter (EnKF) for parameter estimation. This stochastic algorithm partly accounts for the nonlinearities between the input parameters of the semi-empirical ROS model and the fire front position, and is sequentially applied to provide a spatially uniform correction to wind and biomass fuel parameters as observations become available. A wildfire spread simulator combined with an ensemble-based DA algorithm is therefore a promising approach to reduce uncertainties in the forecast position of the fire front and to introduce a paradigm-shift in the wildfire emergency response. In order to reduce the computational cost of the EnKF algorithm, a surrogate model based on a polynomial chaos (PC) expansion is used in place of the forward model FIREFLY in the resulting hybrid PC-EnKF algorithm. The performance of EnKF and PC-EnKF is

assessed on synthetically generated simple configurations of fire spread to provide valuable information and insight on the benefits of the PC-EnKF approach, as well as on a controlled grassland fire experiment. The results indicate that the proposed PC-EnKF algorithm features similar performance to the standard EnKF algorithm, but at a much reduced computational cost. In particular, the re-analysis and forecast skills of DA strongly relate to the spatial and temporal variability of the errors in the ROS model parameters.

## 1 Introduction

Real-time prediction of the direction and speed of a propagating wildfire has been identified as a valuable research objective with direct applications in both fire risk management and fire emergency response (Noonan-Wright et al., 2011). In addition, the perspective of climate change tends to favor extreme drought events and to alter precipitations (Milly et al., 2002; Palmer and Räisänen, 2002; Boé et al., 2009); these conditions dramatically increase the risk for the development of large highly destructive wildfires, commonly known as megafires (Nijhuis, 2012). In this context, accurate predictions of the resulting change in fire regime and intensity cannot only rely on the analysis of past observed wildfire events; the use of a data-driven wildfire spread simulator that

<footnote n="1"><sup>1</sup>CERFACS, 42 avenue Gaspard Coriolis, 31057 Toulouse Cedex 01, France</footnote>

<footnote n="2"><sup>2</sup>SUC/CNRS-URA1875, 42 avenue Gaspard Coriolis, 31057 Toulouse CEDEX 01, France</footnote>

<footnote n="3"><sup>3</sup>Ecole Centrale Paris, Grande voie des vignes, 92295 Châtenay-Malabry, France</footnote>

<footnote n="4"><sup>4</sup>EM2C/CNRS-UPR288, Grande voie des vignes, 92295 Châtenay-Malabry, France</footnote>

<footnote n="5"><sup>5</sup>Institut d'Alembert, Université Pierre et Marie Curie, CNRS-UMR7190, 4 place Jussieu, 75006 Paris, France</footnote>

<footnote n="6"><sup>6</sup>Dept. of Fire Protection Engineering, University of Maryland, College Park, MD 20742, USA</footnote>
takes full advantage of the recent technological advances for geo-referenced front-tracking becomes essential.

Despite our recent progress in computer-based wildfire spread modeling, our ability to accurately simulate the behavior of wildfires remains limited because the underlying dynamics feature complex multi-physics processes occurring at multiple scales (Viegas, 2011). The dynamics of wildfires are determined by interactions between pyrolysis, combustion and flow dynamics, radiation and convection heat transfer, as well as atmospheric dynamics and chemistry. These interactions occur at the following scales: vegetation scales that characterize the biomass fuel; topographical scales that characterize the terrain and vegetation boundary layer; and meteorological micro/meso-scales that characterize atmospheric conditions.

Relevant insight into wildfire dynamics has been obtained in recent years via detailed numerical simulations performed at flame scales (i.e., with a spatial resolution of the order of 1 m). For instance, FIRETEC (Linn et al., 2002) or WFDS (Mell et al., 2007) combine advanced physical modeling and classical methods of computational fluid dynamics (CFD) to accurately describe the combustion-related processes that control the fire behavior (e.g., thermal degradation of biomass fuel, buoyancy-induced flow, combustion, radiation and convection heat transfer). Note that because of the high computational cost, flame-scale CFD is currently restricted to research projects (Linn et al., 2002; Mell et al., 2007; Rochoux, 2014) and is not compatible with operational applications. In contrast, a regional-scale viewpoint (i.e., a viewpoint that considers scales ranging from a few tens of meters up to several kilometers) is adopted in the following: the fire is described as a two-dimensional front that self-propagates normal to itself into unburnt vegetation; the local propagation speed is called the rate of spread (ROS). This viewpoint is the dominant approach used in current operational wildfire spread simulators, see for instance FARSITE (Finney, 1998), FOREFIRE (Filippi et al., 2009, 2013), PROMETHEUS (Tymstra et al., 2010) and PHOENIX RapidFire (Chong et al., 2013). In particular, FARSITE uses a model due to Rothermel (1972) that treats the ROS as a semi-empirical function of biomass fuel properties associated with a pre-defined fuel category (i.e., the vertical thickness of the fuel layer, the fuel moisture content, the fuel particle surface-to-volume ratio, the fuel loading and the fuel particle mass density), topographical properties (i.e., the terrain slope) and meteorological properties (i.e., the wind velocity at mid-flame height). This approach is limited in scope because of the large uncertainties associated with the accuracy of computer models since they do not account for the interaction between the fire and the atmosphere, and since they have a limited domain of validity resulting from a calibration procedure based on experiments (Perry, 1998; Sullivan, 2009; Viegas, 2011; Cruz and Alexander, 2013; Finney et al., 2013). This approach is also limited because of the large uncertainties associated with many of the input param-

eters to the fire problem (Jimenez et al., 2007; Finney et al., 2011).

In order to overcome some of the current limitations of regional-scale wildfire modeling and to build predictive simulations that are compatible with operational framework, the uncertainties in the input data of the ROS semi-empirical model need to be quantified and reduced. The uncertainties inherent in wildfire spread modeling go beyond the limitations of deterministic forecast abilities of the dynamical model (also referred to as the *forward model*) and thus, suggest the use of ensemble forecasts to stochastically characterize the nonlinear response of the front-tracking simulator to variations in the input environmental parameters (D'Andrea et al., 2010; Finney et al., 2011). For instance, Finney et al. (2011) describes an ensemble-based forecasting capability, in which a large number of fire spread scenarios (i.e., the ensemble members) are generated based on a probabilistic uncertainty in the weather conditions and in the moisture content of biomass fuels. Model uncertainties are a combination of epistemic errors that express an imperfect knowledge of the input parameters of the ROS model (that could in theory be removed), and of aleatoric errors that result from natural and unpredictable stochastic variabilities of the physical system (that can be addressed by stochastic models, see for instance Reference by Pagnini and Mentrelli (2014), whose model relies on a stochastic component to represent the transport of firebrands). These uncertainties translate inevitably into errors in the output variables of interest (e.g., time-evolving position of the front, burnt area, maximum value for the ROS). The most classical methodologies for uncertainty quantification in these output variables are random sample-based statistical methods derived from Monte Carlo methodologies. While these methodologies are generic and robust for the simulation of stochastic models, they are however computationally expensive due to the required size of the sample (the computational cost of one realization may be already expensive itself, see Lucor et al., 2007) and each implementation typically requires *ad-hoc* variance reduction techniques (Boyaval, 2012). More efficient sampling methods have been developed to reach a comparable level of accuracy as Monte Carlo-based techniques but with fewer forward model integrations; these sampling methods take advantage of the (possible) regularity of the model response to varying input parameters in order to increase the convergence rate compared to Monte Carlo-based methodologies. In particular, polynomial chaos (PC) non-intrusive techniques issued from spectral-based representations and introduced by Wiener (1938) are very often efficient in terms of precision and cost (Ghanem and Spanos, 1991; Le Maître and Knio, 2010). The key idea is to build a polynomial representation of the forward model response (referred to as the *surrogate model*) to varying input parameters. Once the surrogate model is available, it is possible to benefit from a large sample of realizations (at almost no cost) in order to accurately characterize the model uncertainties. Still, the application
of PC-based sampling techniques to problems of hyperbolic conservation laws remains a challenging task (Desprès et al., 2013).

Recent progress made in airborne remote sensing provides new ways to monitor real-time fire front positions (Wooster et al., 2005, 2013; Riggan and Robert, 2009; Paugam et al., 2013). Unfortunately, these thermal-infrared measurements provide an incomplete description of the fire spread (in particular due to the opacity of the fire-induced thermal plume and/or due to a limited monitoring) and are subject to instrumental errors as well as representativeness errors (i.e., inconsistency between what the sensor can measure and what the computer model can describe). From this perspective, data assimilation (DA) offers a convenient framework for integrating fire sensor observations into a computer model in order to provide optimal estimates of poorly known model parameters and/or model state, and to improve in fine predictions of the fire spread behavior (Mandel et al., 2008; Cowlard et al., 2010; Lautenberger, 2013; Rochoux, 2014). The key idea is that, when used alone, neither measurements nor computer models can provide a reliable and complete description of the real state of the physical system. In the following, the set of model state and/or model parameters to be corrected through DA is gathered in the control vector. The DA algorithm is sequentially applied; each sequence (also referred to as the *assimilation cycle*) is decomposed into two steps: (1) a prediction step, in which the control variables are advanced in time given some uncertainty ranges; and (2) an update step based on the classical Bayes' theorem, in which new observations are considered and the probability density function (PDF) of the control variables is modified consistently with the observations in order to reduce the uncertainties in the model outputs (Gelb, 1974; Tarantola, 1987; Todling and Cohn, 1994; Ide et al., 1997; Kalnay, 2003; Reichle, 2008). The Kalman filter (KF) is the most commonly used sequential DA technique. However, the KF assumes linear dynamics between the control variables and the model outputs as well as a Gaussian statistical distribution for both modeling and observation errors. Extensions of the KF that partly overcome these limitations have been proposed, for instance the extended Kalman filter (EKF) that uses local linearization techniques (Gelb, 1974) or the ensemble Kalman filter (EnKF) that relies on a stochastic description of the model behavior (Evensen, 1994, 2009). An insightful comparison between EKF and EnKF is given within the framework of land DA in Reichle et al. (2002).

In this study, an ensemble-based DA methodology is considered in order to reduce the uncertainties in the ROS model parameters using measurements of the time-evolving location of the fire front. This study is an extension of our previous works presented in Rochoux et al. (2013a, b), in which a prototype data-driven wildfire simulator was developed. The initial prototype featured the following main components: an Eulerian front-tracking solver combined with a model description of the local ROS proposed by Rothermel (1972);

a series of observations of the fire front position; and a cost-effective EKF-based DA algorithm. This prototype was successfully evaluated when applied for estimating the input parameters of the Rothermel-based ROS model (e.g., the fuel moisture content, the fuel particle surface-to-volume ratio, and/or the wind direction and magnitude). However, the EKF algorithm relies on the assumption that the relation between a perturbation in the ROS model parameters and the resulting changes in the fire front position (i.e., the generalized observation operator) can locally be approximated by a linear relation. While the EKF-based studies presented in Rochoux et al. (2013a, b) produced encouraging results and confirmed the value of a DA strategy for improved wildfire spread predictions, the linearity assumption is no longer valid in regional-scale fires, especially when the wind direction and magnitude vary and the vegetation properties are potentially strongly heterogeneous. To better account for nonlinearities in the generalized observation operator, an extension to an EnKF approach was preliminarily explored in Rochoux et al. (2012). This ensemble-based DA approach was originally developed for dynamic state estimation (Evensen, 1994) and has already been used in the field of wildfire modeling for correcting the temperature state variable (Beezley and Mandel, 2008; Mandel et al., 2008, 2011). It was also largely extended to sequential parameter estimation, for instance in the field of hydrology (Durand et al., 2008; Moradkhani et al., 2005). Still, the large number of realizations required by the EnKF algorithm to obtain satisfactory results (Rochoux, 2014) may prove computationally burdensome within an operational framework. This behavior of the EnKF algorithm for parameter estimation is due to four main reasons: (1) the slow convergence rate of the Monte Carlo sampling; (2) the nonlinear interrelation between the control space and the observation space; (3) the complexity of retrieving the specific signature of each control parameter on the resulting distribution of the simulated fire front; and (4) the accumulation of sampling errors along assimilation cycles that can only be addressed by increasing the size of the sample (Li and Xiu, 2008). The required size of the sample significantly increases with the complexity of the physics (multi-parameter estimation) and the model nonlinearities (complex physics), thus emphasizing the need for a reduced-cost EnKF. Efforts have therefore been devoted to designing more efficient EnKF schemes by reducing sampling errors (Szunyogh et al., 2008; Saad, 2007; Li and Xiu, 2008, 2009; Blanchard et al., 2010; Xiu, 2010; Rosić et al., 2013). For this purpose, and following work from Li and Xiu (2009), an EnKF strategy based on a PC approximation (PC-EnKF) is proposed in this paper; the polynomial surrogate model being used during the EnKF prediction step to generate a large number of model simulation trajectories at almost no cost and without loss of accuracy (Birolleau et al., 2014).

In this paper, we present a hybrid PC-EnKF DA algorithm that improves wildfire spread modeling by reducing uncertainty in the vegetation properties used as inputs
of the Rothermel-based ROS model. The objective of this
study is to show the feasibility of this approach for wildfire
spread forecasting under several assumptions, i.e., a mini-
malist treatment of the fire front (idealized as an interface and
consistent with the limited knowledge on the environmen-
tal conditions); a semi-empirical formulation of the ROS;
Gaussianity of the errors on the input parameters of the ROS
model and on the observations; prior values for the control
parameters specified based on user-defined mean and error
standard deviation (STD). In this first part, both the EnKF
and PC-EnKF algorithms are limited to the estimation of spa-
tially uniform parameters of the ROS model due to compu-
tational cost constraints and a lack of high-resolution data
on the environmental conditions. Although it seems appro-
priate to translate the inability of a fire spread model to gen-
erate accurate fire front positions into parameter uncertainty,
other sources of uncertainties such as model structural er-
rors or boundary/initial condition errors also need to be ac-
counted for. For this purpose, in the second part of this se-
ries of two articles (Rochoux et al., 2014), a state estimation
strategy is designed to address anisotropic uncertainties in
wildfire spread as well as to provide observation-informed
initial condition for model integration at future lead times.
Thus, parameter estimation and state estimation are comple-
mentary approaches that are valuable for wildfire behavior
forecasting; it is therefore important to discuss their benefits
and drawbacks for experiments with increasing complexity.

The outline of the paper is as follows. Section 2 presents the available observations of the fire behavior and the wildfire spread model named FIREFLY (i.e., the forward model). The hybrid PC-EnKF algorithm developed for the wildfire application is presented in Sect. 3; in this section, the sequential implementation of the ensemble-based algorithms is also described. Section 4 illustrates how the classical EnKF and the hybrid PC-EnKF allow to properly estimate model parameters on simple test cases, in which the observations are synthetically generated. The performance of the data-driven wildfire spread capability using the reduced-cost approach is demonstrated in a validation test corresponding to a controlled grassland fire experiment.

## 2 Information on wildfires at regional scales: observations and forward model

### 2.1 Observations of the fire front location

#### 2.1.1 Overview of available observations of fire spread

In practice, continental surfaces and vegetation are mainly
observed within the mid- and near-infrared regions of the
electromagnetic spectrum of wavelengths (0.75 to 15 µm). It
is known that for high temperatures as encountered in wild-
fires (varying from 600 K for smoldering to 1200 K in the
flaming zone), the maximum radiant intensity occurs within

the mid-infrared (MIR) region. Thus, current spaceborne and
airborne systems observe wildfires within a narrow wave-
band centered on the 3.9-micron wavelength (Butler et al.,
2004; Wooster et al., 2005, 2013; Paugam et al., 2013; Ro-
choux, 2014), which is both sensitive to flaming and smolder-
ing combustion modes. Beyond fire detection, remote sens-
ing is regarded as a promising approach to provide a quan-
titative description of the fire radiation release to charac-
terize sub-pixel fires (occupying a limited area of the sen-
sor pixel down to 0.1 to 1% of the pixel area) and to esti-
mate fuel consumption as well as smoke emissions (Wooster
et al., 2013). Using spaceborne or airborne platforms, the fire
radiative power (FRP) emissions are detected in the burn-
ing area, while non-active areas remain blank. This infor-
mation is crucial to retrieve the brightness temperature and
thus, to track the time-evolving location of the fire front. For
instance, Paugam et al. (2013) showed that spatiotemporal
variations of the flame front ROS can be accurately retrieved
using FRP analysis on a reduced-scale controlled fire experi-
ment (the final burnt area of the reduced-scale study is about
1000 m<fnref n="2" />); ongoing research aims at extending this FRP anal-
ysis to regional-scale wildfire spread.<fnref n="1" />,<fnref n="2" />

Currently, most spaceborne instruments, including the pi-
oneer generation such as the AVHRR (Advanced Very High
Resolution Radiometer) and the MODIS (MOderate resolu-
tion Imaging Spectroradiometer), offer neither a sufficiently
short revisit period nor a high enough spatial resolution im-
agery for efficient front-tracking at regional scales. While
these objectives no longer seem out of reach for the dual
SPOT-Pléiades constellation,<fnref n="3" /> airborne platforms still seem
the most suitable solution for real-time geo-location of ac-
tive fire contours. Typical examples are the LIVEFIRE sys-
tem (Merlet, 2008; Crombette, 2010) and its US counterpart
FIREMAPPER system deployed since 2004 by the US Forest
Service and the US Department of Interior Bureau of Land
Management (Riggan and Robert, 2009). As a complement,
spaceborne data could be used for validation as well as cali-
bration of models and DA procedures.

2.1.2 Choice of observations for data assimilation

In the present study, we assume that observations of the
fire front position are available and that these observations
can be made at different relevant times with a low mea-
surement error (typically, 0–30 m for the LIVEFIRE sys-
tem). In the following, the observed fire front is represented
as a segmented line using a pre-defined number of equally
spaced markers (i.e., the $N_{fr}^o$ observation markers); the obser-
vation vector noted $y_i^o$ contains the two-dimensional coordi-
nates $(x_i^o, y_i^o)$ of the fire front markers (the subscript $i$ is the

<footnote n="1">¹http://wildfire.geog.kcl.ac.uk/</footnote>

<footnote n="2">²http://gofc-fire.umd.edu/</footnote>

<footnote n="3">³https://directory.eoportal.org/web/eoportal/satellite-missions/
p/pleiades</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Eulerian front-tracking simulator FIREFLY. Left: the fire front is the $c_{fr} = 0.5$ contour line; $\Gamma$ measures the local ROS of the fire along the normal direction to the front $n_{fr}$ (defined by the direction angle of fire propagation $\alpha_{fr}$) given the wind velocity vector $(u_w^*, \alpha_w^*)$. Right: profile of the spatial variations of the progress variable $c$ across the fire front, $(x_i, y_i)$ representing the location of the $i$th fire front marker.</figcaption></figure>

index of a particular marker in the observation vector, with $i = 1, \dots, N_{fr}^o$ observed at the analysis time $t$. The size of the observation vector $y_t^o$ is $2N_{fr}^o$. The coordinates of the fire front markers are assumed to have independent Gaussian-like random errors $\epsilon^o$ with zero mean and with STD $\sigma^o$. Note that this classical assumption of uncorrelated observation errors could be questionable. However, this aspect is out of the scope of this study and is still under active research in the DA field (Brankart et al., 2009; Gorin and Tsyrulnikov, 2011).

Two types of experiments are presented in the following: observation system simulation experiments (OSSE), in which observations are synthetically generated using a reference solution of the FIREFLY fire spread model (called the true evolution) that is modified by random observation errors $\epsilon^o$; and a controlled grassland fire experiment, in which the observations are reconstructed from measured temperature maps and using a definition of the fire front as the 600 K temperature contour line.

## 2.2 The fire spread model (the forward model)

The front-tracking solver, called FIREFLY and formally noted $M$ in the following, simulates the propagation of surface wildfires within the biomass fuel bed and at regional scales, as illustrated in Fig. 1. Note that the present study is limited to flat terrains and problems with complex topography are outside its scope. FIREFLY tracks the time-evolving location of the fire front using the following three components: (1) a sub-model for the ROS noted $\Gamma$; (2) an Eulerian front-tracking solver for the fire front propagation equation; and (3) an isocontour algorithm for the reconstruction of the fire front.

### 2.2.1 The Rothermel-based rate of spread sub-model

#### (a) Original one-dimensional formulation

The ROS sub-model is based on the widely used semi-empirical model due to Rothermel (1972) that describes $\Gamma$ as a function of the local environmental conditions (e.g., vegetation and weather properties). The ROS is derived from the one-dimensional formulation of the energy balance equation per unit volume of the unburnt biomass fuel located ahead of the flame; the physical quantities involved in this energy balance are then parameterized using wind-tunnel experiments. In this formulation, $\Gamma$ [m s$^{-1}$] is expressed as the ratio between the heat flux received by the unburnt vegetation $I_p$ [J m$^{-2}$ s$^{-1}$] and the energy required to ignite the fuel $H_{ig}$ [J m$^{-3}$]. $\Gamma$ reads as follows:

$$ \Gamma = \frac{I_p}{H_{ig}} = \frac{\xi I_r}{\rho_b \chi Q_{ig}} (1 + \Phi_w). \quad (1) $$

$I_p$ is a function of the energy release rate of the combustion $I_r$, of the dimensionless propagating flux ratio $\xi$ (that describes the proportion of energy that is released by the flame and transferred to the vegetation in the non-flaming zone). The wind correction coefficient $\Phi_w$, which was determined for a one-dimensional case corresponding to a head fire configuration, nonlinearly depends on the wind velocity magnitude at mid-flame height $u_w$ such that

$$ \Phi_w \equiv \Phi_w(u_w) = C_w u_w^{B_w} \left( \frac{\beta_v}{\beta_v, \text{opt}} \right)^{-E_w}, \quad (2) $$

with $C_w$, $B_w$ and $E_w$ calibrated parameters depending on the biomass fuel surface-to-volume ratio $\Sigma_v$ [m$^{-1}$], with $\beta_v$ the biomass fuel packing ratio and $\beta_{v, opt} \equiv \beta_{v, opt}(\Sigma_v)$ its optimum value (optimum meaning that $\beta_{v, opt}$ characterizes the optimum arrangement of the biomass fuel particles that produces the most effective mixing between air and fuel gas reactants for the occurrence of combustion). The ignition energy $H_{ig}$ is formulated as $H_{ig} = \rho_b \chi Q_{ig}$, with $Q_{ig}$ [J kg$^{-1}$] the heat of ignition, $\chi$ the dimensionless effective heating number (i.e., amount of fuel effectively involved in the ignition process) and $\rho_b$ [kg m$^{-3}$] the biomass fuel bulk mass density that satisfies $\rho_b = \beta_v \rho_p$ for a porous medium, $\rho_p$ [kg m$^{-3}$] being the biomass fuel particle mass density.

The expression for the local ROS $\Gamma$ due to Rothermel may be written in the following compact form that is equivalent to Eq. (1):

$$ \Gamma \equiv \Gamma(\delta_v, M_v, M_{v, ext}, \Sigma_v, m_v'', \rho_p, \Delta h_c, u_w), \quad (3) $$

where the nomenclature for the input parameters are summarized in Table 1. Note that the fuel loading $m_v''$ [kg m$^{-2}$] satisfies $m_v'' = \rho_b \delta_v / (1 + s_t) = (\beta_v \rho_p) \delta_v / (1 + s_t)$, $s_t$ being the fuel particle total mineral content.
<table><caption>**Table 1.** Input parameters of the Rothermel-based ROS sub-model.</caption><thead><tr><th>Name</th><th>Symbol</th><th>Unit</th></tr></thead><tbody><tr><td>Fuel depth (vertical thickness of the vegetation layer)</td><td>δ<sub>v</sub></td><td>m</td></tr><tr><td>Fuel moisture (mass of water divided by mass of dry vegetation)</td><td>M<sub>v</sub></td><td>%</td></tr><tr><td>Fuel moisture at extinction</td><td>M<sub>v, ext</sub></td><td>%</td></tr><tr><td>Fuel particle surface-to-volume ratio</td><td>Σ<sub>v</sub></td><td>m<sup>-1</sup></td></tr><tr><td>Fuel loading</td><td>m''<sub>v</sub></td><td>kg m<sup>-2</sup></td></tr><tr><td>Fuel particle mass density</td><td>ρ<sub>p</sub></td><td>kg m<sup>-3</sup></td></tr><tr><td>Fuel heat of combustion</td><td>Δh<sub>c</sub></td><td>J kg<sup>-1</sup></td></tr><tr><td>Wind velocity magnitude at mid-flame height (projected onto horizontal plane)</td><td>u<sub>w</sub></td><td>m s<sup>-1</sup></td></tr></tbody></table>

## (b) Extension to two-dimensional surface wildfire spread

The original Rothermel's one-dimensional model is extended to two-dimensional configurations, in order to account for the wind effects on the shape of the fireline, while still maintaining a simple parameterization of the ROS with respect to local environmental conditions. Accounting for wind-induced wildfire spread in FIREFLY is such that when the wind blows in the direction of the fire spread (i.e., a head fire configuration), the wind contribution to the ROS is maximum. On the contrary, the wind contribution to the ROS is zero when the wind blows in the direction opposite to the direction of the fire spread (i.e., a rear fire configuration), meaning that the fire propagates at the value of no-wind ROS on this section of the fire front (i.e., $Φ_w = 0$). On the flanks, the fire front advances faster than in the absence of wind (i.e., $Φ_w > 0$). This implies that the ROS can drastically change along the fireline at a given time. For this purpose, characteristic angles in the horizontal plane $(x, y)$ are defined to represent the direction angle of the wind noted $α_w^*$ and the direction angle of the fire propagation noted $α_{fr}$ (the index fr referring to the front); $α_{fr}$ indicates the outward-pointing normal direction to the fire front noted $n_{fr}$ (see Fig. 1). These angles are defined from the northern direction, namely from the positive y co-ordinates and increasing in the clockwise direction. Since the propagation of wildfires is anisotropic, the normal vector $n_{fr}$ is not uniform along the fireline and is modified over time, with

$$n_{fr} = n_{fr}(x, y, t) = \left( \sin \alpha_{fr}(x, y, t) \cos \alpha_{fr}(x, y, t) \right). \quad (4)$$

Thus, the wind velocity magnitude at mid-flame height $u_w$ (see Table 1) corresponds to the projection of the wind velocity vector $u_w^*$ along the normal direction to the front $n_{fr}$:

$$u_w = u_w(x, y, t) = u_w^* \cdot n_{fr}(x, y, t), \quad (5)$$

with $u_w^*$ defined by its magnitude, $u_w^*$ [m s$^{-1}$], and direction angle, $α_w^*$ [$^\circ$]:

$$u_w^* = \begin{pmatrix} u_w^* \sin \alpha_w^* \\ u_w^* \cos \alpha_w^* \end{pmatrix}. \quad (6)$$

In the following, $u_w^*$ and $α_w^*$ are treated as spatially uniform and time-independent. The projected wind velocity at mid-flame height $u_w = u_w(x, y, t)$ is a time-dependent and spatially varying quantity along the propagating fireline. It is worth noting that the wind contribution $Φ_w$ is forced to a zero-value in FIREFLY when the scalar product $u_w^* \cdot n_{fr}(x, y, t)$ is negative (see Eq. 5) to ensure that the ROS $Γ$ remains positive. This is consistent with the common assumption in the field of fire spread modeling that the fire propagates at least at the no-wind ROS. As for biomass fuel properties, the fuel depth $δ_v = δ_v(x, y)$ is treated as a time-independent, spatially varying quantity; all other ROS model parameters are treated as constant and uniform.

## (c) Sensitivity study of the rate of spread

The identification of which parameters are important to include in the control vector (denoted by $x$) is an essential step towards the application of DA to FIREFLY. The key idea when dealing with parameter estimation is to focus the correction on a reduced set of parameters that have significant uncertainties and to which FIREFLY is the most sensitive.

In order to identify to which input parameters the ROS $Γ$ is the most sensitive among biomass fuel properties and weather conditions, a sensitivity study is carried out with the classical one-dimensional Rothermel's model for short grass. Nominal environmental conditions are as follows: the head fire propagates in presence of a moderate wind $u_w^* = 1$ m s$^{-1}$, and the vegetation is characterized by the moisture content $M_v = 20\%$, the particle surface-to-volume ratio $Σ_v = 11485$ m$^{-1}$, the layer thickness $δ_v = 0.5$ m, and the layer packing ratio $β_v = 0.106\%$. These four parameters are perturbed around these nominal conditions. Note that the moisture at extinction $M_{v, ext} = 30\%$, the fuel particle mass density $ρ_p = 512.6$ kg m$^{-3}$, the effective fuel mineral content $s_e = 1\%$ ($s_t = 5.55\%$), and the heat of combustion $Δh_c = 1.861 \times 10^7$ J kg$^{-1}$ remain constant and correspond to the standard values of the Rothermel's fuel database (Rothermel, 1972).

Figure 2 compares the variability in the ROS $Γ$ when uncertainties are assumed in four parameters, $u_w^*$, $M_v$, $Σ_v$ and $β_v$. It is found that the ROS values are the most sensitive to
<figure><img src="image_4.png" /><figcaption>Figure 2. Sensitivity of the Rothermel-based ROS Γ to environmental parameters; nominal conditions are indicated by vertical lines.</figcaption></figure>

$u_w^*$ As for biomass fuel properties, they feature a wide scatter for $M_v$ and $\Sigma_v$, while $\Gamma$ is less sensitive to $\beta_v$, indicating that a lack of information in $u_w^*$, $M_v$ and $\Sigma_v$ results in a significant uncertainty range in the ROS model predictions. It is also shown that the ROS $\Gamma$ depends nonlinearly on the pair of parameters $M_v$ and $\Sigma_v$; in particular, there is a ROS acceleration when the biomass fuel becomes drier or when the biomass fuel particles become thinner. Note that these nonlinearities are more important when the wind magnitude fluctuates over time or when the fire active area is covered by heterogeneous biomass fuels. This highlights the importance of applying a DA methodology able to handle multiple sources of nonlinearities in the fire spread model.

### 2.2.2 The Eulerian front-tracking solver

An Eulerian front-tracking solver is used to propagate the fire front at the Rothermel-based ROS. FIREFLY adopts a classical approach taken from the premixed combustion literature (Poinsot and Veynante, 2005), in which a reaction progress variable noted $c \equiv c(x, y, t)$ is used as the prognostic variable of the solver and is introduced as a flame marker: $c = 0$ in the unburnt vegetation, $c = 1$ in the burnt vegetation, and the flame front is identified as the contour line $c_{fr} = 0.5$ as illustrated in Fig. 1. In the Eulerian front propagation technique, the progress variable $c$ is calculated as a solution of the following propagation equation:

$$ \frac{\partial c}{\partial t} = -\boldsymbol{\gamma} \cdot \nabla c = \Gamma |\nabla c|, \qquad (7) $$

with $\Gamma = \boldsymbol{\gamma} \cdot \mathbf{n}_{fr}$ the projected ROS given by Eq. (3) and defined along the normal direction to the fire front that satisfies $\mathbf{n}_{fr} = -\nabla c / |\nabla c|$.

Equation (7) is solved using a second-order Runge-Kutta scheme for time-integration and an advection algorithm for spatial discretization based on a second-order total variation diminishing (TVD) scheme combined with a Superbee slope limiter (Rehm and McDermott, 2009; Mallet et al., 2009). Note that FIREFLY requires a two-dimensional field $c(x, y, t-1)$ as initial condition of any time period $[t-1, t]$. This initial condition is constructed such that the transition between $c = 0$ and $c = 1$ is smooth; a tangent hyperbolic function is used to represent this transition.

The validation of the FIREFLY Eulerian front-tracking solver was presented in prior works (Rochoux et al., 2013a; Rochoux, 2014). Model diagnostics were developed to ensure the correct numerical behavior of FIREFLY. These diagnostics were derived from the Kolmogorov-Petrovsky-Piskounov (KPP) analysis valid for uniform fuel conditions (Poinsot and Veynante, 2005) and extended to heterogeneous biomass fuel for application to wildfire spread. They verify that the rate of change of the progress variable $c$ matches the average ROS along the fireline and also that the ROS at the head of the fire is consistent with the Rothermel's 0-D formulation (see Eq. 3). In addition, they also verify that the front thickness, estimated as the average inverse of the maximum gradient of $c$, remains small (i.e., a few mesh step-sizes) and relatively constant over time. In all tests performed to date, these diagnostics have showed the non-diffusive behavior of the numerical scheme underlying FIREFLY, consistently with the physics of the fire spread problem. Further details are provided in Rochoux (2014).

### 2.2.3 Reconstruction of the simulated fire front and comparison with the observed fire front

Once the spatiotemporal variations of the progress reaction $c$ are known, the position of the fire front is extracted using a simple isocontour algorithm such that, formally, the outputs of the FIREFLY model are as follows:

$$ [(x_i, y_i), 1 \le i \le N_{fr}] = \mathcal{M}_{[t-1, t]}(c_{t-1}, \lambda), \qquad (8) $$

where $(x_i, y_i)$ represents the two-dimensional coordinates of the $N_{fr}$ front markers obtained at time $t$ (the index $i$ indicating the marker), where $c_{t-1}$ designates the initial condition (i.e., the spatial distribution of the progress variable $c$ at time $(t-1)$), and where $\lambda$ designates the list of input parameters of the ROS model presented in Table 1, $\lambda = (\delta_v, M_v, M_{v, ext}, \Sigma_v, m_v'', \rho_p, \Delta h_c, u_w)$.

The correction provided by the DA algorithm relies on the comparison between the FIREFLY simulated fire front
<figure><img src="image_5.png" /><figcaption>Figure 3. Construction of the differences between simulated fire front (SFF) and observed fire front (OFF) noted $d_t = [d_{t,1}, \dots, d_{t,N_{\text{fr}}^0}]$. In this illustration, $r = N_{\text{fr}}/N_{\text{fr}}^0 = 4$.</figcaption></figure>

(SFF) described by the $N_{\text{fr}}$ markers (corresponding to a fine-grained discretization of the front) and the observed fire front (OFF) at time $t$. Since observations of the fire front position are likely to be provided with a much coarser resolution and since they may cover only a fraction of the fire front perimeter, the OFF is discretized with a set of $N_{\text{fr}}^0$ markers such that the observation vector $y_t^0$ reads:

$$y_t^0 = [(x_1^0, y_1^0), (x_2^0, y_2^0), \dots, (x_{N_{\text{fr}}^0}, y_{N_{\text{fr}}^0})], \quad (9)$$

with $N_{\text{fr}}^0$ much lower than $N_{\text{fr}}$. In order to compare SFF with OFF, a selection operator $\mathcal{H}$ is introduced. This operator pairs a subset of $N_{\text{fr}}^0$ markers along SFF with the $N_{\text{fr}}^0$ markers along OFF, associating each marker of OFF with its closest neighbor along SFF (see Fig. 3). Preliminary tests reported in Rochoux (2014) have shown that a simple treatment (taking 1 out of every $r$ points) provides reasonable results. Thus, $N_{\text{fr}}^0 = (N_{\text{fr}}/r)$, where $r$ is an integer taking values much larger than 1 that represents the difference in resolution between SFF and OFF. One of the advantages of this representation of the fire fronts is that it provides a local information on the discrepancies between SFF and OFF, and not only a global information such as the difference in the burnt area or in the fireline perimeter. This local information is efficient at representing the anisotropy in wildfire spread.

It is worth noting that the topology of the fire front can be complex in real-world wildfire spread cases, and/or only a section of the fire front can be observed due to the opacity of the fire-induced thermal plume or due to a limited monitoring. Thus, the pairing between simulated markers and observed markers becomes more challenging for complex fire front topology. The generalization of this treatment to complex fire front topology is out of the scope of this study.

It is also worth mentioning that the EnKF and PC-EnKF parameter estimation strategies presented in this paper are valid for any fire spread model; FIREFLY could readily be replaced by any other front-tracking wildfire spread simulator, for instance FARSITE, FOREFIRE, PROMETHEUS or PHOENIX RapidFire.

# 3 Data assimilation algorithm: the polynomial chaos-based ensemble Kalman filter

## 3.1 The standard ensemble Kalman filter

We present here the ensemble Kalman filter (EnKF) algorithm applied, in the context of parameter estimation, for one assimilation cycle between time $(t-1)$ and time $t$.

### 3.1.1 Definition of the control space

The vector $x_t \in \mathbb{R}^n$ corresponds to the control vector that includes the $n$ uncertain parameters to be estimated over the assimilation cycle $[t-1, t]$. This implies that the location of the fire front is not estimated by the EnKF but is indirectly modified by integrating again the fire spread model over the time window $[t-1, t]$ with the newly estimated control parameters. Note that a parameter estimation approach can be considered by itself as an estimation problem and does not need to be combined with a state estimation approach to obtain an optimal EnKF (Pétron et al., 2002; Peters et al., 2005, 2007; Moradkhani et al., 2005; Durand et al., 2008; Ruiz et al., 2013a).

In the present study, the control parameters $x_t$ are assumed global (i.e., spatially uniform) and constant over the time window $[t-1, t]$; they are only modified when moving to the next time window $[t, t+1]$.

### 3.1.2 Generalized observation operator

The generalized observation operator $\mathcal{G}_t$ maps the control space of $x_t$ onto the observation space of $y_t^0$. Within the framework of parameter estimation, $\mathcal{G}_t$ is a composition of the fire spread model $\mathcal{M}_{[t-1,t]}$ (providing the $N_{\text{fr}}$ front marker locations associated with a realization of the control vector $x_t$) with the selection operator $\mathcal{H}_t$ (taking 1 out of every $r = N_{\text{fr}}/N_{\text{fr}}^0$ markers along SFF at time $t$). Formally, $\mathcal{G}_t$ reads:

$$y_t = \mathcal{G}_t(x_t) = \mathcal{H}_t \circ \mathcal{M}_{[t-1,t]}(c_{t-1}, \lambda', x_t), \quad (10)$$

with $y_t$ the location of the $N_{\text{fr}}^0$ fire front markers associated with a set of control parameters $x_t$ at time $t$ (corresponding to the model counterparts of the observed quantities), and with $\lambda'$ the input parameters of the Rothermel-based model that are not included in the control vector $x_t$. In the following, both $x_t$ and $y_t^0$ are considered as random variables.

The observation operator $\mathcal{G}_t$ defined in Eq. (10) is time-dependent since OFF is dynamically evolving: the selection procedure $\mathcal{H}_t$ depends on the location and on the topology of the fire front at a given time, implying that the observation operator is not the same for all members of the ensemble. This formulation is an unusual application of the EnKF algorithm.
<figure><img src="image_6.png" /><figcaption>Figure 4. Flowchart of the EnKF algorithm during the [t − 1 , t ] as- simulation cycle for a parameter estimation approach. Data random- ization (Burgers et al., 1998) is used in the EnKF with ξ o , ( k ) f o- lowing observation error statistics for each member k = 1 , · · · , N$_{e}$ .</figcaption></figure>

3.1.3 Sequential estimation

The EnKF algorithm is sequentially applied over an assimi- lation window [ t − 1 , t ]; each assimilation cycle decomposes into two successive steps for each member of the ensemble indexed by the exponent k as illustrated in Fig. 4:

1. a prediction step (forecast), in which the system is evolved from time ( t − 1) to time t ( t being the next observation time) through an integration of FIREFLY to forecast the fire front position y$_{t}$ given some uncer- tainty ranges in the control vector x$_{t}$ . We note p $^{f}$( x$_{t}$ ) this PDF of the control vector (also called the forecast PDF) at time t . We also note F [ t − 1 , t ] the operator de- scribing the temporal evolution of the control parame- ters from time ( t − 1) to time t , with x$_{t}$ = F [ t − 1 , t $_{]}$( x$_{t}$$_{−}$$_{1}$ ). A temporal evolution of the control vector is introduced here to fit with the classical description of the EnKF al- gorithm: since there is no dynamic model available to describe the evolution of the control parameters, persis- tence forecasting is used to relate the forecast control parameters to the analysis control parameters (Peters et al., 2005). For this purpose, two techniques are re- ported in the literature, inflation on the one hand, a ran- dom walk model on the other hand (West, 1993; Morad- khani et al., 2005; Ruiz et al, 2013b). In this study, the parameter evolution model F [ t − 1 , t ] is artificially set up using a random walk model (see Eq. 23).

2. an update step (analysis), in which new observations are considered at the analysis time t and the forecast PDF of the control parameters is modified consistently with the observations y o t , in order to reduce the uncertainties in the computer model outputs y$_{t}$ . The new PDF, called the analysis and noted p $^{a}$( x$_{t}$ ), is given by the classical

Bayes’ theorem:

$$
p^a(\mathbf{x}_t) \propto p(y_t^o | \mathbf{x}_t) p^f(\mathbf{x}_t), \quad (11)
$$

where the symbol ∝ means proportional to and where p ( y o t | x$_{t}$ ) represents the data likelihood, i.e., the condi- tional PDF of having the observations y o t given the con- trol vector x$_{t}$ .

Based on Bayesian theory, the EnKF algorithm assumes that the errors on the control parameters x$_{t}$ and the errors on the observations y o t are random variables defined by Gaus- sian PDFs with a zero mean value and an error covariance model. Under these assumptions, the forecast PDF may be written as follows:

$$
p^f(\mathbf{x}_t) \propto \exp \left\{ -\frac{1}{2} (\mathbf{x}_t - \mathbf{x}_t^f)^T (\mathbf{P}_t^f)^{-1} (\mathbf{x}_t - \mathbf{x}_t^f) \right\}, \quad (12)
$$

where x f t is the forecast estimate of the control vector, and where P f t ∈ R n × n is the forecast error covariance matrix rep- resenting errors in the ROS model parameters. The data like- lihood may be similarly expressed as follows:

$$
p(y_t^o | \mathbf{x}_t) \propto \exp \left\{ -\frac{1}{2} \mathbf{d}_t^\mathrm{T} \mathbf{R}^{-1} \mathbf{d}_t \right\}, \quad (13)
$$

with R ∈ R 2 N$_{fr}$ × 2 N$_{fr}$ the observation error covariance matrix representing observation errors (assumed constant over time in this study), and with d$_{t}$ the innovation vector of size 2 N o fr corresponding to the differences between SFF and OFF:

$$
d_t = y_t^o - y_t^f = y_t^o - \mathcal{G}_t(x_t^f). \quad (14)
$$

Using the selection procedure (see Fig. 3), d$_{t}$ is simply de- fined as the vector formed by the directed distances between the paired SFF-OFF markers. Note that the statistical mo- ments of d$_{t}$ (e.g., mean and STD) provide a convenient mea- sure of the deviations of model predictions from observa- tions.

Within this framework, the analysis PDF from Eq. (11) is also Gaussian and is written as follows:

$$
p^a(\mathbf{x}_t) \propto \exp \left\{ -\frac{1}{2} (\mathbf{x}_t - \mathbf{x}_t^f)^T (\mathbf{P}_t^f)^{-1} (\mathbf{x}_t - \mathbf{x}_t^f) - \frac{1}{2} \mathbf{d}_t^T \mathbf{R}^{-1} \mathbf{d}_t \right\}, \quad (15)
$$

$$
\propto \exp \left\{ -\frac{1}{2} (\mathbf{x}_t - \mathbf{x}_t^a)^T (\mathbf{P}_t^a)^{-1} (\mathbf{x}_t - \mathbf{x}_t^a) \right\}, \quad (16)
$$

where x a t is the analysis estimate of the control vector, and where P a t ∈ R n × n is the analysis error covariance matrix. Conditional mode estimation searches for the mode of the PDF p $^{a}$( x$_{t}$ ), i.e., the value of the control vector x$_{t}$ that max- imizes the probability to estimate its true value x t $_{t}$. Under Gaussian assumption, this maximum likelihood estimation is
equivalent to a minimization problem:

$$ \max_{\mathbf{x}_t \in \mathbb{R}^n} \mathrm{p}^a(\mathbf{x}_t) \iff \min_{\mathbf{x}_t \in \mathbb{R}^n} \{-\ln[\mathrm{p}^a(\mathbf{x}_t)]\} = \min_{\mathbf{x}_t \in \mathbb{R}^n} \mathcal{J}(\mathbf{x}_t), \quad (17) $$

with $\mathcal{J}$ the cost function of the estimation problem defined as follows:

$$ \mathcal{J}(\mathbf{x}_t) = \frac{1}{2} (\mathbf{x}_t - \mathbf{x}_t^f)^T (\mathbf{P}_t^f)^{-1} (\mathbf{x}_t - \mathbf{x}_t^f) + \frac{1}{2} \mathbf{d}_t^T (\mathbf{R})^{-1} \mathbf{d}_t. \quad (18) $$

The direct minimization of $\mathcal{J}$ leads to the classical KF equations when the generalized observation operator $\mathcal{G}_t$ is linear (denoted by $\mathbf{G}_t$). In the present case, this implies that the fire spread model $\mathcal{M}_{[t-1,t]}$ is linear and that the parameter evolution model $\mathcal{F}_{[t-1,t]}$ is linear (denoted by $\mathbf{F}_{[t-1,t]}$). Using these assumptions, it could be shown that the forecast in the prediction step is obtained via the integration of the following equations:

$$ \mathbf{x}_t^f = \mathbf{F}_{[t-1,t]} \mathbf{x}_{t-1}^a, \quad \mathbf{P}_t^f = \mathbf{F}_{[t-1,t]} \mathbf{P}_{t-1}^f \mathbf{F}_{[t-1,t]}^T, \quad (19) $$

assuming there is no error in the formulation of the parameter evolution model. In this context, the analysis update in Eq. (16) leads to the following equations:

$$ \mathbf{x}_t^a = \mathbf{x}_t^f + \mathbf{K}_t (\mathbf{y}_t^o - \mathbf{G}_t \mathbf{x}_t^f), \quad (20) $$

$$ \mathbf{K}_t = \mathbf{P}_t^f \mathbf{G}_t^T (\mathbf{G}_t \mathbf{P}_t^f \mathbf{G}_t^T + \mathbf{R})^{-1}, \quad (21) $$

$$ \mathbf{P}_t^a = (\mathbf{I}_n - \mathbf{K}_t \mathbf{G}_t) \mathbf{P}_t^f, \quad (22) $$

where $\mathbf{K}_t$ is called the gain matrix. Starting from a prior value of the control parameters (i.e., the forecast $\mathbf{x}_t^f$) and using the observations $\mathbf{y}_t^o$ available at time $t$, the analysis estimate $\mathbf{x}_t^a$ is a feedback information for the fire spread model; $\mathbf{x}_t^a$ is optimal when the variance of its distance to the true value $\mathbf{x}_t^f$ gets to a minimum, meaning, for Gaussian cases, that its PDF is dense around its mean. The expressions in Eqs. (20–22) are the basis of the EKF algorithm used in Rochoux et al. (2013a, b), $\mathbf{F}_{[t-1,t]}$ and $\mathbf{G}_t$ being the tangent linear operators (Jacobian) of $\mathcal{F}_{[t-1,t]}$ and $\mathcal{M}_{[t-1,t]}$ in the vicinity of the control vector $\mathbf{x}_t$, respectively. Thus, in the EKF, a linearized and approximate equation is used for the prediction of errors statistics as well as for the relation between the control space and the observation space.

In contrast, the EnKF algorithm used in this study does not require the explicit use of the linear operators $\mathbf{F}_{[t-1,t]}$ and $\mathbf{G}_t$ in the prediction step. As shown in Fig. 4, the forecast control parameters $\mathbf{x}_t^f$ are stochastically represented at time $t$ based on $N_e$ realizations called the ensemble members

$$ [\mathbf{x}_t^{\mathrm{f},(1)}, \dots, \mathbf{x}_t^{\mathrm{f},(k)}, \dots, \mathbf{x}_t^{\mathrm{f},(N_e)}], $$

with $k$ varying between 1 and $N_e$. These realizations are randomly generated based on mean and error STD according to the user-defined confidence interval for each control parameter over the first assimilation cycle and to previous analysis

results for next assimilation cycles. The temporal evolution of the control parameters is artificially set up using a random walk model so that the $k$-th ensemble member reads:

$$ \mathbf{x}_t^{\mathrm{f},(k)} = \mathcal{F}_{[t-1,t]} (\mathbf{x}_{t-1}^{\mathrm{a},(k)}) = \overline{\mathbf{x}_{t-1}^{\mathrm{a}}} + \mathbf{e}_{t-1}^{(k)}, \quad (23) $$

where $\overline{\mathbf{x}_{t-1}^{\mathrm{a}}}$ is the mean of the posterior estimates obtained at the previous analysis time ($t-1$), and where $\mathbf{e}_{t-1}^{(k)}$ is a randomly generated white noise following a Gaussian distribution of zero mean and given STD (taken equal to the forecast error STD $\sigma^f$ in the following). Thus, the generation of the ensemble of forecast parameters at time $t$ is performed in two steps as in Peters et al. (2005): (1) the mean forecast estimate over the time window $[t, t+1]$ is specified using the mean of the analysis estimates obtained over the previous assimilation cycle $[t-1, t]$; and (2) the ensemble of forecast parameters is obtained by applying a STD to this mean forecast estimate. Additionally, the error STD used in the random walk model remains constant over all assimilation cycles (Peters et al., 2005; Ruiz et al., 2013b). A series of $N_e$ independent forward model integrations up to the analysis time $t$ based on these $N_e$ realizations of the control parameters is performed (starting from the same initial condition at time ($t-1$) that corresponds to the mean of the posterior estimates $\mathbf{x}_{t-1}^{\mathrm{a}}$); this forecast step provides $N_e$ fire front positions at time $t$ corresponding to the model counterparts of the observed quantities and designated as $[\mathbf{y}_t^{\mathrm{f},(1)}, \dots, \mathbf{y}_t^{\mathrm{f},(k)}, \dots, \mathbf{y}_t^{\mathrm{f},(N_e)}]$, with $\mathbf{y}_t^{\mathrm{f},(k)} = \mathcal{G}_t(\mathbf{x}_t^{\mathrm{f},(k)})$ for the $k$th ensemble member.

We note $\mathbf{C}_{xy} \in \mathbb{R}^{n \times 2N_{\mathrm{fr}}^o}$ the matrix that represents the stochastically based relation between the control space (of size $n$) and the observation space (of size $2N_{\mathrm{fr}}^o$); $\mathbf{C}_{xy}$ is expressed as follows:

$$ \mathbf{C}_{xy} = \mathbf{P}_t^f \mathbf{G}_t^T \quad (24) \\ = \sum_{k=1}^{N_e} \frac{(\mathbf{x}_t^{\mathrm{f},(k)} - \overline{\mathbf{x}_t^f}) (\mathcal{G}_t(\mathbf{x}_t^{\mathrm{f},(k)}) - \overline{\mathcal{G}_t(\mathbf{x}_t^f)})^T}{N_e - 1}, $$

where the overline denotes the mean value over the ensemble. Similarly, the symmetric error covariance matrix on the predicted measurements denoted by $\mathbf{C}_{yy} \in \mathbb{R}^{2N_{\mathrm{fr}}^o \times 2N_{\mathrm{fr}}^o}$ is stochastically formulated as follows:

$$ \mathbf{C}_{yy} = \mathbf{G}_t \mathbf{P}_t^f \mathbf{G}_t^T \quad (25) \\ = \sum_{k=1}^{N_e} \frac{(\mathcal{G}_t(\mathbf{x}_t^{\mathrm{f},(k)}) - \overline{\mathcal{G}_t(\mathbf{x}_t^f)}) (\mathcal{G}_t(\mathbf{x}_t^{\mathrm{f},(k)}) - \overline{\mathcal{G}_t(\mathbf{x}_t^f)})^T}{N_e - 1}. $$

This means that the EnKF algorithm approximates the mean and the covariance of the forecast by the mean and the covariance of an ensemble, while still making the assumption that all PDFs are Gaussian. Additionally, the distance between the $k$th prediction $\mathbf{y}_t^{\mathrm{f},(k)} = \mathcal{G}_t(\mathbf{x}_t^{\mathrm{f},(k)})$ and the observation vector $\mathbf{y}_t^o$ is computed according to Burgers et al. (1998), meaning
that an additional noise $\xi^{0,(k)}$ is added to the observation vector to avoid ensemble collapse. Thus, for the $k$th member, the innovation vector $d_t^{(k)}$ reads:

$$d_t^{(k)} = y_t^0 + \xi^{0,(k)} - y_t^{f,(k)}. \quad (26)$$

During the analysis, each ensemble member is updated based on the classical KF formulation presented in Eqs. (20)–(22), with the difference than the generalized observation operator $G_t$ is nonlinear and that the gain matrix $K_t^e$ is now stochastically calculated using Eqs. (24)–(25). The $k$th member analysis satisfies:

$$x_t^{a,(k)} = x_t^{f,(k)} + K_t^e (y_t^0 + \xi^{0,(k)} - G_t(x_t^{f,(k)})), \quad (27)$$

$$K_t^e = C_{xy}(C_{yy} + R)^{-1}. \quad (28)$$

One of the advantages of the EnKF formulation in Eqs. (27)–(28) is that the explicit estimation of the tangent-linear of the observation operator $G_t$ (including the tangent-linear of the fire spread model for parameter estimation) is avoided. This ensemble-based method allows the nonlinearity in the observation operator $G_t$ to be better taken into account than a local estimation $G_t$ achieved for instance through a finite difference scheme as in the EKF (Ros and Borga, 1997; Rochoux et al., 2013a, b). The use of Eqs. (27)–(28) provides an ensemble of posterior estimates at time $t$,

$$[x_t^{a,(1)}, \dots, x_t^{a,(k)}, \dots, x_t^{a,(N_e)}],$$

which is easily used to simulate over the time window $[t-1, t]$ an ensemble of retrospective posterior estimates of the fire front positions $[y_t^{a,(1)}, \dots, y_t^{a,(k)}, \dots, y_t^{a,(N_e)}]$ as well as an ensemble of forecasts of the fire spread beyond time $t$.

Note that in the present study, we assume that observation errors are uncorrelated, i.e., the observation error covariance matrix $R$ is treated as a diagonal matrix, in which each diagonal term is the error variance $(\sigma^0)^2$ associated with the error in the $x$- or $y$-coordinate of the markers along OFF.

As a summary, the main steps of the proposed EnKF algorithm for parameter estimation over the assimilation cycle $[t-1, t]$, are as follows:

1. build an ensemble of forecast control parameters based on Eq. (23), starting from the progress variable field corresponding to the mean analysis field obtained at time $(t-1)$;

2. compute the observation operator through Eq. (10), which includes the FIREFLY fire spread model integration from time $(t-1)$ to time $t$, in order to obtain the model counterparts of the observations at time $t$;

3. apply the Kalman filter update equation at time $t$ for each member of the ensemble based on Eqs. (24)–(28);

<figure><img src="image_7.png" /><figcaption>Figure 5. Flowchart of the PC-EnKF algorithm during the assimilation cycle $[t-1, t]$ decomposed into three steps: (1) construction of the PC expansion of the generalized observation operator; (2) EnKF prediction and update for the assimilation cycle $[t-1, t]$; and (3) parameter evolution to the next assimilation cycle $[t, t+1]$.</figcaption></figure>

4. re-integrate the model Eq. (8) with the analysis parameters over the time period $[t-1, t]$ to obtain the corrected locations of the fire front and the updated progress variable field at time $t$.

To move to the next assimilation cycle $[t, t+1]$, step (1) can be performed again. The integration of the FIREFLY fire spread model starts again from the location of the fire front associated with the mean analysis estimate at time $t$, using the modified control parameters following the random walk model (see Eq. 23).

The proposed parameter estimation algorithm can be regarded as a three-dimensional variational technique with a stochastically based estimation of the error covariance matrices; the three-dimensional variational technique also lacks the dynamic interrelation between analysis and forecast error covariances (Peters et al., 2005; Ruiz et al., 2013b).

## 3.2 Polynomial chaos-based ensemble Kalman filter

In the classical EnKF algorithm, a Monte Carlo sampling is used to generate the forecast members $x_t^f$ and their associated fire front trajectory $y_t^f$. While this provides accurate access to the full statistics of the modeling uncertainties (provided the ensemble $N_e$ is sufficiently large), it involves a large number of forward model integrations (as illustrated in Fig. 4) that becomes time-consuming for regional-scale fire spread problems. To maintain the computational cost of the EnKF algorithm compatible with the objective of wildfire spread forecasting, a numerical strategy based on a polynomial chaos (PC) expansion is introduced; this PC expansion is used in the prediction step of the EnKF algorithm as outlined in Fig. 5. This hybrid EnKF algorithm is denoted by PC-EnKF in the following.
## 3.2.1 General formulation of the surrogate model

The PC-based surrogate model approximates the generalized observation operator $G_t$ at time $t$ and is therefore denoted by $G_{pc,t}$. It is parameterized with respect to the multi-dimensional control vector $x_t^f \in \mathbb{R}^n$ following the forecast PDF $p^f(x_t)$. This random vector may be regarded as a set of second-order random variables (i.e., with finite variance) expressed in terms of a random event $\omega$ such that $x_t^f = x_t^f(\omega)$. It can be projected onto a stochastic space spanned by orthogonal PC functions of independent Gaussian random variables $\zeta(\omega)$ as follows:

$$x_t^f(\omega) = [x_{1,t}^f, x_{2,t}^f, \dots, x_{n,t}^f] = \sum_{q=0}^{\infty} \hat{x}_q \varphi_q(\zeta(\omega)). \quad (29)$$

The simulated positions of the fire front $y_t^f = G_t(x_t^f(\zeta))$ can also be viewed as a random variable and therefore, they can be projected onto a stochastic space spanned by orthogonal PC functions as follows:

$$y_t^f = G_{pc,t}(x_t^f(\zeta)) = \sum_{q=0}^{\infty} \hat{y}_q(t) \varphi_q(\zeta), \quad (30)$$

where $\hat{y}_q = \hat{y}_q(t)$ are time-dependent coefficients, and where $(\varphi_q)_{q=0,\dots,\infty}$ designate the multi-dimensional approximating polynomial functions forming an orthogonal basis with respect to the joint PDF $p^f(x_t) = p^f(x_{1,t}, x_{2,t}, \dots, x_{n,t})$. The choice for the basis functions may depend on the type of random variable functions (Xiu and Karniadakis, 2002). Since the control vector $x_t^f$ is assumed to follow a Gaussian PDF $p^f(x_t)$ within the framework of the EnKF, the surrogate model of the observation operator $G_{pc,t}$ is built upon the basis of the Hermite polynomials (Ghanem and Spanos, 1991). Stated differently, the Hermite polynomials form the optimal basis for random variables following multi-variate Gaussian PDF. Note that the model outputs $y_t^f$ are represented in terms of the same random event $\omega$ as the model inputs $x_t^f$, since the uncertainty in the model outputs is assumed to be mainly due to the uncertainty in the ROS model parameters within the framework of parameter estimation.

In practice, a truncated version of Eq. (30) is used; there are several ways of constructing the approximation space. The most common choice is to constrain the number of terms $N_{pc}$ in the PC expansion by the number of control parameters $n$ and by the maximum order of the polynomial basis $Q_{po}$ such that

$$N_{pc} = \frac{(n + Q_{po})!}{(n! Q_{po}!)}. \quad (31)$$

This choice of $N_{pc}$ ensures that the PC approximation is of highest order $Q_{po}$. Note that $Q_{po}$ is a user-defined quantity that must be chosen carefully according to the model nonlinearity, in order to obtain an accurate representation of the

model outputs $y_t^f$ with a high-order convergence rate. Theoretically, $Q_{po} = 1$ (i.e., only two terms for $n = 1$ corresponding to the mean and STD of the control variable) is enough to approximate exactly a Gaussian random variable. Note also that $N_{pc}$ rapidly grows with $n$ and $Q_{po}$, implying that a balance between accuracy and computational cost must be found. For instance, if $n = 2$ and $Q_{po} = 2$, there are $N_{pc} = 6$ terms retained in the PC expansion. Using this formalism, the surrogate model $G_{pc,t}$ can be formulated as follows:

$$y_t^f \approx G_{pc,t}(x_t^f(\zeta)) = \sum_{q=0}^{N_{pc}} \hat{y}_q(t) \varphi_q(\zeta), \quad (32)$$

where the unknowns are the following time-dependent vectors:

$$\hat{y}_q = \hat{y}_q(t) = [(\hat{x}_1, \hat{y}_1)_q, \dots, (\hat{x}_{N_{fr}^o}, \hat{y}_{N_{fr}^o})_q]_t, \quad (33)$$

$q$ varying between 1 and $N_{pc}$, with $N_{fr}^o$ the number of markers along OFF at time $t$. Note that the size of the $q$th vector $\hat{y}_q$ is $2N_{fr}^o$ (each marker location being represented with both the x- and y-coordinate on the horizontal plane) and thereby, the computation of $(2N_{fr}^o N_{pc})$ coefficients (also referred to as the PC modes) is necessary to build the surrogate model $G_{pc,t}$.

## 3.2.2 Calculation of the polynomial chaos modes

Due to the orthogonality of the PC basis, it can be shown that the $q$th PC coefficients $\hat{y}_q$ are given by the following:

$$\hat{y}_q = \frac{\mathbb{E}[G_{pc,t}(x_t^f) \varphi_q(\zeta)]}{\mathbb{E}[\varphi_q(\zeta)^2]}, \quad (34)$$

where

- $\mathbb{E}[\cdot]$ refers to the expectation operator satisfying $\mathbb{E}[\varphi_q(\zeta) \varphi_l(\zeta)] = 0$ if $q \neq l$, with the following definition for the inner product:

$$\mathbb{E}[\varphi_q(\zeta) \varphi_l(\zeta)] = \int_{\mathbb{R}^n} \varphi_q(\zeta) \varphi_l(\zeta) p(\zeta) d\zeta = \delta_{ql} [\varphi_q^2], \quad (35)$$

with $\delta_{ql}$ the Kronecker delta-function;

- $\mathbb{E}[\varphi_q(\zeta)^2]$ is a normalization factor equal to 1 if the basis is constructed orthonormal;

- $\mathbb{E}[G_{pc,t}(x_t^f) \varphi_q(x_t^f)]$ is computed using a Gauss-Hermite quadrature rule, with $[x_t^{f,(1)}, \dots, x_t^{f,(j)}, \dots, x_t^{f,((N_{quad})^n)}]$ the quadrature roots vector of size $(N_{quad})^n$ constrained by the maximum order of the polynomial basis $Q_{po}$.
such that $2Q_{po} \le 2(N_{quad} - 1)$. Thus, this term is computed as follows:

$$ \mathbb{E}[\mathcal{G}_{\text{pc},t}(\mathbf{x}_t^f) \varphi_q(\zeta)] = \int_{\mathbb{R}^n} \mathcal{G}_t(\mathbf{x}_t^f) \varphi_q(\zeta) \mathrm{d}\zeta \quad (36) \\ \approx \sum_{j=1}^{(N_{quad})^n} \mathcal{G}_t(\mathbf{x}_t^{f,(j)}) \varphi_q(\zeta^{(j)}) w^{(j)}, $$

where $\mathbf{y}_t^{f,(j)} = \mathcal{G}(\mathbf{x}_t^{f,(j)})$ corresponds to the FIREFLY forward model integration evaluated at the jth quadra-ture root $\mathbf{x}_t^{f,(j)}$ with its associated weight $w^{(j)}$, and where $\varphi_q$ is the qth multi-dimensional basis function formulated as tensor products of one-dimensional polynomial functions:

$$ \varphi_q \equiv \varphi_q(\zeta) = \prod_{l=1}^{n} \varphi_{i(l)}^{1D}(\zeta_l), \quad (37) $$

with $\varphi_{i(l)}^{1D}$ the one-dimensional polynomial basis and its multi-index $i(l)$ varying between 0 and $Q_{po}$ to determine the proper term in the multi-variable space.

Based on this formulation, the construction of the surrogate model $\mathcal{G}_{\text{pc},t}$ over the assimilation window $[t-1, t]$ requires a limited number of $(N_{quad})^n$ forward model integrations (see the first step in Fig. 5). The polynomial approximation $\mathcal{G}_{\text{pc},t}$ calculated in Eq. (32) is then used in the prediction step of the EnKF algorithm (instead of the observation operator $\mathcal{G}_t$) to compute the predictions of the time-evolving fire front locations $\mathbf{y}_{\text{pc},t}^f$ for a large number of members $N_e$ (see the second step in Fig. 5). This ensemble of forecasts is used to accurately estimate the covariance matrices $\mathbf{C}_{xy}$ and $\mathbf{C}_{yy}$ that are required in the formulation of the Kalman gain matrix. Thus, the EnKF update can be performed with reliable covariance matrices at a reduced computational cost compared to the standard EnKF algorithm based on a Monte Carlo sampling. This approach leads to analysis estimates of the control parameters $\mathbf{x}_t^{a,(k)}$ and to accurate PDF of the fire front locations $\mathbf{y}_t^{a,(k)}$ ($k=1, \dots, N_e$) using the same surrogate model as for the forecast estimates.

In order to reduce the computational cost of the EnKF algorithm, a surrogate model based on a PC expansion is used in place of the forward model (i.e., the FIREFLY regional-scale wildfire spread model) in the DA procedure. The performance of the resulting PC-EnKF algorithm is assessed on synthetically generated fire spread cases based on preliminary work presented in Rochoux et al. (2012) as well as on the controlled grassland fire experiment.

## 3.3 Numerical implementation

In practice, the EnKF and PC-EnKF ensemble-based DA algorithms were implemented with the fire spread simulator FIREFLY using the OpenPALM dynamic coupling software (Lagarde et al., 2001), co-developed at CERFACS and

ONERA.<fnref n="4" /> OpenPALM allows for the coupling of independent code components with a high-level of modularity in the data exchanges and treatment, while providing a straightforward parallelization environment (Fouilloux and Piacentini, 1999; Buis et al., 2006). In this study, it is used as a task-parallelism manager to handle communications and data exchanges between FIREFLY and the mathematical units required to sequentially apply the EnKF and PC-EnKF algorithms. The PALM-PARASOL functionality in OpenPALM was used to efficiently and independently run the FIREFLY time-integrations in parallel, on the available processors. The master processor of PALM-PARASOL spawns multiple copies of the same computer program (i.e., the slaves), each on one or several processors with a different set of input parameters of the ROS model. Each slave integrates FIREFLY using one realization of the control vector $\mathbf{x}_t$ to provide the associated fire front position $\mathbf{y}_t$, subsequently used for the computation of the covariance matrices $\mathbf{C}_{xy}$ and $\mathbf{C}_{yy}$. As illustrated in Figs. 4 and 5, this integration is performed for the $N_e$ ensemble members using the forward model FIREFLY for the classical EnKF. In contrast, for the PC-EnKF algorithm, a limited number of FIREFLY model integrations $(N_{quad})^n$ is used to build the surrogate model and subsequently, a large number of evaluations of the surrogate model ($N_e$) are computed using the PALM-PARASOL functionality.

# 4 Data assimilation experiments

## 4.1 Convergence of the ensemble-based algorithms

The EnKF and PC-EnKF algorithms are compared on an OSSE experiment, in which the Rothermel-based ROS model of Eq. (3) is reformulated as $\Gamma(x, y) = P \delta_v(x, y)$ with $P$ [$s^{-1}$] a proportionality coefficient and $\delta_v \equiv \delta_v(x, y)$ a spatially varying function that is assumed to be perfectly known. Note that this formulation takes advantage of the proportionality between the ROS $\Gamma$ and the fuel layer thickness $\delta_v$ in the Rothermel's formulation. Thus, the control vector is limited to a single parameter, $x = P$, which encompasses different uncertainties that are not distinguished here.

The fire is ignited at $(x_{ign}, y_{ign}) = (100 \text{ m}, 100 \text{ m})$ as a circular front with a radius of 5 m; it spreads upon a random fuel distribution $\delta_v(x, y)$ over a $200 \text{ m} \times 200 \text{ m}$ domain. Observations (represented using $N_{fr}^o = 20$ front markers) are synthetically generated at 50 s intervals with FIREFLY and a chosen true value $\mathbf{x}^t = P^t = 0.4 \text{ s}^{-1}$. An observation error characterized by the error STD $\sigma^o$ is also introduced. The ensemble of prior values is drawn from a Gaussian distribution centered in $\mathbf{x}^f = 0.2 \text{ s}^{-1}$ with an error STD $\sigma^f = 0.05 \text{ s}^{-1}$ (assumed constant along the assimilation cycles). Note that the true value of the control parameter $\mathbf{x}^t$ is at the tail of the Gaussian PDF associated with the forecast estimates. This case is chosen on

<footnote n="4">⁴http://www.cerfacs.fr/globc/PALM_WEB/</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Convergence of the mean analysis estimates of the pro-
portionality coefficient *P* [s⁻¹] with respect to the number of en-
semble members *N*<sub>e</sub> for a fixed observation error STD σ° = 2 m and
a single assimilation cycle: comparison of the performance between
the EnKF and PC-EnKF algorithms. The orange triangled-dashed
line corresponds to *Q*<sub>po</sub> = 2; and the red circled-dashed line corre-
sponds to *Q*<sub>po</sub> = 4 for the PC-EnKF algorithm. Black squares cor-
respond to the analysis estimates obtained using the standard EnKF.
Vertical error bars correspond to the associated error STD.</figcaption></figure>

purpose, in order to evaluate the ability of the parameter esti-
mation approaches (EnKF and PC-EnKF) to retrieve accurate
values of the control parameter, even though the prior value
is far from the true control parameter and its uncertainty is
high (compared to the observation error STD).

A PC approximation (with a polynomial order $Q_{po} = 4$ and subsequently a quadrature order $N_{quad} = 5$, see Sect. 3.2) is used to build the model response surface $\mathcal{G}_{pc}$ to the control parameter $\mathbf{x} = \mathbf{P}$ corresponding to the forecast ($\mathbf{x}^f, \sigma^f$).

**4.1.1 Sensitivity to sampling errors**

Convergence properties of the EnKF-based analysis esti-
mates are studied in Fig. 6 with respect to the number of
ensemble members *N*<sub>e</sub> for a fixed observation error STD
*σ*° = 2 m and for one assimilation cycle. Since there is no
analytical solution to the problem, the convergence of the
EnKF is assumed to be achieved if the mean value of the con-
trol parameter and its STD remain constant when increasing
*N*<sub>e</sub>. The performance of the PC-EnKF algorithm is compared
to that of the standard EnKF algorithm (black squares) for
different PC polynomial orders, *Q*<sub>po</sub> = 2 (orange triangled-
dashed line) and *Q*<sub>po</sub> = 4 (red circled-dashed line).

Figure 6 shows that in the present configuration, the EnKF
algorithm converges for a minimum of N_e = 48 members
(meaning that FIREFLY is integrated 48 times to produce

48 fire front trajectories associated with each realization of
the control parameter). In particular, below this threshold,
the error bars corresponding to the error STD of the analy-
sis parameter estimates are narrower for both EnKF and PC-
EnKF algorithms. The error STD computed with a low num-
ber of members is therefore not reliable and the ensemble-
based algorithms require a larger sample to accurately repre-
sent the tails of the Gaussian PDF related to the control pa-
rameter *P*. It is shown that the PC-EnKF algorithm provides
a comparable result as the EnKF (in terms of mean and STD)
above *N*<sub>e</sub> = 40 members for a polynomial order *Q*<sub>po</sub> = 4.
However, the results achieved with PC-EnKF are obtained
for a lower number of FIREFLY time-integrations (i.e., 5
FIREFLY model integrations only since *N*<sub>quad</sub> = 5 quadra-
ture points are used to build the model surface response *G*<sub>pc</sub>)
than the standard EnKF, while considering the same number
of members *N*<sub>e</sub> to generate the forecast/analysis estimates.
Thus, the PC-EnKF algorithm provides a solution that re-
produces the converged solution of the EnKF for a compu-
tational cost that is reduced by a factor of at least 8. This
implies that for more complex fire spread cases where more
members are required to track spatial variations in wind and
vegetation conditions, the PC-EnKF algorithm appears as
a promising alternative to obtain accurate simulations of fire
spread at a reasonable computational cost. Additionally, the
PC-EnKF algorithm provides a mean estimate that is less
fluctuating than the EnKF algorithm, with a slightly reduced
scatter for low values of *N*<sub>e</sub>, indicating that the PC-EnKF
strategy requires less ensemble members *N*<sub>e</sub> to reach conver-
gence.

Figure 6 also illustrates the sensitivity of the PC-EnKF-
based analysis to the choice of the PC polynomial order
$Q_{po}$ for a varying number of ensemble members $N_e$. While
$Q_{po} = 2$ (i.e., $N_{quad} = 3$) provides a reasonable approxima-
tion of the mean analysis estimate when considering the
standard EnKF as reference, $Q_{po} = 4$ (i.e., $N_{quad} = 5$) leads
to a more accurate estimate without loss of accuracy. Even
though the fire front marker locations exhibit approximate
Gaussian PDF and in theory $n = 1$ is sufficient to character-
ize their distributions, a high polynomial order is required
in this case. The true value ($P^t = 0.4 \text{ s}^{-1}$) is indeed not in
the zone of high probability occurrence of the forecast esti-
mates ($P^f = 0.2 \text{ s}^{-1}$ with $\sigma^f = 0.05 \text{ s}^{-1}$); the true fire front
locations are at the tail of the forecast PDF, which makes the
estimation of the fire front locations more difficult. This diffi-
culty shows the ability of the PC-EnKF procedure to retrieve
accurate estimates of the fire spread at a low computational
cost and without loss of accuracy, even though the prior in-
formation is very uncertain.

**4.1.2 Example of polynomial chaos-based surface
response**

Figure 7 provides a comparison in the observation space
between the observed fire front and the forecast/analysis
<figure><img src="image_9.png" /></figure>

**Figure 7.** Comparison of fire front locations using the PC-EnKF approach with an EnKF update at 50 s for $\sigma^0 = 2$ m, $Q_{po} = 4$ and $N_e = 1000$; estimation of the proportionality coefficient $P$ [$s^{-1}$]; all fronts correspond to time 50 s. Black crosses correspond to observations; the blue circled-dashed line corresponds to the mean forecast estimate of the fire front and the red squared-solid line corresponds to the analysis counterpart. Horizontal and vertical error bars correspond to the associated error STD along the x- and y-directions, respectively. The location of the fire front marker indexed by the character m is indicated.

estimates obtained through the PC-EnKF algorithm for an observation error STD $\sigma^0 = 2$ m, a PC polynomial order $Q_{po} = 4$ and a number of ensemble members $N_e = 1000$ (a single assimilation cycle is considered). As expected, the analysis estimates (red squared-solid line) provide a more accurate approximation of the observed fire front location (black crosses) than the forecast estimates (blue circled-dashed line).

To offer insight into the main ideas underlying the PC-EnKF algorithm, Fig. 8 illustrates the mapping between the control space and the observation space for one marker of the fireline (its position on the forecast/analysis fireline is indicated in Fig. 7 by the character m). The variations in the x and y coordinates of this marker are represented with respect to variations in the control parameter $P$: black crosses indicate the simulated marker positions associated with the $N_{quad} = 5$ quadrature roots (i.e., FIREFLY model integrations) corresponding to the first step of the PC-EnKF algorithm; and blue circles indicate the forecast estimates obtained through the surrogate model evaluation combined with a Monte Carlo sampling ($N_e = 1000$) corresponding to the second step of the PC-EnKF algorithm (see Fig. 5). These fire front estimates are associated with the forecast control parameter $P^f = 0.20 s^{-1}$ and its error STD $\sigma^f = 0.05 s^{-1}$. In contrast, red squares are produced by the EnKF update applied for any of the 1000 ensemble members, they correspond

<figure><img src="image_10.png" /><figcaption>(a) Forecast estimates</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) Analysis estimates.</figcaption></figure>

**Figure 8.** Model surface response (or surrogate model) of the x and y coordinates of the front marker indexed by m on the fireline (see Fig. 7), with respect to the control parameter $P$ [$s^{-1}$]. Black crosses correspond to quadrature roots (i.e., FIREFLY forward model integrations); blue circles correspond to (a) forecast estimates, and red squares to (b) analysis estimates obtained through the PC-EnKF algorithm at time 50 s. The vertical solid line indicates the true value $P^t = 0.4 s^{-1}$; the vertical dashed lines indicate the mean forecast and analysis estimates of the proportionality coefficient $P$ [$s^{-1}$].

to the analysis estimates related to $P^a = 0.38 s^{-1}$ and $\sigma^a = 0.01 s^{-1}$. The scatter of the ensemble is significantly reduced in the analysis, around the true value $P^t = 0.40 s^{-1}$, highlighting the uncertainty reduction achieved through the ensemble-based DA.

### 4.1.3 Sensitivity to observation errors

For verification purposes on the behavior of the PC-EnKF algorithm, Fig. 9 examines the influence of the observation error on the performance of the EnKF and PC-EnKF algorithms (the EnKF algorithm is used as reference). Statistics (in terms of mean value and STD) of the analysis obtained for $N_e = 48$ members over one assimilation cycle at time $t = 50$ s are presented as a function of the magnitude of the observation errors measured by $\sigma^o$ (up to $\sigma^o = 30$ m); vertical bars give a graphical representation of the magnitude of the STD within the analysis ensemble. The results show the consistency of the PC-EnKF algorithm with the EnKF in retrieving realistic values for the control parameter, even
<figure><img src="image_12.png" /></figure>

**Figure 9.** Mean and STD of the analysis estimates of the proportionality coefficient *P* [s⁻¹] as a function of the observation error STD σ° for a fixed number of members *N*ₑ = 48 and for one assimilation cycle (with an EnKF update at 50s): comparison between EnKF and PC-EnKF. The black solid line corresponds to the true value 0.4s⁻¹; the blue dashed line corresponds to the mean forecast estimate 0.2s⁻¹; and the red circled-dashed line corresponds to the mean analysis estimate obtained using the PC-EnKF algorithm. Black squares correspond to the mean analysis estimates obtained by the standard EnKF. Vertical error bars correspond to the associated error STD.

though the observation error is significant. When the observation error STD σ° is small, the PC-EnKF algorithm successfully drives the analysis ensemble towards the true value of the parameter *P*¹ = 0.4s⁻¹; the resulting analysis exhibits a much reduced scatter by at least a factor 4 in comparison to the forecast STD σ¹ = 0.05s⁻¹. In contrast, when σ° is large, the PC-EnKF algorithm has reduced effects and the analysis ensemble remains close to the forecast ensemble (the analysis STD is similar to the forecast STD σ¹ = 0.05s⁻¹). For intermediate values of σ°, the PC-EnKF algorithm produces optimized analyses lying between forecast and observation; as expected, the more accurate the observations, the more certain the analysis for a given forecast error.

## 4.2 Temporal variability of the parameter error

Sequential application of the EnKF allows for a temporal correction of the parameter *P* for a case in which the time-varying profile of the true parameter *P*¹ was artificially set up between 0.3 and 0.6s⁻¹ over seven assimilation cycles (the true profile is shown in Fig. 10a in black solid line). While the mean value of the forecast estimates is set to 0.2s⁻¹ for the first assimilation cycle, it is set to the mean analysis estimate from the previous assimilation cycle for all further as-

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

**Figure 10.** Sequential EnKF estimation of the coefficient *P* [s⁻¹] over seven assimilation cycles with *N*ₑ = 48 members and σ° = 5m; time-varying true control parameter. The green triangled-dashed-dotted curve corresponds to the free run (without DA); the blue circled-dashed curve corresponds to the mean forecast estimate; the red squared-solid curve corresponds to the mean analysis estimate; and the black solid line corresponds to the true control parameter. (a) Parameter estimates (vertical error bars correspond to the associated error STD). (b) Mean distance to the observed fire front.

simulation cycles. Since there is no explicit dynamic model for the control parameter *P*, it is of primary importance to track the temporal variability of the error in the parameter through a sequential estimation combined with a random walk model in the EnKF prediction step (see Eq. 23). Note that for this experiment, *N*ₑ = 48 members are considered in
the ensemble and a constant observation error STD $\sigma^0 = 5$ m is assumed.

Figure 10a shows the temporal variations of the EnKF estimates along the assimilation cycles. The EnKF solution (red squared-solid curve) provides an optimal mean value of the control parameter, resulting in an ensemble of fire fronts that is coherent with the observation error statistics (see Fig. 10b). In contrast, the model without DA (green triangled-dashed-dotted curve) significantly underestimates the ROS. While being not as accurate as the analysis at the assimilation time, the forecast (blue circled-dashed curve) provides a significant improvement in the prediction of wildfire spread at future lead times compared to the model without DA (i.e., free run). Note that there is a temporal shift between the forecast and analysis estimates in Fig. 10a. The analysis estimates are obtained at the current observation time and thereby, provide the most recently updated information. In contrast, the forecast estimates only contain information up to the previous analysis time:

- if the error in the control parameter does not change, the correction obtained at the previous analysis time is still valid and adapted to track the actual fire front position: the forecast estimates provide reliable scenarios of wildfire spread at future lead times;

- if this error significantly varies in-between two successive analysis times (this is the case of the present DA experiment), the correction is no longer suitable to predict wildfire spread at long lead times. For instance, over the assimilation cycle indexed by 5, the analysis estimates provide a good approximation of the actual fire front location at time $t_4$; however, the forecast obtained at time $t_5$ when starting from the analysis at $t_4$ overestimates the true ROS, meaning that a new observation is required to gain information on the wildfire behavior.

In the present case, Fig. 10b shows that the error in the forecast is systematically higher than that of the error in the analysis in the observation space. This means that the assimilation needs to be renewed according to the temporal variability of the error in the control parameter to ensure a high-level performance of the data-driven simulated forecast.

## 4.3 Application to a controlled grassland fire

The EnKF and PC-EnKF algorithms are applied to a real-world case study, corresponding to a reduced-scale controlled grassland fire (4 m × 4 m), propagating over a flat terrain and occurring under moderate wind conditions (Paugam et al., 2013). These wind conditions are assumed to be uniform and constant, $u_w^* = 1.0 \text{ m s}^{-1}$, blowing into a western direction, $\alpha_w^* = 307^\circ$. The grass is assumed to exhibit a uniform layer thickness, $\delta_v = 0.08 \text{ m}$, and controlled properties, for instance the moisture content and the grass surface-to-volume ratio are (approximately) known, $M_v = 22\%$ (field

<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

**Figure 11.** Model surface response (or surrogate model) of the x and y coordinates of the front marker indexed by m on the fire-line (see Fig. 12) with respect to the control vector $\mathbf{x} = [M_v, \Sigma_v]$ ($n=2$). Black crosses correspond to quadrature roots (FIREFLY integrations). (a) Forecast estimates (blue circles) and (b) analysis estimates (red squares) of the x (top panel) and y coordinates (bottom panel) of the fire front positions are mapped onto the PC-based model surface response.

measure) and $\Sigma_v = 11485 \text{ m}^{-1}$ (values related to short grass taken from the Rothermel's fuel database). The observed fire front locations are extracted from thermal-infrared imaging at 28 s intervals; they are discretized with $N_{\text{fr}}^0 = 40$ front markers with a measurement error estimated to $\sigma^0 = 0.05 \text{ m}$ (based on the spatial resolution of the camera).

The performance of both algorithms is studied over one assimilation cycle [50, 78 s] with an EnKF update at time $t_1 = 78 \text{ s}$ (analysis mode) as well as over a forecast time period [78, 106 s] with an EnKF forecast at time $t_2 = 106 \text{ s}$ (forecast mode); the initial condition at $t_0 = 50 \text{ s}$ is taken as the observed fire front at $t_0$. The control parameters are the fuel moisture $M_v$ and the fuel particle surface-to-volume ratio $\Sigma_v$ such that $\mathbf{x} = [M_v, \Sigma_v]$ with $n=2$. A PC approximation (with a polynomial order $Q_{\text{po}} = 4$ and a quadrature order $N_{\text{quad}} = 5$) is used to build the model response surface to the two control parameters $M_v$ and $\Sigma_v$; this response surface is shown for one particular simulated front marker in Fig. 11a (black crosses). A forecast ensemble of $N_e = 1000$ members (blue circles) is generated at no cost using the PC-based surrogate model assuming uncertainties in $M_v$ and $\Sigma_v$; the forecast estimates of these control parameters are described in Table 2 along with the associated STD. Note that the blue circles are contained within the surface response described by
<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 12.** Comparison between simulated and measured fire front positions for the controlled grassland fire experiment at (a) analysis time and (b) forecast time: black crosses correspond to observations, the blue circled-dashed line corresponds to the mean forecast estimate constructed through the PC-based surrogate model; the red squared-solid line corresponds to the mean analysis estimate obtained by the PC-EnKF procedure applied at time $t_1 = 78$ s. Black squares correspond to the standard EnKF used as reference. The location of the fire front marker indexed by m is indicated.</figcaption></figure>

the black crosses that represent the $(N_{quad})^2 = 25$ FIREFLY forward model integrations performed to build the PC-based surface response. Thus, the PC decomposition properly approximates the observation operator.

### 4.3.1 Analysis mode

In the analysis mode, the forecast ensemble is corrected by assimilating the fire front at time $t_1 = 78$ s. The comparison between the observations (black crosses), the forecast estimates (blue circled-dashed line) and the PC-EnKF-based analysis estimates (red squared-solid line) obtained at time $t_1 = 78$ s are presented in Fig. 12a. The forecast trajectory represents the ensemble mean of the surrogate model simulations obtained without DA (i.e., using the forecast estimates of the control parameters), while the analysis trajectory derives from an EnKF update at $t_1$ using the analysis estimates in the surrogate model integrations.

<table><caption>**Table 2.** PC-EnKF-based experiment for the controlled grassland fire experiment: error statistics (in the parameter space) of the forecast and analysis ensemble estimates for $x = [M_v, \Sigma_v] (n = 2)$. The number of FIREFLY integrations is also presented as indicator of the computational cost.</caption><thead><tr><th>FORECAST</th><th>Cost</th><th>Ens. mean</th><th>Ens. STD</th></tr></thead><tbody><tr><td>PC-EnKF</td><td>25</td><td>15.0%<br>11500 m<sup>-1</sup></td><td>4.0%<br>3000 m<sup>-1</sup></td></tr><tr><td>EnKF</td><td>1000</td><td>15.0%<br>11500 m<sup>-1</sup></td><td>4.0%<br>3000 m<sup>-1</sup></td></tr></tbody><thead><tr><th>ANALYSIS</th><th>Cost</th><th>Ens. mean</th><th>Ens. STD</th></tr></thead><tbody><tr><td>PC-EnKF</td><td>25</td><td>13.8%<br>22583 m<sup>-1</sup></td><td>1.4%<br>1157 m<sup>-1</sup></td></tr><tr><td>EnKF</td><td>1000</td><td>13.5%<br>22345 m<sup>-1</sup></td><td>1.4%<br>1170 m<sup>-1</sup></td></tr></tbody></table>

It is found that the PC-EnKF strategy allows to significantly decrease the distance between the observations and the simulated fronts with a comparable level of accuracy as the standard EnKF algorithm (the PC-EnKF algorithm provides similar analysis mean and STD, see Table 2). As illustrated in Fig. 11b, the uncertainty in the fire front positions is significantly reduced in comparison to the forecast since the STD related to the analysis estimates is much smaller than that of the forecast estimates. This indicates that the PC-EnKF algorithm allows reliable statistical information to be retrieved for only 25 FIREFLY model integrations (in contrast, the standard EnKF algorithm requires 1000 members to correct $n = 2$ control parameters and thereby, correct $2N_{fr}^o = 80$ fire front marker coordinates).

Consistently, Fig. 13 shows that the support of the analysis PDF (see Fig. 13b) is significantly reduced compared to the forecast PDF (see Fig. 13a) for the x and y coordinates of the $N_{fr}^o = 40$ observed front markers. The topology of the PDF along the observed fire front is found to be overall preserved through the EnKF update, implying that the assumption of Gaussian error statistics for the modeling error statistics seems not to deteriorate the performance of the ensemble-based DA algorithms. Some regions of the PDF related to the x coordinates of the front marker locations (nearby $x = 2$ m) are not sensitive to variations in $M_v$ and $\Sigma_v$. These regions correspond to the flank of the fire, meaning that the x coordinates of the surrounding front markers do not vary and the growth of the burning area only induces variations in the y coordinates.

As discussed for the OSSE test cases, the nonlinear response of the observation operator to the control parameters induces a slightly non-Gaussian PDF for the forecast estimates: it is indeed found that the mode of the PDF does not exactly coincide with the mean value. Note that the PDF exhibits a relatively flat tail for decreasing x and increasing y
<figure><img src="image_19.png" /><figcaption>**Figure 13.** Colormap of the PDF of the fire front marker locations (in terms of x and y coordinates) for the controlled grassland fire experiment at the analysis time $t_1 = 78$s. (a) Forecast PDF with respect to the x (left panel) and y (right panel) coordinates of the observed fire front markers. (b) Analysis PDF with respect to the x (left panel) and y (right panel) coordinates of the observed fire front markers.</figcaption></figure>

coordinates of the observed fire front markers: this is due to a sharp ROS acceleration when decreasing the fuel moisture content $M_v$ or alternatively, when increasing the fuel particle surface-to-volume ratio $\Sigma_v$ (see Fig. 2b and c).

### 4.3.2 Forecast mode

In the forecast mode, Fig. 12b compares the fire front position at $t_2 = 106$s obtained using the forecast estimates (without DA) and the analysis estimates derived from a DA update at $t_1 = 78$s. The PC-EnKF algorithm appears to properly represent the forecast trajectory at $t_2 = 106$s in comparison to the standard EnKF. This result illustrates that a PDF sampling based on PC (instead of Monte Carlo in the standard EnKF) can significantly reduce the computational cost of the EnKF prediction/update steps (in terms of number of FIREFLY model integrations that constitute the most time-consuming task in the PC-EnKF algorithm) and thereby, provide accurate error statistics on the inputs and outputs of the wildfire spread model. For instance, 1000 FIREFLY model integrations were used in the EnKF algorithm to accurately represent the error statistics; in contrast, only 25 FIREFLY model integrations were performed in the PC-EnKF algorithm. Thus, the number of FIREFLY model integrations is here divided by a factor of 40.

Additionally, Fig. 12b shows that the errors in the control parameters do not significantly change in-between the two

observation times (i.e., at $t_1 = 78$s and $t_2 = 106$s), meaning that an observation time period of 28 s seems appropriate for applying DA (relatively to the temporal variability of the errors in the control vector $\mathbf{x}$).

While the improved accuracy of EnKF-based data-driven simulations is obtained at the expense of heavy computational cost (in the context of multi-parameter estimation for instance), the PC-EnKF strategy appears as a promising strategy for solving Bayesian filtering problems at a low computational cost that is a requirement of operational frameworks.

## 5 Conclusions

A data assimilation (DA) strategy based on the ensemble Kalman filter (EnKF) with parameter estimation is demonstrated to account for both experimental and modeling uncertainties in wildfire spread modeling and thereby, to provide optimized forecast of wildfire behavior.

The proposed filtering strategy relies on a stand-alone sequential parameter estimation approach (the model state is not included in the control space), in which the control parameters are assumed spatially uniform and constant over the time window over which optimal values are sought for, and in which the observation operator is dynamically evolving to track the actual location of the fire front over time. This strategy was found efficient at reducing uncertainty in the numerical predictions of fire spread for synthetic measurements as well as for a (reduced-scale) controlled grassland fire experiment. It was also found that the nonlinear interrelation between the environmental parameters and the fire front positions induced by the nonlinearities of wildfire spread can be stochastically described over the ensemble members. Since there is no suitable dynamical model for the evolution of the control parameters, a random walk model based on the idea of persistence forecasting and on the assumption of constant error standard deviation is used to relate the forecast to the analysis parameter estimates. Thus, the parameter estimation approach can be regarded as a three-dimensional variational technique with stochastically based estimation of the error covariance matrices. In this context, it was highlighted that the duration of the assimilation cycle is of primary importance in the success of the proposed DA approaches: the assimilation must be renewed according to the temporal variability of the parameter errors, in order to track the actual fire behavior.

In order to reduce computational cost and balance sampling errors due to frequent assimilation cycles, a probabilistic sampling based on polynomial chaos techniques (PC-EnKF) was shown to significantly reduce the computational cost of the EnKF-based parameter estimation approach (by a factor of at least 10 in the present configurations), and thereby to provide access to accurate error statistics on both model inputs and outputs for the formulation of the Kalman gain matrix.
In both EnKF and PC-EnKF DA approaches, the formula-
tion of the observation operator can be regarded as the com-
putation of a front-to-front distance, in which the location
and topology of each front is evolving over time. The under-
lying selection operator that pairs the simulated front markers
with the observed front markers differs between the members
of the ensemble. This is an unusual application of the EnKF
algorithm that has demonstrated very good results to track
coherent features but that still needs more extensive verifi-
cation to ensure that the optimality of the filter is preserved.
Indeed, the computation of the distance between simulated
and observed markers becomes more challenging for com-
plex fire front topology. The generalization of this treatment
in the EnKF algorithm to complex fire front topology will be
revisited in future work and needs to be extended to the oper-
ational context, where wildfires propagate in a highly hetero-
geneous environment and where in-situ, airborne or space-
borne monitoring is currently limited. Projection schemes re-
ported in Rochoux (2014) are expected to provide a valuable
answer to this issue and could be integrated to the proposed
DA algorithms.

Future plans include the extension of the proposed PC-
EnKF approach to cases with spatially varying biomass fuel
conditions and time-varying wind conditions to further val-
idate and optimize this strategy in terms of computational
cost, in order to meet operational requirements. Note that
there is the possibility to change the parameters that are in-
cluded in the control vector between the assimilation cycles,
according to the level of information available (uncertainties
in the surface wind and in the biomass fuel properties are
not time-invariant); automatic sensitivity tests could be per-
formed inline to focus the EnKF update on the most uncertain
parameters over a given time period.

Future plans also include addressing the correlations of
observational errors along the fireline at a given time, simi-
larly to spaceborne data along the pass of polar-orbiting satel-
lites (Brankart et al., 2009; Gorin and Tsyrulnikov, 2011),
in order to improve the performance of both EnKF and PC-
EnKF DA algorithms. In the present study, the assumption
of uncorrelated observation errors is believed not to signifi-
cantly deteriorate the DA results since the procedure of front
selection with a low number of markers (N<sub>fr</sub><sup>o</sup> ≪ N<sub>fr</sub>) can be
regarded as a filtering procedure that tends to reduce the spa-
tial correlations between the observed front markers.

Even though both EnKF and PC-EnKF DA algorithms
used in this study are effective in correcting the input pa-
rameters of the rate of spread (ROS) model, the treatment of
highly anisotropic uncertainties remains an important chal-
lenge for wildfire spread forecasting. Due to the operational
constraints on the computational cost, the parameter estima-
tion approach could be extended to the case of weak spatial
variations of the ROS model parameters. Assuming that the
errors on the parameters vary slowly in time, the DA cor-
rection could then reasonably be used for forecast at future
lead times, thus allowing for mid- to long-term forecast. Still,

there is a need to design a strategy to address spatially dis-
tributed error correlations along the fireline in order to be
able to correct the shape of the time-evolving fire front. There
is also a need to address all possible sources of uncertainty
in the fire spread model, not only in the input parameters of
the ROS model but also in the parameterization of the ROS
that is limited in scope due to a lack of physical modeling
(e.g., steady-state assumption, spotting). The second part of
this series of two articles (Rochoux et al., 2014) is dedicated
to the evaluation of a state estimation approach that is able
to account for both anisotropic uncertainties and modeling
uncertainties. While out of the scope of this series of two ar-
ticles, a proper representation of the model errors could be
performed in the EnKF by introducing a model error covari-
ance matrix (Trémolet, 2007), which could be modeled using
a stochastic model such that proposed by Pagnini and Men-
trelli (2014) for the transport of firebrands.

*Acknowledgements.* The financial support provided by the Agence Nationale de la Recherche under the IDEA<fnref n="5" /> project grant ANR-09-COSI-006 (2010-2013), by the LEFE-MANU grant (INSU-CNRS program, 2011-2013) and by the 2012 CTR Summer Program<fnref n="6" /> (Uncertainty Quantification group) is greatly appreciated. The authors acknowledge Ronan Paugam and Martin Wooster (King's College London) for sharing the experimental data, Florent Duchaine and Thierry Morel (CERFACS) for support on Open-PALM, Albert Simeoni (University of Edinburgh) for clarifications on Rothermel's model as well as Gianluca Iaccarino and Paul Constantine (CTR) for helpful discussions on uncertainty quantification. They also acknowledge Marc Bocquet (CEREA), Eric Blayo (INRIA), and Serge Gratton (CERFACS) for helpful discussions on the properties of the ensemble Kalman filter for parameter estimation and dual parameter-state estimation.

Edited by: J.-B. Filippi

Reviewed by: G. Pagnini and two anonymous referees

References

Beezley, J. D., and Mandel, J.: Morphing Ensemble Kalman filters, *Tellus A*, 60, 131–140, doi:10.1111/j.1600-0870.2007.00275.x, 2008.

Birolleau, A., Poëtte, G., and Lucor, D.: Adaptive Bayesian in-
ference for discontinuous inverse problems, application to hy-
perbolic conservation laws, Commun. Comput. Phys., 16, 1–34,
2014.

Blanchard, E. D., Sandu, A., and Sandu, C.: A polynomial chaos-
based Kalman filter approach for parameter estimation of me-
chanical systems, J. Dyn. Sys., Meas., Control, 132, 061404,
doi:10.1115/1.4002481, 2010.

Boé, J., Terray, L., Martin, E., and Habets, F.: Projected changes
in components of the hydrological cycle in French river basins
during the 21st century, Water Resour. Res., 45, W08426,
doi:10.1029/2008WR007437, 2009.

<footnote n="5">⁵http://anridea.univ-corse.fr</footnote>

<footnote n="6">⁶http://www.stanford.edu/group/ctr/SummerProgram/</footnote>
Bouttier, F. and Courtier, P.: Data Assimilation Concepts and Methods, ECMWF, Meteorological Training Course Lecture Series, March 1999.

Boyaval, S.: A fast Monte-Carlo method with a reduced basis of control variates applied to uncertainty propagation and bayesian estimation, Comput. Methods Appl. Mech. Engrg., 241–244, 190–205, 2012.

Brankart, J.-M., Ubelmann, C., Testut, C.-E., Cosme, E., Brasseur, P., and Verron, J.: Efficient parameterization of the observation error covariance matrix for square root or ensemble Kalman filters. Application to ocean altimetry, Mon. Weather Rev., 138, 1908–1927, 2009.

Buis, S., Piacentini, A., and Declat, D.: PALM: a computational framework for assembling high performance computing applications, Concurrency Computat. Pract. Exper., 18, 247–262, 2006.

Burgers, G., van Leeuwen, P., and Evensen, G: Analysis scheme in the ensemble Kalman filter, Mon. Weather Rev., 126, 1719–1724, 1998.

Butler, B., Cohen, J., Latham, D., Schuette, R., Sopko, P., Shanon, K., Jimenez, D., and Bradshaw, L.: Measurements of radiant emissive power and temperatures in crown fires, Can. J. Forest Res., 34, 1577–1587, 2004.

Chong, D., Tolhurst, K. G., Duff, T. J., and Cirulis, B.: Sensitivity Analysis of PHOENIX RapidFire, Bushfire CRC, University of Melbourne, 2013.

Cowlard, A., Jahn, W., Abecassis-Empis, C., Rein, G., and Torero, J.: Sensor-assisted fire fighting, Fire Technol., 46, 719–741, 2010.

Crombette, P.: Optimisation et poursuite des développements du système LIVEFIRE de géolocalisation automatisée et temps réel de prises de vue aéroportées, Master thesis, Université de Toulouse (France), 2010.

Cruz, M. G. and Alexander, M. E.: Uncertainty with model predictions of surface and crown fire rates of spread, Environ. Modell. Softw., 47, 16–28, 2013.

Daley, R.: Atmospheric data analysis, Cambridge atmospheric and space science series, Cambridge University Press, 1991.

D'Andrea, M., Fiorucci, P., and Holmes, T. P.: A stochastic Forest Fire Model for future land cover scenarios assessment, Nat. Hazards Earth Syst. Sci., 10, 2161–2167, doi:10.5194/nhess-10-2161-2010, 2010.

Després, B., Poëtte, G., and Lucor, D.: Review of robust uncertainty propagation in systems of conservation laws with the entropy closure method, in: Uncertainty Quantification in Computational Fluid Dynamics, edited by: Bijl, H., Lucor, D., Mishra, S., and Schwab, C., Vol. 92 of Lecture Notes in Computational Science and Engineering, Springer-Verlag, 105–149, 2013.

Durand, M., Andreadis, K. M., Alsdorf, D. E., Lettenmaier, D. P., Moller, D., and Wilson, M. D.: Estimation of bathymetric depth and slope from data assimilation of swath altimetry into a hydrodynamic model, Geophys. Res. Lett., 35, L20401, doi:10.1029/2008GL034150, 2008.

Evensen, G.: Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics, J. Geophys. Res., 99, 143–162, 1994.

Evensen, G.: Data Assimilation – The Ensemble Kalman Filter, Springer, 2009.

Filippi, J.-B., Bosseur, F., Mari, C., Lac, C., Le Moigne, P., Cuenot, B., Veynante, D., Cariolle, D., and Balbi, J.-H.: Coupled

atmosphere-wildland fire modelling, J. Adv. Model. Earth Syst., 1, 210–226, 2009.

Filippi, J.-B., Pialat, X., and Clements, C. B.: Assessment of ForeFire/Meso-NH for wildland fire/atmosphere coupled simulation of the FireFlux experiment, Proc. Combust. Inst., 34, 2633–2640, 2013.

Finney, M. A.: FARSITE: Fire Area Simulator – model development and evaluation, Forest Service, US Dept. of Agriculture, Research Paper RMRS-RP-4, 1998.

Finney, M. A., Grenfell, I. C., McHugh, C. W., Seli, R. C., Trethewey, D., Stratton, R. D., and Brittain, S.: A method for ensemble wildland fire simulation, Environ. Model. Assess., 16, 153–167, 2011.

Finney, M. A., Cohen, J. D., McAllister, S. S., and Jolly, W. M.: On the need for a theory of wildland fire spread, Int. J. Wildland Fire, 22, 25–36, 2013.

Fouilloux, A. and Piacentini, A.: The PALM Project: MPMD paradigm for an oceanic data assimilation software, Lect. Notes Comput. Sc., 1685, 1423–1430, 1999.

Gelb, A.: Applied Optimal Estimation, Cambridge Massachusetts MIT Press, 1974.

Ghanem, R. and Spanos, P.: Stochastic Finite Elements, A Spectral Approach, Springer Verlag, 1991.

Gorin, V. and Tsyrulnikov, M.: Estimation of multivariate observation error statistics for AMSU-A data, Mon. Weather Rev., 139, 3765–3780, 2011.

Hirsch, K. G.: Canadian Forest Fire Behavior Prediction (FBP) System: User's guide, Northern Forest Centre, Special Report No. 7, 1996.

Ide, K., Courtier, P., Ghil, M., and Lorenc, A. C.: Unified notation for data assimilation: operational, sequential and variational, J. Meteorol. Soc. Jpn., 75, 181–189, 1997.

Jimenez, E., Hussaini, M. Y., and Goodrick, S. L.: Uncertainty quantification in Rothermel's model using an efficient sampling method, in: The fire environment–innovations, management, and policy; conference proceedings, edited by: Butler, B., W., Cook, W., comps, 26–30 March 2007, Destin, FL, 2007.

Kalnay, E.: Atmospheric Modeling, Data Assimilation and Predictability, Cambridge University Press, 2003.

Lagarde, T., Piacentini, A., and Thual, O.: A new representation of data assimilation methods: the PALM flow charting approach, Q. J. Roy. Meteor. Soc., 127, 189–207, 2001.

Lautenberger, C.: Wildland fire modeling with an Eulerian level-set method and automated calibration, Fire Safety J., 62, 289–298, 2013.

Le Maître, O. P. and Knio, O. M.: Spectral Methods for Uncertainty Quantification with Applications to Computational Fluid Dynamics, Scientific Computation, Springer, 2010.

Li, J. and Xiu, D.: On numerical properties of the ensemble Kalman filter for data assimilation, Comput. Meth. Appl. Math. Eng., 197, 3574–3583, 2008.

Li, J. and Xiu, D.: A generalized polynomial chaos based ensemble Kalman filter with high accuracy, J. Comput. Phys., 228, 5454–5469, 2009.

Linn, R., Reisner, J., Colman, J. J., and Winterkamp, J.: Studying wildfire behavior using FIRETEC, Int. J. Wildland Fire, 11, 233–246, 2002.
Lucor, D., Meyers, J., and Sagaut, P.: Sensitivity analysis of LES to subgrid-scale-model parametric uncertainty using Polynomial Chaos, J. Fluid Mech., 585, 255–279, 2007.

Mallet, V., Keyes, D. E., and Fendell, F. E.: Modeling wildland fire propagation with level-set methods, Comput. Math. Applic., 57, 1089–1101, 2009.

Mandel, J., Bennethum, L. S., Beezley, J. D., Coen, J. L., Douglas, C. C., Minjeong, K., and Vodacek, A.: A wildland fire model with data assimilation, Math. Comput. Simulat., 79, 584–606, 2008.

Mandel, J., Beezley, J. D., and Kochanski, A. K.: Coupled atmosphere-wildland fire modeling with WRF 3.3 and SFIRE 2011, Geosci. Model Dev., 4, 591–610, doi:10.5194/gmd-4-591-2011, 2011.

Mell, W., Jenkins, M. A., Gould, J., and Cheney, P.: A physics-based approach to modeling grassland fires, Int. J. Wildland Fire, 16, 1–22, 2007.

Merlet, N.: Evaluation des possibilités de géolocalisation automatisée et temps réel de prises de vue aéroportées dans le cadre de la lutte contre les feux de forêt, Master thesis, Université de Toulouse (France), 2008.

Milly, P. C., Wetherald, R. T., Dunne, K. A., and Delwort, T. L.: Increasing risk of great floods in a changing climate, Nature, 415, 514–517, 2002.

Moradkhani, H., Sorooshian, S, Gupta, H., and Houser, P.: Dual state-parameter estimation of hydrological models using ensemble Kalman filter, Adv. Water Resour., 28, 135–147, 2005.

Nijhuis, M.: Forest fires: burn out, Nature, 489, 352–354, 2012.

Noble, I. R., Bary, G. A. V., and Gill, A. M.: McArthur's fire danger meters expressed as equations, Aust. J. Ecol., 5, 201–203, 1980.

Noonan-Wright, E. K., Opperman, T. S., Finney, M. A., Zimmerman, G. T., Seli, R. C., Elenz, L. M., Calkin, D. E., and Fiedler, J. R.: Developing the US Wildland Fire Decision Support System, J. Combust., 2011, 168473, 14 pp., doi:10.1155/2011/168473, 2011.

Pagnini, G. and Mentrelli, A.: Modelling wildland fire propagation by tracking random fronts, Nat. Hazards Earth Syst. Sci., 14, 2249–2263, doi:10.5194/nhess-14-2249-2014, 2014.

Palmer, T. N., and Räisänen, J.: Quantifying the risk of extreme seasonal precipitation events in a changing climate, Nature, 415, 512–514, 2002.

Paugam, R., Wooster, M. J., and Roberts, G.: Use of handheld thermal imager data for airborne mapping of fire radiative power and energy and flame front rate of spread, Geosci. Remote Sens., 51, 3385–3399, 2013.

Perry, G. L. W.: Current approaches to modelling the spread of a wildland fire: a review, Prog. Phys. Geogr., 22, 222–245, 1998.

Peters, W., Miller, J. B., Whitaker, J., Denning, A. S., Hirsch, A., Krol, M. C., Zupanski, D., Bruhwiler, L., and Tans, P. P.: An ensemble data assimilation system to estimate CO₂ surface fluxes from atmospheric trace gas observations, J. Geophys. Res., 110, D24304, doi:10.1029/2005JD006157, 2005.

Peters, W., Jacobson, A. R., Sweeney, C., Andrews, A. E., Conway, T. J., Masarie, K., Miller, J. B., Bruhwiler, L., Pétron, G., Hirsch, A. I., Worthy, E. J., van der Werf, G. R., Randerson, J. T., Wennberg, P. O., Krol, M. C., and Tans, P. P.: An atmospheric perspective on North American carbon dioxide exchange: CarbonTracker, Proc. Natl. Acad. Sci., 104, 18925–18930, 2007.

Pétron, G., Granier, C., Khattatov, B., Lamarque, J.-F., Yudin, V., Muller, J.-F., and Gille, J.: Inverse modeling of carbon monoxide surface emissions using Climate Monitoring and Diagnostics Laboratory network observations, J. Geophys. Res., 107, ACH 10-1–ACH 10-23, doi:10.1029/2001JD001305, 2002.

Poinots, T. and Veynante, D.: Theoretical and Numerical Computation, 2nd Edn., R. T. Edwards Inc., 2005.

Rehm, R. G., McDermott, R. J.: Fire front propagation using the level-set method, NIST, Technical Report 1611, 2009.

Reichle, R. H.: Data assimilation methods in the Earth sciences, Adv. Water Resour., 31, 1411–1418, 2008.

Reichle, R. H., Walker, J. P., Koster, R. D., and Houser, P. R.: Extended versus Ensemble Kalman filtering for land data assimilation, J. Hydrometeorol., 3, 728–740, 2002.

Riggan, P. and Robert, G.: Airborne remote sensing of wildland fires, in: Wildland Fires and Air Pollution, Developments in Environmental Science 8, edited by: Bytnerowicz, A., Arbaugh, M., Andersen, C., and Riebau, A., Elsevier, 139–168, 2009.

Rochoux, M. C., Ricci, S., Lucor, D., Cuenot, B., Trouvé, A., and Bart, J.-M.: Towards predictive simulations of wildfire spread using a reduced-cost Ensemble Kalman filter based on Polynomial Chaos approximations, Proceedings of the Summer Program, Center for Turbulence Research, July 2012, NASA AMES, Stanford University, USA, 2012.

Rochoux, M. C., Cuenot, B., Ricci, S., Trouvé, A., Delmotte, B., Massart, S., Paoli, R., and Paugam, R.: Data assimilation applied to combustion, C. R. Mecanique, 341, 266–276, doi:10.1016/j.crme.2012.10.011, 2013a.

Rochoux, M. C., Delmotte, B., Cuenot, B., Ricci, S., and Trouvé, A.: Regional-scale simulations of wildland fire spread informed by real-time flame front observations, Proc. Combust. Inst., 34, 2641–2647, doi:10.1016/j.proci.2012.06.090, 2013b.

Rochoux, M. C., Emery, C., Ricci, S., Cuenot, B., and Trouvé, A.: Towards predictive data-driven simulations of wildfire spread – Part 2: Ensemble Kalman Filter for the state estimation of a front-tracking simulator of wildfire spread, Nat. Hazards Earth Syst. Sci. Discuss., 2, 3769–3820, doi:10.5194/nhessd-2-3769-2014, 2014.

Rochoux, M. C.: Vers une meilleure prévision de la propagation d'incendies de forêt: evaluation de modèles et assimilation de données, Ph.D. Thesis (written in English: Towards a more comprehensive monitoring of wildfire spread Contributions of model evaluation and data assimilation strategies), Ecole Centrale Paris (France), 2014.

Ros, D. D. and Borga, M.: Adaptive use of a conceptual model for real time flood forecasting, Nord. Hydrol., 28, 169–188, 1997.

Rosić, B. V., Kucerová, A., Sýkora, J., Pajonk, O., Litvinenko, A., and Matthies, H. G.: Parameter identification in a probabilistic setting, Eng. Struct., 50, 179–196, 2013.

Rothermel, R. C.: A mathematical model for predicting fire spread in wildland fuels, USDA Forest Service, Research Paper INT-115, Intermountain Forest and Range Experiment, Ogden, UT:40, 1972.

Ruiz, J. J., Pulido, M., and Miyoshi, T.: Estimating model parameters with ensemble-based data assimilation: a review, J. Meteorol. Soc. Jpn., 91, 79–99, 2013a.

Ruiz, J. J., Pulido, M., and Miyoshi, T.: Estimating model parameters with ensemble-based data assimilation: parameter covariance treatment, J. Meteorol. Soc. Jpn., 91, 453–469, 2013b.
Saad, G. A.: Stochastic data assimilation with application to multi-phase flow and health monitoring problems, Ph.D. Thesis, Faculty of the Graduate School, University of Southern California, 2007.

Sullivan, A. L.: Wildland surface fire spread modeling, 1990–2007, 2. Empirical and quasi-empirical models, Int. J. Wildland Fire, 18, 369–386, 2009.

Szunyogh, I., Kostelich, E., Gyarmati, G., Kalnay, E., Hunt, B., Ott, E., Satterfield, E., and Yorke, J.: A local ensemble transform Kalman filter data assimilation system for the NCEP global model, Tellus A, 60, 113–130, 2008.

Talagrand, O.: Assimilation of observations, an introduction, J. Meteorol. Soc. Jpn., 75, 191–209, 1997.

Tarantola, A.: Inverse Problem Theory, Methods for Data Fitting and Parameter Estimation, Elsevier, 1987.

Thirel, G., Martin, E., Mahfouf, J.-F., Massart, S., Ricci, S., and Habets, F.: A past discharges assimilation system for ensemble streamflow forecasts over France – Part 1: Description and validation of the assimilation system, Hydrol. Earth Syst. Sci., 14, 1623–1637, doi:10.5194/hess-14-1623-2010, 2010.

Tippett, M. K., Anderson, J. L., Bishop, C. H., Hamill, T. M., and Whitaker, J. S.: Ensemble Square Root Filters, Mon. Weather Rev., 131, 1485–1490, 2003.

Todling, R. and Cohn, S. E.: Suboptimal schemes for atmospheric data assimilation based on the Kalman Filter, Mon. Weather Rev., 122, 2530–2557, 1994.

Trémolet, Y.: Model error estimation in 4D-Var, Q. J. Roy. Meteorol. Soc., 133, 1267–1280, 2007.

Tymstra, C., Bryce, R. W., Wotton, B. M., Taylor, S. W., and Armitage, O. B.: Development and structure of Prometheus: the Canadian Wildland Fire Growth Simulation Model, Natural Resources Canada, Canadian Forest Service, Northern Forestry Centre, Edmonton, Alberta, Information Report NOR-X-417, 2010.

Viegas, D. X.: Overview of forest fire propagation research, Proc. Int. Assoc. Fire Safety Sci., 10, 95–108, 2011.

West, M.: Approximating posterior distributions by mixture, J. R. Stat. Soc., 55, 409–422, 1993.

Wiener, N.: The homogeneous chaos, Am. J. Math., 60, 897–936, 1938.

Wooster, M. J., Roberts, G., Perry, G., and Kaufman, Y. J.: Retrieval of biomass combustion rates and totals from fire radiative power observations: FRP derivation and calibration relationships between biomass consumption and fire radiative energy release, J. Geophys. Res., 110, D24311, doi:10.1029/2005JD006318, 2005.

Wooster, M., Robert, G., Smith, A., Johnston, J., Freeborn, P., Amici, S., and Hudak, A.: Thermal remote sensing of active vegetation fires and biomass burning events, Volume Remote Sensing and Digital Image Processing 17, Chapter Thermal infrared remote sensing, Springer, 347–390, 2013.

Xiu, D.: Numerical methods for stochastic computations – a spectral method approach, Princeton University Press, Princeton, NJ, 2010.

Xiu, D. and Karniadakis, G.: The Wiener–Askey polynomial chaos for stochastic differential equations, SIAM J. Scient. Comp., 24, 619–644, 2002.