Wildfire smoke-plume rise: a simple energy balance parameterization

The buoyant rise and the resultant vertical distribution of wildfire smoke in the atmosphere have a strong influence on downwind pollutant concentrations at the surface. The amount of smoke injected vs. height is a key input into chemical transport models and smoke modelling frameworks. Due to scarcity of model evaluation data as well as inherent complexity of wildfire plume dynamics, smoke injection height predictions have large uncertainties. In this work we use a coupled fireatmosphere model WRF-SFIRE configured in large eddy simulation (LES) mode to develop a synthetic plume dataset. Using 5 this numerical data, we demonstrate that crosswind integrated smoke injection height for a fire of arbitrary shape and intensity can be modelled with a simple energy balance. We introduce two forms of updraft velocity scales that exhibit a linear dimensionless relationship with the plume vertical penetration distance through daytime convective boundary layers. Lastly, we use LES and prescribed burn data to constrain and evaluate the model. Our results suggest that the proposed simple parameterization of mean plume rise as a function of vertical velocity scale offers reasonable accuracy (30 m errors) at little computational 10 cost.


Introduction
Predictions of surface concentrations of wildfire smoke by regional and global chemical transport models depends on the initial equilibrium height of the smoke plume. Plume rise, which determines this equilibrium height, is widely recognized as an area of uncertainty (Goodrick et al., 2013;Paugam et al., 2016). Traditionally, many operational smoke modelling frameworks 15 relied on plume rise equations originally developed by Briggs (1975) for industrial smokestacks (Larkin et al., 2010;Pavlovic et al., 2016). Yet several studies suggest that this approach may not be appropriate for wildfires (Pavlovic et al., 2016;Heilman et al., 2014;Freitas et al., 2007).
In a recent review of existing plume rise parameterizations, Paugam et al. (2016) highlight three notable models that stand out in literature, as that of Freitas et al. (2007), Sofiev et al. (2012) and Rio et al. (2010). Both Freitas and Rio's approaches 20 use first principles to characterize plume temperature, vertical velocity and entrainment. While the former provides prognostic 1-D equations that can be solved as a stand-alone "offline" model, the latter is implemented as a sub-grid effect within a host chemistry transport model. Notably, both consider an idealized heat source to represent the fire. To initialize the plume at the lower boundary, simplified fire geometry (circular and rectangular for Freita's and Rio's models, respectively) with a uniform heat flux is assumed. Sofiev's semi-empirical approach relies on energy balance and dimensional analysis, while using 25 satellite data to both initialize and constrain the parameterization. Unlike Briggs's equations, all of the above models address wildfire plumes specifically, yet much research is needed to reduce the large uncertainties associated with the model predictions (Mallia et al., 2018). Moreover, it is unclear, whether unreliable predictions should be attributed to the fire input parameters or the plume rise model itself.
One of the central challenges in plume rise model development has been the scarcity of comprehensive model evaluation 30 data (Coen et al., 2012b;Ottmar et al., 2016). To date, information on wildfire smoke emissions and dispersion has largely been derived from two distinct sources: remotely sensed data and prescribed burn campaigns. While increasing numbers of satellite observations contribute to a more complete plume climatology (Val Martin et al., 2010), the data is subject to biases and lacks direct spatiotemporal links to fire behavior (Ichoku et al., 2012). In contrast, field campaigns, such as Prescribed Fire Combustion and Atmospheric Dynamics Research Experiment (RxCADRE) (Ottmar et al., 2016) and Fire and Smoke Model 35 Evaluation Experiment (FASMEE) (Prichard et al., 2019), provide the necessary level of detail for model validation studies.
However, such datasets typically capture a modest range of fire and atmospheric conditions. Our approach, therefore, is to develop a synthetic plume dataset that addresses the limitations of the available observational data. As vast majority of smoke plumes remain in or just above the atmospheric boundary layer (ABL) (Val Martin et al., 2010;Mallia et al., 2018), we use large eddy simulations (LES) to focus on local-and meso-scale plume dynamics. Using 40 a coupled model, we simulate a wide range of fire and atmospheric conditions (Sect. 2). Based on this synthetic LES data (hereafter referred to as "data") we propose a simple energy balance model for predicting plume rise of crosswind integrated (CWI) smoke from a non-uniform fireline (Sect. 3). We use the synthetic plume dataset to constrain and evaluate our plume rise parameterization. We then demonstrate with both numerical and prescribed burn data, that within the range of tested conditions this parameterization offers high speed and accuracy (Sect. 4). Moreover, it provides the means for classifying penetrative vs. 45 non-penetrative plumes, which is key for subsequent dispersion modelling (Sofiev et al., 2012;Val Martin et al., 2012).
The proposed approach is geared toward regional smoke modelling frameworks (e.g. BlueSky and BlueSky Canada). Government agencies, air quality managers and fire response teams depend on these operational tools and their accuracy to issue air quality warnings, evacuation orders and to help mitigate human health impacts. Yet, model evaluation studies suggest that plume rise estimation remains a weak link within smoke modelling systems (Raffuse et al., 2012;Val Martin et al., 2012;50 Chen et al., 2019). Moreover, existing methods struggle to reliably differentiate which plumes remain in the ABL and which penetrate it. The broad goal of the work is, therefore, to address some of these challenges and improve the accuracy of plume rise predictions for regional air quality applications.

Development of a Synthetic Plume Dataset
We devise a synthetic plume data set using a coupled fire-atmosphere model WRF-SFIRE, which combines the well-established 55 Weather and Research Forecasting Model (WRF) with a semi-emperical fire spread algorithm called SFIRE (Mandel et al., 2014;Mallia et al., 2020;Coen et al., 2012a;Kochanski et al., 2013;Clements et al., 2006;Kochanski et al., 2019). The model allows one to explicitly resolve plume dynamics, while parameterizing fuel combustion. One of the primary advantages of using WRF-SFIRE is that it supports two-way coupling between the atmosphere and the fire behavior model, allowing it to capture some of the complex dynamical feedbacks that exist between the fire and the atmosphere (Prichard et al., 2019). Heat 60 and moisture fluxes from the simulated burn provide forcing to the atmosphere, affecting local wind flow and thermodynamics. This in turn influences the modelled fire behavior. The following sections detail the numerical setup of WRF-SFIRE, scope of the dataset, as well as our approach to defining "ground truth" for model evaluation.

Numerical Configuration
WRF-SFIRE was configured in idealized large-eddy resolving mode. Much of our numerical setup was adopted from a case study of a real prescribed burn as detailed in Moisseeva and Stull (2019), to ensure the simulations represent physical conditions backed by model evaluation. Due to high computational demands of LES runs, we focused on the local-and meso-scales, considering only the initial buoyant plume rise of smoke in typical daytime atmospheres. Key parameters varied were ambient wind, fuel category, vertical potential temperature profile and fireline length, denoted as conditions W,F,R and L, respectively (detailed further in Sect. 2.2).

70
Each 10 km x 20 km domain with 40 m horizontal grid spacing was initialized with uniform ambient west wind W and vertical temperature profile R. Depending on the sounding R, the simulations were performed in either a shallow (3000 m) or a deep (5000 m) domain, with 51 or 71 hyperbolically stretched vertical levels, respectively. A constant uniform lower boundary surface thermal flux (tke_heat_flux) in the ambient environment and lateral periodic boundary conditions were imposed to produce a turbulent well-mixed layer. We used full surface initialization (sfc_full_init =.true.), with the lower boundary characteristics 75 set to USGS values for land use most closely matching the Anderson fuel category F (Anderson, 1982). The corresponding surface roughness lengths added various levels of wind shear to each domain to produce a more realistic non-uniform vertical wind profile during spinup of the environment before the fire was initialized in the LES.
Initial convection in the ambient environment was triggered using a perturbed surface temperature field. On average, a nearstationary turbulence spectrum was achieved within the first 30 min of run start. The "restart" file generated at the end of one 80 hour of spinup was used to initialize the main burn simulation, ensuring the fire was ignited in a well-mixed turbulent ABL.
The fire was initialized over a one-minute interval using a straight line of length L. The ignition line was placed one kilometer downwind of the western edge of the domain and centered in the north-south direction (sample illustration of this setup can be found in Appendix A). With a refinement ratio of 10 in each horizontal direction, the fire was simulated on a 4 m sub-grid mesh.

85
The "smoke plume" was modelled with a passive tracer emitted proportionally to the mass and type of fuel burned. The rate of release was controlled by an assigned emission factor representing PM 2.5 for each fuel category, based on values provided by Prichard et al. (2017) Figure 1. Pre-ignition potential temperature profiles (stability condition R). Colors correspond to initial soundings used for model spinup.
The range of ambient winds tested was bound largely by numerical constraints. Due to cyclic boundary conditions, wind 95 speeds higher than 12 ms −1 would require a much larger domain to prevent smoke recirculation. For the lower bound on our wind condition W, we needed to ensure that sufficient wind speed was maintained to propagate the fire. The spread algorithm used within the LES applies a correction factor under low wind speed conditions to prevent the fire from extinguishing itself.
While necessary for numerical reasons this effect is not physical, so winds below 3 ms −1 were excluded from our dataset.
We used 9 different atmospheric profiles (R condition) to initialize the model. We varied the following features for each 100 initialization: initial ABL height (500 m -1600 m) potential temperature lapse rate above inversion (0 K km −1 -20 K km −1 ) initial ABL temperature (290 K -300 K) We tested all fuel categories available within the model (F condition), and varied the length of the fireline (L condition) between 1 and 4 km. Weakly buoyant non-penetrative plumes whose smoke remained within the well-mixed ABL were excluded from the dataset, as their behavior is governed by different physics. A tabulated summary of all combinations included can be found in Appendix B.

110
Note, that varying a single condition while holding the rest constant does not result in a controlled experiment isolating its impact on plume rise. Because WRF-SFIRE incorporates fire-atmosphere coupling, the problem is not well-constrained. For example, by varying fuel type F alone, while holding the rest of test conditions constant, we obtain a set of fires with diverse shapes, sizes, intensities, fireline depths, rates of spread and heat release. This reflects the complexity of non-linear interactions that exist between the fire and the atmosphere. As a result, the parameter space captured within our LES dataset is much greater 115 then the four conditions described in Table 2.

Defining Smoke Injection Height
Given non-stationary fire and atmospheric conditions, determining a consistent definition of an equilibrium smoke injection height is not a trivial task. It requires separating buoyant rise from dispersion, while excluding the effects of initial momentum overshoot and accounting for the advection due to varying ambient and fire-generated winds.

120
A common way of examining vertical distributions of pollutants in the context of air quality is to consider CWI concentrations. This allows to reduce the problem to two dimensions, with plume centerline being defined simply as the CWI concentration maximum at each location downwind of the source. Theoretically, under stationary conditions there exists an equilibrium height, around which the centerline eventually oscillates. In reality, as well as in our LES experiments, neither the ambient nor the fire conditions are stationary. The changing location, shape and intensity of the fire, ABL warming and growth, as well as 125 the development of fire-coupled winds and vorticity continually modify the conditions. As a result, our approach is based on defining a region, where the concentration distribution is quasi-stationary. We consider the last frame of each simulation for this analysis. Using CWI integrated tracer values, we locate the plume centerline ( Fig.   2a). Due to random effects of ABL thermals as well as fluctuations in fire intensity and propagation speed, both centerline height and concentration can vary near the heat source. These oscillations are naturally suppressed in the stable layers above 130 the ABL, as the plume travels downwind and undergoes additional widening and mixing. To obtain the quasi-stationary region for each individual plume, we first calculate the change in tracer concentration along the centerline. We then use a smoothing function to reduce the effect of random turbulent oscillations in both the centerline height and the tracer concentration gradient along the centerline. The downwind region where both of these parameters are not changing rapidly are then then considered quasi-stationary. Additional details of this filtering method are provided in Appendix C.

135
The vertical CWI distribution of tracers are then averaged in the downwind direction over the identified quasi-stationary regions (shaded in grey on the Fig. 2c) to produce a representative downwind distribution for each plume (Fig. 2d). We define the "true" injection height z CL as the mean height of smoothed centerline over the averaging region. The resultant dataset of z CL values is used to constrain and evaluate the proposed smoke injection height parameterization introduced in the following sections.

Smoke Injection Height Model for Penetrative Wildfire Plumes
A common approach to predicting the final equilibrium centerline height of wildfire smoke is to first estimate the initial buoyant energy of the hot rising smoke (Briggs, 1975;Sofiev et al., 2012;Anderson et al., 2011). After the smoke plume entrains surrounding ABL environmental air and cools, the remaining energy is spent doing work to push the cooled smoke plume up into the statically stable capping inversion. The relationship between final and initial energies is often rewritten to show that the potential energy per unit mass (PE) of smoke penetration equals some fraction c 1 of initial heat released from the fire. In kinematic units, the initial heat input has units similar to a kinetic energy per unit mass (KE). The empirical parameter c 1 is usually estimated based on concepts of entrainment into the rising smoke plume (Cushman-Roisin, 2014).
The PE of smoke-plume penetration into the capping inversion can be written as where the penetration distance z of the final equilibrium smoke centerline z CL above reference height z s (near the top of the well-mixed portion of ABL) is The static-stability variable g for the plume-penetration region is where θ CL and θ s are the potential temperatures of the ambient environment at z CL and z s , respectively, and θ CL − θ s = θ .
The KE can be estimated using a velocity scale w f as Traditionally, the bulk potential-temperature difference across the smoke-plume penetration region θ is expected to be relevant for only the PE portion of Eq. (1). However, we found from the LES runs for a wide range of fire and environment conditions that the KE also depends on the same potential temperature difference. This dependence can be expressed in the velocity scale: This velocity scale is related to the fireline intensity parameter I, which is the kinematic heat flux into the atmosphere integrated across the fireline depth (in units of Km 2 s −1 ), and to the mixed-layer depth z i . Note, that I effectively corresponds to the kinematic form of Byram's Fireline Intensity (in units of kWm −1 ).
We speculate that this interesting result is because smoke from a fire does not rise through a passive environment, as is often assumed for Briggs types of plume entrainment models. Instead, the fire and the environment interact in many complex ways.

170
Some of these include: vertical-to-bent-over vortices on the ends of the fire line that rapidly mix environmental air into the buoyant smoke plume; modulation of fire intensity and fire updrafts by translation of ambient thermals across the fire line; plumes of enhanced convergence and updraft along the fire line; mass conservation as descending air beneath the extended smoke plume lowers the local mixed-layer depth; and other factors.
The above can be rearranged into the following form where the dimensionless empirical parameter is C ≈ 1. The factors in square and curly brackets with their corresponding 180 powers have units of time and velocity, respectively. This relationship is plotted in Fig. 3. It provides quite an acceptable fit to the data over a wide range of 140 combinations of fire and atmospheric conditions simulated.
Equation (8) suggests that the relevant length and temperature scales (z , θ ) depend not on the capping inversion strength alone, or on the tropospheric lapse rate above the capping inversion alone, but on the bulk potential-temperature differences across the smoke-plume penetration region, z CL −z s . Eq. (8) is implicit, in that the desired plume centerline equilibrium height 185 z CL appears in both the left and right sides of the equation. The plume centerline height also defines where θ CL is retrieved from the atmospheric sounding; namely, z CL is implicit in both Eq. (7) and (8). However, for any specific fire and environment conditions, values of z CL are easily found by iteration (see Appendix E). Steps to estimating input parameters required for the proposed injection model from the LES data are summarized in Appendix D.
Alternatively, for a small sacrifice in accuracy, we can obtain an explicit solution by considering an idealized version of the 190 atmospheric profile, consisting of an adiabatic mixed layer, entrainment zone and a stable uniformly stratified free atmosphere above (Fig. 4). In such case γ is defined as the overall potential temperature gradient of the free atmosphere and z s as the height corresponding to the intercept of γ and the well mixed portion of the ABL profile. Then, using Eq. (8), z CL can be found explicitly as:

Results
To assess the accuracy of the proposed smoke injection height parameterization (Eq. (8)), we performed two sets of verification studies. The first approach is based on using the synthetic plume dataset to perform model evaluation, bias correction and sensitivity analysis with idealized data. The second portion of this section applies our approach to a case study of a real prescribed burn (RxCADRE 2012).  Fig. 3 are "true" and parameterized smoke injection heights. The former is obtained directly from the LES, as per Sect. 2.3. The latter is determined iteratively using the proposed smoke injection height parameterization (see Appendix E for implementation details).

Shown in
Individual prediction errors do not appear to be a function fireline intensity, as indicated by scatter point color in Fig. 3, or 205 ambient winds (not shown). While overall the model performance is encouraging, the small discrepancy between the unity and regression lines suggests a linear bias. This can be remedied by applying bias correction using regression parameters from the fit shown in Fig. 3. This optimized model produces errors on the order of 20 -30 m, as suggested by the interquartile range shown in Fig. 5d. Model bias will be addressed in further detail in Sect. 5.
Given smooth averaged profiles from the synthetic dataset and excluding condition R8 (adiabatic free atmosphere), the 210 explicit solution using Eq. (9) offers comparable accuracy to the iterative version for both raw and bias corrected datasets (Fig.   6) . We address the limitations of using the explicit approach in Sect. 5.5.

Model Sensitivity
To asses how sensitive the smoke injection model performance is to the particular choice of bias correction parameters, we partition our original plume dataset into training and testing groups through random sampling. We obtain the linear bias 215 correction parameters using training data only (80% of runs). We then apply our bias-corrected iterative solution to the test group (remaining 20% of the runs) and assess model accuracy. Figure 7 summarizes model performance and sensitivity, based

Evaluation with Observations
Next, we apply the proposed model to a real-life case-study. We use observational data from the RxCADRE L2G prescribed burn (Ottmar et al., 2016) and it's numerical simulation (Moisseeva and Stull, 2019). This case study was selected based on (i) comprehensive nature of the observational dataset (ii) penetration of the plume above ABL top and (iii) availability of a completed model validation study, confirming that WRF-SFIRE reasonably captures the smoke plume produced during the 225 burn. Fig. 8 is the strip headfire pattern used to ignite the grass plot. We estimate the burn's input fireline intensity parameter I in two different ways: from raw data collected during the burn as well as from the numerical simulation.

Shown in
The observations-based value I obs is derived from the integral heat flux data obtained from the Highly Instrumented Plots (HIPs) fire behavior package (FBP) sensors (Jimenez and . We use the provided time-integrated values, averaging 230 between all sensors with confirmed fire at the sensor location (as indicated by video footage ). We then We apply our iterative solution (Eq. (8)) to find two z CL estimates based on I obs and I LES , and compare them to the 245 CWI smoke injection height obtained from the LES. The results are shown in Fig. 9. The parameterized injection heights are under-predicted by 20 m and 70 m for LES-and observations-derived I values, respectively.

Plume Classification
In previous sections we apply an energy balance parameterization to predict the mean smoke injection height z CL of a given 250 penetrative plume. For this purpose, only plumes rising above ABL top z i were included in the synthetic plume dataset used to constrain and evaluate the approach (see Table B1). In this section, we step back and consider all performed simulations, to determine whether the same equations can also be used to classify penetrative vs. non-penetrative plumes.
The synthetic dataset described in 2.2 consisted of 140 runs and excluded 7 simulations, where the plume remained trapped in the ABL (see Table B1 and Table 3). We determined this by visual analysis of CWI centerline and smoke fields. The excluded 255 plumes typically exhibited oscillatory or irregular centerline behavior (within the ABL, such as shown in the example in Fig.   B1) with little or no smoke injected above z i . For several combinations of fire and atmospheric conditions, however, making the distinction was challenging. For this reason, we included these "marginally-penetrative" plumes in the dataset.
In real-world applications, classification is a fundamental first step in plume rise parameterization process Sofiev et al. (2012). A viable automated method for categorizing penetrative vs. non-penetrative plumes requires that the distinction be 260 made based on available input parameters, rather then smoke observations (as such are typically not available at the time of making a forecast).
Conveniently, we can use Eq. (8) to obtain a z CL estimate for any combination of input parameters without prior knowledge of plume type. It can, hence, be applied as a classifier by requiring that for a penetrative plume where z i+ denotes the height of the upper edge of the numerical grid box (or ambient atmospheric sounding) containing z i . In other words, this definition ensures that z i and z CL are not in the same vertical model level. If this condition is not satisfied, the plume is assumed to be non-penetrative.
This approach correctly classifies all non-penetrative plumes that had been identified by visual analysis (Table 3). In addition, several plumes exhibiting marginal behavior are also classified by Eq. (11) to be non-penetrative.

270
For the purpose of subsequent dispersion modelling within real-world applications, non-penetrative plumes (i.e. all plumes listed in Table 3) would be assumed to become uniformly mixed in the vertical within a few convective turnover distanced downwind of the fire. Turbulent eddies within the ABL produce a well-mixed layer, resulting in relatively homogeneous vertical distribution of pollutants between the surface and z i . In contrast, for plumes that extend above z i , spanning the ABL, the entrainment layer, and/or the free troposphere, subsequent dispersion is typically handled by trajectory models.

Comparison with Existing Models
The above model evaluation indicates encouraging performance for the proposed smoke injection parameterization (Eq. (8)) at little computational cost. An additional advantage of our method is that it does not require making simplifying assumptions regarding the shape and heat flux distribution of the fire. This allows us to easily apply our model to complex heat sources, such as one produced with the strip head fire ignition pattern during the RxCADRE L2G prescribed burn (Fig. 8).

280
Unlike most existing plume rise parameterizations (Briggs, 1975;Rio et al., 2010;Freitas et al., 2007) we focus on a CWI centerline. Our model can be viewed as a "bulk method", having some common ground with the thermodynamic approach used in the FireWork modelling framework (Anderson et al., 2011;Chen et al., 2019) and the energy balance approach proposed by Sofiev et al. (2012). More specifically, we make no attempt to predict the full evolution of the rising plume centerline velocity or temperature before it reaches its equilibrium height. Rather, we focus on the energy balance of the plume over a "penetration we apply the energy balance approach to a layer well above the surface, starting from a reference height z s close to the top of the ABL.
As noted in Sect. 3, the implicit functional form of our solution (Eq. (8)) can be interpreted as a characteristic timescale multiplied by the characteristic velocity scale w f . By rearranging Eq. (7) and substituting Eq. (8) for z it can be shown that the two expressions for w f are equivalent, namely: The scaling relationship between vertical plume velocity and cubic root of fire heat has been previously established with both Rio's and Freita's models (Rio et al., 2010;Freitas et al., 2007), although our formulation includes different variables inside the radical. While both of our forms for w f and both model formulations (the simplified Eq. (7) and the expanded Eq. (8)) are mathematically equivalent, conversion from one form to another requires raising terms to 6 th power. This results in large 300 prediction errors; hence, for practical applications, the full Eq. (8) should be used.

Dimensionless Relationship
As discussed in Sect. 3, we can obtain an explicit solution for z CL by making additional assumptions about the vertical profile of potential temperature above the ABL. This allows us to reduce our Eq. (9) to a similarity relationship with two dimensionless groups z and H, denoting the left hand side (LHS) and right hand side (RHS) of Eq. (13), respectively. Nondimensional z and

305
H are linearly related, as shown in Fig. 10. The simple relationship suggests that our modelling results could fairly easily be scaled to a wider range of fire and atmospheric conditions, beyond those captured by the synthetic dataset presented in the paper.

310
The raw, non-bias-corrected form of the model suffers from a positive bias for tall plumes, as suggested by Fig. 5c and 6c.
In other words, z CL is overpredicted for plumes injected high above the ABL. We speculate that this is due to the simplifying assumption that most of the cooling, mixing, and dilution occurs below the reference level z s in upper portion of the ABL.
As the distance between z s and z CL increases for tall plumes and as the smoke travels further into the free atmosphere, this assumption becomes increasingly less accurate. Additional radiative cooling and entrainment of ambient air is, therefore, unaccounted for, resulting in over-prediction for z CL .
This issue is largely resolved for our dataset with the applied bias-correction. However, cases with strong shear turbulence and active smoke mixing above the ABL are still likely to be overestimated.

Limitations
The most significant limitation of the proposed smoke injection height parameterization is that it applies only to smoke plumes 320 with no water vapor condensation. Latent heat effects are not considered. Hence, smoke injection level for extreme pyroconvective events (e.g. flammagenitus clouds (WMO, 2017)) will likely be grossly under-predicted with the given formulation.
Therefore, in its current form, our parameterization is unlikely to be suitable for large-scale applications (e.g. global chemical transport models). However, it has the potential to improve regional air quality tools (e.g. BlueSky), since wildfire emissions sources are largely dominated by in-or near-ABL non-condensing smoke plumes (Val Martin et al., 2010.

325
Given the energy-balance formulation of our plume rise parameterization, it may be possible to incorporate latent heat effects by including an extra PE term in Eg.
(1). Similarly to the iterative process for finding a level of neutral buoyancy with Eq. (8) using potential temperature, it may be possible to predict plume condensation level using ambient humidity profile. However, a big obstacle to this development is that, to our knowledge, WRF-SFIRE has not been validated for such conditions. Unlike many existing methods, our parameterization relies on fireline intensity parameter I, rather than average fire heat flux 330 values, as input. While this approach offers an advantage for modelling plumes from complex ignition sources (such as shown in Fig. 8), fireline intensity is difficult to observe in the field.
Another limitation is the inherently implicit form of the full model Eq. (8). While we have not encountered any issues using an iterative solver to find z CL , atypical (or extremely noisy) ambient atmospheric soundings could potentially affect convergence. The explicit form (Eq. (9)) derived using the idealizing ambient sounding (Fig. 4) offers a possible solution for 335 such cases. However, it fails for weakly stable and adiabatic free atmosphere (eg. condition R8 in Fig. 1), as θ s is extrapolated into lower levels of ABL.
Lastly, the model has been developed and tested only for typical daytime atmospheric conditions. We have not assessed model performance for stable night-time atmospheric profiles or in the presence of strong vertical windshear.

340
Plume rise estimation remains one a weak link in our ability to forecast where and how smoke from wildfires travels in the atmosphere. In this study we present a simple parameterization (Eq. (8)) for predicting CWI smoke-plume centerline height from a wildfire of an arbitrary shape and intensity. Our approach is based on energy-balance of the plume over a penetration region. We constrain and evaluate the proposed method using a synthetic LES-derived plume dataset developed for a wide range of fire and atmospheric conditions.

345
Based on the results of cross-evaluation with LES data as well as a real prescribed burn case study, the parameterization offers reasonable accuracy at little computational cost. We demonstrate that the approach can also be applied as a classifier to distinguish penetrative and non-penetrative plumes. This information is key for subsequent dispersion modelling, as plume AMBIENT WIND → Figure A1. Numerical domain setup.
behavior is governed by different physics above and below the ABL. The proposed method can be used as a sand-alone deterministic model or embedded in a host smoke modelling framework.

350
We hope that parameterization presented in this study will be of interest to air-quality researchers to provide a low-cost solution for regional wildfire emissions-modelling applications.
Code availability. S1: WRF-SFIRE sample initialization files (sample_simulation.zip) Red cells in Tab. B2 highlight simulations that were completed, but subsequently excluded from analysis presented in Sect.
3. This was done based on visual inspection of LES fields. There were two possible reasons for exclusion: (i) the plume reached the top of the domain or (ii) the plume appeared to be non-penetrative. In the former case, it's questionable whether the fields are physical, as the plume could potentially be affected by the absorbing layer near domain top, designed to prevent numerical 365 instability. The latter rendered the plume irrelevant for the purpose of analysis presented in Sect. 3. These non-penetrative runs, however, were included for testing the plume classification method presented in Sect. 5.1.

Appendix C: Identifying Quasi-Stationarity
We define the quasi-stationary downwind region for each plume based on two factors: the height of the centerline and tracer concentration gradient along the centerline. Our filter attempts to extract only those portions of the downwind CWI smoke 370 distribution, where both of these factors are changing slowly.
First, we remove the effect of random turbulent oscillations by applying a smoothing function (Savitzky-Golay filter provided by SciPy library with polynomial order set to 3) to both the concentration gradient along the centerline and the centerline height.
We vary the size of the smoothing window as a function of mean ambient wind condition W, such that window_length = max(W · 10 + 1, 51) grid points.

375
The filter then applies the following criteria to extract quasi-stationary regions: smoothed tracer concentration along the plume centerline varies by less then 10% of the maximum concentration gradient smoothed centerline height varies by less then a 100 m the location is downwind of the maximum tracer concentration gradient the location is at least 10 grid points away from the maximum in smoothed and non-smoothed centerline height 380 the location is at least 50 grid points away from the downwind endpoint of the centerline The above thresholds were determined through an informal sensitivity analysis (not shown), based on the filter's ability to effectively identify regions of near-stationary plume centerline height for all simulations in our dataset.

Appendix D: Estimating Model Input Parameters
Summarized in Table D1 are parameters associated with an iterative solution for z CL using Eq. (8). Below is our approach to 385 estimating these parameters from LES data.
As noted above, we consider the problem in crosswind direction. Given a three-dimensional fire of an arbitrary shape (eg. The values of θ s and θ CL are then determined from the pre-ignition sounding for each simulation using the definitions of z s and z CL (as described in Sect. 2.3).
The numerical implementation of our iterative solution using SciPy's fsolve function (scipy.optimize.fsolve) is as follows. We rewrite bias corrected Eq. (8) into an input function toSolve as: where C = 1.005, B 1 = 0.924 and B 2 = 116.417 are bias correction parameters, T 0 is the potential temperature sounding vector, dz is the vertical step and int() is a standard Python function converting the bracketed value into an integer.

405
A possible issue for some solvers is that we are, effectively, iterating over the vertical index of the column vector T 0 corresponding to z CL . As the numerical solver attempts to converge on a solution it may query a non-existent index and fail.
We are able to obtain a fast and consistent performance by ensuring we set z i as the initial guess for z CL and by minimizing the initial step bound option of the solver z CL = f solve(toSolve, z i , f actor = 0.1) (E2)