Role of Horizontal Eddy Diffusivity within the Canopy on Fire Spread

: Wind proﬁle observations are used to estimate turbulent mixing in the atmospheric boundary layer from 1 m up to 300 m height in two locations of pine forests characteristic of the southeast US region, and to 30 m height at one location in the northeast. Basic turbulence characteristics of the boundary layers above and within the canopy were measured near prescribed ﬁres for time periods spanning the burns. Together with theoretical models for the mean horizontal velocity and empirical relations between mean ﬂow and variance, we derive the lateral diffusivity using Taylor’s frozen turbulence hypothesis in the thin surface-fuel layer. This parameter is used in a simple 1D model to predict the spread of surface ﬁres in different wind conditions. Initial assessments of sensitivity of the ﬁre spread rates to the lateral diffusivity are made. The lateral diffusivity with and without ﬁre-induced wind is estimated and associated ﬁre spread rates are explored. Our results support the conceptual framework that eddy dynamics in the fuel layer is set by larger eddies developed in the canopy layer aloft. The presence of ﬁre modiﬁes the wind, hence spread rate, depending on the ﬁre intensity. large horizontal diffusivities. modeled


Introduction
Numerous factors affect the spread of fire in a wildland environment, beginning with the vegetation or fuel present, and the underlying surface topography. Atmospheric conditions, especially humidity and wind speed, are critical factors in the development of fire and rate-of-spread (RoS) calculations. Wind is rarely steady, and fluctuations at many time scales are typically present, with a dominant eddy length and time scale increasing with height away from a smooth wall. Such eddies are responsible for much of the turbulent flux in the boundary layer. Over rough forest canopies of height h, a background mean flow U produces a time scale h/U that describes eddy mixing in neutral stability [1]. The presence of fire within the canopy generates its own circulation which interacts with the background flow in which it is embedded, enhancing turbulence levels [2] and modifying the mean flow structure near the boundary. A dynamical view of fire behavior emphasizes the importance of these interactions and their role in fire spread.
The mean flow structure within the canopy layer is drastically different from that in the conventional atmospheric surface layer. A logarithmic profile for the mean stream-wise velocity is a good approximation for the region above the canopy (referred to as the conventional surface layer), whereas other models such as exponential and hyperbolic functions best fit the observations within the canopy layer (an overview of these models can be found in [3]). A layer directly above the canopy where the vertical profile of horizontal velocity deviates from the logarithmic profile will be referred to as the roughness sublayer (after Harman and Finnigan [4]). Figure 1b provides a sketch of these distinct turbulence layers above and within the canopy.
The properties of the flow in the roughness sublayer resemble a plane mixing layer [5], where the flow is controlled mainly by unique length and velocity scales. In the canopy layer, turbulent motion is more effective [3] with the dominant eddies being height independent [4,6], and with additional eddies introduced by wake flow around branches and stems. Because of this fact, a single length scale alone cannot fully represent the spectrum of multiscale interactions, and a more sophisticated approach would be needed to determine an effective scale [7]. In any case, the vertical transport of horizontal momentum is typically absorbed by the canopy layer through canopy drag. The amount of wind reduction varies dramatically from one type of canopy to another [8] due to structural differences between trees and other vegetation comprising the canopy.
For the case of ground fires in the areas described here, the burning fuel layer-typically pine needles and grass-occupies a layer less than a few tens of centimeters height. This fuel layer may be below the lowest wind measurement level in an observational system, and is often embedded in shrub that, while it might ignite, may not carry the fire. This situation is represented by a fuel layer that itself acts as a thin "fuel canopy" within the larger tree canopy. As such, the fuel layer is expected to show distinct flow dynamics, subject to sub-canopy eddy structures, and strongly damped flow. Here, we focus on the dynamics of the shrub/fuel layer within the canopy layer ( Figure 1). The purpose of this paper is to examine the role of lateral mixing within the canopy and fuel layer in fire spread. Our goal is to estimate mean flow and lateral eddy diffusivity in the canopy and fuel layer as a function of height using observational datasets and theoretical models, and then test to their effects in an idealized 1-D numerical model for fire spread. The lateral diffusion of heat influences the rate of fire spread as it displaces the fire front. We specifically focus on the extrapolation of available wind measurements within the canopy layer into the fuel layer. The goal of the model is to illustrate quantitatively the role of lateral diffusivity in fire spread, in an idealized context. This model is not meant at this stage to simulate the fire spread in a realistic situation.

Data
The data analyzed in this paper was collected during different burn events at the following locations: Pebble Hill Plantation (4/2018, Thomasville, GA, USA; low-intensity prescribed burn; head fire), Tall Timbers Gallien Burn Unit (4/2018, Tallahassee, FL, USA; low-intensity prescribed burn; head fire), The Pinelands National Reserve (3/2011 and 3/2012, New Lisbon, NJ, USA; low-intensity prescribed burn; backing fires). The ignition was by lines in all experiments.
Both the Florida and Georgia burn sites ( Figure 2) used two 1 Hz sonic anemometers placed at 1 m and 3 m above ground, and a ZXLidars ZephIR 300 collecting data from 10 m to 100 m and 300 m with various preset vertical resolutions. Nominal accuracy of the ZephIR data is 0.1 ms −1 and 0.5 • . Both sets of instruments were placed within 100-200 m of the prescribed burn. The burn sites are mainly composed of Pitch pine trees. The air relative humidity during the experiments were 37.1 ± 4.1% and 43.7 ± 1.9% at Pebble Hill Plantation and Tall Timbers Gallien Burn Unit, respectively. The 1-hour forest floor fuel moisture content was 7.9 ± 0.4% at Pebble Hill Plantation, and 7.7 ± 0.2% at Tall Timbers Gallien Burn Unit. The fire front did not pass through the instruments. The purpose of these instruments was to characterize the vertical structure of the surrounding surface atmospheric boundary layer during a low-intensity burn event.
The NJ Pinelands experiments provided data from a 30 m tower at three levels: 3 m, 10 m, and 30 m. Here, we use velocity measurements from sonic anemometers and temperature from Vaisala HMP50 sampled at 10 Hz frequency. The measurement towers were located within the burn site and a fire front passed through the instruments unlike the Florida and Georgia burn sites. In these experiments, Pitch pine along with dense scrub oaks and shrubs dominated each of the plots. The air relative humidity values during the fire periods in 2011 and 2012 experiments were approximately 30% and 20%, respectively. The 1-h forest floor fuel moisture content was 22.6 ± 11.4% in 2011, and 42.7 ± 13.4% in 2012. Further details on experimental design, locations of the towers and measurements techniques can be found in Heilman et al. [9] and Heilman et al. [10] (their Figure 1).

Theoretical Formulation
The recognition of canopy effects on the mean flow profile started around half a century ago (e.g., see an overview in [11]), yet the flow structure deep within the canopy, in the fuel layer, is typically not well documented. We start with a brief review of standard models within and above a canopy followed by a discussion of the flow structure in the lowest surface-fuel layer.

Horizontal Wind
Before proceeding to the overview of the theoretical formulations for the horizontal wind speed, note that we employ a standard Reynolds decomposition where each variable is decomposed into a mean (denoted by overbars) and a fluctuating part (denoted by primes). The mean parameters were computed as 30-min block-averaged values [12]. The along-wind horizontal velocity is u = u + u . We also use a micrometeorological coordinate system where there is no cross wind mean velocity (i.e., v = 0), and the vertical velocity component w is directed opposite to the gravitational force (for a more detailed description, see, e.g., [12]). In the conventional surface layer above the canopy (Figure 1), the mean velocity profile can be approximated as where k = 0.4 is the von Karman constant, u * (h) is the friction velocity at the top of the canopy of height h, d is the displacement height, z 0 is the roughness parameter, and z is height above the ground surface.
Given that the displacement height is estimated to range within 0.6-0.8 of the canopy height [13], we assume d = 0.7h. The friction velocity, in principle, can be computed as u * (h) = (−u w ) 1/2 | z=h , where u w is the eddy momentum flux near the canopy top. Since d (or h) and z 0 are not known a priori, we consider them as free parameters that are estimated by minimizing the residual sum of squares [14] between the theoretical profile (1) and the observations. For the same reason, u * (h) is considered to be a free parameter. Please note that we assume a neutral stratification ignoring thermal effects, and therefore (1) does not include Monin-Obukhov Similarity function [4,15]. In the canopy layer, under the assumption that the pressure gradient force is negligible, the analytical solution could take at least two forms: exponential and hyperbolic [3]. Here we consider the simplest versions of these relations as follows where α 1 and α 2 are empirical coefficients, and u(h) is the mean velocity at the top of canopy (that is determined from the observations given the estimated h from the log-law profile (1)). Solution (2) is derived following Prandtl's mixing-length hypothesis where the vertical momentum diffusivity, K m , depends on background velocity gradient as K m = l 2 ∂u/∂z; whereas the solution (3) assumes K m ∼ u [16]. The exponential dependence fits better in the upper portion of canopy, and the hyperbolic dependence fits better closer to the ground [12]. Recent studies have refined theoretical models of the mean velocity profile in the canopy to allow smoother transitions between the canopy sublayers ( Figure 1b) yielding a continuous function for u(z) from the ground up to the top of the canopy [17] and from the canopy through the roughness sublayer into the log-law layer [4]. Here, we considered only the simplest formulations in the canopy layer; these mean velocity profiles will be used in Section 2.2.2 to estimate the horizontal eddy diffusivity.
To continue the basic description of the vertical structure of wind in the atmospheric boundary layer it is also of interest to examine the local turbulent kinetic energy (TKE) e = 1/2(u 2 + v 2 + w 2 ). A short time-mean TKE is computed using a 1-min averaging window as opposed to a 30-min window, for computing the mean velocity. This is done to resolve higher frequency variability and the vertical propagation of eddy energy.

Horizontal Diffusivity
Although there has been a great deal of work devoted to the understanding of the vertical structure of K m , the vertical structure of the horizontal eddy diffusivity K h in a complex boundary layer structure remains an open question. Here we estimate K h using Taylor's frozen turbulence hypothesis: with T being the integral timescale determined from the velocity auto correlation and angle-brackets denote averaging over the time T. The limitation of this method, however, is that it is only applicable at heights where the observations are available (i.e., above the fuel layer). To estimate K h in the fuel layer, we use a model. The model is K h = Lu , where the length scale of the dominant flow in this layer is L. It is assumed that u (z) = βu(z). A similar scaling relation was derived by Nepf [18] who showed that the ratio of the turbulent velocity scale to the mean velocity scale depends on several factors, including drag and stem geometry in a model of turbulent intensity within emergent vegetation, similar to a fuel layer. We simplify the relation to a constant factor β, based on the observed correlations as discussed in Section 3.3.
Given that the size of the dominant eddies within the canopy layer is approximately height independent [11], the length scale governing the flow within the canopy layer and the fuel layer is assumed to be the height of the canopy L = h. Thus, having a reasonable estimate of u(z) in the lowest 1 m based on the theoretical formulations specified in Section 2.2.1 yields an estimate for K h in the fuel layer The full duration of the Gallien Burn and Pebble Hill data sets comprises roughly 2 h of observations. Given that the fire front never passed the site where observations were taken, the estimated K h from these data reflects the background levels of turbulent mixing. The NJ data sets, on the other hand, have at least 24-h of data spanning the period when the fire fronts were passing the towers. We identified these fire-passing periods based on the spikes in temperature measurements (∼30 min, as shown in Figure 3, marked with "fire"). This 30 min period is used to compute K h during the fire. To estimate the background K h , a 1.5 h period was selected before the fire front reaches the tower (Figure 3, marked with "pre-fire"). Please note that the 30-m tower was placed near the perimeters of the burn blocks where the impact of the fire lines was much less than the impact found at the towers in the interior [9].

Vertical Profile of Horizontal Wind
In the Gallien Burn data sets, the conventional logarithmic model (1) fit for the region above the canopy with h ≈ 24 m, u * (h) ≈ 0.26 ms −1 , and z 0 ≈ 0.1 m. The Gallien Burn wind profile clearly shows three distinct layers: conventional surface layer above ∼ 24 m, the upper canopy layer ∼10-24 m approximated with (2) where α 1 = 0.1 and the lower canopy layer below 10 m approximated with (3) where α 2 = 0.6 ( Figure 4a). For the models (2) and (3) we used the estimated canopy height h to determine u(h) from the observations (this value was allowed to vary by 10% for the best fit).
On the other hand, the Pebble Hill data set shows just two distinct layers: conventional surface layer above h ∼ 40 m and the canopy layer below it (Figure 4b). The best log-law fit suggests u * (h) ∼ 0.27 ms −1 , and z 0 ∼ 0.01 m for this site. The exponential model for this data set was not determined, instead (3) fits the data relatively well, where α 2 = 0.1 and u(h) is computed from the observations.
Finally, in the case of the subset of the Pinelands National Reserve dataset used here, the vertical sampling was limited to a few levels. Given this limited vertical resolution we only applied model (3). Using h = 20 m [19], the best-fit parameters are u(h) = 2.5 ms −1 and α 2 = 1 (Figure 4c).
In Section 3.3, we use the theoretical predictions for u(z) in the lower parts of the surface layer (i.e., below 1 m) to estimate lateral diffusivity in the fuel layer.

Vertical Structure of Turbulent Kinetic Energy
TKE variability for the Pinelands datasets has been thoroughly described by [19] and references therein. The greater vertical extent and upper layer resolution of the Gallien Burn and Pebble Hill LIDAR wind profiler data provides some additional insight to the nature of the turbulence influencing the canopy layer.
The temporal evolution of the TKE at these sites shows strong turbulent mixing events in the upper layers above the canopy and relatively weaker turbulence within the canopy layer ( Figure 5). At the same time, a substantial temporal variability in TKE occurs within the canopy layer itself. In both the Pebble Hill and Gallien sites, the TKE Hovmoeller diagrams show distinct vertically propagating structures ( Figure 5), in which energy reaching the surface layer is somewhat delayed by passage through the canopy. These observations suggest that turbulent momentum transport at the top of the canopy by downward penetration of high-momentum air (e.g.,"sweeps", [20]) drives variability in the surface layer.

Horizontal Diffusivity
An assessment of the horizontal diffusivity from the data following (4) shows that the background K h increases with height (Figure 6a,b). In agreement with previous studies [21], the integral timescale in (4) is approximately constant with the height within the canopy. During the fire period available in the Pineland datasets (Figure 6c), K h increases by an order of magnitude in the 2011 data and by a factor of 2 to 5 in the 2012 data. Such differences can be attributed to the intensity of the fire, with the 2012 fire experiment having a lower intensity compared to 2011 [9]. A composite plot of all sites (Figure 6d) shows that only the more intense fire of the 2011 Pineland burn stands out with relatively higher diffusivity. Diffusivity estimates are all at or above the height of the surface-fuel layer and tend, even in the case of the 2011 burn, to converge within the fuel layer toward values below about 1 m 2 s −1 . The estimates of K h following model (5) fit the observations within the canopy relatively well (Figure 6a,b) and suggest that within the surface-fuel layer, the lateral diffusivity ranges from 0.1 m 2 s −1 to 1 m 2 s −1 . In these computations we used β = 0.24 derived as a correlation coefficient between the standard deviation of u and u(z) (Figure 7).
In the next section we use these K h estimates in the context of a simple 1D model to examine fire spread rates in different wind conditions. To simulate different wind conditions we vary both the mean background wind and, with the horizontal eddy diffusivity, the background turbulent mixing component of heat flux.

Model Setup
A simplified 1D turbulent heat transport model is implemented to estimate the RoS inside a flat, homogeneous, thermally thin fuel bed (Figure 1b). The model represents the gas temperature, T(x, t), of the mixture of air, vaporized fuel, and smoke, which is undergoing combustion when the temperature is above some given threshold ignition value. More sophisticated models have been used to compare RoS directly with laboratory experiments (e.g., [22,23]), but our purpose here is confined to examining the basic dependence on wind speed and revealing the role of diffusion in a highly idealized context.
The gas mixture moves through a porous fuel matrix; the fuel bed loading, F(x, t), is assumed to be homogeneous initially. When the gas within the fuel matrix exceeds the ignition temperature, the fuel is consumed, heat is released, and the fuel depletes according to a simplified combustion equation where the amount of heat released is proportional to the amount of fuel remaining. Eventually the temperature drops below the ignition threshold, and the combustion source turns off. At time t 0 , the temperature at the center of the fuel bed is specified to be above the ignition threshold value in a simple Gaussian distribution, representing the ignition pattern.
The temperature of the gas is governed by the advection-diffusion-reaction equation

∂T ∂t
where u is an advective velocity, Q is the reaction term, and R is a radiation term. The diffusion coefficient, K h , is estimated as described in Section 2.2.2. Temperature has been shifted so that 0 corresponds to an ambient 290 K. The fuel mass is governed by where H is the Heaviside function. Therefore, fuel is consumed at an exponential rate when the temperature is above the ignition temperature T ign . The rate of consumption depends on a burn time, b, which we take from our previous work [24]. For each combustion interval (a k , b k ), we assume a constant upward velocity across the buoyant plume. By conservation of mass, the horizontal velocity is linear. To model eddies at the edge of the plume, we introduce a linear increase in the horizontal velocity leading up to the edge of the combustion region. In the case of a single combustion region, the horizontal velocity is illustrated in Figure 8. Finally, we add a uniform background flow u BG to the plume-induced horizontal velocity. The reaction term is and it corresponds to the amount of heat released by the burning fuel. Please note that this reaction is only active when the temperature is above the ignition temperature. Finally, the radiation of heat from the plume is modeled with a linearized Stephan-Boltzmann equation where λ is tuning parameter with dimensions s −1 . Ranges and units of the parameters for our model are summarized in Table 1.
The governing equations are non-dimensionalized by introducing a characteristic temperature scale T * , length scale x * , time scale t * , fuel mass scale F * , and velocity scale u * = x * /t * . Then, the governing Equations (6) and (7) become ∂T ∂t where each of the variables is now dimensionless. The constant to model radiation is set to λ = 0.2 s −1 for all the simulations.

Model Solution
Equation (6) is solved by applying a Strang splitting [25] to separate the advection term and the diffusion-reaction term. Then, the advection equation is solved with a semi-Lagrangian method while the diffusion and reaction equation is solved with a non-stiff finite difference method. The fuel Equation (7) is discretized with forward Euler. We assume homogeneous Dirichlet boundary conditions, and the simulation is stopped before the boundary conditions interfere with the solution.
The average RoS of the head fire is determined by finding the left-most point where the temperature exceeds the ignition temperature ( Figure 9). This RoS is reported in Figure 10 for a variety of horizontal diffusivity constants and background wind velocities. The results show the expected basic dependence of fire RoS on increasing background wind speed. In addition, stronger heat diffusion drives a greater RoS. At lower velocities and larger horizontal diffusivity, our model captures a backing fire ( Figure 11). The average RoS of the backing fire is reported in Figure 12 for a variety of horizontal diffusivity constants and background wind velocities. As expected, the RoS of a backing fire is largest at small background wind speeds and large horizontal diffusivities. The modeled RoS values agree with field observations reported for grassland fires (e.g., [26]).

Discussion
In many wildland fire situations, the fuel that carries the fire is concentrated in a thin layer close to ground, with a thickness much smaller than the typical canopy height. This fuel layer may be considered a distinct sub-canopy within the larger canopy structure. The mean wind and turbulence level within the fuel layer controls, to some extent, the combustion and spread rate. We analyzed new wind profiler observations of the boundary layer structure above and within a canopy and used canopy flow models to infer the basic mean flow and variance of the wind in the fuel layer.
Horizontal fluctuations of wind receive less attention than vertical fluxes but can play an important role in canopy heat and momentum balances. Strong downdrafts behind the fire line revealed by both numerical simulation Sun et al. [27] and direct measurement Clements et al. [2] are thought to play a key role in fire spread and fire-atmosphere interactions. A strong downdraft near the surface boundary generates, in turn, strong horizontal flow. Heilman et al. [28] noted that the horizontal diffusion of TKE in the fire environment in the Pinelands experiment was of the same order as vertical diffusion. We demonstrate here that the effective heat diffusion due to lateral wind fluctuations, or convective turbulence, also plays an important role in fire spread.
The horizontal extension of flame is an example of lateral heat flux of the hottest gases, and this is well-known to be subject to complex turbulent fluid behavior, even in the absence of background turbulence. Tang et al. [29] (see also [30]) carried out laboratory wind tunnel experiments with stationary burners to study flame intermittency and horizontal extension. In their setup the background turbulence was suppressed to focus on fire-generated effects. Their laboratory results show flame extension variance of l ≈ 10 cm and a frequency of f ≈ 10 Hz (an intermittency time scale of 0.1 s), consistent with an effective "flame" mixing coefficient of f l 2 = 0.1 m 2 s −1 in background mean flows of about 1 ms −1 . This is not much different from observed values in the field, taking into account that the fire spreads in the thin fuel layer at the base of the fire, where similar time scales and somewhat longer flame lengths occur and can raise the effective mixing coefficient.
The roughly linear wind dependence of the laboratory flame extension frequency and extension variance is also consistent with stronger diffusion in higher winds; similarly, simple advective eddy scaling implies a diffusion coefficient of size u l, and the relation shown here between mean wind and eddy velocity implies larger diffusion coefficient with higher mean wind. Our approach places the diffusion coefficient K at the forefront of spread rate consistent with the range of observed values, whether due to intrinsic fire turbulence and flame length or background eddies.
Efforts to directly measure the characteristics of the wind across the thin fuel layer require high resolution instrumentation. One way to achieve this is to use fine-wire thermocouples or static probes to measure high resolution gas temperatures and velocities in wildland fire experiments [31]. At the same time, LIDAR systems for quantifying the fine-scale structure of the canopy and fuels have become powerful tools to investigate fuel properties; we anticipate their value to investigate flow properties as well. High resolution wind measurements would go together with observations of fine-scale fuel structure. Together, these will help to elucidate the relation between flow in the fuel layer and the overlying canopy to determine flow regimes and turbulence intensities in the natural combustion environment.

Conclusions
Motivated by the need to accurately predict the RoS of the fire in a fuel layer, we explore four prescribed burn experiments: two in the southeast US region and two in New Jersey. The analysis of the available horizontal wind speed observations from 1 m up to 300 m allows us to conjecture the vertical profile of horizontal velocity in the fuel layer. We further use Taylor's frozen turbulence hypothesis to estimate the lateral diffusivity in the fuel layer. Our estimates suggest K h = O(1)m 2 s −1 in the fuel layer (the lower 1 m above the ground). In the low-intensity prescribed burn experiments when the fire front was passing the observational towers, the data suggests that the vertical diffusivity still decreases with height, highlighting the importance of eddy dynamics in the fuel layer.
Using a simplified 1D model and applying different wind conditions, we demonstrate that the RoS depends on horizontal diffusivity-larger horizontal diffusivity results in larger RoS. This shows that background turbulent fluctuations in the fuel layer affect the RoS of the low-intensity fires. The applicability of the study is to the burning fuel layer only, not to the layers above where wind and turbulence are stronger. Limitations of the model are that it is 1D and the wind and fuel models are idealized. The main idea of using this simplified model is not to simulate a particular fire but to show that lateral diffusivity plays a fundamental role in spread, at least for zero, low, and moderate winds.