1 Introduction
Every year, millions of hectares of forest are devastated by wildfires. This fact causes dramatic damage to innumerable factors as economy, ecosystem, energy, agriculture, biodiversity, etc. It has been recognized that the recent increase in the fire severity is associated with the strict fire suppression policy, that over last decades has led to significant accumulation of the fuel, which when ignited makes fires difficult to control. In order to reverse this effect, prescribed burns are routinely used as a method of fuel reduction and habitat maintenance [22, 28]. The previous strategy of putting out all wildland fires is becoming replaced by a new approach where the fire is considered as a tool in the land management practice, and some of the fires are allowed to burn under appropriate conditions in order to reduce the fuel load and meet the forest management goals.
Fire management decisions regarding both prescribed burns, as well as wildland fires, are very difficult. They require a careful consideration of potential fire effects under changing weather conditions, values at risk, firefighter safety and air quality impacts of wildfire smoke [31]. In order to help in the fire management practice, a wide range of models and tools has been developed. The typical operational models are generally uncoupled. In these models, elevation data (slope) and fuel characteristics are used together with ambient weather conditions or general weather forecast as input to the rate of spread model, which computes the fire propagation neglecting the impact of the fire itself on local weather conditions (see BehavePlus [1], FARSITE [9] or PROMETHEUS [29]). As computational capabilities increase, a new generation of coupled fireatmosphere models become available for fire managers as management tools. In a coupled fireatmosphere model, weather conditions are computed inline with the fire propagation. This means that the state of the atmosphere is modified by the fire so that the fire spread model is driven by the local micrometeorology modified by the firereleased heat and moisture fluxes. CAWFE [6], WRFSFIRE [15], and FOREFIRE/MesoNH [8], are examples of such models, coupling CFDtype weather models with semiempirical fire spread models. This approach is fundamentally similar to socalled physicsbased models like FIRETEC [12] and WFDS [19], which also use CFD approach to compute the flow near the fire, but focus on flamescale processes in order to directly resolve combustion, and heat transfer within the fuel and between the fire and the atmosphere. As the computational cost of running these models is too high to facilitate their use as forecasting tools, this paper focuses on the aforementioned hybrid approach, where the fire and the atmosphere evolve simultaneously affecting each other, but the fire spread is parameterized as a function of the wind speed and fuel properties, rather than resolved based on the detailed energy balance.
This article describes upcoming data assimilation components for the coupled fireatmosphere model WRFSFIRE [11, 13], which combines a mesoscale numerical weather prediction system, WRF [27], with a surface fire behavior model implemented by a level set method, a fuel moisture model [30], and chemical transport of emissions. The coupling between the models is graphically represented in the diagram in Fig. 1. The fire heat flux modifies the atmospheric state (including local winds), which in turn affects fire progression and the fire heat release. WRFSFIRE has evolved from CAWFE [3, 4]. An earlier version [15] is distributed with the WRF release as WRFFire [5], and it was recently improved by including a highorder accurate levelset method [20].
The coupling between fire and atmosphere makes initialization of a fire from satellite detections and/or fire perimeters particularly challenging. In a coupled numerical fireatmosphere model, the ignition procedure itself affects the atmospheric state (especially local updrafts near the fire line and the near fire winds). Therefore, particular attention is needed during the assimilation process in order to assure that realistic fireinduced atmospheric circulation is established at the time of data assimilation. One possible solution to this problem, assuring consistency between the fire and the atmospheric models, is defining an artificial fire progression history, and using it to replay the fire progression prior to the assimilation time. In this case, the heat release computed from the synthetic fire history is used to spin up the atmospheric model and assure consistency between the assimilated fire and the local micrometeorology generated by the fire itself.
Fire behavior models run on a mesh given by fuel data availability, typically with about 30m resolution and aligned with geographic coordinates. The mesh resolution of satellitebased sensors, such as MODIS and VIIRS, however, is typically 375m1.1km in flightaligned swaths. These sensors provide planetwide coverage of fire detection several times daily, but data may be missing for various reasons and no detection is possible under clouds; such missing pixels in the swath are marked as not available or as a cloud, and distinct from detections of the surface without fire. Because of the missing data, the statistical uncertainty of detections, the uncertainty in the actual locations of active fire pixels, and the mismatch of scales between the fire model and the satellite sensor, direct initialization of the model from satellite fire detection polygons [7] is of limited value at the fuel map scale. Therefore, the satellite data should be used to steer such models in a statistical sense only.
In this study, we propose a new method of fitting fire arrival time to data, which can be used to generate artificial fire history, which can be used to spin up the atmospheric model for the purpose of starting a simulation from a fire perimeter. In combination with detection data likelihood, the new method can be used also to assimilate satellite fire detection data. This new method, unlike position or additive time corrections, respects the dependence of the fire rate of spread on topography, diurnal changes of fuel moisture, winds, as well as spatial fuel heterogeneity.
2 Fire spread model
The state of the fire spread model is the fire arrival time at locations in a rectangular simulation domain . The isoline is then the fire perimeter at time
. The normal vector to the isoline is
. The rate of spread in the normal direction and the fire arrival time at a location on the isoline then satisfy the eikonal equation(1) 
We assume that depends on location (because of different fuel, fuel moisture, and terrain) and time (because of wind and fuel moisture changing with time). Rothermel’s model [24] for 1D fire spread postulates
(2) 
where is the omnidirectional rate of spread, , the wind factor, is a function of wind in the spread direction, and , the slope factor, is a function of the terrain slope. The 1D model was adapted to the spread over 2D landscape by postulating that the wind factor and the slope factor are functions of the components of the wind vector and the terrain gradient in the normal direction. Thus,
(3) 
The fire spread model is coupled to an atmospheric model. The fire emits sensible and latent heat fluxes, which change the state of the atmosphere, and the changing atmospheric conditions in turn impact the fire (Fig. 1). Wind affects the fire directly by the wind factor, and temperature, relative humidity and rain affect the fire through changing fuel moisture.
The fire model is implemented on a rectangular mesh by finite differences. For numerical reasons, the gradient in the eikonal equation (1) needs to be implemented by an upwindingtype method [21], which avoids instabilities caused by breaking causality in fire propagation: for the computation of at a location , only the values from the directions that the fire is coming from should be used, so the methods switch between onesided differences depending on how the solution evolves. Sophisticated methods of upwinding type, such as ENO or fluxlimiters [23], aim to use more accurate central differences and switch to more stable onesided upwind differences only as needed. Unfortunately, the switching causes the numerical gradient of at a mesh node become a nondifferentiable function of the values of at that point and its neighbors. In addition, we have added a penalty term to prevent the creation of local minima. It was observed in [14] that if, in the level set method, a local minimum appears on the boundary, its value keeps decreasing out of control; we have later found out that this can in fact happen anywhere in the presence of spatially highly variable rate of spread, and we have observed a similar effect here during the minimization process.
3 Fitting the fire spread model to data
3.1 Minimal residual formulation
Consider the situation when the two observed fire perimeters and at times are known, and we are interested in the fire progression between the two perimeters. Aside from immediate uses (visualization without jumps, postfire analysis), such interpolation is useful to start the fire simulation from the larger perimeter at time by a spinup of the atmospheric model by the heat fluxes from the interpolated fire arrival time between the fire perimeters; the coupled model can then start from perimeter at time in a consistent state between the fire and the atmosphere. Interpolation between an ignition point and a perimeter can be handled the same way, with the perimeter consisting of just a single point.
In this situation, we solve the eikonal equation (1) only approximately,
(4) 
imposing the given fire perimeters as constraints,
(5) 
We formalize (4) as the minimization problem
(6) 
where is a function such that if and only if , and is the simulation domain. We mostly use the function but other functions, such as have advantages in some situations. There are no boundary conditions imposed on the boundary of .
3.2 Discretization and the constraint matrix
The fire simulation domain is discretized by a logically rectangular grid (aligned approximately with longitude and latitude) and perimeters are given as shape files, i.e., collections of points on the perimeter. We express (5) in the form
(7) 
where is a sparse matrix. Since the points in the shape files do not need to lie on the grid, the rows of are the coefficients of an interpolation from the grid to the points in the shape files, which define the perimeters. We find the coefficients from barycentric interpolation. The rectangles of the grid are split into two triangles each, and, for each triangle, we compute the barycentric coordinates of the points in the shapefile, i.e., the coefficients of the unique linear combination of the vertices of the triangle that equals to the point in the shape file. If all 3 barycentric coordinates are in , we conclude that the point is contained in the triangle, the barycentric coordinates are the sought interpolation coefficients, and they form one row of . For efficiency, most points in the shapefile are excluded up front, based on a comparison of their coordinates with the vertices of the triangle, which is implemented by a fast binary search. When there is more than one point of the shapefile in any triangle, we condense them into a single constraint, obtained by adding the relevant rows of . This way, we avoid over constraining the fire arrival time near the perimeter, which should be avoided for the same reason as limiting the number of constraints in mixed finite elements to avoid locking, cf., e.g., [2].
3.3 Numerical minimization of the residual
To solve (6) numerically, we use a multiscale descent method similar to multigrid, combining line searches in the direction of changes of the value of at a single point, and linear combinations of point values as in [18]. We use bilinear coarse grid functions with the coarse mesh step growing by a factor of 2. See Fig. 6(b) for an example of a coarse grid function with distance between nodes 16 mesh steps on the original, finest level. We start from an initial approximate solution that satisfies the constraint exactly, and project all search directions on the subspace , so that the constraint remains satisfied throughout the iterations.
To find a reasonable initial approximation to the fire arrival time, we solve the quadratic minimization problem
(8) 
where is the normal direction, is the Laplace operator, and is generally noninteger. The reason for choosing is that is the Sobolev seminorm and in 2D, the space is embedded in continuous functions if and only if . Consequently, is not a bound on the value at any particular point, only averages over some area can be controlled. Numerically, when , minimizing with a point constraint, such as an ignition point, results in taking the shape of a sharp funnel at that point (Fig. 5), which becomes thinner as the mesh is refined. That would be definitely undesirable.
The discrete form of (8) is
(9) 
where with a discretization of the Laplace operator with Neumann boundary conditions. To solve (9), we first find a feasible solution , so that , substitute to get
and augmenting the cost fuction, we get that (9) is equivalent to
(10) 
where , is the orthogonal projection on the nullspace of , and is an arbitrary regularization parameter. We solve the minimization problem (10) approximately by preconditioned conjugate gradients for the equivalent symmetric positive definite linear system
(11) 
Since is discretization of the Neumann problem, the preconditioner requires some care. Define as the vector that generates the nullspace of , which consists of the discrete representation of constant functions, and the orthogonal projection on its complement. We use the preconditioner
where is the inverse of on the complement of its nullspace, and recover the solution by . The method only requires access to matrixvector multiplications by and , which are readily implemented by cosine FFT. We only need to solve (11) to low accuracy to get a reasonable starting point for the nonlinear iterations, but the satisfaction of the constraint to rounding precision is important.
4 Assimilation of MODIS and VIIRS fire detections
Data likelihood is the probability of a specific configuration of fire detection and nondetection pixels given the state of the fire. The probability of MODIS Active Fires detection in a particular sensor pixel as a function of the fraction of the area actively burning and the maximum size of contiguous area burning, was estimated in the validation study
[25]using logistic regression. We consider the fraction of the pixel burning and the maximum continuous area burning as a proxy to the fire radiative heat flux in the pixel. The model state is encoded as the fire arrival time at each grid point, and the heat flux can be then computed from the burn model using the fuel properties. Substituting the heat flux into the logistic curve yields a plausible probability of detection for a period starting from the fire arrival time: the probability keeps almost constant while the fire is fresh, and then diminishes.
However, the position uncertainty of the detection is significant, the allowed error is listed in VIIRS specifications [26] as km, and position errors of such magnitude are indeed occasionally observed. Therefore, the probability of detection at the given coordinates of the center of a sensor pixel in fact depends on the fire over a nearby area, with the contributions of fire model cells weighted by , where is the distance of the fire model cell and the nominal center of the sensor pixel, because of the uncertainty where the sensor is actually looking. Assuming that the position errors and the detection errors are independent, we can estimate the contribution of a grid cell to the data likelihood from a combination of the probabilities of detection at the nearby satellite pixels.
Assimilation of data into the fire spread model can be then formulated as an optimization problem to minimize its residual and to maximize the data likelihood. See [10] for further details.
Since the fire model is coupled with an atmosphere model, changing the state of the fire alone makes the state of the coupled model inconsistent. To recover a consistent state, we spin up the atmosphere model from an earlier time, with the modified fire arrival time used instead of the fire arrival time from the fire spread model (Fig. 2). This synthetic fire forcing to the atmospheric model is used to drive atmospheric model [16] and enables establishing fireinduced circulation.
Varying the model state to maximize the data likelihood can also be used to estimate the time and place of ignition as well as other model parameters. The WRFSFIRE [15] model was run on a mesh of varying GPS coordinates and times and the data likelihoods of the relevant Active Fire detection data is evaluated, allowing the most likely place and time of the fire’s ignition to be determined. Fig. 3 shows a visualization of the likelihoods of Active Fire detection data for several hundred ignition points at various times. Work is in progress so that an automated process of determining the most likely time and place of ignition can be initiated from collection of satellite data indicating a wildfire has started in a particular geographic region of interest.
5 Computational experiments
(a)  (b) 
(a)  (b) 
(c)  (d) 
(a)  (b) 
(c)  (d) 
The optimization problem was tested on an idealized case using concentric circles as perimeters in a mesh with nodes. The fire spreads equally in all directions from the center of the mesh. The propagation is set at different rates of spread in different sections (Fig. 4(a)). We also set the fire arrival time at the ignition point and compute the fire arrival time on the two perimeters from the given rate of spread, so in this case there exists an exact solution (Fig. 4(b)).
The constraint matrix was constructed by the method described in Sect. 3.2. The initial approximation of the fire arrival time was then found by solving the quadratic minimization problem described in Sect. 3.3 with . Fig. 5 shows the initial approximation of the fire arrival time imposed by the ignition point and the two concentric circles in our particular case and using different values of from 1 to 1.4. One can see how the unrealistic sharp funnel at the ignition point for disappears with the increasing value of .
Then, we run the multigrid method proposed in Sect. 3.3. The coarsening was done by the ratio of 2. The number of sweeps was linearly increasing with the level. On the coarsest level, the mesh step was 32 and the sweep was done once, the mesh step on the second level was 16 and the sweep was repeated twice, until resolution 1 on the original, finest grid, and sweep repeated 6 times.
Fig. 6c shows the decrease in the cost function with the number of line searches on any level. One can observe that the cost function decreased more in the first cycle and at the beginning of iterations on each level.
The final result after 4 cycles of 6 different resolutions (from 32 to 1 decreasing by powers of two) is shown in Fig. 6(d), which is close to the exact solution.
6 Conclusions
We have presented a new method for fitting data by an approximate solution of a fire spread model. The method was illustrated on an idealized example. Application to a real problem are forthcoming.
7 Acknowledgments
This research was partially supported by grants NSF ICER1664175 and NASA NNX13AH59G, and MINECOSpain under contract TIN201453234C21R. Highperformance computing support at CHPC at the University of Utah and Cheyenne (doi:10.5065/D6RX99HX) at NCAR CISL, sponsored by the NSF, are gratefully acknowledged.
References
 [1] Andrews, P.L.: BehavePlus fire modeling system: past, present, and future. Paper J2.1, 7th Symposium on Fire and Forest Meteorology (2007), http://ams.confex.com/ams/pdfpapers/126669.pdf, retrieved September 2011
 [2] Brezzi, F., Fortin, M.: Mixed and Hybrid Finite Element Methods. SpringerVerlag, New York – Berlin – Heidelberg (1991)
 [3] Clark, T.L., Coen, J., Latham, D.: Description of a coupled atmospherefire model. International Journal of Wildland Fire 13, 49–64 (2004). https://doi.org/10.1071/WF03043
 [4] Coen, J.L.: Simulation of the Big Elk Fire using coupled atmospherefire modeling. International Journal of Wildland Fire 14(1), 49–59 (2005). https://doi.org/10.1071/WF04047
 [5] Coen, J.L., Cameron, M., Michalakes, J., Patton, E.G., Riggan, P.J., Yedinak, K.: WRFFire: Coupled weatherwildland fire modeling with the Weather Research and Forecasting model. J. Appl. Meteor. Climatol. 52, 16–38 (2013). https://doi.org/10.1175/JAMCD12023.1
 [6] Coen, J.L.: Modeling wildland fires: A description of the Coupled AtmosphereWildland Fire Environment model (CAWFE). NCAR Technical Note NCAR/TN500+STR, http://dx.doi.org/10.5065/D6K64G2G (2013)
 [7] Coen, J.L., Schroeder, W.: Use of spatially refined satellite remote sensing fire detection data to initialize and evaluate coupled weatherwildfire growth model simulations. Geophysical Research Letters 40, 1–6 (2013). https://doi.org/10.1002/2013GL057868
 [8] Filippi, J.B., Bosseur, F., Pialat, X., Santoni, P., Strada, S., Mari, C.: Simulation of coupled fire/atmosphere interaction with the MesoNHForeFire models. Journal of Combustion 2011, Article ID 540390 (2011). https://doi.org/10.1155/2011/540390
 [9] Finney, M.A.: FARSITE: Fire Area Simulator – model development and evaluation. Res. Pap. RMRSRP4 Revised 2004, Ogden, UT, USDA Forest Service, Rocky Mountain Research Station (1998). https://doi.org/10.2737/RMRSRP4
 [10] Haley, J., Farguell Caus, A., Mandel, J., Kochanski, A.K., Schranz, S.: Data likelihood of active fires satellite detection and applications to ignition estimation and data assimilation. In: Viegas, D.X. (ed.) VIII International Conference on Forest Fire Research. University of Coimbra Press (2018), submitted
 [11] Kochanski, A.K., Jenkins, M.A., Yedinak, K., Mandel, J., Beezley, J., Lamb, B.: Toward an integrated system for fire, smoke, and air quality simulations. International Journal of Wildland Fire 25, 534–546 (2016). https://doi.org/10.1071/WF14074
 [12] Linn, R., Reisner, J., Colman, J.J., Winterkamp, J.: Studying wildfire behavior using FIRETEC. Int. J. of Wildland Fire 11, 233–246 (2002). https://doi.org/10.1071/WF02007
 [13] Mandel, J., Amram, S., Beezley, J.D., Kelman, G., Kochanski, A.K., Kondratenko, V.Y., Lynn, B.H., Regev, B., Vejmelka, M.: Recent advances and applications of WRFSFIRE. Natural Hazards and Earth System Science 14(10), 2829–2845 (2014). https://doi.org/10.5194/nhess1428292014

[14]
Mandel, J., Beezley, J.D., Coen, J.L., Kim, M.: Data assimilation for wildland fires: Ensemble Kalman filters in coupled atmospheresurface models. IEEE Control Systems Magazine
29(3), 47–65 (June 2009). https://doi.org/10.1109/MCS.2009.932224  [15] Mandel, J., Beezley, J.D., Kochanski, A.K.: Coupled atmospherewildland fire modeling with WRF 3.3 and SFIRE 2011. Geoscientific Model Development 4, 591–610 (2011). https://doi.org/10.5194/gmd45912011
 [16] Mandel, J., Beezley, J.D., Kochanski, A.K., Kondratenko, V.Y., Kim, M.: Assimilation of perimeter data and coupling with fuel moisture in a wildland fire – atmosphere DDDAS. Procedia Computer Science 9, 1100–1109 (2012). https://doi.org/10.1016/j.procs.2012.04.119, Proceedings of ICCS 2012
 [17] Mandel, J., Fournier, A., Haley, J.D., Jenkins, M.A., Kochanski, A.K., Schranz, S., Vejmelka, M., Yen, T.Y.: Assimilation of MODIS and VIIRS satellite active fires detection in a coupled atmospherefire spread model. Poster, 5th Annual International Symposium on Data Assimilation, 1822 July 2016, University of Reading, UK (2016), http://www.isda2016.net/abstracts/posters/MandelAssimilationof.html, retrieved December 2016
 [18] McCormick, S.F., Ruge, J.W.: Unigrid for multigrid simulation. Math. Comp. 41(163), 43–62 (1983). https://doi.org/doi.org/10.2307/2007765
 [19] Mell, W., Jenkins, M.A., Gould, J., Cheney, P.: A physicsbased approach to modelling grassland fires. Intl. J. Wildland Fire 16, 1–22 (2007). https://doi.org/10.1071/WF06002
 [20] MuñozEsparza, D., Kosović, B., Jiménez, P.A., Coen, J.L.: An accurate firespread algorithm in the Weather Research and Forecasting model using the levelset method. Journal of Advances in Modeling Earth Systems (2018). https://doi.org/10.1002/2017MS001108
 [21] Osher, S., Fedkiw, R.: Level Set Methods and Dynamic Implicit Surfaces. Springer, New York (2003)
 [22] Outcalt, K.W., Wade, D.D.: Fuels management reduces tree mortality from wildfires in southeastern United States. Southern Journal of Applied Forestry 28(1), 28–34 (2004)
 [23] Rehm, R.G., McDermott, R.J.: Firefront propagation using the level set method. NIST Technical Note 1611 (March 2009), https://nvlpubs.nist.gov/nistpubs/Legacy/TN/nbstechnicalnote1611.pdf, accessed May 2018
 [24] Rothermel, R.C.: A mathematical model for predicting fire spread in wildland fires. USDA Forest Service Research Paper INT115 (1972), https://www.fs.fed.us/rm/pubs_int/int_rp115.pdf, accessed March 2018
 [25] Schroeder, W., Prins, E., Giglio, L., Csiszar, I., Schmidt, C., Morisette, J., Morton, D.: Validation of GOES and MODIS active fire detection products using ASTER and ETM+ data. Remote Sensing of Environment 112(5), 2711–2726 (2008). https://doi.org/10.1016/j.rse.2008.01.005
 [26] Sei, A.: VIIRS active fires: Fire mask algorithm theoretical basis document (2011), available at https://www.star.nesdis.noaa.gov/jpss/documents/ATBD/D0001M01S01021_JPSS_ATBD_VIIRSActiveFires.pdf, retrieved November 17, 2013
 [27] Skamarock, W.C., Klemp, J.B., Dudhia, J., Gill, D.O., Barker, D.M., Duda, M.G., Huang, X.Y., Wang, W., Powers, J.G.: A description of the Advanced Research WRF version 3. NCAR Technical Note 475 (2008). https://doi.org/10.5065/D68S4MVH
 [28] Stephens, S.L., Ruth, L.W.: Federal forestfire policy in the United States. Ecological Applications 15(2), 532–542 (2005). https://doi.org/doi:10.1890/040545
 [29] Tymstra, C., Bryce, R., Wotton, B., Taylor, S., Armitage, O.: Development and structure of Prometheus: the canadian wildland fire growth simulation model. Information Report NORX147, Northern Forestry Centre, Canadian Forest Service (2010), http://publications.gc.ca/collections/collection_2010/nrcan/Fo1331417eng.pdf, accessed March 2018
 [30] Vejmelka, M., Kochanski, A.K., Mandel, J.: Data assimilation of dead fuel moisture observations from remote automatic weather stations. International Journal of Wildland Fire 25, 558–568 (2016). https://doi.org/10.1071/WF14085
 [31] Yoder, J., Engle, D., Fuhlendorf, S.: Liability, incentives, and prescribed fire for ecosystem management. Frontiers in Ecology and the Environment 2, 361–366 (Jan 2004). https://doi.org/10.1890/15409295(2004)002[0361:LIAPFF]2.0.CO;2
Comments
There are no comments yet.