Inverse modeling of fire emissions constrained by smoke plume transport using HYSPLIT dispersion model and geostationary satellite observations

Smoke forecasts have been challenged by high uncertainty in fire emission estimates. We develop an inverse modeling system, the HYSPLIT-based Emissions Inverse Modeling System for wildfires (or HEIMS-fire), that estimates wildfire emissions from the transport and dispersion of smoke plumes as measured by satellite observations. A cost function quantifies the differences between model predictions and satellite measurements, weighted by their uncertainties. The system then minimizes this cost function by adjusting smoke sources until wildfire smoke emission estimates agree well with satellite observations. Based on HYSPLIT and Geostationary Operational Environmental Satellite (GOES) Aerosol/Smoke Product (GASP), the system resolves smoke source strength as a function of time and vertical level. Using a wildfire event that took place in the southeastern United States during November 2016, we tested the system's performance and its sensitivity to varying configurations of modeling options, including vertical allocation of emissions and spatial and temporal coverage of constraining satellite observations. Compared with currently operational BlueSky emission predictions, emission estimates from this inverse modeling system outperform in both reanalysis (21 out of 21 d; −27 % average root-mean-square-error change) and hindcast modes (29 out of 38 d; −6 % average root-mean-square-error change) compared with satellite observed smoke mass loadings.

1 Introduction

Burning biomass is one of the major factors affecting global air quality (Crutzen and Andreae, 1990). Fire smoke plumes directly emit both particles that can impact cardiopulmonary health and precursors (e.g., NO_x, SO₂, NH₃, and volatile organic carbons, VOCs) (Andreae, 2019) that react to form secondary particulate matter (PM) or other pollutants, such as ozone (Dreessen et al., 2016; Jaffe and Wigder, 2012; Mok et al., 2016; Singh et al., 2012; Valerino et al., 2017). In addition to their impact on air quality, fire emissions influence direct and indirect radiative transfer, aerosol formation, and the formation of cloud condensation nuclei, and they further interact with clouds and, eventually, with the biosphere and climate. Interaction between fire and the climate is an important factor affecting the future direction of the environment (Bowman et al., 2009). While high fire activities are affected by decadal-scale variation of the climate (Carvalho et al., 2011; Flannigan et al., 2005; Spracklen et al., 2009), aerosols released from fires and changed surface albedo due to burned areas are influential, as they disturb the radiative balance in the atmosphere (Liu et al., 2014).

Meeting National Ambient Air Quality Standards (NAAQS) requires US state agencies to understand the primary emission sources for particulate matter. Notwithstanding many states' continuous efforts to control in-state sources of pollutants, it is challenging to account accurately for fire emissions and the out-of-state transport of fire plumes. Due to the huge impact of fires on regional air quality, accurately forecasting their impact is an important public-service task, performed mostly by government agencies. The National Oceanic and Atmospheric Administration (NOAA) Smoke Forecasting System (SFS) was initiated after the large wildfire event in May 1998 (In et al., 2007; Rolph et al., 2009) to predict the movement of smoke from large wildfires (Rolph et al., 2009; Stein et al., 2009). The SFS uses the National Environmental Satellite, Data, and Information Service (NESDIS) Hazard Mapping System (HMS) (Ruminski and Kondragunta, 2006) and the US Forest Service's (USFS) BlueSky framework (Larkin et al., 2009) to detect fires and estimate emissions. The NOAA's HYSPLIT (Stein et al., 2015), a Lagrangian model which is designed to track air parcel trajectories, is then used to calculate transport, dispersion, and deposition of the emitted particulate matter. The SFS provides daily smoke forecasts over the continental United States, Alaska, and Hawaii to provide air quality guidance to the public.

Eulerian systems are also used for smoke forecasting. Smoke emissions from wildfires have been incorporated into the NOAA's National Air Quality Forecast Capability system as real-time, intermittent sources; this system has been forecasting regional air quality for surface ozone and particulate matter concentration since 2015 (Lee et al., 2017). The High-Resolution Rapid Refresh Smoke (HRRR-smoke; https://rapidrefresh.noaa.gov/hrrr/, last access: 29 August 2020) (Ahmadov et al., 2017) system also provides 36 h forecasts for the continental United States using the WRF-Chem modeling system with emissions derived from satellite-measured fire radiative power (FRP). Also, Chen et al. (2019) demonstrated an air quality forecast system over Canada that incorporates near-real-time measurement of biomass burning emissions to forecast smoke plumes from fire events. There are numerous other global smoke forecast systems (Chen et al., 2011; Larkin et al., 2009; Lee et al., 2019; Li et al., 2019b; Pavlovic et al., 2016; Sofiev et al., 2009).

Any improved smoke forecast system must confront several uncertainties, particularly fire (especially wildfire) emission amounts and their allocation spatially and vertically. In general, fire emissions may be estimated in one of two ways: bottom up and top down. Regarding the more traditional bottom-up approaches, fuel consumption is estimated as the product of burned size, pre-burn fuel loading of the fire-affected area, completeness of combustion, and emission factors (Seiler and Crutzen, 1980). Since emission factors are specific to the type of tree ablaze, a completely constructed database is very important in this approach. For example, Wiedinmyer et al. (2006), in estimating emissions from fires in North America, estimated fuel loading based on a combination of satellite and ground-collected data, such as Moderate Resolution Imaging Spectroradiometer (MODIS) thermal anomalies, the Global Land Cover Characteristics data set, the MODIS Vegetation Continuous Fields Product, and emission factors.

Recently, the global coverage of spaceborne instruments has encouraged top-down approaches. A fire's heat signature (that is, FRP) is detectable by satellite, and FRP can be used to estimate the rate of combustion (Giglio et al., 2003; Kaufman et al., 1998). Measurements of FRP from polar-orbiting sensors can both detect active fires and characterize their properties (Freeborn et al., 2009; Jordan et al., 2008; Schroeder et al., 2014). FRP data have been used to quantify biomass consumption, detect the locations of fire emission sources, trace gas and aerosol production (Ellicott et al., 2009; Kaiser et al., 2012; Vermote et al., 2009), and estimate the vertical extension of smoke plumes and other fire emissions (Val Martin et al., 2010). Several global fire emission databases – e.g., Global Fire Emissions Database (GFED), Fire Inventory from NCAR (FINN), Quick Fire Emissions Database (QFED), Global Fire Assimilation System (GFAS), Fire Energetics and Emissions Research (FEER), and Global Biomass Burning Emissions Product (GBBEP) – have been developed using one or both of the bottom-up and top-down approaches (Ichoku and Ellison, 2014; Kaiser et al., 2012; van der Werf et al., 2010; Wiedinmyer et al., 2011; Zhang et al., 2012). Both bottom-up and top-down approaches have their own advantages and limitations. While the bottom-up approach may provide detailed information based on a process- or fuel-specific estimation, it relies on various surveys that require significant time and resources. On the other hand, the top-down approach relies on observations of a few atmospheric variables such as radiation or aerosol optical properties, but it has an advantage from its timely availability and geographical coverage. Both approaches complement each other for better fire emission estimation.

This study extends the current capabilities of the NOAA SFS fire smoke forecast systems, most of which estimate fire emissions using the surface and thermal characteristics of detected fire locations. Transport pathways of smoke plumes are rarely considered in determining emissions strength, vertical extension, and temporal variation. This study aims to develop an inverse modeling system for fire emissions based on a Lagrangian model that can resolve these transport pathways using HYSPLIT simulations and satellite observations. Such an approach has been adapted to inversely estimate various emission sources, including greenhouse gas emissions (Kunik et al., 2019; Nickless et al., 2018; Turnbull et al., 2019), volcanic ash and sulfur dioxide emissions (Boichu et al., 2014; Zidikheri and Lucas, 2020), and radionuclide release from a nuclear power plant incident (Chai et al., 2015; Katata et al., 2015; Li et al., 2019a), but was rarely used in fire emission estimation (e.g., Nikonovas et al., 2017).

The remainder of the paper is structured as follows. Section 2 describes the model and satellite data used to detect fire locations and the transport of fire smoke. Section 3 concerns the methodology and structural design of the inverse modeling system. Results from a case study, sensitivity tests, and comparison with the currently operational system are presented in Sect. 4. Finally, Section 5 summarizes and discusses directions for future work.

2 Data

2.1 HMS and BlueSky

Consistent with the NOAA SFS system, the HMS data are utilized to detect wildfire information. HMS, developed as a tool to identify fires and their smoke emissions over North America in an operational environment, incorporates images from multiple geostationary and polar-orbiting environmental satellites, including the Geostationary Operational Environmental Satellite (GOES)-East/West, Suomi-National Polar-orbiting Partnership (NPP), MODIS, and Advanced Very High Resolution Radiometer (AVHRR) METOP-B, to provide the location and time of detected fires. Automated fire detection algorithms are first employed for each sensor, and then human analysts apply further quality control by examining visible channel imagery for false alarms and missed hotspots (Ruminski et al., 2008; Ruminski and Kondragunta, 2006; Schroeder et al., 2008).

The BlueSky system, developed by the US Forest Service, provides the first guess for fire emission estimation. As a modeling framework, BlueSky links several models of fire information, fuel loading, fire consumption, fire emissions, and smoke dispersion (Larkin et al., 2009; Strand et al., 2012). For the original NOAA SFS system, BlueSky emissions are used as inputs for the dispersion model. In this study, we use the BlueSky emission rate as an initial guess before applying the inverse modeling system.

2.2 GASP

The GOES Aerosol/Smoke Product (GASP) is a retrieval of the aerosol optical depth (AOD) using GOES visible imagery (Kondragunta et al., 2008; Prados et al., 2007). This product is available at 30 min intervals and 4 km × 4 km spatial resolution during the sunlit portion of the day. The Automated Smoke and Tracking Algorithm (ASDTA; https://www.ssd.noaa.gov/PS/FIRE/ASDTA/asdta_west.html, last access: 29 August 2020) detects smoke associated with detected fire source locations. For each pixel, the radiative signatures of an aerosol layer (e.g., dust and smoke) are determined by the scattering and absorption properties of the aerosol. ASDTA also utilizes a pattern recognition technique to isolate smoke aerosols from other type of aerosols, so that it can recognize plumes transported far from fire sources. The GASP product is particularly useful for tracking fast-moving plumes, which polar-orbiting sensors often cannot detect since they provide only one daily image. Since the ASDTA product is a part of the GASP product, in this study, GASP and ASDTA indicate total AOD and smoke AOD, respectively. Hourly data were used for the inverse system.

2.3 HYSPLIT

HYSPLIT computes air parcel trajectories and the dispersion or deposition of atmospheric pollutants (Stein et al., 2015). It has been widely used to simulate pollutant events, including volcanic ash, smoke from wildfires, radioactive nuclei dispersion, and emissions of anthropogenic pollutants. For the inverse modeling system, we used the Transfer Coefficient Matrix (TCM) approach. The unit source calculations give the dispersion factors from the release point for every emission period to each downwind grid location, defining what fraction of emissions are transferred to each location varying as a function of time. This is defined as the TCM (Draxler and Rolph, 2012). The TCM is computed for inert and depositing species and, when quantitative air concentration results are required, the final air concentration is computed in a simple post-processing step that multiplies the TCM by the appropriate emission rate. Results for multiple emission scenarios are easily created and may be used to optimize model results as more measurement data become available.

For the inverse system, 120 h HYSPLIT simulations were conducted daily, starting from 06:00 Z using North American Model 12 km meteorology (NAM12), at each fire source location provided by the HMS fire detection information. A total of 50 000 particles were released for each simulation, and dispersed concentrations were vertically integrated up to 5000 m onto 0.1^∘ spatial grids. Hourly outputs were integrated to match with satellite observational data. HYSPLIT modeling options were configured to be consistent with the SFS system, including options for dry and wet depositions (i.e., 0.8 µm diameter with 2 g cm⁻³ density) (Rolph et al., 2009). In this paper the integrated mass loading of particles from HYSPLIT simulation and satellite products will be compared with each other. For satellite products (e.g., GASP, ASDTA, and MODIS), smoke is converted from AOD using a simple conversion factor (i.e., 1 AOD = 0.25 g m⁻²) which is compatible to 4 m² g⁻¹ mass extinction efficiency (Nikonovas et al., 2017). Although we used a single conversion factor for the study, the actual conversion factors may vary in time and space (i.e., 3.9–5.3 m² g⁻¹) (Chand et al., 2006; Hobbs et al., 1996; Ichoku and Ellison, 2014; Nikonovas et al., 2017; Reid et al., 2005). Therefore, applying more realistic conversion factors and their uncertainties into the system would be another factor in the future improvement of the system.

For HYSPLIT runs, smoke indicates the sum of dispersion simulations (i.e., TCM runs in concentration unit) multiplied by emissions for each source. Since we have integrated dispersion model outputs (in density unit) up to 5000 m height, the results shown in the study are obtained by multiplying the column height (i.e., 5000 m) and are demonstrated as total mass loading for a column (kg m⁻²).

3 Methodology

3.1 Overview

A HYSPLIT inverse system was built and successfully applied to estimate the cesium-137 releases from the Fukushima Daiichi Nuclear Power plant accident in 2011 (Chai et al., 2015). It was then modified to estimate the volcanic ash source strengths, vertical distribution, and temporal variations by assimilating MODIS satellite retrievals of volcanic ash clouds while using ash cloud top height information as well (Chai et al., 2017). It was found that simultaneously assimilating observations at different times produces better hindcasts than only assimilating the most recent observations. In this application, the HYSPLIT-based Emissions Inverse Modeling System for wildfires (HEIMS-fire) is designed to assimilate satellite observations to generate wildfire emission estimates. The smoke plume transport and dispersion, captured by frequent geostationary satellite retrievals, can be used as constraints to obtain the smoke emission estimates. In the system, a cost function quantifies the differences between HYSPLIT model predictions and satellite-observed AOD, weighted by model and observation uncertainties. Minimizing the cost function by adjusting emission rates at different fire locations and at several different release heights thereby provides the fire emission estimates.

3.2 Cost function

Taking a top-down approach, unknown emission terms are obtained by searching for the emissions that provide the model predictions that most closely match the observations. With fire locations mostly identified by the HMS system, unknown emission rates at these specified locations remain undetermined. At each fire location, released smoke can reach different heights under various fuel-loading and meteorological conditions. In addition, emission rates may vary significantly with time. Thus, the unknown elements of the inverse problem are the emission rates q_iktat each wildfire location i at different heights k and time periods t. The cost function F is defined as follows:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}F& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{t=\mathrm{1}}^T\sum _{k=\mathrm{1}}^K\sum _{i=\mathrm{1}}^I{\displaystyle \frac{{\left(q_{ikt}-q_{ikt}^b\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{ikt}^{\mathrm{2}}}}\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{n=\mathrm{1}}^N\sum _{m=\mathrm{1}}^M{\displaystyle \frac{{\left(c_{nm}^{\mathrm{h}}-c_{nm}^{\mathrm{o}}\right)}^{\mathrm{2}}}{{\mathit{\epsilon }}_{nm}^{\mathrm{2}}}}+F_{\mathrm{other}}\end{array},\end{array}$$

where $c_{nm}^{\mathrm{o}}$ is the mth gridded satellite observation (e.g., GASP ASDTA smoke mass loading) at time period n and $c_{nm}^{\mathrm{h}}$ is its HYSPLIT counterpart.

A background term is included to measure the deviation of the emission estimate from its first guess, $q_{ikt}^{\mathrm{b}}$, obtained from the operational BlueSky emission computation. The background term ensures the problem remains well-posed even with the limited observations available in certain circumstances. The background error variance ${\mathit{\sigma }}_{ikt}^{\mathrm{2}}$ measures uncertainties in $q_{ikt}^{\mathrm{b}}$. Pan et al. (2020) compared six global emission estimates and found that the total emission differs by a factor of 3.8. However, emission estimations at specific locations and times can have much larger errors. In addition, the vertical distribution of the smoke emissions is difficult to determine and this adds even more uncertainties to the emission estimates. We chose a large uncertainty for the background term as ${\mathit{\sigma }}_{ikt}=\mathrm{1000}\times q_{ikt}^{\mathrm{b}}+\mathrm{1000}$ kg h⁻¹ at all locations and heights to minimize the adverse impact of inaccurate BlueSky emission estimates. The observational error variances, ${\mathit{\epsilon }}_{nm}^{\mathrm{2}}$, represent uncertainties in both the model and observations, as well as the representative errors. Kondragunta et al. (2008) indicated that GOES aerosol retrievals over land were expected to have uncertainties within 0.15τ±0.05, where τ is the AOD. Paciorek et al. (2008) showed a better performance of GOES aerosol retrievals in eastern US than in western US. Green et al. (2009) demonstrated that GOES AOD correlates best with AERONET in autumn (September to November) than in other seasons. They showed that the RMSE was 0.060 in autumn while the average for all seasons is 0.149. Considering the better performance in the eastern US and in November, AOD uncertainties of 0.10τ±0.06 are assumed in this paper. A slightly larger additive component of the AOD error is chosen to include the effects of the representative errors and model errors that do not vary with the observed AOD values. F_other refers to the other regularized terms that can be included in the cost function. For instance, Chai et al. (2015) has a temporal smoothness penalty term to avoid abrupt changes in the temporal profile of the release rates. While this optimization problem could be solved to obtain optimal emission estimates using many minimization tools, we used the limited-memory Broyden–Fletcher–Goldfarb–Shanno (BFGS) algorithm (Zhu et al., 1997).

3.3 Inverse system

The HEIMS-fire system is designed to conduct a two-step operation: (1) estimation of fire emission using an inverse system and (2) forecast modeling of fire smoke using estimated fire emissions. The inverse system utilizes observations and modeling systems available from multiple agencies. We aim to objectively and optimally estimate wildfire smoke sources' strength, vertical distribution, and temporal variations by assimilating GASP AODs. Figure 1 summarizes the system's incremental stages of data processing, listed below with the required data (and their providing agencies):

1. fire detection from the Hazard Mapping System (NESDIS),

2. HYSPLIT simulations with unit emissions at different locations and release heights,

3. construction of the Transfer Coefficient Matrix with available observations,

4. initial guess for fire emissions (BlueSky, US Forest Service),

5. cost function minimization to estimate smoke emissions,

6. smoke forecast using adjusted smoke emissions.

By minimizing a cost function, the HEIMS-fire system provides adjusted fire emissions that can describe realistic smoke plumes. Results in the following section show that the assimilated smoke plumes agree well with satellite observations. The system requires a first guess as input before it can start to minimize the cost function. Selection of this input is usually critical both for the performance of the minimization calculation and for the final output.

Figure 1 Schematic diagram of the HYSPLIT-based Fire Emission Inverse Modeling System.

Here we also explain the naming conventions for temporal coverage of emission estimation processes and forecasting processes. This inverse system is designed to estimate fire emissions on the target day by analyzing past and present smoke fields and then utilizes them to forecast the future smoke field. The assimilation days (i.e., aday = 0, −1, −2) (see Fig. S1 in the Supplement) indicate the temporal coverage of dispersions and constraining observations. For a target day, 13 November, inversions are conducted using HYSPLIT dispersion simulations and ASDTA observations for 24 h (i.e., aday = 0), 48 h (i.e., aday = −1), 72 h (i.e., aday = −2), and 96 h (i.e., aday = −3). Estimated fire emissions are used to simulate fire smoke for 13 November (i.e.,fday = 0; reanalysis), and the same amount of fire emissions are used in forecast mode for 14 November (i.e., fday = +1) and 15 November (i.e., fday = +2). Application of forecast mode will be further discussed in Sect. 4.5.

4 Results

4.1 Case study

A case study using a November 2016 wildfire event was conducted to test the performance of the HEIMS-fire system. This fire event was a series of wildfires in the southeastern United States in October and November 2016. The US Forest Service reported at least 80 000 acres (∼324 km² burned from 23 October to 9 December 2016. For the case study, we focused on the fire event that occurred in Georgia, South Carolina, North Carolina, and the adjacent states from 10 to 17 November 2016. Figure 2 shows an example of fire smoke detected from a MODIS true-color image and three AOD products from MODIS, GASP, and ASDTA on 10 November 2016. Wildfires in the Appalachians of northern Georgia, western North Carolina, and eastern Tennessee began to produce large smoke plumes moving southeast. Numerous smoke plumes could be seen from active wildfires burning across the region. Changes in the fire events from 8 to 19 November 2016 are also shown as MODIS true-color images in the Supplement (Fig. S2).

Figure 2 Detection of fires over the southeastern region of the United States on 10 November 2016. True-color images from MODIS (left), MODIS AOD (top right), GASP AOD (middle right), and ASDTA AOD (bottom right) are shown. MODIS true-color images and AOD are obtained from https://earthdata.nasa.gov/ (last access: 29 August 2020), and GASP and ASDTA AOD are obtained from NOAA NESDIS.

4.2 Model configuration

The month of November 2016 saw fires nationwide, although the most extensive fires happened in the southeastern US region. We considered four geographic domains in determining fire source inputs, as shown in Fig. 3. Red dots indicate HMS fire detections during November 2016. The results of the sensitivity test using these domains are discussed in Sect. 4.4.

Figure 3 Geographical coverage of case study domains. Red dots indicate HMS fire detections during November 2016. The four labeled domains indicate spatial coverage of fire source inputs for the inverse system used in the sensitivity tests.

The inverse modeling system was tuned using various sensitivity tests. A series of twin experiments was conducted to test the range of uncertainties that comes from the system design. A twin experiment is an idealized modeling test in which we assume that the modeled world adeptly mimics the real world. Using a true solution for the situation, we can test the system's capability to reproduce the true answer. We tested uncertainties of the system across multiple scenarios and four types of potential uncertainty (vertical allocation, temporal coverage, spatial coverage, and impact of observation errors) in these twin experiment cases. Detailed descriptions of the twin experiments and sensitivity test are available in the Supplement.

4.3 Emission estimation

Fire emissions and their vertical distributions for each detected fire location were estimated using the HEIMS-fire system, inversely modeled from ASDTA AOD data as described above. Locations and times of fires detected by HMS were used to initiate HYSPLIT simulations, with emissions released over six layers (100, 500, 1000, 1500, 2000, and 5000 m). On 11 November, 46 fire locations were identified within the assimilation domain (domain 1 in Fig. 3). Thus, the TCM was established based on 276 HYSPLIT simulations (46 fires × 6 release altitudes) and GASP AOD observations. Emissions rates calculated from the BlueSky system were used as an initial condition. Emissions were evenly distributed to all layers used in the system.

Figure 4 Scatterplot comparison between initial and assimilated smoke mass loading using adjusted fire emissions. A 48 h observation (aday = −1) and 6-layer plume release configuration is used (unit: $\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻²).

Minimizing the cost function results in the estimation of fire emissions. Figure 4 shows the agreement between the modeled and observed mass loading from the initial to the adjusted emission estimates. Table 1 presents summary statistics for the changes in reconstructed smoke mass loading from the initial guess to the adjusted emissions. In the end, estimated fire emissions were combined to reconstruct fire smoke plumes. Using adjusted fire emissions, we can reconstruct the integrated smoke columns as a sum of adjusted emissions, q_ikt, applied to each TCM, T_ikt:

$$\begin{array}{}\text{(2)}& c\left(nm\right)=\sum _{ikt}{q}_{ikt}\cdot T_{ikt}(n,m),\end{array}$$

where i, k, and t denote spatial, vertical, and temporal allocation of emission sources and m and n denote the location and time of receptor (i.e., observations), respectively. Figure 5 presents the spatial distribution of reconstructed fire smoke mass loading for the case study in terms of column-integrated density. We applied estimated fire emissions to TCM runs for each detected fire location and vertical release height and then merged them into one hourly concentration field. Reconstructed smoke plumes (i.e., integrated dispersion outputs in 0–5000 m height) show a good agreement with observed smoke (Fig. 5). The first and second columns compare ASDTA and HEIMS smokes for spatially and temporally matching pixels, and the third column shows the full spatial coverage of HEIMS smoke for daytime (around 08:00–17:00 LT, local time) when ASDTA data are available. The 17 November output in Fig. 5 shows how the system responds when observations are limited or missing, while it still provides a robust result by honoring the initial guess information. On 17 November, no ASDTA AOD was provided from the satellite operation. Under 48 h configuration (i.e., aday = −1), the inverse system still produced reasonable outputs using limited observations (16 November) and initial guess emissions (16 and 17 November). This case hints at the importance of both a traditional (e.g., BlueSky emissions) and new inverse system. They complement each other, as one provides the latest data assimilation technique while the other provides prior information and backup stability in a contaminated environment (e.g., excessive cloud cover).

Figure 5 Comparison of observed and reconstructed smoke plumes during 10–17 November 2016. Smoke mass loading from ASDTA (left column) and reconstructed HEIMS smokes for ASDTA-matching (middle column) and for daytime (right column). A 48 h observation (aday = −1) and 6-layer plume release configuration is used.

4.4 Sensitivity tests

Similar to the twin experiments, we conducted a series of sensitivity tests to investigate how the inverse model responds to changes in input data and various configurations of the modeling framework. This will be achieved by focusing on variation in temporal coverage, spatial coverage, and vertical allocation of smoke plumes.

First, we changed the assimilation time windows from 1 d (24 h) to 4 d (96 h). Since the impact of fire emissions easily translates over multiple days, we tested how temporal coverage affects system results. The “1 d” (aday = 0) simulation is run through the inverse model using dispersions and observations for the target day, while the “2 d” simulation uses 2 d (i.e., 48 h) of dispersions and observations (aday = −1). For this test, all observations within the assimilation time windows were selected for the assimilation and the evaluation. The results are shown in Fig. 6a, while the correlation and error statistics are summarized in the top section of Table 2, i.e., (A: 24 h, O: 24 h, E: 24 h), …, (A: 96 h, O: 96 h, E: 96 h). With the exception of 10 and 11 November, in the early stage of the fire event, both the correlation coefficient (R) and normalized root-mean-square error (NRMSE) were improved by the use of more days (i.e., 3 or 4 d) of dispersions and observations for the inverse model. This makes sense because emissions from multi-day fire events spread out and affect concentrations over proceeding days.

Figure 6 Sensitivity tests for (a) temporal coverage, (b) layer selection, and (c) spatial coverage. Temporal coverage of inputs is tested for between 24 h (aday = 0) and 96 h (aday = −3). Selection of layers is tested using two to seven layers among 100, 500, 1000, 1500, 2000, 5000, and 10 000 m. Beginning from two layers at 100 and 500 m, the next higher layer is added for each test. Spatial coverage is tested through domain 1 to 4. The thick blue boxes indicate the intensive fire episode period, 10–17 November 2016.

Table 1 Performance evaluation. Statistics of modeled smoke mass loading ($\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻²) using initial and top-down estimated fire emissions.

Figure 7 Scatterplot comparisons between ASDTA smoke and HEIMS smoke (top row) and ASDTA and SFS smoke (bottom row) for forecast days: fday = 0, +1, +2.

A series of additional simulations were also conducted to test the system's sensitivity to the selection of observations for the assimilation and the evaluation. In these tests, we investigated combinations in assimilation time (“A” in Table 2), observational time (“O”), and evaluation time windows (“E”). Results are also summarized in Table 2. For a fixed assimilation time period (i.e., A: 96 h), using a shorter observational time window resulted in a better result. It is reasonable because we expect a better fitting with a smaller number of data points. However, it can easily be exposed to an overfitting problem if available data for the assimilation are too small.

Table 2 Sensitivity test for temporal coverage. Uses of assimilation (A), observation (O), and evaluation (E) time windows (24 to 96 h) were tested, and statistics, NRMSE, and R were compared. The best performance is marked in bold.

Second, we tested the layers at which fire emissions are initiated in the model. As expected, including more layers results in better statistics, since the transport and dispersion of each smoke plume can vary with the altitude to which their fire emissions are allocated. We tested the model's uncertainties on layers' maximum extension and resolution, with varying selections of two to seven layers at 100, 500, 1000, 1500, 2000, 5000, or 10 000 m. To test the maximum extension, starting from two layers (i.e., with emissions released at 100 and 500 m), we added the next higher layer over six test runs to investigate the effect of maximum extension of smoke plume. Figure 6b shows the results, and error statistics are summarized in Table 3. Including the 5000 m layer resulted in noticeable changes especially, implying the potential benefit of including high-level transport for specific days. Since the 5000 m layer is above typical planetary boundary layer height, emissions injected at this level experience different physical characteristics. Smoke lofted into the free troposphere is less affected by turbulence and scavenging, and easily transports hundreds or thousands of kilometers downwind because of the higher wind speeds. The addition of the 5000 m layer would better represent the potential long-range transport. Smoke plume rise is one of traditionally important questions in smoke modeling, so further research on the topic is warranted. The effects of the layer resolution were also tested. Starting from two layers (i.e., 100 and 5000 m), we added intermediate layers up to six layers and evaluated their performances (Table S3). As expected, including more layers resulted in the better statistics, but its improvement was not significant after four layers.

Table 3 Sensitivity tests for selection of two to seven layers among 100, 500, 1000, 1500, 2000, 5000, and 10 000 m. Beginning from two layers at 100 and 500 m, the next higher layer is added for each test. The best performance is marked in bold.

In the third test, we varied the spatial coverage of input fire information. Although wildfire impacts easily spread by long-range transport, we could not include all the global fire information due to limited computational resources. We therefore tested different spatial domains of fire locations to evaluate what spatial coverage of wildfire detection information is required to estimate fire emissions. Fire sources inside domain 1 through 4 (Fig. 3) were tested in the assimilation constrained by ASDTA AOD inside domain 1. Figure 6c and Table 4 show correlation and error statistics from the sensitivity test of spatial coverage. In most days, we have better results when we include fire emission sources at least within domain 2. It makes sense considering the effects of transported fire plumes from Mississippi and Louisiana (Fig. 3). Maximizing geographical coverage (e.g., domain 4) did not always result in the best performance in our case study. This result, however, should be taken carefully because we do not have strong fire activities outside domain 2 in our study case. Strong long-range transport cases, typically from the northwestern US, Canada, and Alaska, would have bigger impacts.

Table 4 Sensitivity tests for spatial coverage. Use of observational data in terms of spatial availability is tested. Domains 1–4 are shown in Fig. 3. The best performance is marked in bold.

4.5 Hindcast and operation

In this section, we conducted a HEIMS system for hindcast mode, and compared it with operational products from the SFS system. Both SFS and HEIMS use fire detection from HMS for consistency, and HEIMS uses SFS fire emissions for initial guess information. SFS simulates 72 h dispersion of fire smoke for every day in November 2016, which is consistent with fday = 0, +1, +2 of the HEIMS hindcast simulations (as described in Fig. S1).

Notable differences in the configuration of SFS and HEIMS are plume rise estimation, temporal resolution of fire emissions, fire decaying assumption, and meteorology. While SFS computes plume rise using the Briggs' equation (Arya, 1998; Briggs, 1969), which assumes an air parcel's rise is based only on the buoyancy terms, HEIMS determines fire emissions' vertical allocation using an inverse system. At the initial guess, SFS fire emissions are evenly distributed in all layers. Current HEIMS assumes daily emission variation compared to hourly emissions of SFS. Also, SFS assumes 75 % of emissions still happen at the same location the next day; the HEIMS uses 50 % decay assumption after sensitivity tests, which will be discussed in the next section. For HEIMS simulation, we used aday = −1 (2 d temporal coverage) for the simulations shown. On the other hand, HEIMS would benefit from a better meteorology. Although both systems use the NAM12 forecast meteorology, HEIMS hindcast used the first 24 h portion of everyday forecast cycle, and SFS used 72 h forecast.

For forecasting days, smoke is estimated as the summation of impact from previous days and new emissions on the target days. For example, smoke at fday = +2 can be reconstructed as

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{S}_{\mathrm{fday}=+\mathrm{2}}& =q_{f=\mathrm{0}}\cdot {\mathrm{TCM}}_{f=\mathrm{0}}+q_{f=\mathrm{0}}\cdot p\cdot {\mathrm{TCM}}_{f=+\mathrm{1}}\\ & +q_{f=\mathrm{0}}\cdot p^{\mathrm{2}}\cdot {\mathrm{TCM}}_{f=+\mathrm{2}}\end{array},\end{array}$$

where q and p denote emissions and persistency rate, respectively. The persistency rate, p, assumes the change of future day emissions. Its role will be discussed in the next section.

Figures 7 and 8 demonstrate simulated fire smoke by SFS and HEIMS on 11 November and 2 d forecasts (hindcasts for HEIMS) for 12 and 13 November. Both systems reproduced the smoke well in their general patterns and intensity, as shown in ASDTA AOD and MODIS true-color image (Fig. 8).

Figure 8 Spatial distributions of observed and forecasted fire smoke plumes on 11–13 November 2016: true-color image from MODIS (top row), ASDTA smoke (second row, converted from AOD), HEIMS smoke hindcast (third row), and SFS smoke forecast (bottom row, from operation) are shown. MODIS true-color images are obtained from https://earthdata.nasa.gov/ (last access: 29 August 2020).

As expected, the HEIMS shows better agreement at fday = 0, as the fire emissions were assimilated on the day. For fday = +1 (i.e., 12 November), the HEIMS shows better agreement in RMSE and mean bias (RMSE = $\mathrm{58.1}\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻² and bias = $-\mathrm{22.4}\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻², compared with RMSE = $\mathrm{63.0}\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻² and bias = $-\mathrm{42.2}\times {\mathrm{10}}^{-\mathrm{6}}$ kg m⁻²) while SFS has better slope. For fday = +2, HEIMS is better in mean bias but worse in RMSE and R. Table 5 summarizes RMSE statistics from HEIMS and SFS for each day of November 2016. In most of the days, HEIMS posts better statistics compared with SFS, implying the potential benefit of system improvement by adding an additional observational constraint. For the comparisons of HEIMS hindcast and SFS operational simulations, HEIMS system shows better performance in both hindcast days ($\mathrm{16}/\mathrm{19}=\mathrm{84}$ % on fday = +1 and $\mathrm{13}/\mathrm{19}=\mathrm{68}$ % on fday = +2).

Table 5 Performance statistics and RMSEs for HEIMS and SFS smoke mass loading for 0–2 forecast days (e.g., fday = 0, +1, +2). For each forecast day, the better performance from the two systems is marked in bold. Statistics for days without observations or without operational outputs are not available (unit: 10⁻⁶ kg m⁻²).

Change of fire activity is also a problem for both systems. If there is considerable change of fire activity for fday = +1 and +2, the forecast will result in worse performance. If fire activities increase, the simulated smoke from HEIMS will be underestimated, and if fire activities decrease, the HEIMS system will overestimate the impact of smoke. Therefore, information of next day fire activity or fire duration will be important for an accurate fire smoke forecast system, which will be discussed further in the next section.

4.6 Persistency of fire activity

The selection of the persistent rate of daily fire emissions (i.e., persistency = 1 − decaying_rate) and its importance to the smoke forecast system's performance are discussed here. In our current systems, both SFS and HEIMS, we use a simple assumption of fire emission change for next day. Figure 9 shows how the HEIMS responds to the selection of a persistent rate for forecast days fday = +1 and +2. We applied five different persistent rates ranging from 0 % to 100 %; persistency = 0 % assumes no new fire occurs, and persistency = 100 % assumes the same amount of fire emissions released at the same location from the previous day. For the top row of Fig. 9, simulated smokes in fday = +1 and +2 solely originate from fires in fday = 0. On the other hand, persistency = 100 % simulation demonstrates accumulated impacts (target day and previous days), showing denser smokes estimated compared with persistency = 0 %.

Figure 9 Sensitivity test with varying persistent rates. For p=75 %, we assume that 75 % of fday = 0 emissions remain in fday = +1 and 56.25 % (0.75²) emissions remain in fday = +2.

An implication from these comparisons is the importance of persistent rate selection. Indeed, better smoke forecasting may require improvement via two separate steps. The first one is to estimate today's emissions, which can be improved by better assimilation techniques than those we introduce in this paper. The second issue is to predict fire activity, which is more related to the studies of fire behavior. In more detail, we need to predict how long existing fires persist, and also to predict the occurrence of new fires, which may pose the greatest difficulty for daily operational systems. Without a better understanding and modeling of fire behavior, the current system has to rely on the empirical solution. For our case study, choosing a persistent rate of 50 % d⁻¹ produced the best result, but it warrants further study with a long-term data set being used in an operational system. Prediction of wildfire consistency based on the change of meteorological conditions, such as the Fire Weather Index (FWI, https://cwfis.cfs.nrcan.gc.ca/background/summary/fwi, last access: 29 August 2020), will be a good indicator for the change of fire emission. Without this kind of fire behavior model, the fire smoke forecast system could be limited.

5 Summary and discussion

Accurate estimation of emissions from wildfire sources is critical to improving the performance of air quality forecast systems. Wildfire emissions may be estimated based on fire detection information from the surface (bottom up) or instead based on the intensity of radiance measured from space (top down). This study extends the top-down approach by applying an additional constraint, i.e., transported smoke plume recorded by geostationary satellites. We developed an inverse modeling system to estimate wildfire smoke emissions over North America using NOAA's HYSPLIT and GOES Aerosol/Smoke products. This HEIMS-fire resolves the strength of smoke sources as a function of time and vertical level. The system adjusts estimated wildfire smoke emissions until they agree well with satellite observations.

We conducted numerous sensitivity tests, varying the temporal, vertical, and spatial coverage of the input data sets used to initiate the inverse system. Results are mostly consistent with general expectations based on the characteristics and behavior of fire events. As transport from previous days can impact large areas, including multiple days of observations to constrain fire emissions yields statistically better results. Including more vertical layers also leads to better results; for example, including the 5000 m layer resulted in the best improvement. Spatial coverage was tested in terms of four different domains, and while this particular test presented no solid conclusion, adding more information in general yielded better results, as expected. It also should be noted that the uncertainties of the emission estimation and the smoke forecasts thereafter are not quantified in this study. An ensemble of HYSPLIT predictions using different meteorological inputs will be used to estimate the uncertainties of the results in the future.

For operational purposes, including an additional constraint that extends current smoke forecast systems to use smoke plume transport has clear advantages. Future study could improve this approach in several respects. First, the conversion of AOD to smoke mass loading is simply empirical; secondary formation of PM is not considered. Omission of chemical reaction models is the basic characteristic of trajectory- or dispersion-based models compared with Eulerian, full chemistry models. Applying estimated emissions to a chemistry dispersion model could improve results. Second, the system is highly dependent on the quality of constraining observations. Use of the latest satellite instruments could further improve results. Third, we have not yet included surface observations into the inverse system. Utilizing both surface and more columnar observations from other satellite systems will improve the model performance. Fourth, we only used the target day fire emissions for the smoke forecast. Since fire smoke lasts several days, including the previous days' emissions will enhance the background effect.

This study aimed to improve the operational smoke forecast by providing accurate fire emission inputs. Unfortunately, the GASP product was discontinued in early 2018. However, the concept of minimizing a cost function based on satellite observations remains robust and can be applied to other data sets. In particular, we plan to apply the GOES-R Advanced Baseline Imager (ABI) product to constrain the extension of fire smoke.