An evaluation of empirical and statistically based smoke plume injection height parametrisations used within air quality models

Keywords: air quality model, CALIOP, ceilometer, MicroPulse scanning lidar, plume rise, prescribed burns, remote sensing, wildfire.

Worldwide, wildland fires (prescribed and natural) emit large quantities of harmful gas- and particulate-phase pollutants into the atmosphere (Andreae and Merlet 2001; Van der Werf et al. 2010). Globally, fire-related air pollutants are estimated to cause up to 330 000 fire-attributable deaths year⁻¹ (Johnston et al. 2012; Lelieveld et al. 2015). Across the contiguous United States (USA), the economic cost of wildland fire on human health has been estimated to be ~US$100 billion and cause ~8500 premature deaths year⁻¹ (e.g. Rappold et al. 2017; Fann et al. 2018). To minimise economic losses and impacts on human health related to smoke inhalation and smoke forecast, additional assessments of wildland fire emissions are necessary, because modelled emission estimates and meteorological dispersion contain errors and bias due to data limitations (Goodrick et al. 2013; Wilkins et al. 2018; Liu et al. 2019).

Liu et al. (2019) listed knowledge gaps and suggested areas needing wildland fire modelling research. An overarching theme was the vertical distribution of a smoke plume, known as the plume rise parametrisation. Plume rise is controlled by fire and meteorology, which include the energy released by the fire, the size of the combustion zone, fuel composition and ambient atmospheric conditions (Labonne et al. 2007; Achtemeier et al. 2011; Paugam et al. 2016; Walter et al. 2016). The height of the plume is commonly referred to as smoke plume injection height (PIH) (e.g. Paugam et al. 2016). Inaccurate representation of PIH can seriously compromise air quality model performance (e.g. Achtemeier et al. 2011; Mallia et al. 2018; Wilkins et al. 2018; Liu et al. 2019).

Our understanding of wildland smoke behaviour and plume evolution has vastly improved with the use of satellite data (Colarco et al. 2004; Al-Saadi et al. 2008; Soja et al. 2009, 2012; Ichoku et al. 2012; Val Martin et al. 2018; Sokolik et al. 2019), ground-based lidar and Doppler radar (Charland and Clements 2013; Clements and Oliphant 2014; Clements et al. 2006, 2007, 2016, 2018; McCarthy et al. 2018). Remote sensing data have improved the ability to evaluate plume rise models but there are still inherent limitations, such as with satellite overpass times – some occur late morning before the plume has had a chance to fully develop (Giglio et al. 2010). These remote sensing analyses suggest that the plume rise parametrisation problem in models is often linked to the interaction between the smoke plume and the capping inversion at the top of the planetary boundary layer (PBL) (Val Martin et al. 2018; Wilkins et al. 2020). The likelihood for a smoke plume to penetrate the PBL top varies by fire size, fuel loading and atmospheric conditions and can range from a negligible fraction to full penetration (Kahn et al. 2008; Val Martin et al. 2010; Tosca et al. 2011). It is imperative to model the PBL and stable layer (STL) height correctly, because both the PBL and STL can often act as capping atmospheric layers constraining PIH. For smoke to escape the PBL or STL, the heat from a fire must generate sufficient buoyancy to penetrate these layers. Furthermore, the ability of a smoke plume to penetrate the PBL or STL can be highly dependent on the current atmospheric conditions under which a plume is generated, because they can strongly cap or assist in upward motion of air. Whether a plume remains within or penetrates the PBL or STL strongly influences pollution exposure, duration, transport and chemical residence times (e.g. Westphal and Toon 1991; Wotawa and Trainer 2000; Hyer et al. 2007; Wilkins et al. 2020).

A widely used plume rise approach in deterministic air quality models (AQMs) is Briggs (1975), but its accuracy and performance have been questioned. It has been reported that the Briggs approach both underestimates plume rise for small fires and overestimates plume rise for large fires (Wilkins et al. 2018). Gordon et al. (2018) reported a 50% underestimate in stack heights. Other plume height parameterisations exist (e.g. Freitas et al. 2007; Rio et al. 2010; Sofiev et al. 2012); however, implementing these more complex schemes has not always improved model performance (e.g. Freitas et al. 2006; Kahn et al. 2007; Leung et al. 2007; Val Martin et al. 2012; Paugam et al. 2016). It can be difficult to determine the exact reason for the model uncertainty as it could be related to the plume rise algorithm or the inputs driving the plume rise, i.e. heat flux. Some of these discrepancies are likely related to historic model design, previous state of knowledge and input limitations. Historically, regional and global AQMs did not account for smoke from small fires (<500 ha) that penetrated the PBL (Mims et al. 2010; Zhou et al. 2018; Baker et al. 2019), latent heat releases often seen in fire-generated smoke-infused pyrocumulonimbus thunderstorms (Peterson et al. 2015, 2018), or fires that generated multiple cores and were detrained in vertically distinct layers of the atmosphere (Kahn et al. 2007; Achtemeier et al. 2011; Liu et al. 2013; Val Martin et al. 2010).

The present study implemented empirical and statistically based plume rise algorithms in an AQM for two conditions: a large wildfire (>6000 ha day⁻¹) and a series of small prescribed fires (<500 ha day⁻¹). Specific plume rise parameterisations were chosen by their computational cost to implement in an AQM, ease of implementation and availability of variables. A quantitative assessment of plume rise was then conducted and compared with satellite- and ground-based observation data. The goal of this study was to seek improvement of wildland smoke plume injection heights in the US Environmental Protection Agency (EPA)’s Community Multiscale Air Quality (CMAQ) modelling system.

Observational data were obtained from the California Rim Fire in August 2013, representative of a large wildfire (Fig. 1) and the Konza Prairie Biological Station prescribed burn experiment in the Flint Hills, Kansas, in March 2017, representative of a series of small grassland burns (Fig. 2). Table 1 summarises the field-specific information for each of these burns.

Fig. 1.  Burned areas from the 2013 California Rim Fire: Burn 1 (21 August, yellow shading), Burn 2 (24 August, pink shading), and the full burned area perimeter (red line).

Fig. 2.  Burned areas from the March 2017 prescribed fire experiment in Konza Prairie, Kansas. Burn units 3, 4 and 5 are shaded yellow, orange and pink, respectively. Also shown are the locations of ground-based Hexagon Miniaturised Micro Pulse Lidar (MiniMPL) (open white squares) and two Vaisala Model CL-51 ceilometers at the Konza Prairie Biological Station (open black squares). For Burns 3 and 4 (16 Mar 17), the MiniMPL and an additional CL-51 to the one placed at Headquarters was deployed. For all other days, only the CL-51 at Headquarters was used.

Table 1.  Summary of burn information for the California Rim Fire (Burns 1, 2), which burned in August 2013 (from Peterson et al. 2015), and the Konza Prairie Biological Station prescribed burn experiment (Burns 3, 4, 5), which burned in the Flint Hills, Kansas, in March 2017 (from Whitehill et al. 2019)

The Rim Fire was an intense wildfire, burning 104 131 ha, that had an active flaming stage from late August to early September 2013 in the central Sierra Nevada Mountains, California. The Rim Fire contained a few major spread events (periods of rapid growth in burned area) generating plumes that reached well above the PBL as determined by satellite retrievals and aircraft imagery (Peterson et al. 2015). The first spread event (Burn 1), between 21 and 23 August 2013, burned 36 206 ha (~35% of the total) and the second spread event (Burn 2), between 25 and 26 August 2013, burned 12 067 ha (~12% of the total). We analysed the first spread event (12 068 ha day⁻¹) and the day before the second spread event (6033 ha day⁻¹). The ecosystem fuels were mixed forest dominated by a coniferous overstorey and shrubs. Fuel loading was estimated at 1-km resolution using the United States Forest Service Fuel Characteristic Classification System database (McKenzie et al. 2007; Ottmar et al. 2007; Larkin et al. 2009).

During the Konza Prairie prescribed burn, 13 total field units were burned (<4 to 205 ha); however, owing to the size of those individual units, only six fields were chosen for this analysis (threshold >30 ha day⁻¹). Four of those six fields were combined into one fire for modelling purposes because they were close in proximity, duration and burn start times. Each of the resulting three Konza Prairie burns (Burns 3, 4 and 5) occurred from late morning to early afternoon on 16 (Burns 3, 4) or 20 March 2017 (Burn 5). The fuel loading was composed of 95% big bluestem grass, switchgrass and Indian grass. The fields were irregularly shaped and followed natural terrain features with roads used as fire breaks (Whitehill et al. 2019).

Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP) is an instrument on the Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observations (CALIPSO) satellite (Omar et al. 2009; Winker et al. 2009; http://www-calipso.larc.nasa.gov/). CALIOP provides high vertical resolution (30–60 m) 532 nm and 1064 nm (km⁻¹ sr⁻¹) attenuated backscatter data, which can detect both thick and optically thin smoke layers in the atmosphere. Standard products include a Vertical Feature Mask, which distinguishes clean air, clouds, stratosphere, surface, subsurface and aerosols, and an Aerosol Subtype product (Fig. 3; additional tracks are shown in Fig. S1, Supplementary material).

Fig. 3.  CALIPSO track overlaid on a MODIS Visible image, showing both the vertical and horizontal extent of smoke and clouds. MODIS Terra and Aqua fire detection data (red dots) highlight the Rim Fire burning in California on 25 August 2013. The visible smoke was transported across California, Nevada, Oregon and Idaho before intersecting with the CALIPSO track over Montana and Canada. The CALIOP smoke-aerosol vertical-profile data (black and brown) extend from the surface to ~5 km. The vertical extent and height of the clouds are evident to the south and north of the CALIOP swath, and the horizontal extent of the clouds is evident in the visible MODIS image. CALIOP data provide the vertical properties that inform the horizontal view to provide a complete representation of aerosol transport and the atmosphere.

For the Rim Fire, a satellite-based PIH product was produced by extracting the smoke-filled aerosol parcels from the CALIOP data. Figs 3 and 4 highlight the vertical and the horizontal transport of Rim Fire smoke across several states, where the smoke transected the CALIPSO track. In this case, the smoke extended from ~5 km to the surface (Fig. 3). The CALIOP-based plume-detrainment height data used in this work were produced by:

1. overlaying CALIPSO tracks on the NOAA Hazard Mapping System (HMS) smoke product to ensure that the aerosols were smoke (Fig. 4; see Figs S1 and S2 in the Supplementary material for additional CALIPSO tracks);

2. extracting CALIOP aerosol data, which were used to initialise the Langley Trajectory Model (LaTM) (Pierce et al. 2003, 2009);

3. the smoke-laden aerosols were then transported backwards in three-dimensional space and time (15-min time steps) until they horizontally coincided with a fire within 20 km (Fig. 5). Fires were determined using daily Moderate Resolution Imaging Spectroradiometer (MODIS) active-fire detection data (Figs 3 and 4; Giglio et al. 2003; Soja et al. 2012; Thomas et al. 2017; Baker et al. 2018).

Fig. 4.  CALIPSO track overpass path on 25 August 2013 (green dots: see Fig. 3) overlaid on a Hazard Mapping System Fire and Smoke (HMS) product with MODIS Terra and Aqua fire detection data (red dots). Map shows the site of the 2013 California Rim Fire. The aerosols backtracked to the Rim Fire were transported across California, Nevada, Oregon and Idaho before intersecting with the CALIPSO track over Montana and Canada (yellow dots).

Fig. 5.  2013 California Rim Fire. Maximum vertical allocation of smoke-plume tops of Burn 1 (21 August) and Burn 2 (24 August), plotted with the planetary boundary layer (PBL) and fire radiative power (FRP). Smoke plume detrainment was derived using Cloud-Aerosol Lidar with Orthogonal Polarisation (CALIOP) data. Each colour represents a distinct CALIPSO overpass. For example, the 21 August smoke plume was captured by nine distinct CALIOP overpasses, and the 24 August smoke plume was captured by seven CALIOP overpasses. Each CALIOP overpass is named by the date of smoke/CALIOP data coincident, model initialised GMT (Greenwich mean time) time, and the location of the smoke plume in the CALIOP overpass. For example, the smoke represented in purple-blue from 21 August coincided with a CALIOP overpass on 22 August at 41–44°N latitude; back trajectories were initialised to connect the smoke to the fire at 1200 UTC.

NASA’s Modern-Era Retrospective Analysis for Research and Applications, Version 2 (MERRA-2) meteorological data were used to drive LaTM simulations (Gelaro et al. 2017). For details, see Supplementary material S1, CALIOP-derived retrievals.

For the prescribed burns at Konza Prairie, we used a ground-based scanning customised miniaturised version of the Hexagon Micro Pulse Lidar (Spinhirne 1993; Spinhirne et al. 1995), hereafter referred to as Mini Micro Pulse Lidar (MiniMPL), following the procedures described by Kovalev et al. (2005) and Charland and Clements (2013). The MiniMPL used a combination of plan position indicator (PPI) and range height indicator (RHI) scans to provide near-range atmospheric lidar retrievals of plume tops (Welton and Campbell 2002). For the three Konza Prairie burns (Burns 3, 4, 5), Fig. 6 shows the raw MiniMPL retrievals for Burns 3 and 4 that were converted into PIH (Fig. 7; see Fig. S5 in the Supplementary material for Burn 5 results).

Fig. 6.  2017 Konza Prairie Biological Station prescribed burn experiment showing, for Burns 3 and 4 (16 March), smoke-plume range height (m above sea level) and backscatter (km² µs⁻¹ µJ⁻¹), derived using a ground-based Hexagon Miniaturised Micro Pulse Lidar (MiniMPL). The lidar starting height was 440 m above sea level for the plume height measurements. For Burn 5 (20 March), see Fig. S5.

Fig. 7.  2017 Konza Prairie Biological Station prescribed burn experiment showing, for Burns 3 and 4 (16 March) and Burn 5 (20 March), maximum observed smoke plume tops (m above sea level), derived using a ground-based Hexagon Miniaturised Micro Pulse Lidar (MiniMPL). The lidar starting height was 440 m above sea level for the plume height measurements.

To estimate the PBL height during each of the three Konza Prairie burns, two Vaisala Model CL-51 ceilometers (Münkel et al. 2007; McKendry et al. 2009, 2010) were also deployed (e.g. Tsaknakis et al. 2011; Liu et al. 2013; Clements and Oliphant 2014). Operational details for the CL-51 and MiniMPL are summarised in Supplementary material S2 (Operational details for the CL-51 and MiniMPL) and Table S1.

Prior to evaluating alternative plume rise algorithms against the Briggs method, each model formulation was compared with independently measured or derived PBL heights. For the Rim Fire, we used the MERRA-2 PBL, and for the Konza Prairie prescribed fire, we took the PBL height measurements from the ceilometer and lidar data that were on site (Table 1). We then compared independent and modelled PBLs to determine model error and bias (Section S3, Planetary boundary layer (PBL) analysis).

To better understand model estimated predictions of PIH, we employed the US EPA’s CMAQ system (www.epa.gov/cmaq) with the standard option for plume rise, which uses the Sparse Matrix Operational Kernel for Emissions (SMOKE) as the emissions processing system (Pouliot et al. 2005). The CMAQ standard plume rise algorithm is detailed in sections S4 (Smoke model details) and S5 (Weather and Research Forecasting (WRF) model (WRF 12 km and 4 km)).

For the Rim Fire, the model was initialised using fire activity data (Sullivan et al. 2008) extracted from the BlueSky Framework using SMARTFIRE2 (Raffuse et al. 2009, 2012). For the Konza Prairie burns, we used field-specific information provided by Konza Prairie’s Biological Station research staff (Table 1). Parametrisation details for each plume rise method are listed in Table 2 (plume rise model parameterisation), Table 3 (algorithm configurations) and Table S3 (model equations). Specific plume rise parameterisations were chosen by their computational cost to implement in an AQM, ease of implementation and availability of parametrisation variables calculated by CMAQ. A major criticism of the standard CMAQ-Briggs PIH calculated on a 12-km grid (hereafter BASE12) is that it is based on experimental data from non-fire plumes (i.e. stack point sources). Another critique is that the buoyancy flux calculation may be inappropriate for large fires (Sofiev et al. 2012; Zhou et al. 2018). To assess these criticisms, we evaluated the CMAQ BASE12 model sensitivity to alternative PIH algorithms, model grid spacing and temporal profiles.

Table 2.  Plume rise model parameterisation

Table 3.  Plume rise model differences

To evaluate the BASE12 PIH algorithm sensitivity, we tested two models (Table 2, Table S3). The first model was an empirical energy balance parameterisation (like convective cloud formulations) designed for fires (Sofiev et al. 2012), which is implemented at 12- and 4-km grids (hereafter SOFIEV12KMHR and SOFIEV4KMHR). Sofiev is designed to use Fire Radiative Power (FRP); when available, we used the satellite-based measurements, and when not available, it was derived (FRP_calc = heat flux × 0.1 × area burned) using the model heat flux and area burned, assuming that radiative energy was 10% of the total fire heat energy (Wooster et al. 2005; Val Martin et al. 2012). The alternative model was underpinned by a physical rationale based on fire smoke’s general tendency to pool under stable layers (Kahn et al. 2007; Val Martin et al. 2010). This formulation assumes that a smoke plume will reach the PBL and not rise above that layer. The PBL is not a static parameter; therefore, we added 500 m to the PBL maximum value in the model for each hour simulated. This better represented the maximum PIH potential just reaching above the PBL, based on the hourly time-averaged thickness of the model layer used (hereafter PBL50012KM). See section S6 (Sofiev algorithm explanation) and Table S3 (energy balance equations).

To evaluate the BASE12 spatial profile sensitivity, a grid spacing of 12 km was compared with the Briggs algorithm at 4 km (hereafter BRIGGS4KMHR); the only difference was the meteorological inputs used at those respective grid resolutions. To evaluate the BASE12 temporal profile sensitivity, as suggested by Zhou et al. (2018), we adjusted the temporal profiles from the national-scale modelling standard 12-h profile (0600–1800 Local Standard Time) to a 3- or 4-h profile (from fire start time) on the 12-km grid (hereafter BRIGGS12KMHR).

Table 1 summarises plume rise information for both the 2013 Rim wildfire and 2017 Konza Prairie prescribed fire. Burns 1 and 2 of the Rim Fire (Fig. 1) occurred on 21 and 24 August 2013, respectively, during drought conditions of warm temperatures with low relative humidity (RH < 10%). These conditions caused the fires to spread quickly in winds of 3–10 m s⁻¹ (Peterson et al. 2015). The burned area was greater for Burn 2, but the daily fire spread rate for Burn 1 (503 ha h⁻¹) was double that of Burn 2 (251 ha h⁻¹). At Konza Prairie (Fig. 2), Burns 3 and 4 were conducted on 16 March, and Burn 5 occurred in mid-morning on 20 March (1600 Coordinated Universal Time (UTC)) through early afternoon (2200 UTC). Burns were conducted under the influence of an upper-level trough in a post-frontal air mass and a boundary layer range of 1–2 km by mid-afternoon (2100 UTC). Surface conditions contained a weakening high-pressure system (15 March: 1036.0 hPa; 20 March: 1010.5 hPa), southerly winds (2–5 m s⁻¹) and moderate RH (>20%). The large-scale patterns of wind, temperature and pressure affecting eastern Kansas during each burn were consistent except for temperature and cloud cover. During the burns, the temperature ranged from –4° to 33°C, with low RH of <15%, and light wind speeds of <5.5 m s⁻¹ (Whitehill et al. 2019). The vegetation composition was consistent among burns, with a fuel loading range of 530–630 g m⁻² and an average burn rate of 32 ha h⁻¹. Details on the meteorological conditions during the Rim Fire and the Konza Prairie prescribed burns are in Peterson et al. (2015) and Whitehill et al. (2019), respectively. The PBL comparison of the observations with model STL showed that the model performed better for the Konza Prairie prescribed grassland burns (mean bias (MB) of ±400 m)) than for the Rim wildfire (MB of +1000 m; see S3, Planetary boundary layer (PBL) analysis).

The 2013 California Rim Fire has been described as a massive wildfire, as the data here support. As of 2019, the Rim Fire was California’s fifth largest wildfire. This was also a heavily actioned fire, which was fully contained 9 weeks after the fire ignition. CALIOP plume detrainment data are shown in Figs 3–5 and Table 4. Fig. 3 shows one of several CALIPSO curtains (see Fig. S1 for additional curtains and Fig. S2 for all tracks used), where smoke aerosol data were extracted to develop the daily smoke-plume detainment height products (Fig. 5). Data were extracted from nine CALIPSO curtains to develop the daily smoke evolution for 21 August 2013 (Fig. 5), and seven curtains were extracted for 24 August 2013. Smoke-filled air parcels were injected and detrained into and above the PBL throughout the day on both 21 and 24 August, with plume heights increasing as the afternoon RH decreased and temperatures increased (Fig. 5).

Table 4.  Comparison of the percentage of smoke below and above the Planetary Boundary Layer (PBL) for the 2013 California Rim Fire (CALIOP data) and for the larger MISR analysis (Val Martin et al. 2010)
For the Rim Fire, most of the smoke was injected and detrained above the PBL

The maximum FRP was much higher on 21 August (3025 MW; Fig. 5, Table 4), compared with that on 24 August (max. 1734 MW). The plume was injected over 1 km higher on 21 August from ~1500 local time throughout the day. On all days of the Rim Fire, 83% of the smoke was injected and detrained above the PBL (Table 4); 21 and 24 August were particularly active, with 92% of the smoke detrained above the PBL. The Rim Fire was intense, with much more smoke lofted higher in the atmosphere than is typical. This result contrasts with the larger-scale Multi-angle Imaging SpectroRadiometer (MISR) plume-height analysis (Val Martin et al. 2010), which showed that most smoke (88–96%) remained in the PBL, highlighting the differences in extreme fires compared with ‘normal’ fires. As a result, the standard 12-h assumption was retained.

The scanning MiniMPL measurements (Fig. 6, Fig. S5) provide evidence that the largest smoke concentrations during the Konza burns were in the advancing portions of the fire. The lidar retrievals showed that smoke from these prairie fires <500 ha tended to pool in the lower PBL (Fig. 7, Table 1), with 17% entering the free troposphere. A comparison of the maximum PIH with the associated maximum backscatter concentrations indicates a contrasting pattern to the standard vertical fire emissions profile, which suggests that some of the smoke reached a higher elevation, but most pooled lower in the PBL. For example, with nearly triple the burn rate, Burn 5 had plumes nearly two to four times higher than Burns 3 and 4 (Fig. 7). For the individual fields <200 ha, smoke generally reached the top of the highest observed layer within the PBL, with some individual plume cores penetrating the PBL. This finding verifies the capping potential of the PBL for prairie fires under these prevailing weather conditions.

The diurnal plume height and intensity from the Konza Prairie burns (Fig. 7) provides additional evidence that small fires (<500 ha) can penetrate the PBL, even if plumes are generated with short burn times (<1–4 h). Plume heights were measured at 600–4100 m, while the PBL range was 750–1270 m. These prescribed fires were consistent with other fires where the meteorological conditions constrained plume rise (e.g. capping inversion), but the burn rate played a lesser role (e.g. Burn 4 plume tops were, on average, 500 m higher than those of Burn 3, with a burn rate difference of only 32 ha h⁻¹). The fuel loading for these fires showed no clear connection to PIH. For example, Burn 4 contained a fuel loading of 942 kg ha⁻¹ (~17%) more than Burn 3, but the resulting plume tops remained relatively similar.

On 20 March 2017, several fields totalling 405 ha burned concurrently in a small area relative to the AQM grid spacing; for this reason, we treated these as one field, Burn 5 (Fig. S5). Burn 5 produced several plumes with distinct cores that were injected 2000–3000 m above the PBL. Plumes that penetrated the PBL continued to rise to higher stable layers (3190 m at 1500–1700 UTC, 3750 m at 1800–2000 UTC). Still, a significant amount of the smoke remained just above the PBL (1300 m), which is consistent with smoke detraining at multiple levels in the atmosphere. Compared with Burn 5, the backscatter intensity of Burns 3 and 4 demonstrated that, as expected, plumes with a lower vertical extent had higher pollution concentrations near the surface.

The plume rise evaluation was carried out on an hourly basis, for each model layer, and for each of the five burns we evaluated (Table 1). For the Rim Fire (Burns 1, 2), we used hourly averaged observational data from satellite-based retrievals from CALIOP (Fig. 8: maximum plume tops; Figs S3, S4: mean and minimum heights) to compare with each model formulation. For the Konza Prairie burns (Burns 3, 4, 5), we used hourly averaged ground-based retrievals from the MiniMPL (Figs 9–11).

Fig. 8.  2013 California Rim Fire plume tops, derived from the CALIOP data (grey), CMAQ base model (red) with the planetary boundary layer (PBL, black dotted line), Sofiev method (blue), and PBL + 500 m approach (green). See Tables 2 and 3 for details on each simulation.

Fig. 9.  Konza Burn 3 (16 March 2017): comparison of smoke plume tops for observed (grey), base model (black), and six model alternatives (red). See Tables 2 and 3 for model simulation configuration specifics.

Fig. 10.  Konza Burn 4 (16 March 2017): comparison of smoke plume tops for observed (grey), base model (black), and six model alternatives (red). See Tables 2 and 3 for model simulation configuration specifics.

Fig. 11.  Konza Burn 5 (20 March 2017); comparison of smoke plume tops for observed (grey), base model (black), and six model alternatives (red). See Tables 2 and 3 for model simulation configuration specifics.

For the Rim Fire, the CMAQ-Briggs (BASE12) simulation overall demonstrated the ability to capture the mean PIH range observed (Fig. 8, lower end of the shaded area). The BASE12 simulation successfully captured plume ability from the Rim Fire to penetrate the PBL. The model–observations analysis with BASE12 consisted of an MB of less than ±5–20% compared with the alternative methods. The base simulation compared well with the 3000–5000-m mean plume height value range reported in Peterson et al. (2015). Compared with the observed and other PIH algorithms, BASE12 contained an inherent low bias of 1000 m during the day and a corresponding inherent high bias at night. The night-time plume height appeared to be mostly sporadic and only sometimes within the range of measurements. Similarly, Sofiev et al. (2013) stated that night-time plume height could not be substantiated because of the limited number of observations. These findings suggest that a high bias in particulate matter with a diameter smaller than 2.5 μm (PM_2.5), as discussed in Wilkins et al. (2018), could have been skewed by a night-time high bias in PM_2.5 and not an overall daily model PM_2.5 bias for wildfires.

Comparing the BASE12 run with the SOFIEV12KM12HR method provides additional evidence that BASE12 may underpredict PIH during the daytime for big spread events (12 068 ha day⁻¹) and overpredict PIH for moderate-spread events (6033 ha day⁻¹). The physical rationale method (PBL50012KM) determined that, during the day, PIH showed similar heights to CMAQ-Briggs, when the plume was constrained within the PBL (24 August). For times that the plume was not constrained within the PBL, the accuracy was limited by the height of the PBL. Our simulations did not couple meteorology and fire feedbacks as a system; therefore, we did not expect to capture fire emissions impacts to the PBL locally. This could include, but not be limited to, the impacts from fire heat fluxes directly related to the fire and potential smoke radiative impacts on the PBL (e.g. Kochanski et al. 2019). If the system is not coupled, then this rationale may not be reasonable, because this method relies on plumes not escaping the PBL (e.g. Burn 2).

For the prescribed burns (Konza Prairie; Figs 9–11), the BASE12 simulation indicated an overall low PIH bias of 600–2000 m. According to the algorithm differences, the PBL50012KM performed better for the prescribed grassland burns than for the forested wildfire cases (Rim Fire). This difference suggests that the smaller grassland fires (<500 ha) had less impact on the boundary layer physics and contained fewer plumes that escaped the PBL (10–50%). Our results for PBL50012KM are consistent with conclusions from previous studies that when smoke is injected into the free troposphere, it tends to accumulate within atmospheric layers of relative stability aloft (Kahn et al. 2007; Val Martin et al. 2010).

The Sofiev formulation most closely matched these observations. Therefore, we compared the performance of the Sofiev method with the CMAQ-Briggs (BASE12) method. This comparison provided more evidence that the BASE12 model is biased low (Burn 4: BASE12 MB –215.5 m, Sofiev MB –449.2; Burn 5: BASE12 MB –2975.1 m, Sofiev MB –3586.2 m). Lowering the grid spacing from 12 to 4 km provided little improvement in PIH estimates, only slightly increasing the error and bias. This small increase was likely because the meteorological input did not have many differences between the simulations (see S3, Planetary boundary layer (PBL) analysis).

Lastly, using the updated temporal profile, 3–4 h compared with 12 h, the simulation switched the overall PIH estimation bias from negative to positive (e.g. for Burn 3, BASE12 MB –272.4 m; BRIGGS4KM MB 174.8 m). Furthermore, the model’s ability to simulate the overall maximum PIH improved by 200–1000 m. The difference in model performance for estimating PIH for Burns 3 and 4, compared with Burn 5, was because of model design and stability layers. However, CMAQ-Briggs is not designed to capture sharp rising plumes from small fires and rapid changes in the PBL, which typically cap vertical motion.

Each model formulation (Table 2) was compared individually with the BASE12 simulation for the maximum PIH, the change in plume top compared with BASE12, and the percentage of plume change above PBL (Fig. 12). All simulations of plume top heights remained below 5000 m for the Rim wildfires and below 2000 m for the prescribed grassland burns at Konza. For the Rim Fire plume heights, a model-to-model comparison showed a 1500–2500-m high bias for both moderate and high-intensity burn days. For the prescribed grassland burns, there was a 50–1000-m low bias across all cases. For the wildfire cases, the base case placed >60% of the plume above the boundary layer, while the other simulation placed 30–55% of the plume above the PBL. For grassland prescribed fires <500 ha, the base case placed <10% of the plume above the PBL, while the alternative model formulations placed 10–50% of the plume above the boundary layer, depending on area burned (e.g. Burn 5) and meteorological situation (e.g. Burn 4). The PBL50012KM simulation placed 50% of the plume above the PBL owing to the way the model was formulated. Therefore, after removing the PBL50012KM simulation, the actual average indicated that, for small fires, 10–35% of the smoke reached above the PBL.

Fig. 12.  Modelled plume information is shown for each of the five burns selected at the 2013 California Rim Fire (Burns 1, 2) and the 2017 Konza Prairie prescribed burn (Burns 3, 4, 5). Orange: BASE12, grey: BRIGGS12KMHR, yellow: SOFIEV12KMHR, light blue: PBL50012KM, green: BRIGGS4KMHR, dark blue: SOFIEV4KMHR, brown: BRIGGS4KM12HR. Top: Plume tops (above ground level) for each simulation. Middle: change in plume tops compared with BASE12. Bottom: percentage of the plume above the planetary boundary layer (PBL).

We evaluated the model bias against model formulation changes in plume rise algorithm, grid spacing and temporal allocation (Fig. 13). Overall, the BASE12 model had a consistent low bias of 10–3500 m for the Rim Fire burns. For estimating PIH, the standard base case model performed better (MB –2200 m, RMSE 3000 m) than all alternatives (MB –3500 m, RMSE 3500 m).

Fig. 13.  Modelled plume information is shown for each of the five burns selected at the 2013 California Rim Fire (Burns 1, 2) and the 2017 Konza Prairie prescribed burn (Burns 3, 4, 5). Orange: BASE12, grey: BRIGGS12KMHR, yellow: SOFIEV12KMHR, light blue: PBL50012KM, green: BRIGGS4KMHR, dark blue: SOFIEV4KMHR, brown: BRIGGS4KM12HR. Top: mean bias (MB) for each simulation. Bottom: root mean square error (RMSE) change in plume tops compared with the observations.

However, the base case model typically provided PIH estimates with the largest errors and bias for the Konza Prairie prescribed grassland burns. For fires <80 ha, BASE12 with 3–4-h implementation performed best. The model parameterised with Sofiev performed the best for all fires <500 ha, but by a relatively small margin (<400 m). The base model showed the most improvement when changing the temporal allocation (20–1000 m). There was negligible difference between the grid resolution choices (<300 m) owing to the model vertical grid sizes (20 to +500 m, increasing with height). A change lower than 500 m would not be significant computationally unless that model layer was near the PBL. More importantly, an error larger than 500 m could have placed smoke in the incorrect model vertical layer, leading to increased uncertainty in downwind transport (e.g. Burn 5).

The vertical extent of a wildland fire smoke plume is determined by its classification, detection and the fire inventory and modelling system. How a fire is classified or input into a model will determine how, or even if, a plume rise algorithm will be applied to a given fire. Generally, in most AQMs (e.g. CMAQ), only one numerical method is used to simulate plume rise for wildfires and prescribed fires. In some cases, smaller fires will not have adequate plume rise associated with their emissions (Zhou et al. 2018). With a limited sample size of fires, our results suggest that there are benefits to model performance by treating a fire event smoke plume rise with multiple formulations selected based on the event’s size, duration and type (wild or prescribed fire).

For large wildland fires, night-time v. daytime plume rise proved to be an area of concern (Sofiev et al. 2013). The night-time high bias in BASE12-estimated plume height might have been due to missing night-time fire characteristics (e.g. model physics, meteorology, intensity, detection). The Burn 2 simulation in the Rim Fire revealed a higher plume height bias during the daytime and a low night-time bias comparable with Sofiev’s method (potentially due to the heat flux value from BlueSky near zero). Furthermore, the heat flux input from the BlueSky framework, although loosely correlated with satellite-derived FRP, tends to underestimate injection height (Kahn et al. 2007; Val Martin et al. 2010, 2012). This may occur if (i) satellite pixels of the order of 1 km² are only partly filled with a fire (Giglio et al. 2006); (ii) there is overlying smoke opacity (Kahn et al. 2008); or (iii) fire elements have non-unit emissivity, such as smouldering fractions (Val Martin et al. 2018). The suggested solution has been for modellers to multiply heat flux by a factor of 5 or more to match true plume buoyancy (Kahn et al. 2007; Ichoku and Ellison 2014). However, applying this adjustment here did not improve either simulation.

Another concern is that lingering smoke may be captured by observations but missed or ventilated out in model simulations. A full analysis of night-time plume rise for this study was not possible, because no fire observations were taken at night. Another difference between observation and models related to large wildland fires concerned the area burned. For example, plumes from Burns 1 and 2 produced significantly different average observed plume heights (4000–6000-m differences), but model simulations remained similar, in the range of 2000–3000 m. To improve the validation method for modelling, the use of the Soja et al. (2009, 2012) satellite-based method was helpful, because the plume tops derived were within 500–1000 m of those Peterson et al. (2015) and Saide et al. (2015) reported from the same fire.

For the Konza Prairie small prescribed grassland fires (<500 ha), selecting the proper temporal allocation of emissions and heat fluxes was even more critical than the actual algorithm. For the Rim Fire or larger fires (>5000 to 10 000 ha per day) that do not follow the typical diurnal cycle within an AQMs’ temporal profile, we suggest that it can be better to obtain hourly burned area from, e.g. Geostationary Operational Environmental Satellite (GOES), to estimate fire emissions. Model performance here improved with temporal allocation based on hourly observations but there is potential for greater improvement with the use of a proper treatment of sub-hourly fluctuations. When emissions were allocated based on size of fire and active burning phase, model performance improved significantly (+1000 m). Moreover, there was a connection between model plume heights and stability layers, consistent with earlier studies (e.g. Val Martin et al. 2018). Smoke pooled in stability layers. The stability layer in which a plume would be capped was determined by the strength of that stability layer compared with the energy or lift of an individual plume or buoyancy, a measure of stability calculated as buoyancy frequency. Burn 5 exhibits evidence for the importance of these layers, where plumes that penetrated the PBL continued vertically until reaching the next stable layer at ~3000 m.

Evaluating ways to improve the performance of contemporary AQMs should focus on providing numerically efficient algorithm choices and readily available supplementary data streams. To advance our understanding of smoke dispersion events, the model user should be allowed to select (i) plume rise algorithms based on fire type, and (ii) the time constraints and emission allocation windows. Based on this study, we suggest the following modelling improvements: (i) provide the option to use 4–6 h emissions allocation windows based on the fuel loading and field size, instead of 11 h from the detection time; (ii) if the fire start time is unknown, use the time of fire detection and place the time profile before and after the detection hour (e.g. for a fire detected at 1100 UTC, time profile would be 0900–1300 UTC); (iii) provide for the use of region-specific information such as weather, fuels and burn practices, because these results were highly reflective of the Flint Hills environment and associated fuels (grassland fields) of the Konza Prairie prescribed burn; (iv) allow for variable burn rates, because we found that the optimum estimates for intense to mildly intense fire spread rates were 20–35 ha h⁻¹ for the Flint Hills region, in contrast to 250–500 ha h⁻¹ for the Rim Fire.

Smoke plume injection height (PIH) is an important predictor of how smoke is transported and dispersed downwind of a wildland fire. Most air quality models rely on plume rise methods to determine the vertical allocation of emissions. In this study, we compared several methods with observed PIH. We used the CMAQ modelling system to investigate the impacts of model grid spacing, emissions temporal profile and three plume rise algorithms for five burn events (two wildfire and three prescribed).

Although approaches more advanced than the Briggs algorithm offer the potential to incorporate complex features of wildland fires, our results indicate that the Briggs algorithm performed comparably when provided improved inputs (MB of less than ±5–20% and RMSE lower than 1000 m compared with the alternatives). Our results indicate that the standard model formulation for plume rise has a high bias for large fires (MB range of 1000 to 3000 m). This high bias could be due to a high night-time bias, since the model had a pervasive low daytime bias. Predictions of PIH rely on correctly determining the PBL height. For prescribed grassland burns, the maximum PIH using all approaches was largely underpredicted. However, the bias and mean error, MB of 200 to 600 m and RMSE of 600 to 2000 m, were improved with a more resolved temporal profile (3 to 4 h compared with 12 h). Lastly, the assumption that a rising smoke plume will tend to rest near an STL is valid, but the onus is on the model’s ability to accurately capture those layers. For the large wildfire case (Rim Fire), Briggs placed >60% of the plume above the boundary layer, compared with 83% from observations and 30–55% from alternative models. For the Konza Prairie grassland prescribed fires <500 ha, the Briggs model placed <10% of the plume above the PBL, compared with 17% from observations and 10–35% from alternatives. Thus, the Briggs model performs better if not comparably for wildfires compared with the alternatives while for prescribed fires, Briggs requires adjustments to improve performance.

Based on these findings, we suggest the following modifications to the current air quality modelling system:

1. Compute plume rise for small fires (<500 ha). Many current models simply inject smoke from small fires into the boundary layer or the lowest model layer.

2. Assume a temporal profile that more closely matches the active burn period of a prescribed fire. Many models currently assume 12 or 24 h, but this tends to dilute the emissions and heat intensity of these fires.

3. Assume a fire-specific temporal profile, and if information is not available, apply one of two selectable options:

   1. Take the detection time and generate a temporal profile (e.g. if fire is detected by MODIS or GOES-16/17 at 1100 UTC, time profile applied if using 4 h will be 0900–1300 UTC for the burn).

   2. Use the burn rate or the regional average estimate of area burned (e.g. 35 ha h⁻¹ for the Konza Prairie prescribed burn).

Information on wild and prescribed fire is limited. Therefore, a lack of data persists as an inherent limitation of our study, which only considered the use of a small sample size of fires. But the resulting smoke emissions and impacts are substantial and clear enough to extrapolate relevant findings for future recommendations. Given the findings of this study, we suggest providing a temporal allocation tailored to the types of fires presented. To help with this matter, we urge and suggest that forest agencies, fire land managers and others collect information or create a survey of burn practices by region, season and type of biomass, as indicated in Section 7 of the US EPA National Emissions Inventory documentation (US EPA 2018). Fire plume rise modelling can be improved with the incorporation of space-based retrievals (Soja et al. 2012; Val Martin et al. 2018; Sokolik et al. 2019), databases for plume heights that are not typically modelled (e.g. intense pyrocumulonimbus: Lareau and Clements 2016; Peterson et al. 2017; Wilkins et al. 2020), ground-based measurements (Clements et al. 2018) and fuel consumption rate assessments (van Leeuwen et al. 2014). We recommend that future research explore the potential benefits of a hybrid approach combining multiple algorithms that can be used interchangeably based on conditions encountered by the model. For example, Wilkins et al. (2020) demonstrated a potential for increase in local ozone by 10–80 ppbv downwind of major biomass-burning sources. Plume rise algorithms are not limited to those presented here; in the future, we seek to investigate other methodologies, e.g. the 1-D model of Freitas et al. (2007), in order to evaluate them against observations. Lastly, we suggest an enhanced analysis of the impact of meteorological inputs on plume rise models, which could improve information available to prescribed burn decision-makers (e.g. burn windows, plume dispersion height and direction). These suggested improvements could help enable near-real time model predictions for fire modelling (Marsha and Larkin 2019; Shankar et al. 2019).

The authors declare that they have no conflicts of interest.

This research did not receive any specific funding. We are greatly appreciative of NASA funding for the CALIPSO work from the CloudSat and CALIPSO Science Team under grant number NNX16AM30G.