Approaches to Modeling Bed Drag in Pine Forest Litter for Wildland Fire Applications

Modeling flow in vegetative fuel beds is a key component in any detailed physics-based tool for simulating wildland fire dynamics. Current approaches for drag modeling, particularly those employed in multiphase computational fluid dynamics (CFD) models, tend to take a relatively simple form and have been applied to a wide range of fuel structures. The suitability of these approaches has not been rigorously tested for conditions which may be encountered in a wildland fire context. Here, we focus on beds of Pinus rigida needle litter and undertake a two-part study to quantify the drag and evaluate the capabilities of a multiphase large eddy simulation CFD model, the NIST Fire Dynamics Simulator. In the first part, bed drag was measured in a wind tunnel under a range of conditions. The results were fit to a Forchheimer model, and the bed permeability was quantified. A traditional approach employed in the multiphase formulation was compared to the parameterized Forchheimer equation and was found to over-predict the drag by a factor of 1.2–2.5. In the second part, the development of a velocity profile above and within a discrete fuel layer was measured. Using the Forchheimer equation obtained in the first part of the study, the CFD model was able to replicate a qualitatively consistent velocity profile development. Within the fuel bed, the model appeared to under-predict the velocity magnitudes, which may be the result of unresolved pore-scale flow dynamics.


Introduction
A strong understanding of the underlying processes that drive wildland fires is necessary for providing robust, science-based solutions to the challenges they presented to society. Computational Fluid Dynamics (CFD) modeling of wildland fire dynamics offers one valuable avenue for improving this understanding. Studies using CFD tools have shown promise in bringing new insights to a variety of topics related to the dynamics of wildland fires, or at the very least have highlighted areas that require continued scientific investigation. Examples include the importance of fuel structure across a range of scales (Dahale et al. 2013;Pimont et al. 2011), the importance of flame structure (Linn et al. 2012;Frangieh et al. 2020), and the role of char oxidation (El Houssami et al. 2016;Perez-Ramirez et al. 2017). As these models evolve, it is important to ensure that the underlying submodels are formulated appropriately and are tested rigorously against relevant experimental data (Mell et al. 2010;Morvan 2011).
Currently, a common approach to incorporating the influence of wildland fuels (vegetation) into more established CFD frameworks is known as the multiphase formulation. A number of specific CFD models can be found today, but they generally trace their origin to the work of Mell et al. (2010), Grishin (1997) and Larini et al. (1998). In this approach, vegetation is treated as an assemblage of subgrid-scale particles, typically assumed to be thermally-thin, and their influence is spatially averaged into a porous media. Therefore, processes which fundamentally occur at the scale of individual particles must be submodeled, generally employing a combination of simplifying assumptions and empiricism, a number of which have yet to be rigorously verified.
Two such processes which are of interest are the momentum transfer through fluid-solid drag interactions and the heat transfer through the fluid-solid convective interactions. These processes are linked and are fundamental to modeling fire behavior. They will control the thermal degradation of the particles (and thus their ignition) and the flow of oxygen to both gas-phase (flaming) and surface (smoldering) combustion.
In this paper, we focus on the submodeling of drag forces. Vegetation drag has been studied in a number of different contexts (e.g., Lai 1955;Rudnicki et al. 2004;Nepf 2012;Fehrmann et al. 2017); however, in this particular case we are interested in flow through relatively densely packed layers of dead foliage, representative of forest litter. Understanding this flow is important for quantifying the behavior of surface fires (which are dominated by spread in this fuel layer). These can represent the early stages of highintensity wildfire behavior, but are also found in the more low-intensity fires typically employed by forest managers in prescribed burning activities (Clark et al. 2014). Further, a detailed investigation of flow in these types of fuels is intended as a complimentary effort to ongoing research into the characteristics of fire behavior in identical fuels (Mueller et al. 2018;Campbell-Lochrie et al. 2021).
In order to test the ability of the multiphase approach to reproduce flow in pine forest litter, two sets of experiments were carried out. In the first case, drag was quantified using a wind tunnel to measure pressure drop through fully packed beds of pine needles. The bulk density (and thus bed porosity) and velocity were varied to isolate the influence of these variables. This was compared to numerical simulations of the same configuration using a traditional drag formulation based on past studies, and an alternative formulation is proposed.
In the second case, the flow profile above and within a layer of forest litter was measured. This was carried out by filling only the lower half of the wind tunnel test section. Velocities were recorded at a number of locations, and the effects of varying bulk density and upstream velocity were again investigated. Complementary numerical simulations were carried out, but in this case the drag force formulation obtained in the first set of experiments was used as an input for the model. In this way, it was possible to confirm whether the formulation derived from simple pressure drop measurements is applicable to more complex scenarios.

Particle Formulation
In the multiphase approach, as implemented for simulating wildland fire behavior, drag within the fuel bed is represented by adding a locally averaged body force term, ⟨F d ⟩ V , to the momentum equation in volumes for which vegetation is present. Typically, this assumes that drag originates from individual particles in cross-flow (i.e., bluff bodies), and the bulk volumetric force is obtained by summing the contributions of all particles, n, within a control volume, V: where is the surface-to-volume ratio of a single particle, is the solid volume fraction in the control volume, c d is the particle drag coefficient, c s is a shape factor, is the fluid density, and u is the fluid velocity. In this case, the angular brackets, ⟨ ⟩ V , represent a superficial average over the control volume and will be subsequently dropped for simplicity. The shape factor originates from the projected area ( A p ) and is determined based on the assumed particle geometry and orientation. For cylindrical particles, oriented perpendicular to the flow, this becomes c s = 1∕ . Note that a shape factor using a spherical assumption has also been employed for models of fire spread in various vegetative fuel beds (e.g., Morvan and Dupuy 2001;Dahale et al. 2013;Perez-Ramirez et al. 2017). However, assuming a cylindrical particle, the empirical drag coefficient in Eq. 1 has the following velocity dependence (Schlichting and Gersten 2016): where Re p = 4u∕ is the particle Reynolds number, where is the fluid kinematic viscosity. This approach takes advantage of the fact that empirical drag coefficients for simple geometries have long been well-established. Further, and are typically known for simple wildland fuel layers (or at least may be estimated). Thus, the model does not require a-priori tuning.
In reality, particle orientation is likely to be more complex and thus complicates the shape factor determination. The model also does not consider more dense particle packing, such as with forest litter layers. In this case, wake interactions or the formation of interconnected pore structures (with associated tortuosity) may also need to be considered in drag modeling (Kyan et al. 1970). Vegetation flexibility is also ignored. This is known to play a role in reducing the velocity dependence of drag for shrubs (Cao et al. 2012) and trees (Rudnicki et al. 2004;Vollsinger et al. 2005), which streamline in strong flows. In packed beds, it has been suggested that deflection of fibers may contribute to drag through elastic storage of energy (Kyan et al. 1970), but this depends on particle properties and may only play a role in high-flow conditions (Macdonald et al. 1979). Thus flexibility is likely only important for litter layers insofar as determining structure of dense packing. For these reasons, it is necessary to test the performance of this approach. (1)

Packed Bed Formulation
While this work is focused on packed beds of forest litter, we can seek alternative approaches in other applications related to flow in porous media (found in hydrology, process engineering, etc.). In this case, the modified Reynolds number, Re m = u∕(1 − ) , may be used to evaluate the flow regime, where is the bed porosity ( = 1 − ) (Holdich 2002). This ratio is similar to the particle Reynolds number ( Re p ); however, that quantity does not consider the effect of interacting flow around multiple particles. For the same particle size and flow conditions, decreasing porosity would decrease Re m (increase the importance of viscous forces), but Re p would remain the same, making the former more suitable to classify flow in packed beds. For low modified Reynolds numbers, the pressure loss in porous media is dominated by viscous dissipation and depends linearly on velocity. This is the well-known Darcy's Law: Here pressure loss (Pa m −1 ) is written as a force per-unit-volume ( F d , N m −3 ), and k is the Darcy permeability. When the modified Reynolds number becomes greater than 2, inertial effects must also be considered (Holdich 2002). These can be incorporated with an additional quadratic term, resulting in the Forchheimer equation: where Y is an empirical coefficient, and k is the permeability. Note that permeability here is sometimes considered to be separate from the Darcy permeability in Eq. 3, as fitting these equations to the respective flow regimes may yield different values; however, differences tend to be negligible (Skjetne et al. 1999;Balhoff and Wheeler 2009;Chai et al. 2010). The physical source of the quadratic term has been debated (flow separation, kinetic energy changes due to tortuosity, etc.) but has been demonstrated through both theoretical and empirical analysis (Macdonald et al. 1979;Li and Engler 2001).
Dividing the second term by the first term on the right-hand side of Eq. 4 yields the socalled Forchheimer number (Zeng and Grigg 2006): which gives the ratio of inertial to viscous contribution to the drag. The case of strong inertial dominance behaves like a bluff body and strong viscous dominance collapses to a pure Darcyian flow. It has been suggested that the flow can be approximated as Darcyian only if the inertial contribution is less than 10% ( Fo < 0.11) (Zeng and Grigg 2006).
An alternative, non-dimensional form of the Forchheimer equation can also be written: Tanino and Nepf (2008) tested emergent cylinders (as a surrogate for aquatic vegetation) over a range of porosities and found that this equation predicted the drag behavior in the range 30 < Re p < 700 . They proposed empirical values for 0 and 1 , which showed a dependence on . However, given the empiricism, it is unclear how well the values can be extrapolated to other bed structures. It is apparent that applying this approach to modeling the drag in a fire behavior model requires a-priori knowledge of permeability and the coefficient Y (from which 0 and 1 can also be derived, or vice versa). This requires experiments such as the ones performed in this study, as there is currently no well-established universal relationship between permeability and more easily observed parameters (i.e., , or ). Therefore, the approach of Eq. 1 (which doesn't require such knowledge) has been more typically applied in modeling fire spread in wildland fuels, and we return to the question of which approach more accurately describes drag in litter fuel beds.
As mentioned previously, flow through packed beds has been investigated in a number of contexts, extending well beyond wildland fires. In particular, Forchheimer-like equations have been used to describe pressure drop in a wide range of geological material (sand, gravel, etc.). Take for example, the often-used Ergun equation (Ergun 1952), which takes the same functional form as Eq. 4, but uses an estimate of permeability based on particle and bed properties. These estimates are empirically based and have been found to span a wide range (e.g., Macdonald et al. 1979). Further, the parameters of interest for forest litter in fire-like conditions (high porosities compared to other granular material, and typically high Reynolds numbers) are not well-covered by past research, and the results cannot be confidently applied without confirmation of applicability through wind-tunnel tests.
While far fewer in number than other topics related to porous media, there have been some studies quantifying the pressure drop through similar beds of forest litter material. These include investigations on different components of broadleaf litter, at varying states of decomposition (Wang et al. 2019), as well as pine needles (Simeoni et al. 2012;Santoni et al. 2014;Ventura et al. 2002;Fehrmann et al. 2017). However, results have been conflicting with regard to the flow regime (Darcy or Forchheimer) and have yet to be extended to a modeling framework or used in the investigation of more complex flow scenarios, as is the aim here.

Experimental Approach
The experiments were carried out in a stainless steel wind tunnel, designed to generate heated flows characteristic of conditions encountered in wildland fires. The general design of this apparatus, named the Wildfire Exposure and Ignition Response Device (WEIRD), is shown in Fig. 1. Air is drawn from the surrounding environment and an associated pressure increase supplied by a blower (Leister Silence). This air is then forced through a process heater (Leister LHS 61L) capable of elevating the temperature up to 650 • C. For the tests described in this report, only ambient temperatures are used. Air is then passed through a divergence section, a settling chamber with a honeycomb straightener, and a contraction section which was designed following a 3rd-order polynomial matched with a 7th-order polynomial for the outlet (Morel 1975;Su 1991). An enclosed test section, with dimensions of 750 mm × 150 mm × 100 mm, is attached to the contraction outlet.
Within the test section, porous fuel beds were comprised of dead pine needles of the species Pinus rigida (common name: pitch pine). These have been the focus of a number of previous studies on the mechanisms of fire spread in forest litter (e.g., El Houssami et al. 2016; Campbell-Lochrie et al. 2021) and the particle properties have been extensively characterized (Thomas 2016). The surface-to-volume ratio for a single needle was taken to be = 4661 m −1 (giving an equivalent cylindrical diameter of d = 4∕ = 0.86 mm), and the effective material density (assuming a moisture content of 10%) was e = 650 kg m −3 . Figure 2 shows the typical appearance of individual particles. Note that pine needles are clasped together, in species-specific quantities, on one end by a thin structure known in botany as a fascicle. In the case of Pinus rigida, there are 3 needles per fascicle; however, most species have between 2 and 5 needles per fascicle (Harlow and Harrar 1941).
For the first set of experiments, a 700 mm length of the test section was filled completely with pine needles and static pressure was measured at four locations along the bed (see Fig. 3a). For this basic configuration, the measured pressure drop can be attributed to the effect of the porous media (whether in the Darcy or Forchheimer regime), and thus can be equated to the drag force per-unit-volume ( F d ) used to represent the effect in the multiphase model: Static pressures were recorded with digital pressure transducers (Sensirion SDP800-125). These cover a range of ± 125 Pa, with span accuracy of 3% and a zero-point accuracy of 0.08 Pa. The superficial velocity (i.e., the bulk velocity neglecting local variations within bed pores) was measured upstream at the tunnel centerline with a ø 7.9 mm S-type pitot tube (Dwyer 160S-18), connected to the same type of pressure transducer. The combination of pitot tube and transducer was calibrated in an open tunnel for the range of velocities investigated.
For a given packing of the test section, pressure drop was recorded at eight different superficial velocities, nominally at 0.25 m s −1 intervals between 0.25 and 2.0 ms −1 . At a given velocity, pressure data were recorded at 10 Hz for 120 s to ensure a steady state was captured.
The pine needles were distributed in the test section in a quasi-random manner, while attempting to maintain macro-scale homogeneity in the test section volume. The total mass of vegetation within this volume was adjusted to achieve a range of bulk densities ( b ) from 10 to 60 kg m −3 . The consistency of the packing method was improved by sub-dividing the mass of fuel before packing, so that an equal portions of fuel were distributed in the upper and lower sections of the bed, and this was supplemented with careful visual inspection. For b <30 kg m −3 , a loose fine-wire steel supporting mesh was inserted at the mid-plane height of the test section in order to help maintain a quasiuniform distribution of fuel across the full 100 mm depth. This range was selected both to cover the range of beds which could reasonably be constructed in the test section with homogeneous structure and to overlap with conditions studied previously (e.g., Fehrmann et al. 2017). Data on realistic values for Pinus rigida litter density are limited, though values on the order of 10-40 kg m −3 can be inferred from measurements of loading and height Gallagher et al. 2017).
Four test repetitions were carried out at each bulk density (with a complete re-packing of the test section) in order to capture any variability in bed packing technique.  Fehrmann et al. (2017) identified the fact that the arrangement and orientation of the fuel elements can have an impact on permeability, with manufactured litter layers behaving differently from natural ones. Thus, it is acknowledged that in-situ pine litter may yield slightly different results, but this technique for creating the bed structure closely matched that which has been used for flame spread experiments in the same fuel (Campbell-Lochrie et al. 2021), and so the results are still highly relevant.
For the second set of experiments, the test section was filled only to a height of 50 mm, over a length of 300 mm (see Fig. 3b). Tests were run with superficial velocities at 0.5 m s −1 increments between 0.5 and 2.0 m s −1 . The local velocity was measured at 33 points, distributed along five different vertical transects in the tunnel centerline. The measurements were made with a ø 8 mm hotwire anemometer (Kimo CTV 210-R), which has an accuracy of 3% ±0.03 m s −1 up to 3 m s −1 , and 3% ±0.1 m s −1 beyond this. The probe was moved to each point sequentially and left to record for 60 s at 10 Hz. In this case, the bulk densities tested were 20, 40, and 60 kg m −3 , and these packings were repeated twice (the reduced repetition was considered acceptable given the highly consistent results found in the first set of tests, as discussed in Sect. 4.1.1). In all cases, a loose fine-wire steel mesh was placed at the mid-plane height to help create a well-defined upper surface to the fuel bed.

Numerical Approach
CFD simulations, mirroring the experimental configurations, were carried out using the National Institute of Standards and Technology (NIST) Fire Dynamics Simulator (FDS), version 6.7.3 (McGrattan et al. 2019b). This employs a Large Eddy Simulation (LES) approach to solving a low-Mach number formulation of the Navier-Stokes equations. In this work, subgrid-scale turbulence was estimated with Deardorff's eddy viscosity model (Pope 2000). This choice was guided primarily by the fact that it is the default approach in FDS, and the focus here was not on evaluating turbulence models. FDS also includes additional numerical models for incorporating gas-phase combustion, thermal radiation, solid fuel degradation, particle transport, etc.. Subgrid-scale vegetation is incorporated through source and sink terms in the gas-phase conservation equations. As such, the variation in quantities such as velocity at the scale of individual particles or pores is not directly resolved. Further, a continuum approach is used, where the full volume is considered available for gas phase flow (that of the solid phase is ignored). This implies superficial rather than intrinsic averaging. The focus of this work is on the momentum sink term which represents the bulk effect of vegetation drag on the gas phase ( F d ).
Further details of the particular CFD framework are outside of the scope of this discussion, but the specifics of the model formulation can be found in the technical documentation (McGrattan et al. 2019a). FDS was selected primarily as it has the potential to simulate fire dynamics in wildland fuels over a range of scales, including in similar fuel beds (e.g., Mueller et al. 2018). Therefore, its ability to model the drag at this scale is of particular interest.
Simulations of the experiments were conducted using three-dimensional domains with the same dimensions indicated in Fig. 3. In the case of the pressure drop simulations, the numerical grid resolution was x = y = z = 10 mm. For the layer flow simulations, the resolution was x = y = z = 5 mm. The upwind yz-boundary was set to maintain a fixed flow, using either the target superficial velocities for the pressure drop cases or the experimental values from the upwind-most velocity profile for the layer flow cases. The downwind yz-boundary was set to a zero pressure gradient (open outflow), and the lateral boundaries were set to solid surfaces.
For solid boundaries, FDS employs a wall model to approximate the near-wall velocity gradient when the boundary layer is unresolved. The walls were treated as smooth, with the model following the work of Werner and Wengle (1993). This approach determines the tangential velocity component in the first off-wall cell based on whether it is located in the viscous sublayer (linear profile) or not (logarithmic profile). A preliminary investigation indicated that implementing a wall roughness to modify the log law region (following Pope (2000)) did not significantly alter the results, with the fluid mechanics being dominated by the presence of the fuel bed. Further details of the wall model can be found in McGrattan et al. (2019b).
Vegetation, in the form of a bulk drag force, was incorporated in all grid cells corresponding to the shaded volume in Fig. 3. The one exception was the case of b = 60 kg m −3 in the layer flow configuration. The high density of needles used in the experiment pushed up against the wire mesh at the upper surface of the fuel bed, resulting in a bed depth of 55 mm. Thus, this more realistic value was used in the corresponding simulations.

Pressure Drop Measurements
Example pressure drop measurements, showing four repeats of single bulk density condition, are given in Fig. 4. Measurements were steady over the full 120 s duration of a given test, and the average relative standard deviation for a single measurement is 4% ( n = 768 ), with higher flow conditions (higher pressure readings) resulting in even lower variability. The example in Fig. 4 shows that individual tests had a linear pressure drop, as anticipated x is the distance from the upwind edge of the fuel bed by Eq. 7. Across all conditions and repeats, average r 2 for linear fit on a single test is above 0.99 ( n = 192 ). Further, the accuracy of the pressure measurements is confirmed by the fact that linear fits show pressure returning to the background value as the flow exits the fuel bed ( x = 700 mm). Finally, the results are shown to be highly repeatable. For example, a single linear fit of all the data points in Fig. 4 still gives r 2 > 0.99 . This indicates that local variability in bulk density and particle orientation due to individual re-packing of the beds did not significantly bias the results.
The time-average pressure drop values, as determined from the slope of a linear fit as in Fig. 4, are shown for all experiments in Fig. 5. The data have been fit with a second-order polynomial, which is forced to zero at u = 0 , giving the same form as Eq. 4. The results again demonstrate high repeatability for given bulk densities (the 90% prediction intervals for the polynomial fit are shown as a shaded area). Deviation from the fit can most likely be attributed to the upstream measurement of superficial velocity, which had more noise than the static pressure measurements. A Forchheimer-type formulation gives a good representation of the velocity dependence on bed-drag for these pine needle layers under the range of Reynolds numbers considered. This can be seen in the small error values shown in Table 1 as compared with the range of pressure drops measured.

Flow Regime
The permeability, k, and modified Forchheimer coefficient, Y, for the beds in question are given in Table 1. These were derived from the polynomial fit for Eq. 4. Combining them through Eq. 5, the range of Forchheimer numbers can be plotted in order to evaluate the flow regime (Fig. 6). The region where viscous effects are stronger than inertial effects ( Fo < 1 ) only corresponds to a very small portion of the conditions encompassed by the parameter space of this work. Indeed, for high Reynolds numbers and high porosities, the inertial effects approach a dominance by a factor of 10. Figure 6 also highlights the fact that no immediately clear relationship was observed between the drag behavior and porosity. The differences between curves are nonuniform and non-monotonic (e.g., the = 0.923 case has the strongest viscous effects despite having the second lowest porosity). This may be influenced by uncertainty in the second-order fit (see the increased prediction intervals in Fig. 5 for higher bulk density), but may also have some more physical basis in the way flexible pine needles arrange themselves in a bed as they become increasingly compact. Quantifying the latter influence is beyond the scope of this work; however, the fact that different curves are obtained at all demonstrates that the modified Reynolds number does not fully encompass the different behaviors in these beds.
If the curves in Fig. 6 are extended to the proposed threshold for Darcian flow of Re m = 2 , it is found that Fo varies between 0.01 and 0.03. These values are below the alternatively suggested threshold of Fo = 0.11 (Zeng and Grigg 2006), indicating that the criteria of up to 10% inertial contribution is less strict than the modified Reynolds number. Put  . 6 Ratio of viscous to inertial effects for the parameter space covered in the experiments. The relationships are determined from Eq. 5 and Table 1 another way, extending to Fo = 0.11 in Fig. 6 would suggest that Darcyian flow extends to higher Reynolds numbers in our case ( Re m = 9 to 20). Nevertheless, under no test condition was it found that the contribution of inertial effects was small enough to be neglected (minimum Fo = 0.46 found for = 0.923).
The quadratic behavior qualitatively agrees with the work of Erić et al. (2011) for straw, Fehrmann et al. (2017) for pine forest litter, and Wang et al. (2019) for broad leaf forest litter. Conversely, both the study by Simeoni et al. (2012) and the continuation of that work by Santoni et al. (2014) found a linear relationship, and thus suggested that the drag should be modeled as a purely Darcyian regime. Their tests were all reported to be in the range of Re m > 100 , which is above the theoretical threshold for which inertial effects become important. However, this study focused on flow vertically through a column of needles. This may have influenced needle orientation relative to the flow, which has been shown to affect permeability (Fehrmann et al. 2017). Further, reference velocity was recorded downstream, and the influence of local velocity variability at the outlet of pore spaces on this measurement is unknown. Given the small number of studies focused on similar materials, there is a need for future work to explore the mechanistic relationship between bed structure and flow regime.

Permeability
The magnitudes of permeability obtained here are between 1 × 10 −5 m 2 and 7 × 10 −7 m 2 . These are generally high when compared to those observed for artificially constructed layers of Pinus radiata needles by Fehrmann et al. (2017). For example, they reported a permeability of 4.53 × 10 −7 m 2 for = 0.96 , only 16% of the value obtained in this study (inferred from the data in Table 1). However, they also showed a strong sensitivity of permeability to porosity above values of 0.92, suggesting that experimental uncertainty may play a role in the discrepancy. Indeed, for a = 0.95 , they reported an increase in permeability to 1.59 × 10 −6 m 2 , despite the slight decrease in porosity. This is 81% of the estimated value for an equivalent porosity in this study. Santoni et al. (2014), while finding that flow in their case was purely in a Darcyian regime, also report values around an order of magnitude lower. For example, for Pinus laricio needles (which have a surface-to-volume ratio similar to Pinus rigida), permeability was 2.65 × 10 −7 m 2 at a porosity of = 0.95 . This is 14% of the value obtained here. These differences may be traced to alternative experimental designs, as discussed in Sect. 4.1.2).
The relationship between permeability and bed porosity is summarized in Fig. 7. Some suggested approaches exist for approximating this relationship. One commonly found in porous media literature, which has been previously applied to forest litter, is to use the Carman-Kozeny equation (e.g., Santoni et al. 2014;Wang et al. 2019): which gives the permeability as a function of porosity, element surface-to-volume ratio, and a constant, k k which accounts for particle shape and pore tortuosity. Here, we obtain a good fit with k k = 9 , but tend to over-predict permeability at high porosity and underpredict at low (Fig. 7). It should be noted that this formulation was developed for Re m < 2 ; however, permeability, as determined for the viscous term in Eq. 4, is generally assumed to take the same functional form when considering higher Reynolds numbers (Ergun 1952;Holdich 2002). It has also been shown in a numerical studies that the difference in permeability between the Darcy and Forchheimer permeabilities is small, and so the Carman-Kozeny equation may produce a reasonable estimate across flow regimes (Balhoff and Wheeler 2009;Chai et al. 2010).
Following the approach of Tanino and Nepf (2008), a similar formulation for permeability is found by combining Eqs. 4 and 6: It was suggested that 0 = 25 for a porosity of ∼0.9 (but found to increase with decreasing porosity). An average value of 0 = 6.8 was found by fitting the data in this study to Eq. 6, and the resulting permeability prediction is shown in Fig. 7. In this case, the fit is worse than the Carmen-Kozeny equation, and the effect of porosity on permeability for pine needle beds is under-represented.
Equation 9 shows an inverse linear dependence of permeability on (recall that = 1 − ). On the other hand, 3 ∕(1 − ) 2 in Eq. 8 approaches 1∕ 2 (an inverse square dependence) for close to 1, such as in this study. The experimental data suggest that this stronger dependence is more appropriate, particularly in the mid-range of porosities investigated here. Indeed, while examining lower Reynolds-number cases, Kyan et al. (1970) found k k ≈ 9 for = 0.95 in fibrous media, but that this parameter increases significantly close to = 1 . This can help explain the over-prediction of Eq. 8 at high porosity when a constant k k is used. At low porosity, the bed structure may change considerably due to flexible arrangement of the needles as they are packed more densely, again suggesting that a constant k k may only apply to a narrow range of porosity.
A final noteworthy point is the inverse dependence of these permeability formulations on surface-to-volume ratio. Increasing characteristic element size while maintaining porosity should thus increase permeability. This can help explain differences in permeability between this study and those with finer needles (e.g., Fehrmann et al. 2017). Likewise, the fact that needles are not independent but grow in groups on a single fascicle may be important. For example, the Pinus laricio studied by Santoni et al. (2014) tends to have two needles per fascicle, while Pinus rigida has three. Total bed porosity does not account for this local structuring, and species which have greater clustering

Fig. 7
Dependence of permeability on bed porosity. × 's represent values obtained from polynomial fit of all experimental data and error bars are derived from the 90% confidence interval on the fit may have larger pore spaces (similar to decreasing while preserving ), resulting in increased overall permeability.

Pressure Drop Simulations
The results of the numerical simulations of the first set of experiments are shown in Fig. 5 and compared to the theoretical predictions of Eqs. 1 and 7. This comparison is presented to verify that the modeling framework, which includes a full numerical solution of the conservation equations in three dimensions, approximated boundary conditions, etc., collapses to the simple theory in this case. This is confirmed by the fact that the maximum relative error between any one simulation and the theory is 1%. Of greater significance is the clear tendency for the formulation of Eq. 1 to over-predict the experimentally determined drag force. This over-prediction ranges from roughly 1.2 to 2.5 times the experimental data. Given that the initially assumed shape factor of c s = 1∕ does not capture the random and dense arrangement of particles, a poor fit is not unexpected. Further, Combining Eqs. 1 and 2 (recognizing that 15 < Re p < 115 ) shows that the bluff body approach gives a velocity dependence of the following form: This results in a slight discrepancy from the quadratic velocity dependence, which was found to fit the experimental data well.
Nevertheless, an alternative approach to modeling each scenario with the permeability coefficients determined in Table 1 is to propose a new representative shape factor. This was done by minimizing the relative error between the experiments and Eq. 1. A value of c s = 0.16 produces both the lowest mean and cumulative relative error magnitude.
The predictive performance of this optimized value is shown Fig. 8. The mean relative error magnitude is 12% when comparing all test cases (Reynolds numbers and porosities) (10) F d = C 1 u + C 2 u 1.8 . The fact that the error is also a function of shows that the approach fails to fully capture the drag dependence on bed structure. If the bed is assumed to behave as a collection of independent particles, their average projected area relative to the flow increases with bed packing, resulting in the decrease in relative error, from positive to negative, as a function of porosity. The notable exception is the non-monotonic behavior at the lowest porosities. These features point to anisotropic bed characteristics which depend on the packing. The simple approach of a single shape factor ignores this by nature, but such complexities may be important for future studies.

Velocity Profile Development
Having characterized the bed drag, it is necessary to evaluate the ability of the numerical model to capture the development of a flow profile within and a above a pine litter bed of finite height (Fig. 3b). For this investigation, Eq. 1 is replaced in the model with the more appropriate Eq. 4, employing appropriate values from Table 1. An example comparison of the velocity profiles for the experiment and simulation, at discrete locations along the test section centerline, for an upstream velocity of u 0 = 1 m s −1 , and for the three different bulk densities, is shown in Fig. 9. In the context of this discussion, the bed height (5 cm) is referred to as h.
In the case of the experiment, the upstream velocity profile ( x∕h = −1.5 ) is shown to be uniform along height. The flow was also laminar, with a maximum turbulence intensity of 11% at 0.5 m s −1 , and 4% for the higher velocities. As the flow approaches the leading edge of the fuel layer ( x∕h = −0.5 ), a marginal adjustment can be observed, with decreasing velocity below layer height and increasing above, particularly for increasing bulk density. Just downwind of the leading edge ( x∕h = 0.5 ), the flow profile becomes very irregular as a function of height. This can be attributed to the adjustment of the flow and the fact that the layer does not actually have a well-defined upwind face due to its highly porous nature. Further downstream, the flow establishes into a classic profile for flow over a porous canopy (such as a forest) (Finnigan 2000). A clear influence of bulk density (and correspondingly, porosity) can be identified, particularly at x∕h = 4.5 , with higher intra-bed velocities for low bulk density and vice versa above the layer.
The simulations show a very similar development of the velocity profiles. It is expected that upstream-most profiles match, as the experimental data at this point were used to define the numerical boundary condition. As the simulated flow establishes within and above the fuel bed (e.g., x∕h = 4.5 ), it adopts the same qualitative profile as the experiments. An exception is the region of flow adjustment just downwind of the leading edge ( x∕h = 0.5 ). As opposed to the experiment, the profiles appear more established at this point (there is a clear relationship between local velocity and bulk density), though they are still developing in the downstream direction. This can be traced to the fact that, unlike the experiments, the simulations have well-defined edges to the fuel layer, with vegetation becoming suddenly and uniformly present for all grid cells located at x∕h > 0 and y∕h < 1. Figure 11 gives a more direct comparison between the experimental and numerical data at x∕h = 4.5 . In this case, both experimental repeats are included, and it can be seen that the profiles show consistency. This confirms again that local variability in bed structure (due to repacking) did not dominate the results. The simulations also appear to match the quantitative experimental data relatively well. This is particularly true for in the space above the fuel bed ( z∕h ≥ 1 ) for b = 40 kg m −3 , for which the mean relative error in the predicted velocity is 7%. For b = 20 kg m −3 there is more of a discrepancy, which may be linked to the lack of a clearly defined leading edge in the experiments, and the agreement was improved by selecting a sampling location further upstream in the simulations (suggesting the flow established more quickly in the numerical case). For b = 60 kg m −3 , the discrepancy just above the bed may be linked to uncertainty in the true layer height for this dense packing, as discussed in Sect. 3.2. Still, the mean relative error for the upper three measurement points is only 9%.

Intra-Bed Velocity
While the simulations appear to perform well above the fuel bed (e.g., Fig. 11b), the systematic tendency for under-prediction within the fuel bed is of interest. This is demonstrated, across the range of velocities, for an example bulk density of 40 kg m −3 in Fig. 10. By normalizing the velocity values by the mean intra-bed velocity ( u b ), it can be seen that the experiments show peak values above the bed that ranges between 5 and 7 times the intra-bed velocity. The simulations, on the other hand, have peak values between 8 and 19 times the intra-bed value. It is worth noting that, despite the differences in scale, in both cases higher upstream velocity magnitudes tend toward an equalizing ratio of above-to in-bed velocity.
The upstream velocity profiles match very well (again, this is expected, as the experimental measurements at x∕h = −1.5 were used to set the upstream numerical boundary condition), and therefore, the amount of mass entering the domain is the same in both cases. Mass conservation dictates that the same mass flow rate should be found at the yz-plane at x∕h = 4.5 , but the discrepancies in the normalized profiles indicate otherwise. Direct evaluation of the differences in mass flow is difficult, as the shape of the flow profile close to the walls is important for obtaining the integrated volume flow. Unfortunately, it was not possible to resolve this in detail with the hotwire probe used. It was possible, however, to confirm the conservation of mass in the simulations, and it must also be conserved in the experiments. Therefore, the source of the intra-bed velocity discrepancy should be considered further. A possible explanation for the larger measured velocities within the fuel bed is the fact that the presence of the vegetation reduces the area of the section available for flow at the downstream location. This is not accounted for in the simulation, which applies a superficial volume average and ignores the volume of the solid phase, from the perspective of the gas phase. Using the bulk porosity values, the bulk interstitial velocity should be higher ( ∼10% in the 60 kg m −3 case) to achieve the same mass flow rate in the reduced available area (Holdich 2002). This fails to explain the close to 100% difference between the experimental and numerical values for intra-bed velocity. However, bulk porosity does not fully describe the characteristics of the bed, and it cannot be assumed that all free area implied by the porosity is available for flow. There will also be a distribution of velocity, as local boundary layers will be created individual pores, so maximal values will tend to increase beyond the ∼ 10% estimated by assuming plug flow (rather than say, Poiseuille flow).
Some research has focused on resolving the distribution of interstitial velocity in porous media, using non-invasive measurement techniques (Huang et al. 2008;Datta et al. 2013;Holzner et al. 2015). For example, Huang et al. (2008) used refractive index matching to apply image tracking techniques for measuring flow in a packed bed of spheres. They found a distribution of interstitial velocities, with maximum values nearly five times the average flow rate. These studies cover a different range of Reynolds numbers and porosity than is relevant for this work, but qualitatively similar behavior should be expected. The fact that the hotwire measurements are uniformly higher, and generally insensitive to repeat packing of the bed may be attributed to the probe design. The measurement wire is housed in a hole which passes through the ø 8 mm support rod, thus creating a local void space which is repeatable in all scenarios. To gain greater insight into vegetative fuel beds, future efforts should explore alternative, ideally non-intrusive, measurement techniques.

Canopy Analogy
In order to evaluate and contextualize the observed shapes of the velocity profiles, a comparison can be made to past work on flow through terrestrial and submerged canopies of vegetation. A parameter typically encountered in such studies is the leaf area index, a, which gives proportion of solid frontal area in a canopy volume (Kaimal and Finnigan 1994), analogous to c s . In terrestrial canopies, the shape factor is often taken as c s = 1∕2 , while Nepf (2012) suggested a value of c s = 1∕ for submerged canopies of cylindrical vegetation (as was derived for Eq. 1). Here, we take the newly proposed average shape factor of c s = 0.16 . Integrating a along canopy height gives the ratio of frontal area to bed area, and for an assumed height-invariant structure this is equal to ah. The fuel layers here have ah = 1.2 to 3.8. Belcher et al. (2008) suggested two regimes of canopy, sparse and dense: with the latter occurring for ah >> 0.1 . Thus, this study is within the dense regime. For such canopies, an inflection point appears in the velocity profile around the canopy top and Kelvin-Helmholtz vortices are produced by shear, similar to a mixinglayer (Raupach et al. 1996). In this case, we observe the inflection, but there appears to be little turbulence generated (this can be seen by the small standard deviations in Fig. 9, implying low turbulence intensity). Nepf et al. (2007) suggested the length scale of the vortices generated at the canopy top should scale as ≈ 0.23∕c da . Taking c d ≈ 2 , from Eq. 2, gives ≈ 2 to 5 mm. Thus, the vortices approach the size of individual elements (d), at which point a limit of ≈ 2d has been suggested . Consequently, the turbulence generated will be of a very small scale and likely poorly captured by the measurement system here. Moving at 1 m s −1 , a 2 mm vortex would pass a sensor in just 2 ms. Furthermore, the shear layer will potentially only be fully developed for x∕h > 10 ( Dupont et al. 2011), though it has been suggested this distance may scale with canopy density (Belcher et al. 2008;Ghisalberti and Nepf 2009). Therefore, vortices closer to the canopy edge will be even smaller.
For dense canopies, the length scale of shear layer turbulence penetration into the canopy (taken as ) can be used to define an upper and lower region. For ah > 1 , as here, the entire canopy can be considered in the lower region. Due to this small penetration, it is isolated from the flow above-which can be described by a typical logarithmic profile (e.g., Finnigan 2000) with a displacement height equal to the canopy height (Nepf 2012). In the lower canopy, the velocity profile will depend only on the density profile. This agrees well with the fact that no velocity gradient is observed below h, given the height-invariant density. Figure 10 demonstrates that, for a given overall density, the profiles within the bed collapse to the same shape when scaled against a reference velocity.

Perspectives
In order to contextualize the findings of this study, it is worth recalling the intended application space-the numerical simulation of fire dynamics in packed beds of vegetation, representative of forest litter layers. The first half of the study found that an often-used approach for simulating the drag in so-called multiphase models of wildland fire, which assumes a summation of idealized bluff body drag effects, tends to over-represent pressure drop through beds of Pinus rigida. This over-representation can be generalized by a factor of roughly 1.2-2.5 times. As discussed in Sect. 1, the drag will have important implications for fire spread modeling. The effect on flow rates will influence both the temperature of fuel elements-and thus the potential to achieve different states of thermal degradation (e.g., pyrolysis)-and the supply of oxygen to gas-and solid-phase combustion reactionsand thus the rate at which these occur, and the efficiency or 'completeness' of the reaction. Therefore, the determination of an appropriate drag model for these beds is a key step for modeling the fire dynamics.
It was found that a simple correction, by way of using c s = 0.16 , reduced the average error magnitude in drag force to 12%. While still large for some applications, this may be considered an acceptable improvement, given the high degree of uncertainty associated with predicting wildland fire behavior (Cruz and Alexander 2013;Hoffman et al. 2016). It is important to note that this proposed coefficient is centered around the conditions studied. There is a negative error bias for low porosities and positive for high. The apparent cause is a missing representation of the changes which arise in element arrangement as porosity is modified. These findings support the overall importance of element orientation, as discussed previously (and highlighted by Fehrmann et al. 2017). Future work should focus on systematically evaluating the effects of varying orientation, as well as developing a better understanding of the bed structures which may be found in realistic forest conditions.
The use of the Forchheimer equation was explored as an alternative approach. This was found to yield more accurate results, as it is derived directly by fitting the experimental observations. However, given the range of behavior observed for pine needle beds in this and other studies, we cannot yet confidently model the permeability without fitting (e.g., with Eqs. 8-9) to an accompanying experimental study. As this is not feasible, or desirable, for the full range of fuel layers which may be encountered in wildland fire scenarios, future work must focus on illuminating the fundamental drivers of permeability so that drag may be reasonably modeled for a wide range of realistic beds.
The second half of the study demonstrates that a drag formulation determined from simple pressure drop measurements, input into a tool commonly used for modeling fire dynamics, is able to reproduce more complicated scenarios, such as flow over the edge of a porous fuel bed. This serves as an independent confirmation of the results from the first half of the study. However, despite the qualitative agreement between the experiments and the expected theory, the apparent discrepancy in intra-bed velocities between the experiments and simulations highlights an interesting unknown. The experiments suggest that the structure of the bed may lead to significantly elevated flows in pore spaces (and correspondingly low flows close to the surface of individual needles or in densely packed regions). The open question is then whether the potentially coarser description of the flow provided by the model is adequate for simulating problems related to combustion and fire spread. In terms of combustion, it is the overall mass flow rate which matters, and so it may be possible to neglect the highly local details of velocity variation within pore spaces. In terms of convection, the local velocity magnitude may be important, and a compensating factor may have to be introduced into the numerical model in order to capture the bulk behavior. This topic is the focus of ongoing research efforts.

Conclusions
The drag produced in packed beds of pine needles, representative of forest litter, was characterized through pressure drop measurements. A basic approach used for modeling this drag in studies of fire spread was assessed against the data. The characterized drag was then used to refine the modeling approach, and this was tested against velocity data gathered for flow within and above such a packed bed. The key findings were as follows: • The drag within beds of Pinus rigida needles, with bulk densities between 10 and 60 kg m −3 , for superficial velocities between 0.25 and 2.0 m s −1 , was found to be largely dominated by inertial effects, and well-described by the Forchheimer equation (Eq. 4) when parameterized using Table 1; • A classic approach for simulating this drag in multiphase models of wildland fire spread, assuming a simple arrangement of cylinders in cross-flow, was found to overpredict the drag by a factor of 1.2-2.5; • Implementing the Forchheimer equation, using parameters determined for these beds, enabled the model to successfully reproduce overall characteristics of velocity development above a fuel layer • An apparent under-prediction by the numerical model of the interstitial velocity measurements highlighted the fact that pore-scale velocity distributions are still poorly understood for this media and should be investigated in future research.
These new insights serve to improve our current ability to both understand and model the dynamics of fire spread in forest litter.