Coupled atmosphere-wildland fire modeling with WRF-Fire

We describe the physical model, numerical algorithms, and software structure of WRF-Fire. WRF-Fire consists of a fire-spread model, implemented by the level-set method, coupled with the Weather Research and Forecasting model. In every time step, the fire model inputs the surface wind, which drives the fire, and outputs the heat flux from the fire into the atmosphere, which in turn influences the atmosphere. The level-set method allows submesh representation of the burning region and flexible implementation of various ignition modes. WRF-Fire is distributed as a part of WRF and it uses the WRF parallel infrastructure for parallel computing.


Introduction
Wildland fires impact the lives of millions of people and cause major damage every year worldwide, yet they are a natural part of the cycle of nature. Better tools for modeling wildland fire behavior are important for managing fire suppression, planning controlled burns to reduce the fuels, as well as to help assess fire danger. Fire models range from tools based on Rothermel (1972) fire spread rate formulas, such as BehavePlus (Andrews, 2007) and FARSITE (Finney, 1998), suitable for operational forecasting, to sophisticated 3-D computational fluid dynamics and combustion simulations suitable for research and reanalysis, such as FIRETEC (Linn et al., 2002) and WFDS (Mell et al., 2007). BehavePlus, the PCbased successor of the calculator-based BEHAVE, determines the fire spread rate at a single point from fuel and environmental data; FARSITE uses the fire spread rate to provide a 2-D simulation on a PC; while FIRETEC and WFDS require a parallel supercomputer and run much slower than real time.
Wildland fire is a complicated multiscale process, from the flame reaction zone on milimeter scale to the synoptic weather scale of hundreds of kilometers. Since direct numerical simulation of wildland fire is computationally intractable and detailed data are not available anyway, compromises in the choice of processes to be modeled, approximations, and parametrizations are essential. Fortunately, a practically important range of wildland fire behavior can be captured by the coupling of a mesoscale weather model with a simple 2-D fire spread model (Clark et al., 1996a,b). Weather has a major influence on wildfire behavior; in particular, wind plays a dominant role in the fire spread. Conversely, the fire influences the atmosphere through the heat and vapor fluxes from burning hydrocarbons and evaporation of fuel moisture. Fire heat output has a major effect on the atmosphere; the buoyancy created by the heat from the fire can cause tornadic strength winds, and the air motion and moisture from the fire can affect the atmosphere also away from the fire. It is well known that a large fire "creates its own weather". The correct wildland fire shape and progress result from the two-way interaction between the fire and the atmosphere (Clark et al., 1996a(Clark et al., ,b, 2004Coen, 2005). WRF-Fire (Mandel et al., 2009) combines the Weather Research and Forecasting Model (WRF) with the ARW dynamical core (Skamarock et al., 2008) with a semi-empirical fire spread model. It is intended to be faster than real time in order to deliver a prediction.
WRF-Fire has grown out of the NCAR's CAWFE code (Clark et al., 1996a(Clark et al., ,b, 2004Coen, 2005). CAWFE consists of the Clark-Hall mesoscale atmospheric model, coupled with a tracer-based fire spread model. Although the Clark-Hall model has many good properties, it is a legacy serial code, not supported, and difficult to modify or use with real data, while WRF is a parallel supported community code routinely used for real runs. See Coen and Patton (2010) for a further discussion of their relative merits in the wildland fire appli-2 cation. WRF-Fire was started by Patton and Coen (2004), who proposed a combination of WRF with the tracer-based model from CAWFE, formulated a road map, and made the important observation that the innermost domain of the weather code, which interacts directly with the fire model, needs to run in the Large Eddy Simulation (LES) mode. Patton ported the Fortran 77-based fire module to Fortran 90 and developed the initial serial coupled code. However, instead of using the existing tracer-based CAWFE code, the fire module in WRF-Fire was developed based on the level-set method (Osher and Fedkiw, 2003). One of the reasons was that the representation of the fire region by the level-set function was thought to be more flexible than the representation of the burning region in CAWFE by four tracers in each cell of the fire mesh. In particular, the level-set function can be manipulated more easily than tracers for the purpose of data assimilation. Insertion of the heat fluxes, while fundamentally the same as in CAWFE, had to be redone for WRF variables already in Patton's initial code. Thus, only the code for the calculation of the fire spread rate and the heat fluxes remained from CAWFE. While WRF-Fire takes advantage of the experience accumulated with CAWFE, WRF is quite different from the Clark-Hall atmospheric model and the fireline propagation algorithm is also different. Thus, it needs to be demonstrated that WRF-Fire can deliver similar results as CAWFE, and WRF-Fire needs to be validated against real fires (Sect. 10). The level-set method was used for a surface fire spread model in Mallet et al. (2009). Filippi et al. (2009) coupled the atmospheric model Meso-nh with fire propagation by tracers. Tiger (Mazzoleni and Giannino, 2010) uses a 2-D combusion model based on reactionconvection-diffusion equations and a convection model to emulate the effect of the fire on the wind. FIRESTAR (Morvan and Dupuy, 2004) is a physically accurate wildland fire model in two dimensions, one horizontal and one vertical. UU LES-Fire (Sun et al., 2009) couples the University of Utah's Large Eddy Simulation code with the tracer-based code from CAWFE. See the survey by Sullivan (2009) for a number of other models. WRF-Fire was briefly treated as one of the topics in Mandel et al. (2009). The purpose of this paper is to describe the fire module and the coupling with WRF in the current WRF-Fire code in sufficient detail, yet understandable to a reader not familiar with WRF. In addition, the advances since the paper Mandel et al. (2009) was written in 2007 include new, practically important ignition schemes (Sect. 4.4), vertical interpolation of the wind in the boundary layer dependent on land-use (Sect. 6), parallel computing (Sect. 7), input of real data (Sect. 9), and validations in progress on real fires (Sect. 10). This paper also contains reproducible descriptions of the physical model (Sect. 2), the required WRF settings (Sect. 5), and the coupling of WRF with the fire module (Sect. 6).
WRF-Fire is public domain software and it has been distributed as a part of the WRF source code at wrf-model.org since version 3.2, released in April 2010 (Dudhia, 2010). The 3 released version is updated periodically and supported by NCAR. The current development version of WRF-Fire with the latest features and bug fixes and additional visualization tools, guides, and diagnostic utilities, are available directly from the developers at openwfm.org. This article describes WRF-Fire as it is scheduled to be included in WRF 3.3, to be released in March or April 2011. WRF-Fire user's guide is available as a part of the WRF user's guide (Wang et al., 2010), to be updated with the release.
New features since WRF version 3.2 include new ignition models, vertical interpolation of the wind from logarithmic profile, fetching high-resolution geogrid data, terrain gradient interpolation, and optional input of fuel map, land use map, and high-resolution topography in ideal runs.

Physical fire model and fuels
The physical model consists of functions specifying the fire spread rate and the heat fluxes, and it is essentially the same as a subset of CAWFE (Clark et al., 2004;Coen, 2005). The spread rate calculation is in turn based on BEHAVE (Rothermel, 1972;Andrews, 2007). It is described here in more detail for the sake of reproducibility and to point out the (minor) differences.

Fuel properties
Fuel is characterized by the quantities listed in Table 1, which are given at every point of the fire mesh. To simplify the specification of fuel properties, fuels are given as one of 13 Anderson (1982) categories, which are preset vectors of values of the fuel properties. These values are specified in an input text file (namelist.fire), and they can be modified by the user. The user can also specify completely new, custom fuel categories.

Fire spread rate
The fire model is posed in the horizontal (x, y) plane the Earth surface is projected on. The semi-empirical approach to fire propagation used here assumes that the fire spread rate is given by the modified Rothermel (1972) formula where R 0 is the spread rate in the absence of wind, φ W is the wind factor, and φ S is the slope factor. The components of (1) computed from the fuel properties (Table 1), the wind speed U , and the terrain slope tanφ following the equations in Table 2. See Rothermel (1972) for further details, derivation, and justifications. 4 Chaparral is a special fuel in that the spread rate depends only on wind speed. For chaparral, (1) is replaced by (Coen et al., 2001, Eq. 1) S = min{6,max{0.03333,1.2974U }}. (2) The only differences from Rothermel (1972) are the subtraction of the moisture from the fuel load in the computation rather than up front, limiting the slope and the windspeed, the special chaparral spread rate from CAWFE (2), and the explicit reduction of wind from 6.1 m height to midflame height, following Baughman and Albini (1980).
In either case, the spread rate can be written as where S 0 , R 0 , b, c, d, e are the fuel-dependent coefficients that represent the spread rate internally. These coefficients are stored for every grid point. At a point on the fireline, denote by n the outside normal to the fire region, U the wind vector, and z the terrain height. The normal component of the wind vector, U = U · n, and the normal component of the terrain gradient, tanφ = ∇z · n, are used to determine the spread rate, which is interpreted as the spread rate in the normal direction n.

Fuel burned and heat released
Each location starts with fuel fraction F = 1. Once the fuel is ignited at a time t i , the fuel fraction decreases exponentially, where t is the time, t i is the ignition time, F 0 is the initial amount of fuel, and T f is the fuel burn time, i.e., the number of seconds for the fuel to burn down to 1/e ≈ 0.3689 of the original quantity. Since by definition of the fuel weight w (Table 1), the fuel burns down to 0.6 of the original quantity in 600 s when w = 1000, we have which gives The input coefficient w is used in WRF-Fire rather than T f for compatibility with existing fuel models and literature. The average sensible heat flux density released in time interval (t,t + ∆t) is computed as and the average latent heat (i.e., moisture) flux density is given by where 0.56 is the estimated mass ratio of the water output from the combustion to the dry fuel, and L = 2.5 × 10 6 J kg −1 is the specific latent heat of condensation of water at 0 • C, used for nominal conversion of moisture to heat. This computation is from CAWFE.

Domain, grids, and nodes
The atmospheric model operates on a logically quadrilateral 3-D grid on the Earth surface, and uses a sequence of horizontally nested grids, called domains (Kalnay, 2003). Only the innermost (the finest) atmospheric domain is coupled with the fire model; see also Sect. 8. Scalar variables in the atmospheric model are located at the centers of the 3D grid cells, while the wind vector components are at a staggered grid at the midpoints of the cell faces. The fire model operates on a refined fire mesh called the subgrid (Fig. 1), and all of its variables are all represented by their values at the centers of the cells of this fire subgrid.
4 Mathematical core of the fire model

Fire propagation by the level-set method
The model maintains a level-set function ψ, the time of ignition t i , and the fuel fraction F . Denote a point on the surface by x = (x,y). The burning region at time t is represented by a level-set function ψ = ψ(x,t) as the set of all points x such that ψ(x,t) ≤ 0. There is no fire at x if ψ(x,t) > 0. The fireline is the set of all points x such that ψ(x,t) = 0. Since on the fireline, the tangential component of the gradient ∇ψ is zero, the outside normal vector at the fireline is n = ∇ψ ∇ψ .
Now consider a point x(t) that moves with the fireline. Then the fire spread rate S at x in the direction of the normal n is and, from the definition of the fireline, ψ(x(t),t) = 0. By the chain rule and substituting from (7) and (8), we have So, the evolution of the level-set function is governed by the partial differential equation called the level-set equation (Osher and Fedkiw, 2003). The spread rate S is evaluated from (3) for all x, not just on the fireline. Since S ≥ 0, the level-set function does not increase with time, and the fire area cannot decrease, which also helps with numerical stability by eliminating oscillations of the level-set function ψ in time.
The level-set equation is discretized on a rectangular grid with spacing ( x, y), called the fire grid. The level-set function ψ and the ignition time t i are represented by their values at the centers of the fire grid cells. This is consistent with the fuel data given in the center of each cell also.
To advance the fire region in time, we use Heun's method (Runge-Kutta method of order 2), The right-hand side F is a discretization of the term −S ∇ψ with upwinding and artificial viscosity, where n = ∇ψ/ ∇ψ is computed by finite central differences and ∇ψ = ∇ x ψ,∇ y ψ is the upwinded finite difference approximation of ∇ψ by Godunov's method (Osher and Fedkiw, 2003, p. 58), where ∇ + x ψ and ∇ − x ψ are the right and left one-sided finite differences and similarly for ∇ + y ψ and ∇ − y ψ. Further, in (12), ε is scale-free artificial viscosity (ε = 0.4 here), and x + similar term for y is the five-point Laplacian of ψ scaled so that the artificial viscosity is proportional to the mesh step, A numerically stable scheme with upwinding, such as (13), is required to compute the term ∇ψ in the level set equation (10). However, in our tests, the gradient by standard central differences, worked better in the computation of the normal vector n by (7), which is used to evaluate the normal component of the wind and the slope in (3).
Before computing the finite differences up to the boundary, the level-set function is extrapolated to one layer of nodes beyond the boundary. However, the extrapolation is not allowed to decrease the value of the level-set function to less than the value at either of the points it is extrapolated from. For example, when (i,j) is the last node in the domain in the direction x, the extrapolation is used, and similarly in the other cases. This is needed to avoid numerical instabilities at the boundary. Otherwise, a decrease in ψ at a boundary node, which may happen with nonhomogeneous fuels in real data, is amplified by the extrapolation, and ψ keeps decreasing at that boundary node in every time step until it becomes negative, starting a spurious fire.
The model does not support fire crossing the boundary of the domain. When ψ < 0 is detected near the boundary, the simulation terminates. This is not a limitation in practice, because the fire should be well inside the domain anyway for a proper response of the atmosphere.

Computation of the ignition time
The ignition time t i in the strip that the fire has moved over in one time step is computed by linear interpolation from the level-set function. Suppose that the point x is not burning at time t but is burning at time t + t, that is, ψ(x,t) > 0 and ψ(x,t + t) ≤ 0. The ignition time at x satisfies ψ(x,t i (x)) = 0. Approximating ψ by a linear function in time, we have and we take

Computation of the fuel fraction
The fuel fraction is approximated over each fire mesh cell C by integrating (4) over the fire region. Hence, the fuel fraction remaining in cell C at time t is given by Once the fuel fraction is known, the heat fluxes are computed from (5) and (6). This scheme has the advantage that the total heat released in the atmosphere over time is exact, regardless of approximations in the computation of the integral (15). Our objective in the numerical evaluation of (15) is a method that is second order accurate when the whole cell is on fire, exact when no part of the cell C is on fire (namely, returning the value one), and provides a natural transition between these two cases. Just like standard schemes in numerical analysis can be derived from the requirement that they are exact for all polynomials up to a given degree, the guiding principle here is that the scheme should be exact in as many special cases as possible. Then we expect that the scheme should work well overall.
While the fuel burn time T f can be interpolated as constant over the whole cell, the levelset function ψ and the ignition time t i must be interpolated more accurately to allow a submesh representation of the burning area and a gradual release of the heat as the fireline moves over the cell. In addition, we need the fuel fraction computed over each mesh cell, because the heat fluxes in the mesh cells are summed up to give the heat flux in an atmospheric cell. Our solution is to split each cell into 4 subcells C j , interpolate to the corners of the subcells, and add the integrals, cf., Fig. 2. The level-set function ψ is interpolated bilinearly to the vertices of the subcells C j , and the burn time T f is constant on each C j , given by its value at the fire grid nodes. However, to interpolate the ignition time t i we first define t i outside of the fire region and on the fireline by This allows us to omit the condition ψ ≤ 0 in the definition of the integration domains in (16) and integrate on the whole cells, respective subcells, only. Then, we interpolate t i bilinearly to the vertices of the subcells C j and correct the resulting values by applying the compatibility condition (17). To compute the integral over a subcell C j , we first estimate the fraction of the subcell that is burning, by where x k are the the corners of the subcell C j . This approximation is exact when no part of the subcell C j , is on fire, that is, all ψ(x k ) ≥ 0 and at least one ψ(x k ) > 0; the whole C j 10 is on fire, that is, all ψ(x k ) ≤ 0 and at least one ψ(x k ) < 0; or the values ψ(x k ) define a linear function and the fireline crosses the subcell diagonally or it is aligned with one of the coordinate directions. Next, replace t i (x k ) by t when ψ(x k ) > 0 (i.e., the node x k is not on fire), and compute the approximate fraction of the fuel burned as This calculation is accurate asymptotically when the fuel burns slowly and the approximation β of the burning area is exact.

Ignition
Typically, a fire starts from a horizontal extent much smaller than the fire mesh cell size, and both point and line ignition need to be supported. The previous ignition mechanism (Mandel et al., 2009) ignited everything within a given distance from the ignition line at once. This distance was required to be at least one or two mesh steps, so that the initial fire is visible on the fire mesh, and the fire propagation algorithm from Sect. 4.1 can catch on. This caused an unrealistically large initial heat flux and the fire started too fast. The current ignition scheme achieves submesh resolution and zero-size ignition. A small initial fire is superimposed on the regular propagation mechanism, which then takes over. Drip-torch ignition is implemented as a collection of short ignition segments that grows at one end every time step. Multiple ignition segments are also supported.
The model is initialized with no fire by choosing the level-set function ψ(x,t 0 )=const > 0. Consider an initial fire that starts at time t g on a segment a,b and propagates in all directions with an initial spread rate S g until the distance r g is reached. At the beginning of every time step t such that we construct the level-set function of the initial fire, and replace the level-set function of the model by For a drip-torch ignition starting from point a at time t g at velocity v until time t h , the ignition line at time t is the segment a,a + v(min{t,t h } − t g ), and (20) becomes followed again by (21), at the beginning of every time step begining at time t such that The ignition time of newly ignited nodes is set to the arrival time of the fire at the spread rate S g from the nearest point on the ignition segment.

Atmospheric model
We summarize some background information about WRF-ARW from Skamarock et al. (2008), to the extent needed to understand the coupling with the fire module. The model is formulated in terms of the hydrostatic pressure vertical coordinate η, scaled and shifted so that η = 1 at the Earth surface and η = 0 at the top of the domain. The governing equations are a system of partial differential equations of the form where R contains also the advection terms, and Φ = (U,V,W,φ ,Θ,µ ,Q m ). The fundamental WRF variables are µ = µ(x,y), the hydrostatic component of the pressure differential of dry air between the surface and the top of the domain, written in perturbation form µ = µ + µ , where µ is a reference value in hydrostatic balance; U = µu, where u = u(x,y,η) is the Cartesian component of the wind velocity in the x-direction, and similarly V and W ; Θ = µθ, where θ = θ(x,y,η) is the potential temperature; φ = φ(x,y,η) = φ + φ is the geopotential; and Q m = µq m is the moisture content of the air. The variables in the state Φ evolved by (22) are called prognostic variables. Other variables computed from them, such as the hydrostatic pressure p, the thermodynamic temperature T , and the height z, are called diagnostic variables. The variables that contain µ are called coupled. The value of the right-hand side R(Φ) is called tendency. See Skamarock et al. (2008, p. 7-13) for details and the form of R.
The system (22) is discretized in time by the explicit 3rd order Runge-Kutta method where the differential operator R is discretized by finite differences and the tendencies from physics packages, such as the fire module, are updated only the third Runge-Kutta step (Skamarock et al., 2008, p. 16). In order to avoid small time steps, the tendency in the third Runge-Kutta step also includes the effect of substeps to integrate acoustic modes.
6 Coupling of the fire and the atmospheric models The terrain gradient is computed from the terrain height at the best available resolution and interpolated to the fire mesh in preprocessing. Interpolating the height and then computing the gradient would cause jumps in the gradient, which affect fire propagation, unless highorder interpolation is used.
In each time step of the atmospheric model, the fire module is called from the third step of the Runge-Kutta method. First the wind is interpolated to a given height z f above the terrain (currently, 6.1 m following BEHAVE), assuming the logarithmic wind profile where z is the height above the terrain and z 0 is the roughness height. For a fixed horizontal location, denote by z 1 , z 2 ,... the heights of the centers of the atmospheric mesh cells; these are computed from the geopotential φ, which is a part of the solution. The horizontal wind component u(z f ) under the u-points (Fig. 1) is then found by vertical log-linear interpolation, that is, u(z f ) is found by 1-D piecewise linear interpolation of the values u(z 0 ) = 0, u(z 1 ), u(z 2 ),... at lnz 0 , lnz 1 , lnz 2 ,... to lnz f . If z f ≤ z 0 , we set u(z f ) = 0.
The v component of the wind is interpolated vertically in the same way. Each horizontal wind component u, v is then interpolated separately to the cell centers of the fire subgrid by bilinear interpolation. The fire model then makes one time step: 1. If there are any active ignitions, the level-set function is updated and the ignition times of any newly ignited nodes are set following Sect. 4.4. 13 2. The numerical scheme (11)-(13) for the level set equation (10) is advanced to the next time step.
3. The time of ignition set for any any nodes that were ignited during the time step, from (14).
4. The fuel fraction is updated following Sect. 4.3.
5. The sensible and latent heat flux densities are computed from (5) and (6) in each fire model cell.
6. The resulting heat flux densities are averaged over the fire cells that make up one atmosphere model cell, and inserted into the atmospheric model, which then completes its own time step.
The heat fluxes from the fire are inserted into the atmospheric model as forcing terms in the differential equations of the atmospheric model into a layer above the surface, with assumed exponential decay with altitude. Such scheme is needed because WRF does not support flux boundary conditions. This is code originally due to Clark et al. (1996a,b) and it was rewritten for WRF variables in Patton and Coen (2004). The sensible heat flux is inserted as an additional source term to the equation for the potential temperature θ, equal to the vertical divergence of the heat flux, where R Θ (Φ) is the component of the tendency in the atmospheric model equations (22), σ is the specific heat of the air, (x,y,z) is the density, and z ext is the heat extinction depth, given as parameter fire ext grnd in namelist.input. The latent heat flux is inserted similarly into the tendency of the vapor concentration q m by where L is the specific latent heat of the air. Cf. Clark et al. (1996a, Eqs. 10, 12, 13, 18).
7 Software structure 7.1 Parallel structure Parallel computing imposes a significant constraint on user programming technique. WRF parallel infracture (Michalakes, 2000) divides the domain horizontally into patches. Each patch executes in a separate MPI process and it may be further divided into tiles, which execute in separate OpenMP threads (Fig. 3). Communication between the tiles is accomplished by exiting the OpenMP parallel loop over the tiles. The fire grid tiles are colocated with the atmospheric grid tiles (Fig. 1). The patches are declared in memory with larger bounds than the patch size, and communication between the patches is accomplished by HALO calls (actually, includes of generated code), which update a layer of array entries beyond the patch boundary from other patches. The fire module computational code itself is designed to be tile-callable as required by the WRF coding conventions (WRF Working Group 2, 2007). Tile-callable code updates array values on a single tile, assuming that it can safely read data from a layer of several array entries beyond the tile boundary. The communication (OpenMP loops or HALO calls) happens outside; this means that every time when communication is needed, tile-callable code must exit, and then the next stage can resume on the next call (Fig. 4). The fire module code executes in 6 stages interleaved with communication, 3 stages for initialization and 3 stages in every time step.

Software layers
The fire module software is organized in several isolated layers (Fig. 5). The driver layer contains all exchange of data between the tiles in parallel execution. The rest of the code is tile-callable. The driver layer calls the interpolation and other coupling between the fire and the atmospheric grids, and the fire code itself. The atmospheric physics layer mediates the insertion of the fire fluxes into the atmosphere, as described in Sect. 2. Only these two layers depend on WRF; the rest of the fire module can be used as a standalone code, independent on WRF. The utility layer contains interpolation and other service code, such as stubs to control access to WRF infrastructure, so that WRF calls can be easily emulated in the standalone code. The model layer is the entry point to the fire module. The core layer is the engine of the fire model, described in Sect. 4. The fire physics layer evaluates the fire spread rate and heat fluxes from fuel properties. One of the goals of the design is that the only components that will need to be modified when the fire module is connected to another atmospheric model in future are the driver layer, the atmospheric physics layer, and the WRF stubs in the utility layer.
8 Recommended WRF settings

Domains and nesting
WRF-Fire may be run in both "ideal" and "real" modes, which require slightly different setups. In both cases, the model requires a set of data defining model initialization (wrfinput). In the real cases, boundary conditions in a form of wrfbdy files must be also provided, and both types of files are created by real.exe preprocessor from the WRF Preprocessing System (WPS). These files contain not only meteorological and topographical data but also fire related information, such as the fuel type map and high-resolution topography on the fire mesh. Since the WRF-Fire initialization for the real cases does not differ from the one for the regular WRF, all physical and dynamical options available in the regular WRF are also available in WRF-Fire. Therefore, the same general rules apply to the configuration of WRF-Fire as to the configuration of the regular WRF. However, one should keep in mind that resolutions of the finest domains in fire simulations are usually significantly higher than in weather forecasting applications. This has two consequences in terms of the proper WRF-Fire setup. First, if the resolution of any of the inner domains is less than 100 m, this domain should be actually resolved in the large eddy simulation (LES) mode, without the boundary layer parameterizations. At this resolution, the model should be able to resolve the most energetic eddies responsible for mixing within the boundary layer, so the boundary layer parameterization in this case is not needed. Second, since in the nested mode, vertical levels are common for all domains, the height of the first model level selected for the most outer (parent) domain, defines also the level of the first model layer for all inner (child) domains, even if their horizontal resolutions are an order of magnitude smaller. The fact that the vertical model resolution is the same for all domains significantly limits the minimum height above the ground of the first model level. This in turn is crucial for the fire model, which uses the wind speed interpolated to 6.1 m above the ground. Therefore, in the cases when the first model level must be relatively high above the ground it is recommended to perform downscaling using the ndown.exe program, being a part of the WRF distribution. In this case the outer domains are run separately without the fire, and then based on the output from this simulation, ndown.exe creates a set of new initial and boundary condition files (wrfinput and wrfbdy) for the separate simulation from the innermost domain(s). This allows for a new setup of vertical levels for the innermost domains, and selecting proper physical options for them.

Large Eddy Simulation and surface properties
To enable the high-resolution simulation in Large Eddy Simulation (LES) mode, user should first disable the boundary layer parameterization (bl pbl physics=0). The LES mode requires the proper surface fluxes in order work properly. We recommend the option isfflx=1, which makes WRF use a surface model to compute the surface fluxes. Other options with constant heat fluxes and drag are not well suited for fire simulations. Out of all surface exchange parameterizations only the classic Monin-Obukhov theory (sf sfclay physi is recommended for the LES cases. This option assures a proper computation of surface transfer coefficients that are used together with the surface properties (provided by the surface model) for computation of the surface fluxes of the momentum, heat and moisture. The surface model itself computes properties of the surface, but does not compute the surface exchange coefficients, which are needed for computation of the surface fluxes. Hence, in order to compute them, the surface properties must be provided by a surface model, which is enabled by choosing a non-zero sf surface physics. The subgrid scale parameterization used by the WRF in LES mode is defined by the km opt parameter, which should be set to 2 (TKE closure), or 1 (Smagorinsky scheme). In real cases, real.exe automatically provides proper initialization for the selected land surface model, and all other components. In idealized cases, users have an option of the basic surface initialization, intended to be used without the surface model, or the full surface initialization (sfc full init=1). One should keep in mind that without the full surface initialization, there is no direct way to define surface properties such as temperature or roughness. For idealized cases with the full surface initialization, the surface scheme utilizes a table containing records of land-use categories and corresponding surface properties like roughness length, heat capacity, etc. All these properties are defined in a text file LANDUSE.TBL, which may be edited by the user. Therefore, setting up the land-use category is enough to provide all static surface properties. The basic parameters required by this surface model like land use index, surface air temperature and soil temperature, may be defined directly in namelist.input by the variables sfc lu index, sfc tsk, and sfc tmn if they are intended to be the same over the whole domain. If they are not spatially uniform, they may be read in from external files if fire read lu, fire read tsk, or fire read tmn are set to true. For details about the input data for real cases, see Sect. 9.

Fire subgrid refinement ratios
The fire mesh needs to be about 10 times finer than the atmospheric mesh to allow for gradual heat release into the atmosphere, even if fuel and topography data may not be available at such fine resolution. The fire mesh refinement in the x and y direction (sr x and sr y) must be defined in the domain section of namelist.input. Since these refinement factors define dimensions of the fire-related variables, they must be selected before execution of real.exe, which generates the WRF input files. Any change in atmospheric to fire grid ratios requires re-running real.exe and creating new input files. The atmospheric mesh step should be about 60 m or less for proper feedback of the wind on the fire line; larger mesh step can result in incorrect fire spread rates and atmospheric behavior (Clark et al., 1996a, p. 887).

Time step
In real WRF-Fire simulations performed in multi-domain configurations the time step requirements for the outer domains (run without fire) do not differ from general meteorological cases. The recommended time step of 6 times the horizontal grid spacing (in km) may be used as a starting point. However, for the finest domains run with fire simulations, the time step in most cases must be significantly smaller. For domains with low vertical resolution and simple topography, the horizontal mesh step is crucial for numerical stability, since the horizontal velocity is greater than the vertical one. In fire simulations with high vertical resolution, the vertical velocity induced by fire may violate the CFL condition. Therefore, it is advisable to use a vertically stretched grid, with finer resolution at the surface (where updraft velocities are not that high) and lower resolution at higher levels where stronger updrafts are expected. This allows for having the first model level relatively close to the ground, yet with vertical spacing aloft big enough to handle strong convective updrafts without violating the CFL condition.
In real cases, the pressure levels may be defined directly in the namelist.input file. In ideal WRF-Fire runs, there is now an option stretch hyp, which turns on hyperbolic grid stretching. The grid refinement may be adjusted using the z grd scale namelist variable. One should keep in mind that running the WRF-Fire simulations with high-resolution topography in most cases limits the maximum numerically stable time step. Steep terrain often induces high vertical velocities that may violate the CFL condition. Therefore, these cases usually require significantly smaller time steps than similar simulations run with low-resolution, smooth topography.

Data input
A WRF ideal run is used for simulations on artificial data. An additional executable, ideal.exe, is run first to create the WRF input. A different ideal.exe is built for each ideal case, and the user is expected to modify the source of such ideal case to run custom experiments. The ideal run for fire supports optional input of gridded arrays for land properties, such as terrain height, roughness height, and terrain height. This allows to run simulations which go beyond what would normally be considered an ideal run and simplifies custom data input; the simulation of the FireFlux experiment (Sect. 10) was done in this way.
A WRF real run is used for prediction and analysis of natural events. For a real run, a user must supply data for the initial and boundary conditions for the WRF simulation. The WRF Preprocessing System (WPS) (Wang et al., 2010, Chapter 3) contains a number of utilities useful for preparing standard atmospheric and surface datasets for input into WRF.
1. Geogrid creates the surface mesh from a specified geographic projection and interpolates static surface data onto the mesh. It supports several interpolation methods as well as data smoothing and creation of gradient fields. Geogrid reads data in a tiled binary format described by a text file and writes to a NetCDF file for each nested mesh. All data required for atmospheric simulations up to 30 arcseconds resolution globally are provided by NCAR.
2. Ungrib extracts atmospheric data from standard GRIB files and writes to a simple binary format. Ungrib does not do any interpolation; it only searches through a number of files for necessary variables within the time window of the simulation. Data for ungrib must be obtained by the user. Several free sources of atmospheric GRIB data are available online from production weather simulation.
3. Metgrid reads the output from geogrid and ungrib and produces a series of NetCDF files read by WRF's real.exe binary. The geogrid output is copied directly into each of these files, while the ungrib output is interpolated horizontally on to the computational mesh.
The metgrid files produced by WPS are portable and relatively compact so they can be transferred to a computer cluster for the simulation's execution. From this point, the real.exe program in WRF handles the vertical interpolation of atmospheric fields and all processing for the creation of WRF's initial (wrfinput) and boundary (wrfbdy) files.
WPS has been extended with the ability to produce data defined on the refined surface meshes used by WRF-Fire (Sect. 8); however, it is not possible to distribute high resolution, global fields as is done in the standard dataset. Instead, the user must download any necessary high resolution fields and convert them into geogrid's binary format for each simulation. WRF-Fire is distributed with an additional utility, convert geotiff.x, which can perform this conversion from any GeoTIFF file. This utility is written entirely in C and depends only on the GeoTIFF library.
For a WRF-Fire simulation, it is only strictly necessary to download one additional dataset for input into geogrid. This dataset contains fuel behavior categories and is stored 19 in the variable NFUEL CAT. For simulations within the United States, this data can be obtained in GeoTIFF format from the USGS at http://www.landfire.gov. WRF-Fire uses an additional variable for topography, ZSF, which is allowed to be different from the topography used used by the atmospheric code defined by HGT. This is useful because a high resolution WRF simulation generally requires the topography to be highly smoothed in preprocessing for numerical stability. The fire code can benefit from a rougher topography for more accurate fire spread computations.
Once the static data is converted into the geogrid binary format, the GEOGRID.TBL should be edited to inform geogrid of the location of each supplementary dataset. WRF-Fire expects two variables to be created on the refined subgrid (NFUEL CAT and ZSF), this is indicated by the line subgrid=yes; all other variables will be defined on the standard atmospheric grid.
For atmospheric data, it is best to use the highest resolution dataset available to initialize a WRF-Fire simulation to capture as much of the local conditions near the fire as possible. Generally, publicly available atmospheric data is limited to around 10 km resolution. As a consequence, one should create several nested grids, each with a 3 to 1 refinement ratio, and a long spin-up prior to ignition in order to recreate local conditions. Preliminary results indicate that assimilation of data from weather stations or satellite radiances may be required for an accurate simulation .
10 Computational simulations Kim (2011) has verified that the level-set method in the fire module advects the fire shape correctly, on some of the same examples that were used to verify the tracer code in CAWFE (Clark et al., 2004).
A number of successful simulations with WRF-Fire now exist. Jenkins et al. (2010) have demonstrated fireline fingering behavior for a sufficiently long fireline (Figs. 8, 9) on an ideal example, with similar results as in Clark et al. (1996a,b). Kochanski et al. (2010) have demonstrated the validity of WRF-Fire on a simulation of the Clements et al. (2007) FireFlux grass fire experiment and obtained good agreement with data (Figs. 6, 7). Dobrinkova et al. (2010) simulated a fire in Bulgarian mountains using real meteorological and geographical data, and ideal fuel data. Beezley et al. (2010) simulated the 2010 Meadow Creek fire in Colorado mountains using real data from online sources. Topography (Fig. 10) at up to 3 m horizontal resolution was obtained from the National Elevation Dataset (NED, http://ned.usgs.gov) and fire fuel datasets from Landfire (http://landfire.cr.usgs.gov) at up to 10 m resolution. Six nested domains were required to scale the simulation down from the atmospheric initialization (32 km) to the fire grid 20 resolution (10 m). Cloud physics was enabled in domains 1-3. The fire subgrid refinement ratio was 10 times on the finest domain to capture fire surface variables and for a gradual release of the heat flux near the fireline. Realistic fire and atmosphere behavior was obtained (Figs. 11,12).

Discussion
11.1 Additional features WRF-Fire does not yet support canopy fire, although canopy fire colocated with ground fire is contained in CAWFE. The reason was the desire to keep the code as simple as possible early on and add features only as they can be verified and validated. The support for canopy fire will be added in future. Adding smoke from the fire to WRF is also under consideration. A list of desired features and a record of the progress of the development are maintained at http://www.openwfm.org/wiki/WRF-Fire development notes.

Atmosphere
Rothermel's spread model (1) assumes wind as if the fire was not there. In practice, the wind was measured away from the fire. In a coupled model, however, the feedback on the fire is from the wind that is influenced by the fire. Clark et al. (2004) noted that the horizontal wind right above the fireline may even be zero, and proposed to take the wind from a specified distance behind the fireline. Also, the strong heat flux from fire disturbs the logarithmic wind profile, and the rate of spread as a function of wind at a specific altitude may not be a good approximation; rather, the fire spread may depend more strongly on the complete wind profile  and on turbulence (Sun et al., 2009).
The assumption of horizontal homogeneity in the Monin-Obukhov similarity theory is not satisfied here; the horizontal dimension of the active part of fire is not orders of magnitude larger than the boundary layer height as required, and it may be in fact smaller. Another indication that the Monin-Obukhov theory may not apply for fires is a strong drop in the heat transfer in the case of strong temperature gradients, shown in our preliminary tests.
Horizontal wind could be interpolated vertically to different heights for different fuels like in CAWFE model, which takes the wind from different mesh levels for different fuels. However, here we follow a classical approach of Rothermel (1972) and Baughman and Albini (1980), where the wind speed is evaluated at the common 6.1m height, and then converted to the mid-flame height using the fuel-specific wind correction factors.
Very strong vertical components of the wind caused by the fire result in the need for short time steps to avoid violation of the vertical CFL condition (Sect. 8.4). It would be interesting to couple the fire module also with the Non-hydrostatic Mesoscale Model (NMM) core 21 of WRF, which is implicit in the vertical direction (Janjic et al., 2005), and it may perform better in the presence of strong convection (Litta and Mohanty, 2008). The ARW core is semi-implicit in the vertical direction in the vertical wind component and the geopotential.

Fire
The more recent Scott and Burgan (2005) fuel categories are more detailed than Anderson (1982) categories, they are supported by BehavePlus, and fuel maps using them are available from Landfire. But instead of describing additional categories in namelist.fire, it may be more useful to support the import of fuel files from BehavePlus, which is also well suited for editing and diagnosing fuel models. More accurate fuel models (Albini et al., 1995;Clark et al., 1996a), including those in BehavePlus, consider fuels to be mixtures of components with different burn times, which results in a different heat release curve. While the spread rate of established fire in the simulation of the FireFlux experiment was reasonably close, the simulated fire still arrived at the observation towers too soon , because it started too quickly. A better parametrization of the ignition process seems to be in order. The fire spread in the Meadow Creek fire simulation was also too fast, but for a different reason. It is well known that the actual spread rates of wildland fires tend to be lower than the spread rates in simulations, which are derived from laboratory experiments. This effect might be attributed to irregularities on scales not captured by the simulation (Finney, 1998, p. 34), including granularity of the fuel supply not reflected in the data. Refining the semi-empirical model from detailed numerical simulations and parametrizing complex fire behavior are suggested important research areas.
The computation of the heat fluxes in (5) and (6) does not take into account the evaporation of moisture present in the fuel, only the production of water by burning of hydrocarbons. This error is typically just few %, however, which is small in comparison with other uncertainties. The fuel models should be dynamic (with variable fuel moisture) as in BehavePlus. Coen (2005) added an explicit diurnal cycle for the moisture into CAWFE. Here, moisture content could be coupled with existing WRF land surface models, which could take into account air humidity and precipitation. The radiative and convective parts of the sensible heat flux should be treated differently. The release of surface heat and moisture into the atmosphere are already present in WRF soil models. Their scale, however, is different from the powerful heat release from a fire.

Numerical methods
In a numerical implementation, the level-set method is global, unlike tracers, which move locally. In spite of the fact that the level-set equation determines the fire spread locally from the spread rate at the fireline, the behavior of the fireline depends slightly on the wind, the fuel, and the level set function in certain other locations from previous time steps, because of the discretization errors and the artificial diffusion. This nonlocal behavior has not been practically significant, however.
The fuel fraction calculation (19) can have significant error in the fire subgrid cells near the fireline, which will to some degree average out over the atmospheric mesh cells. Rigorous error analysis will be done elsewhere. We are currently testing an alternative method which is always second order in the sense that it is exact when the time from ignition and the level-set function are linear in space. The alternative method is more computationally expensive, but, on the other hand, it might allow to decrease the subgrid refinement ratio; with large meshes, it is possible to run against 32 bit integer limits.

Data assimilation
Data assimilation for wildland fires is an area of great interest. Methodologies for a reactiondiffusion model were proposed based on the ensemble Kalman filter (EnKF) and the particle filter (Mandel et al., 2004). Unfortunately, statistical perturbations can cause spurious fires, which do not dissipate. Combination of the EnKF with Tikhonov regularization alleviates the problem somewhat (Johns and Mandel, 2008;Mandel et al., 2009), but the resulting method is still not robust enough. A new method, called morphing EnKF and based on combined amplitude and displacement correction (Beezley and Mandel, 2008), was shown to work with WRF-Fire (Mandel et al., 2009), and it is under continued development (Mandel et al., , 2011. We are not aware of any work elsewhere on data assimilation for a coupled fire-atmosphere model. Particle filters were proposed for discrete cell-based fire models (Bianchini et al., 2006;Gu et al., 2009), using fitness functions involving the area burned rather than intensities of physical variables.
Starting the model from a known fire perimeter is important for many potential users. This can be understood as a data assimilation problem, but we are considering a simpler method for this particular case: prescribe the fire history up to the time of the given perimeter to allow the atmospheric conditions to evolve, then allow the coupled model take over. Tools to produce such artificial fire history are being developed. Possibly the simplest alternative is an interpolation from a given ignition point and time to the given perimeter. A more complex version would run the fire model (without atmosphere) backwards in time and attempt to find the ignition point automatically. The latter approach could be also interesting for forensic purposes.

23
We have described the coupled atmosphere-fire model WRF-Fire. The software is publicly available and it supports both ideal and real runs. Visualization and diagnostic utilities are available. Currently, the model is suitable for research and education purposes. Validation is in progress. Table 1. Fuel properties. The notation is from Rothermel (1972) except as indicated. The identifiers are as used in WRF-Fire and CAWFE. In the input files, some quantities are given in English units per Rothermel (1972); see Wang et al. (2010, p. A-5).
symbol description identifier a wind adjustment factor (Baughman and Albini, 1980) windrf from 6.1 m to midflame length w fuel weight (i.e., burn time) (s) 40% decrease of fuel in 10 min for w = 1000 weight w total fuel load (kg m −2 ) fgi δ m fuel depth (m) fueldepthm σ fuel particle surface-area-to-volume ratio (1/m) savr M x moisture content of extinction (1) fuelmce ρ P ovendry fuel particle density (kg m −3 ) fueldens S T fuel particle total mineral content (1) st S E fuel particle effective mineral content (1) se h fuel heat contents of dry fuel (J kg −1 ) cmbcnst M f fuel particle moisture content (1) fuelmc g Table 2. Computation of the fire spread rate factors in (1) from the fuel properties (Table 1), the wind speed U at 6.1 m, and the terrain slope tanφ. All equations are from Rothermel (1972) unless otherwise indicated. All input quantities are first converted from metric to English units (BTU-lb-ft-min) to avoid changing the numerous constants in the Rothermel (1972) computations. Further, following CAWFE, the wind is limited to between 0 and 30 m s −1 and the slope is limited to nonnegative values.   (16) over the cell C with the center a 3 is computed as the sum of integrals over the subcells C 1 ,...,C 4 . While the values of ψ and t i are known at a 3 = x 3 , they need to be interpolated to the remaining corners x 1 , x 2 , x 4 of the subcell C 1 from their values at the points a 1 ,...,a 4 . The computational domain is divided into disjoint rectangular patches. Each patch is updated by a single MPI process (distributed memory parallelism), and the process may read arary data in a strip around the patch, called halo region. The communication between the patches is by halo calls to the RSL parallel infrastructure (Michalakes, 2000), which update the halo regions by the values from the neighboring patches. Each patch may be divided into tiles, which execute in separate OpenMP threads (shared memory parallelism). Following WRF coding conventions (WRF Working Group 2, 2007), computational kernels execute in a single tile. They may read array values from a strip beyond the tile boundary but no explicit communication is allowed. 3-D arrays are divided into patches and tiles in the horizontal plane, cf., Fig. 1.