Effects of High-Density Gradients on Wildland Fire Behavior in Coupled Atmosphere-Fire Simulations

Abstract

Coupled atmosphere-fire modeling is recognized as a relevant approach for the representation of the interaction between a wildland fire and local meteorology at landscape scales. The atmospheric model component used in the coupled system is based on several approximations, which are adopted for computational efficiency or physical processes representation, including the widely used anelastic approximation. The validity domain of the anelastic approximation may be questioned in the context of high-resolution wildland fire modeling due to the large fire-induced heat releases near the surface. This study aims to study this question with the MesoNH anelastic model coupled with the Blaze fire model. A compressible version of the MesoNH-Blaze coupled model has been developed for comparison with the anelastic system. The FireFlux I experimental fire is used for this comparative study conducted at a 10-m and a 25-m horizontal atmospheric resolution. Results show significant anelastic/compressible differences at a 10-m resolution on the physical processes occurring near the fire with higher horizontal velocities and the presence of gravity waves downstream of the fire. This is in addition to the fire plume with realistic larger vertical velocities. Differences at a 25-m resolution are much smaller in all evaluated processes. The compressible system only enriches the physics underlying fire-atmosphere interactions at a very high resolution, which means that the anelastic approximation remains relevant for large-scale coupled atmosphere-fire simulations, considering the significant economy concerning numerical costs.

Introduction

Wildfires result from complex physical, chemical and biological processes interacting over a wide range of spatial and temporal scales (Gollner et al., 2015). Understanding the fundamental processes driving wildfire behavior is a key point for predicting fire spread across the landscape and the induced atmospheric dynamics, as well as for the anticipation of the human and environmental impacts of extreme wildfire events (Tedim et al., 2018). During a wildfire event, active flaming areas release very large amounts of heat into the atmosphere, which induces the development of a buoyant plume, creating air entrainment effects (fire-induced flow) toward the active flaming areas that can enhance fire propagation and more generally modify wildfire behavior. They are subject to significant temperature and air density gradients and these are important to capture so as to correctly predict the interactions between a wildland fire and the atmosphere.

Coupled atmosphere-fire models have been developed in order to represent the fire-induced flow at landscape-to-meteorological scales, among which MesoNH-ForeFire (Filippi et al., 2009, 2013, 2018), WRF-SFIRE (Kochanski et al., 2013; Mandel et al., 2011), WRF-FIRE (Coen et al., 2013; Mandel et al., 2009; Muñoz-Esparza et al., 2018), and more recently MesoNH-Blaze (Costes et al., 2021). These coupled models aim to solve the flow dynamics near active flaming areas without representing the physical processes at the flame scale. They represent the fire as a moving front separating the burnt area from the unburnt vegetation using a parameterization of the front speed, or rate of spread (ROS), which takes near-surface wind conditions, biomass fuel moisture content and properties, and terrain topography as inputs (Sullivan, 2009). The atmospheric model imposes the surface wind field as input to the fire spread model, while the fire spread model imposes latent and sensible heat fluxes as surface boundary conditions to the atmospheric model. This two-way coupling provides access to the temporal and spatial variability of the near-surface wind conditions which are perturbed by the development of the wildland fire.

Coupled atmosphere-fire models are subject to several uncertainties, for instance uncertainties in the fire model inputs (Allaire et al., 2020) and model uncertainties associated with the simplified physics of fire parameterizations (Cruz et al., 2018). A recent study has shown the relative importance of the atmospheric variability on the simulated wildland fire behavior and the need to move to an ensemble modeling approach (Costes et al., 2021). It is essential to explore to which extent atmospheric model uncertainties affect the simulated wildland fire behavior. This present work focuses on the uncertainties associated with the atmospheric model formulation, that is, on the choice of the atmospheric model equations to solve in order to properly capture the atmosphere-fire interactions in the vicinity of a wildland fire.

In atmospheric models, it is well known that the fast propagation of acoustic waves, without meteorological interest, sets a very strong computational constraint on the time step in the Euler equation system. Mainly, three different approaches are used to avoid this problem: (a) the use of two different time steps for acoustic and non-acoustic processes; (b) the implicit treatment of certain terms, which usually involves a linearization around a mean atmospheric state; and (c) the “anelastic” approximation. In the third approach, acoustic waves are eliminated from the equations through the use of a constant density profile instead of the actual fluid density in the continuity equation and in the momentum equation, all this not counting the buoyancy term. The fluid is therefore formally considered as pseudo-incompressible, and the pressure is deduced from the solution of an elliptic equation. Some atmospheric models (Lafore et al., 1998; Prusa et al., 2008; Smolarkiewicz & Margolin, 1997) use the anelastic approximation with different formulations (e.g., Durran, 1989; Lipps & Hemler, 1982).

As with any modeling hypothesis, the anelastic approximation is subject to a validity domain. This refers to the modeled processes or to the size of the computational domain for example, Most studies on the validity domain of the anelastic system focus on strong convection cases (Kurowski et al., 2014, 2016; Smith & Bannon, 2008) or horizontally very large domains (Davies et al., 2003; Kurowski et al., 2015). In the framework of coupled atmosphere-fire simulations at decametric resolution, the validity of the anelastic approximation is also questionable. First, the assumption of small perturbations for state variables as potential temperatures (Bannon, 1996; Durran, 1989) can prove inaccurate in the vicinity of the fire front, where temperature fluctuations can reach the same order of magnitude as the reference state. Second, temperature gradients can be up to four orders of magnitude higher near the fire front, which can be a challenge for the anelastic system (Klein et al., 2010). Third, the strong wind divergence constraint of the anelastic system can affect the fire-induced wind and therefore the coupled model ability to represent wildland fire dynamics. It is therefore of primary importance to assess the validity of the anelastic approximation to simulate landscape-scale wildland fire behavior.

A fully compressible atmospheric model coupled to a fire model is required as a comparison for the anelastic version of the coupled model. Several non-hydrostatic atmospheric models have already solved the compressible version of the Navier-Stokes equations, among which (Japan Meteorological Agency Non-Hydrostatic Model, Saito et al., 2006), WRF (Weather Research and Forecasting model, Skamarock et al., 2008), ECMWF IFS (Kühnlein et al., 2019), and (Application of the Research to Operations at MEsoscale, Seity et al., 2011). It is worth mentioning that only WRF is currently used in the context of wildland fire behavior modeling (Mandel et al., 2009, 2011). In the present work, a fully compressible version of the MesoNH atmospheric model has been developed and coupled with the Blaze fire model in a similar way to the anelastic version used in Costes et al. (2021). A comparison of the anelastic and compressible versions of MesoNH-Blaze is carried out at decametric atmospheric resolution for the FireFlux I experimental fire (Clements et al., 2006, 2007). This comparison is of great interest first for the validation of the compressible version implementation and second for the evaluation of the possible benefits of the compressible version on the simulated wildland fire behavior.

The outline of this paper is as follows. The MesoNH atmospheric model is presented in Section 2, with details on the equations of the compressible version. The Blaze fire model is presented in Section 3. Section 4 presents the FireFlux I experiment (referred to as FireFlux in the following) and the associated numerical settings of the MesoNH-Blaze coupled model. A comparison of the anelastic and compressible results with respect to the observations is provided in Section 5. A more detailed analysis of the physical processes reproduced with the anelastic and compressible versions of MesoNH-Blaze at two horizontal resolutions (10 and 25 m) is provided in Section 6. The conclusion and outlook are presented in Section 7.

The MesoNH Atmospheric Model

MesoNH (Lac et al., 2018; Lafore et al., 1998) is a non-hydrostatic anelastic atmospheric model developed by both CNRM (Météo-France/CNRS) and Laboratoire d’Aérologie (Université de Toulouse/CNRS). It is used to simulate meso-scale (a few kilometers to less than a thousand kilometers) up to micro-scale (less than addeda hundred meters) atmospheric flows (Aumond et al., 2013; Bergot et al., 2015; Filippi et al., 2018; Auguste et al., 2019, among others).

Standard Anelastic Version of MesoNH

The Durran (1989) pseudo-incompressible system, referred to as anelastic in the following, is the sound-proof formulation adopted in MesoNH and will be briefly presented hereafter. It assumes a constant dry density profile called the dry reference density such as ρ_d(x, y, z, t) ≡ ρ_d,ref(z). Thus, the continuity equation may be written as: [Image Omitted. See PDF]

The atmosphere is formally incompressible and an elliptic equation is required to diagnose the pressure. This continuity equation formulation sets a strong constraint on the wind field. The momentum conservation equation for the dry air fraction is: [Image Omitted. See PDF]with, from left to right, the Eulerian time variation and the momentum advection, the pressure gradient force $${\boldsymbol{\mathcal{F}}}_{{\Pi}}$$, the buoyancy force $${\boldsymbol{\mathcal{F}}}_{B}$$, the Coriolis force and the other sources of momentum $$\boldsymbol{\mathcal{F}}$$ associated with physical processes such as subgrid-scale transport (i.e., turbulence, convection, vegetation drag). The only term in the anelastic system where the density variations are considered is the buoyancy term $${\boldsymbol{\mathcal{F}}}_{B}$$. These density variations are accounted for using the virtual potential temperature perturbation: [Image Omitted. See PDF]

For simplicity purpose, hydrometeor mass and volume are neglected in this model presentation. The only non-zero mixing ratio is the water vapor mixing ratio r_v, defined by the ratio of the mass of water vapor to the mass of dry air in a volume of air. It means that no phase change is considered in the following equations. The thermodynamic equation can be written as: [Image Omitted. See PDF]with, from left to right, the Eulerian time evolution and advection of the energy variable, the moist correction in the absence of phase change, as well as the radiative, diffusive and fire sensible heat flux effects contained in $$\mathcal{H}$$. The state equation used in the Durran (1989) system is not linearized: [Image Omitted. See PDF]where the total density ρ is retrieved from: [Image Omitted. See PDF]

New Compressible Version of MesoNH

Moving the atmospheric dynamic equations to a compressible formulation requires fundamental changes in the dynamic kernel of MesoNH since air density is no longer constant in time and uniform in space.

Because this study only aims to focus on the relevance of using compressible equations in order to evaluate the role of density gradients for wildland fire spread, no effort was undertaken here on CPU efficiency. Further numerical developments would be necessary to optimize the time-step for general research atmospheric applications. So hereafter only the changes in the physical equations are described.

Prognostic Variables

A new prognostic variable for air density is added to the set of equations. For MesoNH, dry air density and potential temperature have been chosen respectively as mass and energy variables. First, dry air density is selected over pressure because it is a conservative variable. Pressure is then determined using the state equation, which remains the same as in the anelastic system (Equation 5). Second, potential temperature is chosen as the energy variable over temperature for its conservative property in adiabatic processes. Moreover, as potential temperature is already used in the anelastic system, it is natural to rely on the same energy variable in the compressible system.

Mass Conservation

Dry air density ρ_d (x, y, z, t) is now used as a prognostic variable in the continuity equation: [Image Omitted. See PDF]

Dry air density is preferred to total density ρ in order to avoid source terms, stemming from phase change for example, Assuming dry air pressure is close to total pressure, total density can be retrieved from dry air density as follows: [Image Omitted. See PDF]

The equation used to describe the conservation of water vapor mixing ratio is derived from the anelastic system with replacing dry reference density by dry prognostic density: [Image Omitted. See PDF]with $${\mathcal{Q}}_{v}$$ the volume source term accounting for subgrid-scale transport and fire latent heat flux.

Momentum Conservation

The reference density used in the anelastic version of the momentum conservation equation (Equation 2) is replaced in the compressible version by dry density. The terms corresponding to the rapidly evolving terms related to acoustic waves, such as the pressure term $${\boldsymbol{\mathcal{F}}}_{{\Pi}}$$ and the buoyancy term $${\boldsymbol{\mathcal{F}}}_{B}$$, are modified. To maintain consistency between the model equations, the pressure gradient term is written using the Exner function: [Image Omitted. See PDF]

The buoyancy term is also directly calculated since density is now a model prognostic variable. To remain consistent with the pressure term formulation (Equation 10), the total density ρ is replaced by an Exner function Π derived from the state equation (Equation 5): [Image Omitted. See PDF]

This formulation allows the strict compensation of the effects of pressure and buoyancy in the framework of the hydrostatic equilibrium, which is written in the compressible version as follows: [Image Omitted. See PDF]

The anelastic set of equations imposes a strong constraint on wind divergence, which no longer exists in the compressible version. As a consequence, it is necessary to add a viscosity term $${\boldsymbol{\mathcal{F}}}_{\mu }$$ in order to efficiently diffuse the acoustic waves produced within the computational domain. For this purpose, a very high dynamic viscosity is imposed, the objective of which is to decrease the rapid production of momentum related to acoustic waves: [Image Omitted. See PDF]where μ is the dynamic viscosity and λ is the volume viscosity. Based on these elements, momentum conservation translates into the following equation: [Image Omitted. See PDF]

Energy Conservation

Some simplifying assumptions used in the anelastic approximation for energy conservation no longer apply in the more general compressible framework. In particular, pressure variations are not limited to the vertical direction in the moisture correction term in the absence of phase changes. Energy conservation translates into the following equation: [Image Omitted. See PDF]

This new equation features a temporal variation in pressure, which requires a two-stage numerical resolution method (Appendix A).

It is worth mentioning that before considering the wildland fire behavior application, the compressible system of equations has been validated against two canonical test cases (Appendix B).

The Blaze Fire Model

Blaze is a fire model embedded in the MesoNH code and described in detail in Costes et al. (2021). The purpose of Blaze is to represent the fire spread and its effects on the atmospheric flow. For this purpose, Blaze is composed of two main components: (a) a fire front propagation model, and (b) a surface heat flux model. At each time step, the near-surface wind computed in MesoNH affects the rate of fire spread calculated in Blaze. Blaze then computes the fire heat fluxes that are given as source terms in the energy and water vapor balance equations of MesoNH.

In BLAZE, the fire spread model is based on an Eulerian level-set front-tracking numerical approach where the ROS is estimated using Balbi's parameterization (Balbi et al., 2009; Santoni et al., 2011). The mesh used by Blaze corresponds to a refined surface grid of MesoNH. The Hamilton-Jacobi equation is solved by using third-order Runge-Kutta and third-order weighted essentially non oscillatory (WENO) schemes. Heat flux parameterizations are based on the arrival time and residence time of the fire. Sensible and latent heat fluxes decay exponentially during the flaming time period. Smoldering sensible heat flux can be accounted for as a constant flux after the flaming time period.

MesoNH and Blaze can be coupled either in a one-way or a two-way mode in order to analyze the response of the different components of the coupled model in a separate way or when fully coupled. In a one-way mode, the fire-induced heat flux feedback to the atmosphere can be deactivated to focus solely on the fire spread forced by the atmospheric model. A fire replay mode is also available for forcing the atmosphere with an arrival time map that can be derived from Blaze or from observations. This study focuses on the two-way coupled mode, where the fire spread and the heat fluxes are computed at each time step, in order to study both the fire spread and the fire-induced wind.

Application to the FireFlux I Experimental Fire

The main objective of the present study is to provide a detailed comparison of the anelastic and compressible versions of MesoNH-Blaze for the FireFlux I experimental fire. The first objective is to assess the behavior of the compressible version of MesoNH for a fire-induced flow. Measurements of the FireFlux I experiment are used for that purpose and detailed analysis is presented in Section 5. Once the confidence level in the compressible version is sufficient, a phenomenological analysis is conducted in Section 6 so as to compare the anelastic and compressible versions of MesoNH more deeply.

This section provides a description of the FireFlux I experiment and the numerical configuration of the coupled MesoNH-Blaze model used in Sections 5 and 6.

Experiment Summary

The FireFlux I experimental fire conducted on February 2006 in Texas corresponds to a tall grass burn plot of about 30 ha (Clements et al., 2007). This experimental fire aimed to document the near-surface fire-induced flow. For this purpose, it was instrumented with two towers equipped with thermocouples and anemometers at different heights. This experimental data set was already used to validate WRF-SFIRE (Kochanski et al., 2013), MesoNH-ForeFire (Filippi et al., 2013), and MesoNH-Blaze (Costes et al., 2021). Observations at the instrumented towers were used here to validate the compressible version of MesoNH-Blaze. A detailed analysis of the fire-induced flow near the surface and of the fire spread was carried out in areas where the density gradients were the strongest so as to evaluate their possible effects on the simulated wildland fire behavior.

An overview of the experimental field is shown in Figure 1. The ignition point is on the North side of the plot. Two firefighters carried out the ignition process by simultaneously lighting the western ignition line (170 m lit in 153 s) and the eastern ignition line (215 m lit in 163 s). The total duration of the burn was about 15 min. During the fire, the wind was blowing from the north-east with a limited magnitude (below 10 m s⁻¹) so the southward fire spread was at a moderate speed (1.6 m s⁻¹). The two measurement towers (main and small towers) were 43 and 10 m high, respectively. Experimental observations of turbulence and thermodynamics were reported in Clements et al. (2006, 2007, 2008).

Figure 1. Schematic of the FireFlux experimental field. The red dot represents the ignition point. The blue dots represent the main and small tower positions. The green area corresponds to the burn lot up until the small tower. Reprinted from: Subgrid-scale fire front reconstruction for ensemble coupled atmosphere-fire simulations of the FireFlux I experiment, Aurélien Costes, Mélanie C. Rochoux, Christine Lac, Valéry Masson, Fire Safety Journal, Copyright (2021), with permission from Elsevier [OR APPLICABLE SOCIETY COPYRIGHT OWNER].

Numerical Configuration

The computational domain of MesoNH-Blaze is 4 km by 7.5 km in the horizontal directions, and 1 km in the vertical direction. The main study is conducted at a 10-m horizontal resolution and a 4-m vertical resolution (within the first 50 m) for the atmospheric model, and at a 5-m resolution for the fire model. The atmospheric model resolution is selected so as to be consistent with previous studies (Costes et al., 2021; Filippi et al., 2013; Kochanski et al., 2013). The choice of the fire model resolution results from a mesh convergence analysis conducted in Costes et al. (2021). The Blaze fire model is used with Balbi's spread-rate parameterization (Santoni et al., 2011) and the exponential with smoldering heat flux parameterization. In order to capture high gradients, the numerical schemes within MesoNH are the fifth-order WENO (WENO5) in space and the five-stage third-order explicit Runge-Kutta scheme (RK53) in time for the wind, and the definite positive piecewise parabolic method (PPM) for scalars.

An ensemble of fifteen MesoNH-Blaze simulations is considered for both anelastic and compressible versions to account for the impact of atmospheric variability on the simulated fire propagation. Each ensemble member corresponds to a different ignition time and thereby a different turbulence initialization.

These numerical settings are identical to those described in Costes et al. (2021) except for the following two points:

1. The ignition temperature T_i used in the Blaze fire model is modified here to better match the observed ROS for both anelastic and compressible ensembles at 10-m horizontal resolution, and in order to have the same fire parameters to carry out the anelastic/compressible comparison. T_i is changed from 590 to 610 K.

2. The near-surface winds supplied by MesoNH to Blaze are temporally filtered for the compressible version only. The strong and rapid temporal variations of the near-surface winds are indeed inconsistent with the assumptions underlying the spread-rate parameterization. For this purpose, a temporal filtering based on the exponentially weighted moving average method, with a time constant of 20 s, is applied to the near-surface wind before computing the rate of fire spread.

In the following, the anelastic/compressible comparison is structured in two sections: (a) a comparison to the FireFlux I observations in terms of fire-induced flow and ROS is carried out in Section 5, and (b) a phenomenological analysis of the anelastic/compressible results is provided in Section 6 in terms of fire-induced structures and energy spectrum, with an additional impact study of the horizontal atmospheric resolution.

Comparison to FireFlux I Observations Main Tower Fire-Induced Flow

Figure 2 shows the time series of air temperature at different heights at the main tower for the anelastic and compressible ensembles. Available measurements are also presented for comparative purpose; they are filtered by a Hodrick-Prescott filter (Hodrick & Prescott, 1997) to extract the trend from the observed time series (the variability is represented by the 5-s moving standard deviation). The time axis corresponds to the time since the fire passed at the main tower (denoted by t^a). Near the surface, that is, at 2 and 10 m, the compressible and the anelastic results show a temperature peak that is almost identical. The compressible ensemble shows an increase in temperature however ahead of the fire front passage (around −25 s), which means that there is a stronger horizontal advection of hot air. At higher altitudes, that is, at 28 and 42 m, the peak temperature of the compressible ensemble arrives earlier than for the anelastic ensemble, indicating a more tilted plume in the compressible simulations, which seems more consistent with the measurements. The temperature peak is also less intense. The ensemble spread related to atmospheric variability is slightly smaller in the compressible results than those of the anelastic ones. For both systems, the ensemble spread increases with altitude.

Figure 2. Temporal evolution of air temperature at different heights of the main tower: (a) 42 m, (b) 28 m, (c) 10 m, and (d) 2 m, for anelastic (blue) and compressible (red) ensembles. The ensemble mean is indicated in thick solid lines. The ensemble standard deviation is indicated in colored areas. Available measurements are shown in gray and black solid lines.

Figure 3 shows the horizontal wind speed time series at the main tower. At 2 m and 10 m, the fire-induced wind is greater in the compressible simulations than in the anelastic simulations (around 2 ms⁻¹ higher) and deviates more from the measurements. The stronger wind near the surface with the compressible version explains the higher advection of hot air noticed in Figure 2d. The differences between the anelastic and compressible results are larger than the atmospheric variability due to turbulence. This means that the high wind divergence constraint of the anelastic version is a strong flow assumption. It also seems that the ensemble variability due to turbulence is lower in the compressible simulations than in the anelastic simulations.

Figure 3. Temporal evolution of horizontal wind speed at different heights of the main tower: (a) 42 m, (b) 28 m, (c) 10 m, and (d) 2 m, for anelastic (blue) and compressible (red) ensembles. The ensemble mean is indicated in thick solid lines. The ensemble standard deviation is indicated in colored shaded areas. Available measurements are shown in gray and black solid lines (black dots represent failed sensor moments).

The vertical wind time series at the main tower is presented in Figure 4. At all heights, the compressible ensemble is closer to the measurements. In particular, the compressible simulations feature a better time positioning of the main updraft as well as a correct representation of the subsidence induced by the fire front passage. The maximum value of the vertical wind speed is lower in the compressible simulations than in the anelastic simulations. The updraft simulated with the compressible model is more consistent with the measurements.

Figure 4. Temporal evolution of vertical wind speed at different heights of the main tower: (a) 42 m, (b) 28 m, (c) 10 m, and (d) 2 m, for anelastic (blue) and compressible (red) ensembles. The ensemble mean is indicated in thick solid lines. The ensemble standard deviation is indicated in shaded colored areas. Available measurements are shown in gray and black solid lines (black dots represent failed sensor moments).

Overall, the compressible version of the coupled MesoNH-Blaze model tends to better reproduce the fire-induced flow at the main tower of the FireFlux I experiment (except for the 2-m and 10-m horizontal wind). The near-surface plume dynamics seems better represented with the compressible version due to less approximations on the buoyancy force. It appears also less sensitive to atmospheric inflow turbulence variability.

Rate of Fire Spread

For each ensemble member, the fire front position at a given time t can be derived from the arrival time matrix t^a (x, y). If all simulated fire fronts within an ensemble cross the same point (x, y), one can compute the averaged arrival time of the ensemble $$\overline{{t}_{k}^{a}}$$ as the average of the arrival times of all members: [Image Omitted. See PDF]where N_e = 15 is the ensemble size. The mean fire front plotted in Figure 5 corresponds to the mean arrival time matrix $$\overline{{t}^{a}(x,y)}$$. To obtain a fire front spread, the quartiles Q₁ and Q₃ of the propagation are determined. The interquartile range provides information about the variability of the spread within the ensemble.

Figure 5. Time-evolving fire front position for the anelastic (blue) and the compressible (red) ensembles 1, 2, 5, and 7 min after ignition. The mean fire front derived from the ensemble-averaged arrival time matrix is represented using a thick solid line. The ensemble variability is represented by the interquartile range in thin dashed lines. Statistics for the rate of spread between the two instrumented towers are also provided (in terms of mean and standard deviation for the simulated results). The ignition point is represented by the red dot. The + symbol represents the main tower. The × symbol represents the small tower.

Results presented in Figure 5 indicate that the compressible ensemble is significantly faster in terms of ROS than the anelastic ensemble. This is consistent with a stronger horizontal wind speed near the surface for the compressible simulations (Figure 3d). The mean fire front obtained with the compressible ensemble is found to be more extended with a more pronounced peak shape than that of the anelastic one. The difference in ROS between compressible and anelastic ensembles is larger at the fire head than on the flanks. This can be explained by two factors. First, on the eastern flank, the wind is nearly perpendicular to the fire front normal, meaning that a wind speed increase has less influence on the projected wind speed used by the spread-rate parameterization. This effect explains the more pronounced peak shaped front in the compressible ensemble. Second, the fire spread is forced at the top boundary of the burn plot by the drip torch ignition process. Lastly, the effect of the atmospheric variability on the fire is not as great on the flanks due to this delay in the free propagation regime.

The variability in the fire front propagation becomes larger as the fire propagates but also depends on the angle between the fire front propagating direction and the wind direction. Indeed, the closer this angle is to 90°, the larger the variability. In such a situation, a small change in the wind direction or in the fire front propagation direction can significantly modify the ROS.

The average ROS between the two instrumented towers was measured experimentally (1.61 ms⁻¹). The same quantity can be estimated for the compressible ensemble (ensemble mean value equal to 1.77 ms⁻¹) and the anelastic ensemble (ensemble mean value equal to 1.41 ms⁻¹). The statistics of the simulated ROS are reported in Figure 5. The anelastic ensemble features a lower ROS than in the observations. This result was expected by increasing the ignition temperature (compared to the validation study presented in Costes et al. (2021)). Contrarily, the compressible ensemble shows a ROS that is higher than in the observations. The relative differences between the ensemble-averaged ROS and the observed ROS for the anelastic and compressible ensembles are −12.2% and +9.8%, respectively.

In this section, results have shown that accounting for the high-density gradients has a significant impact on the fire spread in terms of fire front speed and geometry. In the present coupled model configuration (ROS and flux parameterizations), the compressible version of MesoNH-Blaze leads to a fire front propagation that is faster than the actual propagation observed during the FireFlux I experiment, whilst the anelastic version is lower. After this validation step, differences in terms of processes are analyzed.

Phenomenology Analysis

Beyond observations at the main tower, it is of primary interest to analyze the fire-induced flow in a more thorough way. High heat fluxes released by the fire induce structures in different directions and of different sizes. This section investigates these structures and some of their characteristics. For this purpose, it is necessary to consider the ensemble and the flow dependency on atmospheric turbulence. If it is easy to aggregate time series and spectrum, it is not as easy for 2-D or 3-D physical fields. For simplicity purpose, only one member of each anelastic/compressible ensemble is selected here for comparison. The selected member is that closest to the ensemble mean (in terms of correlation and mean error on 2-m temperature and horizontal wind at the main tower).

Horizontal and Vertical Fire-Induced Structures

The objective here is to compare the fire-induced flow throughout horizontal and vertical cross-sections, and to study the 3-D structures induced by the fire heat fluxes.

Figures 6a and 6b represent the horizontal cross-section of the 2-m air temperature for the anelastic (Figure 6a) and compressible (Figure 6b) simulations after 200 s. The compressible simulation provides a slightly lower maximum temperature than the anelastic simulation, but the warm area (where T_2m > 20°C) is more widespread downstream of the fire front. This is indeed referred to as downstream from a flow perspective, that is, following the streamlines. These streamlines tend to converge downstream of the fire front for the anelastic simulation, resulting in a limited temperature advection to the unburnt area. In the compressible simulation, streamline oscillations can be observed along the eastern flank of the fire front. Such oscillations are driven by small structures which develop when the wind and the fire front normal are orthogonal. The wind at this location follows the fire front line and the air heats up due to the intense heat fluxes. This leads to a wind convergence, itself inducing an updraft. Due to the highly unsteady nature of the flow (Finney et al., 2015; Frangieh et al., 2018), these convergences are intermittent and so are the induced updrafts.

Figure 6. Horizontal cross-section of 2-m air temperature (top panels) and 2-m horizontal wind speed (bottom panels) 200 s after ignition for a member of the (a, c) anelastic and (b, d) compressible ensembles. Streamlines represent the 2-m horizontal wind speed. The thick black curve represents the fire front simulated by Blaze. The thin black lines represent the burn plot boundaries. The ignition point is represented by the red dot. The + symbol represents the main tower. The × symbol represents the small tower. The dotted rectangle represents the domain used to calculate kinetic energy spectra (Figure 12).

Horizontal cross-sections of the 2-m horizontal wind speed (Figures 6c and 6d) show that some areas are subject to a significant wind acceleration (up to 20 ms⁻¹) in the compressible case. This phenomenon occurs at regular time intervals during the fire front propagation. In contrast, the anelastic case features a zone of low wind speed, downstream of the fire front, due to wind convergence (which is noticeable by looking at the streamlines) as an effect of the strong constraint on the wind divergence.

The horizontal divergence, $$\delta =\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}$$, is shown in Figures 7a and 7b. In the compressible case, gravity waves are found downstream of the fire front, causing this intermittent character by the alternate presence of maxima and minima. This is due to the shift of the propagation regime toward a wind-driven regime. These waves are noticeable in the vertical wind (Figures 7c and 7d). Gravity waves are clearly identified on the vertical density anomaly cross-sections (Figure 8b). The section plane of Figure 8 is indicated as a green line in Figures 7c and 7d. Note that the streamlines oscillate beyond the fire front in the compressible case. This type of oscillation can be found in wind-driven configurations (Frangieh et al., 2018; Ghaderi et al., 2021).

Figure 7. Same legend as Figure 6 for the 2-m horizontal wind divergence (top panels) and the 4-m vertical wind (bottom panels).

Figure 8. Vertical cross-section of the ρ′ density anomaly, 200 s after ignition, according to Figure 6 cross-sectional plane, for one member of the (a) anelastic and (b) compressible ensembles. Streamlines correspond to the wind. The T = 60°C isotherm is shown in a solid gray line.

It is also interesting to compare the flow on a vertical cross-section within the thermo-convective plume formed by the fire. This is done 500 s after ignition. For this purpose, the direction vector is determined $${\underline{e}}_{s}$$ of the slice plane in the direction of the mean flow over the first 50 m of the atmosphere. The point of origin of the slice plane is chosen so that the slice plane passes through an area of significant updraft. The updraft regions are shown on a horizontal cross-section of the 28-m vertical wind speed (Figure 9).

Figure 9. Same legend as Figure 6 for the 28-m vertical wind 500 s after ignition. The dotted rectangle represents the domain used to compute probability density distributions (Figure 11).

The main updraft is located above the fire front at an angle determined by the downstream entrainment of the warm air by the mean wind (Figure 10). The very hot air at the surface accelerates due to buoyancy. As density decreases with altitude, the rising warm air eventually comes into equilibrium with the air parcels around it, but continues to rise by inertia (Figures 8a and 8b). Above the density equilibrium altitude, buoyancy tends to slow down the rising air (because the density anomaly becomes positive), leading to a decrease in the vertical wind speed and reaching a zero speed value at an altitude about 1000-m. At this altitude, the air coming from the surface is denser than the surrounding air (positive density anomaly) and will therefore descend. This is the reason a subsidence is observed downstream of the plume top (Figures 10c and 10d). Air will descend once more to reach the density equilibrium altitude, which here corresponds to the maximum altitude of the boundary layer (about 600 m). The simulated fire plume obtained with the compressible version is subject to higher vertical velocities than in the anelastic case, and rises at a higher altitude.

Figure 10. Vertical cross-section of the density anomaly ρ′ (top panels) and of the vertical wind w (bottom panels) 500 s after ignition, according to Figure 9 cross-sectional plane, for one member of the (a, c) anelastic and (b, d) compressible ensembles. Streamlines correspond to the wind. The T ⩽ 60°C region is shown in gray (top panels) and in orange (bottom panels).

Plume updrafts and downdrafts statistics are represented in Figure 11 in terms of probability density distribution profiles of vertical velocity (computed within the dashed domain represented in Figure 6 to capture the entire plume). The compressible simulation gives much stronger updrafts of around 450 m in height, and much stronger downdrafts near the plume top. Relaxing the anelastic constraint has an important effect therefore on the ability of the coupled model to represent highly buoyant flows and the compensatory subsidences.

Figure 11. Vertical profiles of probability density distributions of the vertical velocity for the (a) anelastic and (b) compressible simulations 500 s after ignition computed in the dashed domain of Figure 9.

The fire-induced dynamical structure indicates the significant impact of the horizontal density gradients and buoyancy with the compressible version of MesoNH-Blaze at 10-m horizontal resolution. The fire front downstream region differs significantly in compressible simulations with the triggering of gravity waves which seem realistic by analogy with gravity waves produced by strong cloud convection (Lac et al., 2002, among others).

Kinetic Energy Spectra

The compressible version does not impose any constraint on the wind divergence unlike that of the anelastic one. The compressible simulations should be able from a theoretical viewpoint to form finer structures, especially in areas where there is a strong horizontal density gradient. In order to test this assumption, an interest in computing and analyzing the kinetic energy spectra for ensembles is required. Energy spectra are computed using Ricard et al. (2013)'s method for meridional and vertical wind components near the ground, over a rectangle around the burning area represented by the dashed box in Figure 6, for each ensemble member (Figure 12). In order to compare the situation with and without fire, kinetic energy spectra associated with the turbulence generation spin-up are also calculated at a 2-min time interval.

Figure 12. Kinetic energy spectra 200 s after ignition for the (a) 2-m meridian wind and (b) 4-m vertical wind, computed in the dashed rectangle of Figure 6. Blue curves correspond to the anelastic ensemble, red curves to the compressible ensemble, and gray and black curves to the FireFlux I case prior to ignition for anelastic and compressible versions, respectively. Solid lines correspond to the ensemble mean, and shaded areas to the ensemble standard deviation. The theoretical −5/3 slope is indicated in panel (a).

Without fire (black and gray curves), no significant differences between the anelastic and compressible versions can be noted. This allows a verification, once again, the coherence of the compressible system with the anelastic system without convection. With the fire, both anelastic and compressible ensembles show an increase in energy at all scales compared to the no-fire configuration. The near-surface meridian wind energy (Figure 12a) for the compressible ensemble is between four and 10 times larger at small scales (<70 m) than the anelastic ensemble. Relaxing the anelastic constraint promotes the generation of fine scale structures carrying energy. The compressible ensemble also shows much more variability than the anelastic ensemble at finer scales. With the 4-m vertical wind spectrum (Figure 12b), gravity waves already identified in the horizontal cross-sections (Figure 7) are clearly visible. An energy peak is noticeable between 50-m and 80-m wavelength for the compressible ensemble. This corresponds to the typical period of the gravity waves generated by the fire forcing.

It is also notable that the effective resolution here is consistent between both systems, around 7Δx – 8Δx, and in agreement with Lac et al. (2018) for the fifth-order WENO wind advection scheme. As a reminder, the effective resolution corresponds to the abscissa of the intersection between the −5/3 slope line estimating the energy cascade (the inertial zone) at the largest scales (>300 m), and the steepest slope line estimating the diffusive zone at the finest scales (<100 m).

Impact of Horizontal Resolution on High-Density Gradient Effects

We have shown that the effects induced by the compressible version of MesoNH-Blaze are pronounced at a 10-m horizontal resolution because of the strong horizontal density gradients. When the resolution is reduced, the gradients will clearly decrease, because of the coarser atmospheric cells on one hand, and because of the reduced intensity of the heat fluxes injected by Blaze on the other. Averaging the heat fluxes from the fire cell in a coarser atmospheric cell implies a temporal spread of the energy injection and therefore a drop in their instantaneous value. To further evaluate the impact of the horizontal resolution, the same coupled simulations as in the previous sections are run at a 25-m resolution and compared to the 10-m resolution results.

Figure 13 shows the horizontal density gradients for the anelastic ensemble (Figures 13a and 13b) and for the compressible ensemble (Figures 13c and 13d) at a 25-m and a 10-m resolution after 200 s. In the anelastic cases, air density is deduced from the state equation. The differences between the anelastic and compressible simulations are clearly visible at a 10-m resolution with a very strong gradient in the anelastic case and the presence of small ondulatory structures in the compressible case. At a 25-m resolution, density gradients are about three times smaller than at a 10-m resolution. Differences between the anelastic and compressible simulations are significantly reduced, and gravity waves are no longer present in the compressible version. The flow is also relatively similar between the two configurations.

Figure 13. Horizontal cross-section of 2-m horizontal density gradient ∇x,yρ for 10-m atmospheric resolution (top) and 25-m atmospheric resolution (bottom) 200 s after ignition for one member of the (a, c) anelastic and (b, d) compressible ensembles.

Time series of the 28-m vertical wind, 2-m horizontal wind speed, and 2-m temperature for the anelastic and compressible ensembles at 25-m resolution (Figure 14) show very small differences between both systems compared to the 10-m resolution results. The difference between the two ensemble means is less than the atmospheric variability.

Figure 14. Evolution of (a) 28-m vertical wind speed, (b) 2-m horizontal wind speed, (c) 2-m air temperature for the anelastic and compressible ensembles at 25-m horizontal resolution.

The ensemble-mean ROS at 25-m resolution is 1.37 ms⁻¹ for the anelastic ensemble and 1.59 ms⁻¹ for the compressible ensemble. The slightly stronger horizontal wind in the compressible version (Figure 14b) leads to a faster fire front propagation than in the anelastic version, but it remains slower than the ROS in the compressible simulations obtained at 10-m resolution.

This study has shown significant differences between the anelastic and compressible systems for horizontal density gradients higher than 10⁻² kg m⁻⁴. For gradients on the order of 10⁻³ kg m⁻⁴, some differences exist but they are more tenuous. It would be interesting to push toward even lower gradients to observe whether or not there is a convergence of the two systems of equations. This would mean that the thermal expansion effects are negligible. The validity domain of the anelastic approximation in the framework of the simulation of fires can not be represented by a strict limit which manifests itself rather as a range of resolutions in which the anelastic solution drifts progressively away from the compressible solution.

Conclusion

The focus of this article is to assess the validity of the anelastic approximation of the atmospheric model component in the context of wildland fire behavior simulation, where strong horizontal density gradients can occur and challenge a coupled atmosphere-fire model. For this purpose, a detailed comparison of the anelastic and compressible equation systems has been carried out for the FireFlux I grassland experimental fire simulated by the MesoNH-Blaze coupled model.

To achieve this goal, a compressible version of the MesoNH atmospheric model has been developed and validated on canonical cases. The compressible system has been applied to the FireFlux I fire spread case with MesoNH coupled to the Blaze fire model. Ensembles of 15 members have been produced in order to account for the intrinsic variability due to atmospheric turbulence.

Several key differences between the anelastic and compressible systems have been identified at a 10-m resolution, schematically represented in Figure 15. Relaxing the anelastic constraint on the wind divergence leads to faster horizontal wind speed and the development of gravity waves. Additionally, fine-scale eddies are much more present in the compressible simulations. The stronger horizontal surface wind implies a 25% faster fire spread compared to the anelastic simulations, and thereby a more pronounced peak shape for the fire fronts simulated with the compressible system. The plume dynamics are also affected. In the first decameters, the stronger horizontal wind associated with a weaker updraft gives a more tilted plume in the compressible simulations. Higher in the atmosphere, the compressible simulations feature stronger updrafts, as well as stronger downdrafts near the plume top as a compensatory contribution. All these dynamical characteristics seem to be realistic despite the absence of observations formally attesting to them. Better consideration of thermal expansion effects through the compressible system of equations seems to be valuable for the stimulation of wildland fire behavior at metric resolution.

Figure 15. Scheme highlighting differences in the fire-atmosphere interactions for the anelastic and compressible versions of the coupled fire-atmosphere model in case of high horizontal density gradients.

Using a coarser atmospheric grid (25 m instead of 10 m) induces lower horizontal density gradients (about 10⁻³ kg m⁻⁴ instead of 10⁻² kg m⁻⁴ at a 10-m resolution). For this order of magnitude of horizontal density gradients, differences between the anelastic and compressible systems are limited. This result suggests that the compressible formulation loses its added value as the horizontal atmosphere model resolution becomes coarser and the horizontal density gradients become lower. The validity domain of the anelastic approximation in the framework of the simulation of fires can not be represented by a strict limit which manifests itself rather as a range of resolutions in which the anelastic solution drifts progressively away from the compressible solution. This study has shown that the 10-m resolution solution belongs to this range of scales, whilest the 25-m one does not seem to. The actual validity domain of the anelastic approximation depends on the expected level of accordance between the anelastic and compressible responses according to a given set of criteria. In the framework of coupled atmosphere-fire modeling, the 25-m resolution anelastic simulations give a satisfactory level of accordance to compressible simulations. This is not the case for the 10-m resolution simulations.

A more thorough study could be conducted to define more precisely the limits of the anelastic approximation validity domain. Evaluating horizontal density gradients may be an interesting proxy for the assessment of the choice of the atmospheric formulation that a given wildland fire simulation requires.

While this study focused on the physical aspects of the compressible solution, future work is necessary to make the compressible system more efficient from a numerical perspective. As it is classically performed in fully compressible atmospheric models, the compressible system could be made less numerically expensive by using a time-splitting method, which allows two different time steps for the fast and slow terms of the system of equations (Klemp et al., 2007). This approach has been explored but the robustness still needs to be validated before being applied to the context of coupled atmosphere-fire simulations.

For future real case studies of large-scale wildland fires with the MesoNH-Blaze coupled model, the horizontal resolution of the atmospheric model will be limited to a few dozen meters due to the size of the computational domain and also due to the need to run ensembles to characterize atmospheric variability. In this context, the anelastic system of equations seems to be an efficient way to simulate the atmospheric dynamics associated with the development of a wildland fire.

A AppendixPredictor Corrector Scheme for the Energy Equation

This section will focus on certain numerical aspects of the compressible system implementation in the MesoNH existing algorithmic environment. The compressible version of MesoNH does not need to redefine the discretization of the derivative operators compared to the anelastic version. In particular, the high-order numerical schemes (of order three to five) used for the advection terms can be reused for the compressible version. Several schemes are available, including the third-order and fifth-order WENO (Liu et al., 1994; Lunet et al., 2017), and the fourth-order centered scheme (RKC4). The WENO schemes are implicitly diffusive and do not require additional numerical diffusion (unlike the fourth-order centered scheme). Nevertheless, they allow proper dealings with problems subject to strong gradients. These schemes are recommended for the use of compressible MesoNH because the added value of the compressible version is in the strong gradient areas precisely. The advection of scalar quantities is handled by a PPM scheme (Colella & Woodward, 1984). In general, these spatial operator discretization schemes are associated with high-order temporal discretization schemes such as the Runge-Kutta methods, which are iterative temporal integration methods. Again, the same methods can be used for compressible and anelastic formulations.

A predictor/corrector type method is applied in the compressible system to compute the potential temperature θ. The energy conservation equation (Equation 15) introduces a time derivative of pressure, which is a diagnostic variable. To understand the numerical scheme used to solve this equation, one rewrites the energy conservation equation (Equation 15) as follows: [Image Omitted. See PDF]

The purpose of the method is to obtain the potential temperature θ at time (t + Δt), denoted n + 1, from the quantities at time t, denoted n. The previous calculation steps allow the computing of $${\rho }_{d}^{n+1}$$ and $${r}_{v}^{n+1}$$ (and thus ρⁿ⁺¹). Using a forward-in-time Euler scheme and considering $${\left({\rho }_{d}\theta \right)}^{n}={\rho }_{d}^{n}{\theta }^{n}$$, one can write: [Image Omitted. See PDF]

To evaluate this expression, we need Πⁿ⁺¹, which is a diagnostic quantity depending on θⁿ⁺¹ through the state equation. We will therefore use intermediate variables, noted *, which will be used to estimate Πⁿ⁺¹ (Π* ≈ Πⁿ⁺¹). One starts by calculating an estimate of the potential temperature without the term of pressure temporal variation: [Image Omitted. See PDF]which allows to estimate the intermediate pressure: [Image Omitted. See PDF]

This estimator is then used to calculate the time derivative of pressure and to compute the potential temperature θⁿ⁺¹: [Image Omitted. See PDF]

The diagnosis of the pressure Πⁿ⁺¹ can then be carried out: [Image Omitted. See PDF]

This numerical method allows the potential temperature to be calculated without making any assumption about a preferred direction, in space or time, of the pressure variation in the moisture correction term in the absence of phase change.

From an algorithmic perspective, the computation steps for buoyancy, via the thermodynamic variables (Equation 3), and the pressure solver are removed and replaced by a routine computing both pressure and buoyancy forces (Equations 10 and 11).

B AppendixValidation of the Compressible Version of MesoNH

This section describes two 2-D canonical test cases used to validate the compressible version of MesoNH: (a) a linear orographic wave, and (b) a propagating pressure wave. The first test case considers an air flow over a non-flat terrain; the objective is to validate the compressible solution for a stationary flow through comparison with an analytical solution and to the anelastic simulation. The second test case considers a density perturbation that is not pressure-balanced and that induces the formation of a spherical pressure wave; the objective is to validate the propagating speed of the pressure wave against theoretical considerations and to highlight the potential benefits of a compressible system of equations. It is worth noting that a simulation of the atmosphere at rest has been performed on a preliminary basis to verify that the prognostic variables of the compressible version remain at hydrostatic equilibrium.

B1. Linear Orographic Wave

This test case corresponds to a theoretical atmospheric case combining synoptic wind and terrain, for which an analytical solution exists. The goal of this test is to assess the validity of the compressible version of MesoNH for typical atmospheric flow. The passage of a stably stratified current over a terrain heterogeneity induces a vertical displacement of the fluid parcels, which may oscillate from their equilibrium state and lead to a gravity wave propagation. This gravity wave is primarily governed by the spatial extension and the shape of the terrain heterogeneity. The wave is hydrostatic when the horizontal extent of the terrain heterogeneity is very large compared to its vertical extent. In this case, the wave propagation is vertical and the vertical wind speed has a reference solution (Alaka, 1960; Scorer, 1949).

The computational domain corresponds to a 2-D rectangle of L_x = 800 km by L_z = 30 km, with a terrain heterogeneity following an Agnesi curve $$h(x)={h}_{m}\left(\frac{{a}^{2}}{{x}^{2}\,+\,{a}^{2}}\right)$$, with h_m = 1 m and a = 10 km. The surface air temperature is 250 K, the surface pressure is 10⁵ Pa, and the mean flow velocity is U = 20 m ⁻¹. The atmosphere is considered to be stable with a Brunt-Väisälä frequency N = 0.2 ⁻¹. The spatial resolution is Δx = 500 m and Δz = 250 m. The numerical schemes are the fifth-order WENO (WENO5) in space and the five-stage third-order Runge-Kutta (RK53) scheme in time for the wind, and the positive definite PPM for scalar variables. The anelastic time step is 10 s, and the compressible time step is 0.2 s. The numerical configuration is similar to the one used in Lunet et al. (2017).

Both anelastic and compressible solutions show a very good agreement with the reference solution (Figure B1). There are some very small discrepancies between the numerical solutions and the reference solution, but the compressible and anelastic results are identical. They are also consistent with WRF's results in the same configuration (Doyle & Skamarock, 2005).

Figure B1. Vertical cross-section of the vertical wind speed for the linear orographic wave test case. Dashed line represents the analytical solution. Colors represent the MesoNH solutions in its (a) anelastic version, and (b) compressible version after 10 hr of simulation.

B2. Propagating Pressure Wave

The goal of this test is to assess the validity of the compressible formulation for a flow with compressible effects. A density anomaly ρ′ is considered so that $$\rho =\bar{\rho }+{\rho }^{\prime }$$, where $$\bar{\rho }$$ is the hydrostatic density. From the equation system, this density anomaly satisfies the following acoustic propagation equation: [Image Omitted. See PDF]with $${c}^{2}=\frac{\gamma P}{\rho }=\gamma RT$$ the speed of sound. The compressible version of MesoNH is expected to correctly compute such pressure waves.

The size of the computational domain is L_x = 16 km and L_z = 10 km. Initially, the atmosphere is considered at hydrostatic equilibrium and without moisture. The surface air temperature is 300 K, and the surface pressure is 10⁵ Pa. The atmosphere is considered to be stable, $$\frac{\partial \theta }{\partial z} > 0$$, with a Brunt-Väisälä frequency $$N=\sqrt{\frac{g}{\theta }\frac{\partial \theta }{\partial z}}=0.02$$ s⁻¹. To generate a pressure wave, an initial density perturbation $$\tilde{\rho }$$ of amplitude ρ_a and characteristic size h_a is added as follows: [Image Omitted. See PDF]where ρ_a = 0.01 kg m⁻³ and h_a = 300 m.

The spatial resolution is uniform with Δx = Δz = 50 m. The time step is fixed with Δt = 0.05 s. The selected numerical schemes are the centered fourth-order scheme in space and the centered fourth-order Runge-Kutta scheme (RKC4) in time for wind transport, and the monotonic PPM in space and the forward-in-time temporal scheme for meteorological and scalar variables (Lac et al., 2018). No numerical diffusion is added.

Figure B2a represents the time-evolving density anomaly ρ′, which is computed as the difference between the density at time t and the density at hydrostatic equilibrium. Results show that the perturbation is initially at the domain center. It then propagates quasi-isotropically.

Figure B2. (a) Time-evolving vertical cross-sections of the ρ′ density anomaly related to the propagating pressure wave. (b) Pressure wave position as a function of time. Solid lines represent the simulated wave position determined from the local pressure maximum for the upward and downward waves. Symbols represent the position of a theoretical wave propagating at the sound speed c.

The propagating speed of the pressure wave is compared with the local sound speed c of a theoretical wave. To construct this theoretical wave, a recursive procedure is used starting from the initial time t₀ at which the initial density perturbation forms two pressure perturbations (an ascending wave and a descending wave): [Image Omitted. See PDF]where z_th,up(t) is the position of the theoretical ascending wave, z_up(t) the position of the simulated ascending wave, $$c\left({z}_{th,up}(t)\right)=\sqrt{\frac{\gamma P\left({z}_{th,up}(t),t\right)}{\rho \left({z}_{th,up}(t),t\right)}}$$ the sound speed calculated from the atmosphere properties at the considered altitude, and Δτ is an arbitrary time step. The same procedure can be applied to the descending wave. Figure B2b shows that the simulated wave propagates at the correct speed, both in the upward and downward directions. The differences between the theoretical and simulated descending waves after 14 s are due to the simulated wave interaction with the ground surface.

NomenclatureAbbreviation
• HWS
• Horizontal wind speed
• PPM
• Piecewise parabolic method
• RK
• Runge-kutta
• ROS
• Rate of spread
• WENO
• Weighted essentially non oscillatory

Superscript
• n
• Time index

Subscript
• d
• Dry atmosphere
• ref
• Reference state
• v
• Water vapor

Symbol
• c
• Sound speed
• C_pd
• Dry air calorific capacity at constant pressure
• C_ph
• Moist air calorific capacity at constant pressure
• C_vd
• Dry air calorific capacity at constant volume
• g
• Gravitational constant
• P
• Pressure
• P₀₀
• Reference pressure at surface level
• r_v
• Water vapor mixing ratio
• $$\mathcal{R}$$
• ROS
• R
• Gas constant
• t
• Time
• T
• Temperature
• T_i
• Ignition temperature
• t^a
• Fire arrival time
• U
• Wind vector
• w
• Vertical wind speed
• Δt
• Time step
• Δx
• Atmospheric horizontal mesh size
• Π
• Exner function
• ρ
• Moist air density
• ρ_d
• Dry air density
• θ
• Potential temperature
• ∇
• Gradient operator
• Ω
• Earth angular velocity

Acknowledgments

The Authors gratefully acknowledge support from ANR (ANR-16-CE04-0006, FIRECASTER). The Authors gratefully acknowledge Jean-Baptiste Filippi (CNRS) and Ronan Paugam (UPC) for helpful discussions on coupled atmosphere-fire modeling and wildland fire modeling. They acknowledge Thomas Burgot (CNRM/Météo-France) for the stimulating discussions on anelastic and compressible models. They also acknowledge Juan Escobar (LAERO) for support on MesoNH.

Data Availability Statement

MesoNH-BLAZE (v5.4.4) is available at mesonh.aero.obs-mip.fr/mesonh54. Data used for this paper are available at Costes (2022).