Zonal-Based Emission Source Term Model for Predicting Particulate Emission Factors in Wildfire Simulations

A physics/chemistry-based numerical model for predicting the emission of fine particles from wildfires is proposed. This model implements the fundamental mechanisms of soot formation in a combustion environment: soot nucleation, surface growth, agglomeration, oxidation, and particle fragmentation. These mechanisms occur on a scale too fine for the discretization of most wildfire models, which need to simulate landscape-scale dynamics. As a result this model implements a zonal approach, where the computed soot particle distribution is partitioned into process zones within a single resolved grid cell. These process zones include: an inception zone (for nucleation), a heating zone (for coagulation, surface growth, and fragmentation), a reaction zone (for oxidation), and a quenched zone (for atmospheric processes). Governing mechanisms are applied to the appropriate zones to predict total particle growth and emission. The proposed model is implemented into HIGRAD/FIRETEC, a physics-based wildfire simulation code which couples interactions between fire, fuels, atmosphere, and topography on a landscape scale. Fire simulations among grasslands and conifer forests are performed and compared against experimental data for emission factors.


Introduction
Particulate formation and emission is an increasingly important area of research to the fire sciences and management. A major constituent of smoke emissions is soot particles which are believed to be the bulk of emitted fine particles (less than 2.5 microns in diameter). To predict amount and effects of emitted soot particles, there are two important stages to model: source and downwind regimes. The source regime consists of the origin of soot particles, in particular predicting how much soot is emitted directly from the flames. The downwind regime is the behavior of soot particles in smoke plumes, this includes transformation of particles (particle coagulation, species' condensation to particle surfaces, etc), transport of * Correspondence should be addressed to: Alexander J. Josephson, E-mail: alexanderj@lanl.gov Fire Technology, 57, 943-971, 2021 Ó 2020 The Author(s) Manufactured in The United States https://doi.org/10.1007/s10694-020-01024-7 particles (advection, deposition, etc), and interactions between particles and the atmosphere (light absorption/scattering, ice nucleation, etc).
Whereas there has been decades of modeling focused on the smoke particulate processes occurring in smoke plumes or the downwind regime [1][2][3][4][5], modeling of smoke particulate processes in the source regime is still rudimentary at best. This modeling is predominantly in the form of experimentally-derived emission factors, mass of emitted particles over mass of consumed fuel, to determine mass of soot emitted from a fire. In principle, when evaluating an individual fire, researchers would correlate an emission factor (found in an inventory of emission factors) to aspects of the fire environment, such as fuel type, as well as ambient and fire-induced local conditions. To this regard, there has been an extensive effort to expand the existing emission factor inventory [6][7][8][9][10], but the fundamental challenge with this approach is that there are so many variables affecting emissions, which include not only the conditions present before the fire but those created by the fire, that it is impossible to match all fires exactly to an inventory no matter how extensive. To this regard, the current work seeks to develop a first-pass model that predicts a particle emission factor of fire based on the fundamental chemical and physical mechanisms that govern soot formation and emission, which to the authors' knowledge has never been produced before.
These fundamental mechanisms for soot formation and emission have been researched by the combustion community for decades. And while there have been some detailed and sophisticated models developed to predict soot formation [11][12][13][14][15][16][17], the vast majority of work has been focused on soot formation in gaseous flames. Unfortunately, fires regularly occur in solid-fuel systems and the early stages of soot formation vary greatly in these systems from gaseous fuel systems [18]. With regard to solid fuel systems, there are very few models developed, only a few intended for coal systems [19,20] and even fewer for biomass systems [21] where we'd expect most fires to be occurring. All the above cited soot models require a high resolution (millimeters to centimeters) to implement as the models require local values of temperature and chemical species concentrations resolved on a combustion scale.
There have been a plethora of numerical models for wildland fire behavior [22][23][24][25]; however, soot formation and emission in wildland fires has not been a focus of any of these models. Rather, focus has always been on fire behavior particularly rates of spread on landscape scales. Soot particulate emission has been modeled in some wildfire models, such as Porterie et al. [26], however, interest has been more on soot particles' role in radiative heat transfer and fire spread, thus, modeled soot emissions have not been compared to experimental data.
One of the challenges in the numerical study of soot emissions in wildland fires comes from the scale differences between soot formation process and fire behavior. Wildfire simulation models intended for simulation of fires at landscape scales (HIGRAD/FIRETEC [27], WRF-Fire [28], WFDS [29], etc.), resolve flow dynamics at much larger scales (meters to kilometers) to capture fire-atmosphere coupling with limited computing resources, which translates to the inability to resolve details of the processes controlling soot-formation phenomena [30]. This work seeks to overcome this challenge for modeling soot emissions from wildland fires for use in the context where processes are resolved at meter scales.
We, the authors, do not claim the model produced in this work would predict particulate emission from wildland fires better than measured smoke yields in conditions similar to conditions of existing emission factor inventories; but instead serve to associate soot production values with simulated fire environment metrics (cell-wise conditions available in the model calculations). With the developed smoke model we are able to 1) simulate scenarios (both past and future) where no data is available and still obtain an effective estimation of smoke emissions and 2) perform fast/economical simulations to inform relations and theories (such as used for smoke management scenarios) which can then be validated through a reduced number of more expensive but now targeted experiments.

Model Development
The basis of this model comes from a previously developed computationally inexpensive soot formation model [31] for use in simulating large solid-complex fuel systems at the combustion scale. That model evaluates three variables: number density of tars, bulk mass density of soot particles, and particle number density. Each of these variables are influenced by soot formation phenomena, and that model evaluates the most common of these phenomena: soot nucleation, surface reactions, and coagulation.
Unfortunately, the previously developed model was designed to solve the mechanisms of soot formation at the micro-combustion scale, on the order of millimeters to centimeters depending on the system. As stated in the introduction, most developed wildfire simulation softwares are not designed to resolve combustion kinetics (such as soot formation) at this scale as the computational economics of such a resolution would prohibit landscape scale domains. For instance, HIGRAD/FIRETEC, the wildfire simulation software used in this work, executes simulations on domains of 100s to 1000s of meters, with resolutions on the order of 2 m at the ground level. The emissions model of this work was developed to account for the aggregated net meter-scale impacts of processes occurring at the soot formation scales without resolving their details.
To accommodate the scale differences between the soot source processes and those that can be resolved in landscape-scale simulations, we apply a zonal approach. This approach partitions individual flames into a series of smaller zones based on distinguishing mechanisms affecting particle characteristics. This proposed zonal model partitions the processes of soot formation into four zones within a simulated grid cell. These zones include the inception, heating, reaction, and quenched zones. The generation and emission of soot particles is computed as a cumulation of the processes occurring in each zone. One key assumption, is that at each given timestep the mechanisms in each zone are in a pseudo-steady state. This means that we assume the flame to be saturated with soot particles at every time-step and we are computing the portion of those particles that are emitted at the given time-step based on current flame structure and chemistry.
In the following sections, the partitioning of soot generation into the processes occurring in different zones based on flame characteristics and the fundamentals of soot formation applied within each zone will be presented along with the overall approach for determining a local particle emission factor given these processes.

Flame Characteristics
To determine the extent of each zone within a localized combustion environment and thus the extent of soot formation within each zone, characteristics of the flame are evaluated.
An average flame length in each cell (d F ) is determined using a correlation developed by Heskestad [32], Hesk À 1:02 ð1Þ D eq is the diameter equivalence of the consumed fuel bed if all the fuel was in a circular pattern. For the constant N Hesk : c p is the specific heat of air, T 1 is the ambient temperature, DH c is the fuel's heat of combustion, Z st is a stoichiometric coefficient, and _ Q is the heat release rate. This correlation was originally created by Heskestad to predict the luminous height of a turbulent diffusion flame but was reconfirmed by Newman and Wieczorek [33] to give accurate predictions of the chemical flame height in surface fires. It is the chemical flame height which is important to soot formation and particulate emission. This average flame length encompasses the inception, heating, and reaction zones.
The reaction layer thickness is determined using the correlation presented by Bilger [34], where D is the average molecular diffusivity of gaseous species through the reaction layer, c is the turbulent strain rate (defined as the turbulent dissipation divided by the kinematic viscosity), k is the chemical rate of oxidation of fuel, m 1 and m 2 are stoichiometric mass coefficients, and R b is a scalar defining the diffusive coupling of reacting species defined by Bilger [34]. With locally resolved gas velocities taken from the parent software, HIGRAD/ FIRETEC in the case of this study, residence times in both the flame and in the reaction layer, rt F and rt R , are computed as simply the gas velocity multiplied by the average length or thickness (either flame or reaction zone).
The last needed flame characteristic is an entrainment of air in the flame structure. In concept, when fuel is pyrolyzed there is some air that is sucked into the root of the flame. This air is entraped in the structure of the flame and shortens the average flame length. With this entrained air, small amounts of oxidizer are introduced within the structure of the flame itself. To compute the entrainment (g) we compare the above computed average flame length (d F ) against an ideal maximum flame length of a laminar diffusive flame (d Ã ) with an equivalent fuel consumption rate [35], where Q f is the volumetric flow (m 3 s À1 ) of volatile gases coming off the fuel, D 1 is the diffusivity (m 2 s À1 ) of oxygen in the surrounding air, and T f /T H are the average temperatures (K) of the burning fuel and flame, respectively.

Inception Zone
In the inception zone, biomass is pyrolyzed and from the products of this pyrolysis soot particles are nucleated. Pyrolysis is a thermochemical conversion of biomass into combustible products that include tars. Tars are large hydrocarbons which condense if cooled to standard temperature and pressure. These tars typically contain some sort of aromatic ring within their structure [36] and are the primary soot precursors from fires burning a solid fuel [21]. To find the fraction of volatiles that are tar, as well as the size of that tar, we executed ensembles of CPDbio (Coal-Percolation model for Devolatilization, Biomass adaptation), a detailed network devolatilization model designed to predict the products of biomass pyrolysis [37], using thousands of simulations to sample over various pyrolysis conditions expected in a wildfire scenario. These conditions consisted of temperatures ranging from 800 K to 2500 K and pressures ranging from 0.5 atm to 1.5 atm. From these thousands of executions, we were able to create an easy-to-evaluate relationship between parent-fuel and product tar under wildfire conditions, a relationship referred to from now on as the 'sooting potential'.
Biomass is made up of four basic components, cellulose (y cell ), hemicellulose (y hemi ), lignin (y lig ), and extractives; each of which has a different tendency to pyrolyze and produce tars at different yields and sizes. CPDbio predicts the pyrolysis behavior of cellulose, hemicellulose, and lignin independently, with hemicellulose and lignin further broken up into softwood and hardwood types. Extractives vary greatly species-to-species with some extractives possibly producing large amounts of tar and others producing almost no tar [30]. In this work, extractives are lumped in with the other species as a source of tar, because at this point in model development we are seeking to create a generalized model for particulate emissions for all fuel species which is currently impossible with the variation in extractives, without more species' specific information. In the future, when more information becomes available, we hope to refine this capability.
The results of the performed simulations are summarized in Table 1. This table shows the mean and standard deviation of tar mass yield, as a fraction of the pyrolyzed parent fuel, and the averaged size of the tar molecules of the ensemble of CPDbio simulations. For the most part, there was very little variance in these values with the exception of the tar mass from softwood lignin. To form the sooting potential we averaged softwood and hardwood components together thus forming a simple correlation, for estimating the tar yield (y tar ) and average tar mass (m tar ) pyrolyzed from the fuel.
Using this sooting potential we are able to find the bulk density of tar produced in the inception zone where q Vol represents the observed density of volatiles during pyrolysis (approximated as 0:45 kg m À 3). The (1 À g) factor is present because no tar is expected to be contributed from entrained air which was defined in Eq. (4). Equation 8 only computes the tar number density in the inception zone assuming the tar species are unreactive. In reality, these tar species are highly reactive and subject to thermal cracking and soot nucleation at very fast rates. To find the actual tar concentration, we consider three chemical processes of tar formation and consumption, tar inception (r TI ), thermal cracking (r TC ), and soot nucleation (r SN ) derived in the earlier work [31] with details available in Appendix A; here we will just note that the dependency of each term on N tar . If we assume tar concentrations to be in a pseudo-steady state (dN tar ¼ 0) then we can solve for the actual tar number density using the quadratic formula, giving us the real tar concentration within the inception zone. Using this solution for the actual tar number density, we can compute rates of tar cracking (r TC ) and soot nucleation (r SN ). These rates are used to compute the number concentration of soot particles exiting the inception zone and entering the heating zone, The second fraction of the above equation is the ratio of tar molecules nucleating to soot particles versus being cracked back into lighter gases. The factor 4 comes from an approximation of four tar molecules needed to coagulate to form an incipient soot particle stable enough to not thermally crack at flame temperatures [38]. To compute the mass density of these particles, we simply multiply the number of incipient particles by the weight of those particles, N SH and M SH are the soot number and bulk mass densities exiting the inception zone and entering the heating zone.

Heating Zone
The heating zone of a flame makes up the bulk volume of a flame. In this zone, volatile gases undergo a large number of transformations as gases heat up and mix with entrained air. From a soot formation perspective, there are three major types of reactions occurring which we model: particle coagulation, surface reactions, and fragmentation. By solving the instantaneous rate of coagulation (See Appendix B) analytically, assuming a free-molecular regime, and integrating across the particle residence time in the heating zone, % rt F , we are able to compute the total extent of particle coagulation within the heating zone where e is a Van der Waals enhancement factor of 2.2, q soot is the density of soot particles (1850 kg m À3 ), and k B is the Boltzman constant (1.3806E23 J K À1 ). We use rt F , the average residence time of the flame, as it is assumed the total volume of the inception and reaction zones are negligible compared to the heating zone. Surface reactions within the heating zone include surface growth, via the HACA (Hydrogen Abstraction followed by Carbon Addition) mechanism (k HACA ) [39], gasification (k gas ) [40], and oxidation (k oxi ) [40]. Computing rates of these reactions, requires temperatures for the heating zone (T H ) which can be computed from a calibrated correlation shown in Table 2. Chemical partial pressures of particular gaseous species, also needed for computing reaction rates, are interpreted from a chemistry surrogate table provided in the supplementary information. These tables have averaged concentration values of gaseous species in the heating and reaction zones tabulated by set values for air entrainment and the mass fraction of O 2 in the surrounding air. For details of how the temperature surrogate correlations of Table 2 and the chemistry surrogate tables of the supplementary information were constructed please refer to Appendix C.
By solving the instantaneous rate of surface reactions (See Appendix B), and integrating across rt F we are able to compute the cumulative change in bulk soot mass density across the heating zone In order to analytically solve the instantaneous rate of surface reactions and integrate, we have to approximate a constant number density of particles (N S ). In this zone, changes in N S are only governed by coagulation, and instantaneous rates of coagulation have a squared dependency on N S implying rates of change will decrease exponentially as particles coagulate. Because of this, we expect averaged values of N S to be much closer to N 0 SR than N SH . Therefore, we use N 0 SR as an approximate constant for N S with which we can obtain Eq. (14).
We've found through repeated simulations that surface growth mechanisms (through the HACA mechanism) plays a minimal role in particle mass concentrations in flames from solid complex fuels which release high quantities of tars, that in turn, dominate the mass production of soot particles. However, in the heating zone there is a non-negligible amount of oxidation, and gasification, which occur due to air entrainment. This oxidation also causes particle fragmentation which we model using a scheme developed by Sirignano [41] which is again solved analytically using the updated values of N SR and M SR computed from coagulation and surface reactions where m C is the molecular mass of carbon, v is the density of activation sites on the soot particle surface (1:7 m À2 ), and d p is the diameter of primary particles which make up a soot aggregate, assumed to be 50 nm in this work. N SR and M SR are the particle number and bulk mass densities exiting the heating zone and entering the reaction zone.

Reaction Zone
The reaction zone is a thin outer rim of the flame where oxidation reactions dominate and the heat of a flame is produced. In this layer, oxidation and gasification rates are computed as before but with a reaction temperature from Table 2 and chemical partial pressures interpolated from the supplemental chemistry surrogate tables.
As the reaction zone is thin compared to the heating zone, we assume there is inadequate time for significant amounts of particle coagulation to occur, and thus the overall particle number density does not change through the reaction zone Oxidation occurs heavily in the zone and affects the bulk mass density greatly Above, S SR is the effective bulk surface area of particles available for oxidation, and is computed by a correlation developed by Balthasar and Frenklach [42] m 0 is the mass of the basic incipient soot particle, which in this model is 4m tar from Sect. 2.2. d is a shape factor for soot particles describing the morphology of soot particles on a range from 2/3 to 1, where 2/3 minimizes the surface area to volume ratio of particles (spherical) and 1 maximizes that ratio [42]. In this study a surrogate model was developed to predict the shape factor of particles based on particle size, Above, d 2=3 is the average diameter of particles (in nm) if they were spherical and is computed For the derivation of Eq. (19), see Appendix D. N QR and M QR are the particle number and bulk mass densities exiting the reaction zone and entering the quenched zone.

Quenched Zone
Once particles pass into the quenched zone, they are emitted to the atmosphere.
To compute the mass of particle emissions (M PE ) we compute an emission factor and multiply this by the mass of raw fuel consumed by fire. The emission factor is a multiplication of two terms EF ¼ŷ tarŷsoot : ð21Þ y tar , was defined in Sect. 2.2, andŷ soot is a mass fraction of the tar, which nucleates to soot and survives to be emitted.ŷ soot is computed simply as a ratio of the mass density of tar in the inception zone to the predicted mass density of soot emitted to the quenched zonê To obtain the number of particles emitted to the atmosphere, simply multiply the total mass emission by the ratio of particle number to mass computed entering the quenched zone M PE and N PE are the total mass and particle number emitted from a fire to the surrounding environment. For use in a CFD model, these values should be divided by the volume of computations cells to which these particles are emitted to obtain intrinsic values of bulk particle number and mass densities which can transported in any CFD scheme.

Model Implementation
While not part of the model itself, this section describes how the model presented in the previous sections was implemented into HIGRAD/FIRETEC to inform the reader of aspects and challenges to its implementation which may help researchers with future implementations into other platforms. HIGRAD is the atmospheric CFD model [43] and FIRETEC is the fire physics model [44] used in this study. In implementing the zonal-based emission model we've rewritten HIGRAD to transport (with advection and diffusion considerations) two additional variables: N, the local bulk number density of particles, and M, the local bulk mass density of particles. Future works include the evolution of particles as they continue to collide and transform within a fire's plume, but for now we're focused on the source regime. FIRETEC uses a kinetic-mixing model [44] for predicting a rate of fuel consumption at a high temporal resolution (typically on the order of a millisecond) in each cell where a portion of the solid fuel's temperature is above the critical threshold (600 K). In cells where fuel is being consumed, this proposed model is implemented on a cell-to-cell basis.
Within an ignited cell, the total amount of fuel (m fuel ) consumed within a given timestep is used along with a material density (q fuel ) and sizescale (ss fuel ), typically 500 kgm À3 and 0:5 mm respectively, of the consumed fuel, to reconstruct the diameter equivalence used in Eq. (1), Using this D eq and N Hesk from Eq. (2), we compute the average flame length in each cell where active burning occurs. This flame length represents the extent of the heating zone and when divided by the local gas velocity resolved in HIGRAD/FIRETEC gives us rt F which is all we truly need with regard to the extent of the heating zone. It is possible, and very likely, that the length of the flame extends past the boundaries of the computational cell. In FIRETEC, burning is done locally, meaning the conversion of fuel to products of combustion happens within the same cell as the fuel consumption with no consideration of the flame envelope. Likewise for now, particles are emitted locally within the cell of burning rather than at a distance of the flame length.
With the rate of fuel consumption translated to a chemical rate of oxidation (simply consumption rate over mass consumed) and an estimated flame temperature of 1250 K along with other models for turbulent dissipation [45], Eq. (3) is used to compute the reaction zone length, again translated to rt R by dividing the local gas velocity. This gives us the extent of the reaction zone local to this cell. For now, the average molecular diffusivity of gaseous species through the reaction layer is assumed to be 6.25E-4 (m 2 s À1 ) derived from the one-dimensional flame simulations used to produce the chemistry surrogate tables in Appendix C.
Local entrainment is computed using Eqs. (4) and (5) with a stoichiometric coefficient of 0.1775 (switch grass) and an estimated flame temperature of 1250 K. With this entrainment value we are able to update the heating and reaction zone temperatures using Table 2 and the chemical species' partial pressures using the tables constructed in Appendix C and provided in the supplementary materials.
Sections 2.2-2.5 accurately describe the calculation of the emission factor within a given cell at a single time step. Once this local emission factor is calculated, it is multiplied by the consumption rate of raw fuel, not including consumed char, to provide the source rate term for M, bulk mass density of particles, for transport by HIGRAD. Source rate term for N, the bulk number density, is computed through Eq. (23).

Results and Discussion
To begin validation of the proposed model, a series of simulations were carried out and compared to data collected by Urbanski et al. [6]. In that work, Urbanski et al. reported near-field average emission factors for several species and dozens of fires. Fires were classified by ecosystem, where a rough Gaussian distribution of emission factors was computed for each ecosystem. From these data, emission factors for PM 2:5 from two ecosystems were extracted for comparison with the proposed model.
The two fire ecosystems used from Urbanski et al.'s work were: (1) grass/shrublands of mid-western and southeastern United States and (2) interior west mountain conifer forests of United States and southwestern Canada. For both of these fire ecosystems, a representative fuel map was obtained and three simulations were performed containing logarithmic wind profiles with 2, 6, and 10 m s À1 velocities 10 m above the ground.
The simulations were executed using the HIGRAD/FIRETEC CFD software [44], with a domain size of 800 Â 400 m discretized into 2 Â 2 m horizontal cells and 1.5 m vertical resolution near the ground with a parabolic vertical stretching component. A 100 m fireline was initialized 100 m into the simulation domain. All simulations were performed on flat topography.

Grasslands
The default fuel map for a HIGRAD/FIRETEC simulation is a flat heavily-fueled grassland field with a uniform fuel distribution of 1 kg m À3 at 5% moisture content and a height of 0.7 m [44]. This default setting was used to compare against Urbanski et al.'s grass/shrublands ecosystem. Figure 1 shows a visualization of the HIGRAD/FIRETEC linefire simulations in a uniform grassland. The density of fuel is shown in green and is initially uniform, elevated ( > 500 K) gas temperatures are shown in red, and the bulk mass of emitted particles (M) are shown on a grayscale from 0 (white) to 1E-5 (black) kg m À3 . As can be seen in the figure, this fire is relatively low intensity and thus emitted particles remain relatively near the ground. In addition, due to the mass of fuel consumed and low intensity flames, quantities of particulate emissions are generally low. Figure 2 shows the comparisons between model-predicted emission factors and experimental data, where the black dashed line represents a normal distribution of expected emission factors reconstructed from the mean and standard deviation of field data collected by Urbanski et al. for the grass/shrubland ecosystem. The vertical lines represent the cumulative emission factors from the low, mid, and high wind simulations up to 340 s of simulation time. This image shows that the proposed model, or more specifically its implementation in HIGRAD/FIRETEC, consistently over-predicts the particulate emissions factor for the grassland fires, however all simulations fell within (or close to) a single standard deviations of the field datas' mean, showing an agreement in trends and acceptable value agreements for a model to which there is no equivalent at this point.
An interesting trend arises in the simulation data. Figure 2, shows that emission factor and wind velocity are correlated but not linearly proportional. At the low wind case, 2 m s À1 , we get the highest predicted emission factor (16.5), and as velocity increases to 6 m s À1 the predicted emission factor (13.6) also decreases. However, at the highest velocity, 10 m s À1 , the predicted emission factor (10.9) continues to drop but less so then before. This behavior is due to interactions between wind and fire structure.
Hosseini et al. [9] noted a higher emission factor from a heading fire line than from a backing fire line, this is generally extrapolated to mean that higher intensity fire tends to have higher emission factors. Higher wind velocities tend to create more intense fires, thus we'd expect to see increasing emission factor. However, higher winds also increase the entrainment of air, and thus oxidizing species, into the fire structure. This entrainment leads to greater rates of particle oxidation and thus lower the emission factor. The nonmonotonically decreasing emission factor with wind speed illustrated in Fig. 2 shows just a small part of this play between entrainment and fire intensity as driven by wind velocity, and becomes more pronounced in canopy fires.

Conifer Forests
While the grasslands provide a base comparison between the proposed model and experimental data, wildland fire commonly occur in other fuels. For a broader picture of how this model does with various fuel types, simulations were performed using a simulated ponderosa pine forest based on data collected as a part of a fire/fuel surrogate study near Flagstaff, Arizona, which was led by Carleton Edminster of the Rocky Mountain Research Station of the United States Forest Service. Ponderosa pines reflecting the measured tree inventory were placed within the computational domain along with a moderate groundcover of grass and litter distributed according to the proximity of trees following a procedure described by Reference [46]. Figure 3 shows a visualization of the HIGRAD/FIRETEC simulations for a linefire in this ponderosa forest. The three images show the smoke plume at 140 s of simulation time for three simulations, from left to right: 2, 6, and 10 m s À1 winds at 10 m above the ground. Figure 4 is similar to Fig. 2 in that it gives a comparison between field/experimental data and the HIGRAD/FIRETEC simulations. The figure shows that in a forest case the HIGRAD/FIRETEC implementation of this model matched field data very well, with all three simulations within (or close to) a single standard deviation of the mean field data. Consistent with the field data, emission factors from the conifer forest fires are generally higher than the grassland fires most likely due to the higher intensity of the fires. As with the grassland simulation, Fig. 4 shows a snapshot of the interplay between fire intensity and air entrainment by the order of wind velocities and emission factors.
Comparing Fig. 2 to Fig. 4, we see that the proposed model responded well to the difference in ecosystems with an increase in emission factor as we moved from the grassland to the conifer forest. Perhaps a primary reason for this increase in emission factor comes from the canopy fire. Canopy fires have greater access to winds which extend out the flames causing longer flames, compared to ground fires, in which particles have longer times to coagulate and grow larger, resulting in less oxidation in the reaction zone and therefore higher emission factors. However, the grassland fires did not respond as much to changes in wind velocity as the conifer forest fires.
In the conifer forest, the lowest emission factor (12.7) was with the lowest wind, 2 m s À1 , which might be expected according to Hosseini et al.'s observations. With an increase in wind, 6 m s À1 , came an increase in emission factor (22.1), again not surprising as the increasing wind led to a high intensity fire which became a consistent canopy fire whereas the low-wind fire only occasionally canopied. How- ever, increasing the wind speed more, 10 m s À1 , brought a decrease to the emission factor (16.6) due to the increased entrainment of air into the structure of the canopy flames.

Conclusions and Final Comments
A physics and chemistry-based model for predicting a particulate emission factor and particle sizes from wildland fires was developed. This model implements mechanisms developed in a previous model but applies a zonal approach to adapt to a coarse grid. By partitioning the hypothetical combustion zone into the various flame regions characteristic of all flames, this model is able to apply the mechanisms of soot formation in a pseudo-steady state to predict the overall emission of particles, both in mass and size. Through this partitioning, we are able to partially apply decades of detailed model development for soot formation [13][14][15][47][48][49] to the particulate emission source problem of fire behavior. Simulations were performed using HIGRAD/FIRETEC and compared against experimental data. Emission factors from simulations were compared against field data from fires in different ecosystems. The implemented model was able to capture trends in soot emission between different ecosystems, and the resulting predicted emission factors were within 1 standard deviation of field data and thus within bounds of uncertainty. However, the simulation-predicted emission factors were on the fringes of the field data's uncertainty bonds indicating room for improvements. The presentation of this zonal model and global comparisons of emission factors to field experiments represents a first-pass approach to tackling the problem of predicting particulate emission factors using a physics and chemistry-based model. While continued validation and model refinement are certainly needed to improve emission predictions, this model, or an earlier version of it, has already been used to explore difference in emission profiles from fires under critical (low moisture fuels and high winds) and marginal (high moisture fuels and low winds) conditions [50]. That previous work also contained a sensitivity analysis of the proposed model, revealing this model's strong dependencies on local oxygen concentrations and combustion rates of reaction, with weaker dependencies on turbulent dissipation and gas density. That performed sensitivity analysis led to model improvements reflected in this current work but further analysis should be continued.
To the authors' knowledge, the proposed model is the first of its type, a physics and chemistry based model for predicting fine particulate emissions in a wildfire scenario. While we acknowledge the model certainly has deficiencies and should undergo further verification and validation, we hope that this work would serve as a starting point for researchers to expand and refine in this important aspect of smoke modeling.
where k B is Boltzmann's constant, and d tar , the effective diameter of the tar, can be computed using a geometric relationship [53] Although the rate of the second coalesence reaction should theoretically vary from the first as the diamers are now twice as large, we assume this rate to sufficiently small that we can approximate computed reaction to be the same.

Appendix B: Heating and Reaction Zone Chemical Rates
This appendix defines the rates of soot formation mechanisms used in the heating and reaction zones of this work.

Particle Coagulation
Soot coagulation occurs through the collision of two spherical soot particles which stick and form one larger sphere conserving mass. The rate of particle coagulation, assuming a free-molecular regime, is given by the frequency of particle collisions [31] dN soot dt Here b T represents a frequency of collision between tars and e is a steric factor, the van der Waals enhancement factor, with a value of 2.2 [52]. where d soot is the effective diameter of the soot particles

Surface Reactions
Surface reactions affect only the mass of particles in the following way, r SS is the rate of surface reactions (kg=m 2 s) for which we consider two types of surface reactions: surface growth, through the hydrogen-abstraction-carbon-addition mechanism (HACA), and consumption through oxidation and gasification both of which are given in detail by Josephson et al. [31].

Particle Fragmentation
Particle fragmentation occurs when the weak point of an aggregate, where two primary particles are attached, is attacked via oxidation, or gasification, thus breaking the particle into two. We model particle fragmentation using a model developed by Sirignano et al. [41], where m C is the molecular weight of carbon, v is the site density of activation sites on the surface of the soot (1.7E19 1=m 2 ), and n p is the number of primary particles within the soot aggregate, where d p is the assumed diameter of a primary particle, in this study 50 nm. In Sirignano et al.'s work Eq. (34) was n p À 1 instead of n p , however we make the assumption that n p > > 1 which simplifies the math considerably and allows us to analytically solve 34.

Appendix C: Derivation of the Temperature Surrogate Models and Chemistry Surrogate Tables
This appendix depicts how the surrogate tables for predicting species concentrations were created. These tables are a result from an executed series of one-dimensional flame simulations with different air entrainment and surrounding oxygen depletion values. A one-dimensional flame is simulated by linearly varying the mixture-fraction (n) from 1 to stoichiometric. Mixture fraction is an expression of the fraction of mass coming from the fuel stream m f , as opposed to the oxidizer stream m o . In this scenario, the fuel stream is the mass of fuel volatiles resulting from pyrolysis plus the mass of air entrained at the root of the flame. Air entrainment values were varied from 0 (no air) to 1 (all air). The oxidizer stream is surrounding air with possible oxygen depletion. Oxygen mass fractions in the surrounding air were varied from 0 (all N 2 no O 2 ) to 0.233 (ambient air). A mixture-fraction/temperature mapping was extracted from a DNS simulation performed by Lignell et al. [54]. These DNS simulations were of an evolving nonpremixed ethylene jet flame; as ethylene is often used as a surrogate fuel for biomass volatiles [55], we presumed these DNS results to be viable for this application. The extracted mapping was used to determine the temperature profile of each one-dimensional flame with respect to the linear mixture-fraction profile.
For the fuel of biomass we used a dry switch grass species with an elemental analysis of 52.3% carbon, 6.3% hydrogen, 40.5% oxygen, and 0.5% nitrogen [56]. This switch grass species is balanced stoichiometrically into parts CO, CO 2 , CH 4 , and C 2 H 4 . Using the GRI 3.0 mechanism [57], an optimized detailed chemical reaction mechanism designed for natural gas flames, and Cantera [57], an objectoriented software toolkit for chemical kinetics, thermodynamics, and transport processes, we equilibrated each mixture fraction to the mapped temperature and an atmospheric pressure. This allowed us to estimate chemical concentrations at each mixture fraction for each of our species of interest mapped within the mixture-fraction space: OH, O 2 , CO 2 , H 2 O, C 2 H 2 , H 2 , and H.
To translate the mixture-fraction space to the zonal-based system used by this model, the heating zone is considered to be 1:0 > n > 2n st and the reaction zone is considered to be 2n st > n > 0:5n st , where n st is the stoichiometric mixture fraction. In the case of no entrainment, where there is no air in fuel stream, n st ¼ 0:175.
The average concentration of each species of interest is computed for each air entrainment and O 2 mass fraction value giving the data points which make up the chemistry surrogate tables which are available in the supplementary material. Temperatures did not depend on the air's O 2 mass fraction, only on the air entrainment value. Thus we were able to fit a simple 4th-order polynomial, to the temperature data. In this surrogate, Eq. (37), g is the entrainment value on a scale from 0 to 1, and T i is the predicted temperature of zone i. The resulting calibrations are shown in Table 2.  Figures 5 and 6 show examples of a single contour from the chemistry surrogate tables (where O 2 mass fraction was ambient) and the fitted temperature correlations for the heating and reaction zones respectively.

Appendix D: Derivation of the Shape Factor Surrogate Model
Balthasar and Frenklach [42] derived the concept of the shape factor for soot particles to compute an effective surface area of particles at which surface reactions may take place. The shape factor of a particle (hdi) describes the morphology of a particle on a scale from 2/3 to 1, where 2/3 minimizes the surface area to volume ratio of particles (spherical) and 1 maximizes that ratio (an aggregate of connected but not overlapping incipient particles). This model has been used by a variety of researchers [21,42,58] in many of simulations, and while shape factor does depend on the formation history of particles [58] when comparing the three cited references we found that there was a very close correlation between absolute particle size and shape factor. Using data from Balthasar and Frenklach's work we were able to fit a logarithmic function to estimate shape factor hdi ¼ a lnðd p þ bÞ þ c ð38Þ to the extracted data. This fit resulted in Eq. (19) and the fit to extracted data can be seen in Fig. 7 (R 2 ¼ 99:92). To verify this fitted function, the shape factor requires that an incipient particle be spherical, hdi = 2/3, with this model a particle of 10 nm, which is approximately the size of our incipient particle, has a shape factor of 0.6612, a little low but not unreasonable. On the other end, a particle can never have a shape factor above 1. This function crosses the x axis at 77,000 implying a particle must have an equivalent diameter of 77 lm which is much larger than anything we could see in these simulations.

ELECTRONIC SUPPLEMENTARY MATERIAL
The online version of this article (https://doi.org/10.1007/s10694-020-01024-7) contains supplementary material, which is available to authorized users. . Shape factor data extracted from Balthasar and Frenklach [42] and fitted with Eq. (19) (R 2 ¼ 99:92). On the x-axis is the diameter of particles if they were spherical, essentially a mass of particles, and the y axis is the shape factor of the particles.