Wildfire burn severity and emissions inventory: an example implementation over California

Wildfire severity is a key indicator of both direct ecosystem impacts and indirect emissions impacts that affect air quality, climate, and public health far beyond the spatial footprint of the flames. Comprehensive, accurate inventories of severity and emissions are essential for assessing these impacts and setting appropriate fire management and health care preparedness strategies, as is the ability to project emissions for future wildfires. The frequency of large wildfires and the magnitude of their impacts have increased in recent decades, fueling concerns about decreased air quality. To improve the availability of accurate fire severity and emissions estimates, we developed the wildfire burn severity and emissions inventory (WBSE). WBSE is a retrospective spatial burn severity and emissions inventory at 30 m resolution for event-based assessment and 500 m resolution for daily emissions calculation. We applied the WBSE framework to calculate burn severity and emissions for historically observed large wildfires (>404 hectares (ha)) that burned during 1984–2020 in the state of California, U.S., a substantially more extended period than existing inventories. We assigned the day of burning and daily emissions for each fire during 2002–2020. The framework described here can also be applied to estimate severity for smaller wildfires and can also be used to estimate emissions for fires simulated in California for future climate and land-use scenarios. The WBSE framework implemented in R and Google Earth Engine can provide quick estimates once a desired fire perimeter is available. The framework developed here could also easily be applied to other regions with user-modified vegetation, fuel data, and emission factors.


Introduction
Though some wildfires are essential to the Earth system and have beneficial effects on ecosystems that have evolved with fire, human influences have altered wildfire risk via increasing ignitions, climate change, and cumulative impacts of land management and use (Bowman et al 2009, Stephens et al 2013, Balch et al 2017, Abatzoglou et al 2018, McLauchlan et al 2020. In addition to ecosystem effects, wildfires directly impact people through property losses, injuries, mortality, and evacuations, while smoke and its health impacts are how the largest numbers of people are affected by wildfires through the regional transport of smoke in the atmosphere (Finlay et al 2012, Reid et al 2016, Bowman et al 2017, Thomas et al 2017, Burke et al 2021. The two most commonly used criteria for analyzing fire effects are burned area and burn severity (Key and Benson 2006, Kolden et al 2012, Meng and Zhao 2017. The degree to which a fire-disrupted ecosystem has changed is referred to as burn severity (Eidenshink et al 2007, Keeley andSyphard 2019). Because of differences in wind, topography, fuel conditions, and other factors, burned areas typically consist of complex landscape mosaics of low, moderate, and high burn severity (Turner et al 1994, Perry et al 2011, Birch et al 2015, Prichard et al 2020. Thus, landscape-scale quantification of burn severity is vital to identify the immediate and long-term influence on fire regimes, air quality, emissions, and other impacts. Passive spectral reflectance satellites and field observation have been used to determine burn severity. The change in Normalized Burn Ratio (NBR) between pre-and post-fire satellite images (dNBR) has become the primary index for remote sensing of severity. In the United States, the Monitoring Trends in Burn Severity Project (MTBS) (Eidenshink et al 2007), which is derived from long-running Landsat earth observing sensors, mapped burn severity for large fires dating back to 1984. However, the MTBS products have limitations, with sources of potential errors (e.g. subjective and not field validated severity classification thresholds) and a 2 year delay in including new fires in the database, that make the them less than ideal for research . The composite burn index (CBI) field method tabulates severity within 30 m plots by averaging measurements of several rating factors on a continuous scale ranging from 0.0 (unburned) to 3.0 (high severity) (Key and Benson 2006). On-the-ground burn severity was assessed following the CBI protocol only for 234 fires that burned in the contiguous United States (CONUS) during 1996-2018 (Picotte et al 2019).
Emission inventories are a critical input for a variety of modeling applications that seek to understand how emissions affect air quality, climate, and human health (Larkin et al 2014, Jaffe et al 2020, Pouliot et al 2020. Several global biomass burning emissions inventories have been produced and updated in recent decades. 'Bottom-up' inventories such as the Global Fire Emissions Database (GFED; van der Werf et al , 2006van der Werf et al , 2010van der Werf et al , 2017 (French et al 2014) provides tools for assessing wildland fire emissions in area-of-interest or for an identified fire event. It also provides pre-calculated emissions data on state or country level on monthly/yearly scale with the earliest records back to 1984. Comparisons of fire emissions inventories have revealed variability in inventory results, with each having advantages and disadvantages (Larkin et al 2014, Faulstich et al 2022. Despite all efforts, there is an urgent need for burn severity and emissions inventories at the scale of actionable management (e.g. fuel reduction treatments) and over longer temporal periods.
Here, we introduce the WBSE and present a detailed description of the inventory model, initial results using California historical large wildfires as an example, comparison with other inventory estimates, and a discussion of uncertainties. The primary objective of this manuscript is to describe the model, and our comparison with other emissions inventories is limited in scope. Similarly, the historical reconstructions of wildfire severity maps and emissions are analyzed in more detail in a companion paper. This study benefits from the recently available regression models developed at national and regional scales transforming the dNBR/NBR to CBI to consistently map burn severity in a transparent manner (Picotte et al 2021). Following Wiedinmyer et al (2011), we use modified algorithms and updated fuel emissions and consumption factors drawing on recently published field and remote sensing data to better capture emissions from more extreme fires. Larger, severe wildfires in California increasingly impact ecosystems, human health and life, communities, and infrastructure due to expanding wildland-urban interfaces, fuel accumulation from a century of fire exclusion policies, and increasing aridity, which increases the amount of biomass available to burn. (Abatzoglou and Williams 2016, Westerling 2016, Williams et al 2019, Goss et al 2020, Goodwin et al 2021. The comprehensive, long-term event and daily emissions records described here could be used to study health effects of wildfire smoke, either by combining them with transport modeling to model air quality and estimate exposures, or by incorporating them into statistical models predicting health impacts as a direct function of estimated emissions. These data will also facilitate analyses of changing emissions impacts on the carbon cycle over the last three decades. High resolution severity and emissions raster maps are generated for each fire event to support further spatial analysis. While the emissions calculated for California with WBSE are not a substitute for real-time daily emissions estimates, it is designed to extend the estimated emissions record back to 1984 with a finer spatial resolution and provide more up-to-date estimates on emissions factors reflecting information from California's recent extreme fires.

Methods
The WBSE provides estimates of 30 m resolution burn severity, and emissions of CO 2 , CO, CH 4 , nonmethane organic compounds (NMOC), SO 2 , NH 3 , NO, NO 2 , nitrogen oxides (NO x = NO + NO 2 ), PM 2.5 , OC, and BC. We implemented WBSE for California large wildfires on a per-fire event scale since 1984 and also a daily scale since 2002. The inventory implementation steps, input datasets, and output data are summarized in figure 1. Emissions of all species are calculated as a function of area burned, fuel loading, the fraction of vegetation burned based on burn severity, and an emissions factor specific to each vegetation type using the following equation modified from the FINN model : Where E i is the mass of emission species i emitted, A v,s and cr v,s are the area burned and fuel consumption rate of general vegetation class v (grass, shrub, forest <5500 ft, forest 5500-7500 ft, forest >7500 ft) under severity class s (low, moderate, and high severity), F v is the fuel loading of vegetation class v, and ef v,i is the emission factor of emission species i for vegetation v.
The FINN model was developed to estimate global emissions from open burning at 1 km resolution. It has been used in various modeling studies for modeling atmospheric chemistry and air quality (Jiang et al 2012, Val Martin et al 2013, Brey et al 2018. The primary differences of emission calculation methods between WBSE and FINN are that WBSE calculates emissions per fire event while FINN provides global emission estimates in real time. Other differences include (a) obtaining the area burned by fire severity class from a new CBI-based burned area product (based on fire perimeters published by MTBS and California Department of Forestry and Fire Protection (CAL FIRE)) while the Moderate Resolution Imaging Spectroradiometer (MODIS) Thermal Anomalies Product is used in FINN; (b) assigning vegetation type from the LANDFIRE existing vegetation type (EVT) product to general vegetation classes with an updated crosswalk from Hurteau et al (2014) while the MODIS land cover product is used in FINN; (c) calculating fuel loading and consumption based on burn severity for each vegetation class while fuel consumption is based on a function of tree cover and fuel loading for generic land cover classifications for the global regions in FINN; (d) Using region specific emissions factors for California and Western U.S. ecosystems to calculate emissions for the example implementation over California; and (e) assigning the day of burning based on MODIS and VIIRS data for each large wildfire burned since 2002. The inventory has a spatial resolution of 30 m for per-fire emission estimates and 500 m for daily emissions rather than the 1 km resolution in FINN. The area burned in low, moderate, and high burn severity is calculated from Landsat-derived indices (NBR and dNBR) using a hierarchy of regression models developed by Picotte et al (2021). We recognize that other fire severity metrics have been developed (e.g. RdNBR and RBR; Miller andThode 2007 &Parks et al 2014, respectively), but we opted to use NBR/dNBR because of the regionally specific models developed by Picotte et al (2021) that used NBR/dNBR. We do not have a quantitative assessment of the uncertainty; however, we assign a factor of 2 to the uncertainty due to all of the inputs to the model following Wiedinmyer et al (2011). A matrix table of 30 m pixel level PM 2.5 emission uncertainty for five general vegetation classes due to burn severity classification error is given in supplementary table 3.
Fire records for California from 1984 to 2019 were retrieved from MTBS (https://mtbs.gov/viewer/ index.html) via interactive viewer on 8 May 2021, resulting in a dataset with a total of 1623 wildfires. We also acquired fire perimeters for 74 large wildfires in 2020 from CAL FIRE (https://frap.fire.ca.gov/ frap-projects/fire-perimeters/) and calculated dNBR for each 2020 fire using the dNBR calculation tool with Google Earth Engine (GEE) (figure 2). This process first selects either initial assessment or extended assessment for each fire. The initial assessment utilizes Landsat images acquired immediately after a fire to capture first-order fire effects. The extended assessment uses images obtained during the growing season following the fire to identify delayed first-order effects and dominant second-order effects (Eidenshink et al 2007). We utilized LANDFIRE Biophysical Settings (BPS) to determine which assessment type to apply for each fire burned in 2020. After Picotte et al (2021), we used extended assessment if the majority of general vegetation groups within the fire perimeter are forests, while initial assessment is used when the majority of general vegetation groups are grassland/shrubland. By contrast, MTBS uses extended assessment for forest and shrubland types. We did not delineate grasslands into burn severity categories. Instead, we classified them as burned ('grass burn') because of difficulties in assessing vegetation change. Post-fire images for extended assessment were selected during the next peak of the green season (June-September) using the mean compositing approach suggested by Parks et al (2018). Composite post-fire images acquired immediately within two months after the fire containment dates were used for the initial assessment. Composite pre-fire images for extended and initial assessments were acquired with the matching periods from the preceding year. The dNBR images were produced by quantifying the spectral difference between composite pre-fire and post-fire Landsat scenes.
Burn severity for a given fire can have either a linear or non-linear relationship to dNBR/NBR (Picotte et al 2021). So to make an equivalent burn severity product across all fires, we calculated the unitless, continuous CBI variable from dNBR/NBR values using the linear and Sigmoid B regression models developed for the CONUS by Picotte et al (2021). We first applied the MTBS or CAL FIRE perimeter shapefile for each wildfire event to crop the BPS raster file. The BPS data are categorized into one of 12 vegetation groups. For CBI calculations, we further grouped 'Conifer' , 'Hardwood' , 'Hardwood-Conifer' , 'Shrubland' , and 'Riparian' into the 'forest/shrub' classification; 'Savanna' , 'Sparse' , and 'Grassland' into the 'grass' classification; and the rest 'Barren-Rock/Sand/Clay' , 'Open Water' , and 'Perennial Ice/Snow' into the 'unchanged' group. We applied the decision tree framework (figure 3) for each pixel categorized as forest/shrub to convert the dNBR/NBR value into a CBI value. CBI values less than or equal to 0 were arbitrarily assigned a value of 0.001 to be classified into the 'unburned' group.
Modeled CBI values larger than 3 were assigned the value 3 to be classified into the 'high severity' group. For the grass group, we set a CBI value of 3.5 to  indicate that this pixel was a grass burn. CBI values were then classified following thresholds modified based on Crotteau et al (2014) into six severity classes: unburned, low severity, moderate severity, high severity, grass burn, and non-processing area.
Fuel categories were assigned from LANDFIRE EVT products. EVT represents the current distribution of the terrestrial ecological systems classification (Rollins 2009), and the vegetation layers were updated every couple of years to account for landscape disturbances. The most appropriate version was applied for each fire event depending on the burn year-EVT data version LF 1.0.5, 1.2.0, 1.3.0, 1.4.0, 2.0.0 were applied to wildfires burned before the year of 2002, 2002-2010, 2011-2012, 2013-2014, and after 2014, respectively. For emissions calculations, EVT data were then categorized into five general vegetation categories: grass, shrub, forest under 5500 feet (1676 m), forest between 5500-7500 feet (1676-2286 m), and forest above 7500 feet (2286 m), updated for California ecosystems (supplementary table 1). Fuel consumption was determined following Hurteau et al (2014) assigning fuel loading and consumption values for each severity class for the five general vegetation categories based on the First Order Fire Effects Model v5 (Reinhardt et al 1997) (supplementary  table 2).
The LF EVT classes were assigned to generic vegetation classes of FINN so that emission factors could be applied. Emission factors for greenhouse gases, particulate matter, and reactive trace gases were updated with recent data for each general vegetation class using results from recent field campaigns and studies specific for California ecosystems and Western U.S. ecosystems (see table 1).
To assign the day of burning for individual pixels, NASA fire information for resource management system (FIRMS) active fire products from MODIS (Collection 6) within 750 m of the fire perimeter shapefiles supplied by MTBS or CAL FIRE were selected for interpolation to account for detections that might be outside the boundary due to detection radius. VIIRS 375 m data, when available since 2012, was added to complement MODIS data with improved performance to assign burn dates using the fire progression raster tool (figure 4). We filtered the MODIS/VIIRS detection points to the date range of interest and created a 500 m buffer around each point. Points were then converted to circle polygons to represent each point's detection extent properly. The average date was selected as the proper date in regions of overlapping buffers. We then calculated daily emissions and assigned them to the centroids of the aggregated daily progression polygons. Versions were also implemented using the earliest or latest date of detection points to assign the day of burning (results not reported here, but maps are available).

Results and discussion
Burn severity and emissions were calculated for 1697 large wildfire events in California during 1984-2020 at a 30 m resolution scale using the WBSE framework described above. Day of burning and daily

Burn severities
Over the 1984-2020 period, the average annual area burned (excluding obscured area) was 0.23 million ha with substantial interannual variability (figure 5). Obscured area includes data gaps associated with clouds, smoke, or the Landsat 7 ETM+ scan-line corrector failure (Key 2006), affecting 0.3% of the area mapped within historical burn perimeters. The maximum annual area burned was 1.6 million ha in 2020, about 220 times the area burned in the minimum year of 1991. Of the approximately 8.5 million ha mapped by this project, about 63% burned at moderate or high severity (figure 6). As MTBS has been previously established as a rough estimate of burn severity rather than a rigorous calculation for a specified need (e.g. emissions calculations), a direct comparison of the total area burned by severity category between MTBS and WBSE is helpful to contextualize the more rigorous approach applied here. MTBS analysts subjectively classify NBR and dNBR into six classes (unburned to low,  low, moderate, high, increased greenness, and nonprocessing area mask) . Here, we empirically classified modeled CBI values into four classes (unchanged, low, moderate, and high) and two categories of 'grass burn' and 'non-processing area' using an objective, standardized procedure. A pixel-level classification comparison (excluding nonprocessing area) indicates that ∼45% of the pixels in both datasets were classified as the same burn severity (cells highlighted in blue in table 2); however, WBSE burn severity was generally higher than MTBS (cells highlighted in pink in table 3). We found that 22% percent of pixels classified as low burn severity by MTBS were classified as moderate burn severity by WBSE. By contrast, less than 2% of pixels were classified as having higher burn severity by MTBS than WBSE (cells highlighted in green in table 3). The remaining pixels were mainly classified as low burn severity by WBSE but unburned to low burn severity by MTBS.

Emissions estimates
Summary statistics of trace gases and particulate species emitted by wildfires in California during 1984-2020 are shown in table 3. The annual average of CO 2 and PM 2.5 emitted are 18 Tg and 108 Gg, respectively, with a significant interannual variation. Emissions from wildfires varied significantly and depended on fire size, fire severity, and vegetation characteristics. The largest wildfire was the 2020 August Complex, which burned 417 898 ha in six counties in the Coast Range of Northern California and produced the maximum amount of all types of emissions. This single fire emitted approximately 2.5 times the annual average wildfire emissions for all of California for the past 37 years. We compared our emissions estimates of PM 2.5 to five global fire emissions inventories (GFEDv4s, FINNv1.5, GFASv1.2, QFEDv2.5r, FEERv1.0-G1.2) that are commonly used in atmospheric modeling simulations by extracting monthly and annual emissions for California from 2003 to 2020 using the FIRECAM online tool (Liu et al 2020). We note that GFEDv4s values assessed here for 2017-2020 emissions are preliminary. Annual emissions were also extracted from California Air Resource Board (CARB 2021) for 2000-2020 and WFEIS for 1984-2020. The estimates agree reasonably well with other emissions inventories, although they may vary by as much as a factor of 3 annually across inventories (figure 7).

Daily emissions
We calculated daily emissions for a total of 841 fires since 2002 when MODIS data were available. For 2002-2020, 6% percent of the original FIRMS fire points were removed annually due to low detection confidence. Seventy-four fire events were removed from the daily emissions database due to less than ten MODIS/VIIRS fire detections falling within the WBSE fire perimeter; the daily emissions from these events are considered uncertain given the low number of fire detections (after Parks 2014). These 74 fires were also small fire events, and their total emissions only accounted for less than 1% of the emissions during 2002-2020. The average daily emissions over the summer during 2002-2020 are higher compared to the average amount (0.47 Gg 8 ) of daily PM 2.5 emissions from all other sectors in California over 2002-2017 (figure 8). The extremes that peak around the middle of August and September are primarily contributed by the 2020 wildfires (figure 9).

Limitations and uncertainties
The results of this research should be considered with some limitations in mind, as uncertainties are associated with several aspects of the calculation process. Uncertainties related to the FINN model, including fire location, timing and area burned, vegetation identification, fuel loading, and fuel consumption, and emission factors, are described in detail by Wiedinmyer et al (2011). The modified FINN model used here reduces uncertainties related to fire numbers and sizes by using documented fire records rather than satellite detections. However, fire perimeters used in this calculation for wildfires during 1984-2019 were generated by on-screen interpretation and delineation of dNBR images by MTBS analysts. We used MTBS perimeters to better compare the burn severity difference between WBSE and MTBS. The results may change slightly when using other documented fire perimeter data; for example, we also calculated emissions since 1984 using perimeters from  CAL FIRE (not shown here but available). Uncertainty also arises from misidentification of land cover: (a) the burn severity classification is based on the LF BPS dataset, which characterizes vegetation systems that may have been dominant on the landscape before Euro-American settlement, but may not be representative of current vegetation; (b) fuel loading is calculated based on LF EVT, which represents the current distribution of the terrestrial ecological systems classification, but this dataset is infrequently updated and may not accurately represent the conditions at the time of burning. Fuel loading, consumption, and emission rates were constant values for general vegetation categories and severity classes, which may not reproduce the full heterogeneity of land cover. The day of burning was assigned with an average date in regions of overlapping buffers in this analysis, and the results change slightly if the earliest date or latest date is used as the proper date. In addition, around 5% of emissions were not assigned to a valid date since no fire points were detected by satellites. This may lead to a slight underestimate of daily emissions and thus the derived impact on air quality and public health. The open-source, open data framework presented here can be further refined as needed to address any of these limitations. Likewise, Picotte et al (2021) 's severity algorithms implemented here are designed to be robust to limitations in the available plot-level CBI data and can be extended to incorporate new plot-level data as they become available.

Conclusions
The WBSE framework implemented in R and GEE produces quick estimates of burn severity and emissions from wildfires at a 30 m spatial resolution and day of burning and daily emissions at 500 m spatial resolution. The burn severity estimates have advantages over MTBS products in that: (a) burn severity is categorized consistently based on field-validated equations, rather than subjective thresholds for each fire event; (b) an initial assessment can be estimated as soon as the desired fire perimeter is available; (c) they enable consistent annual estimates for wildfires burned in California during 1984-2020, which facilitates objectively quantifying changes over time and across fires. Further, daily emissions estimates can be calculated for periods when daily fire activity data are available (beginning in 2002). The daily emission estimates are currently being incorporated into chemical transport models to evaluate the impact of wildfire smoke on downwind air quality and population exposure. This work facilitates climate vulnerability assessments that include the effects of smoke on public health and provides estimates of greenhouse gas emissions during a fire. The publicly available WBSE fire severity history and emissions inventory offers a valuable, open-source resource for the fire research community, readily updated after each fire season, with a framework that allows for iterative improvements as more data become available. Future planned improvements include automated fire perimeter mapping of observed fires and incorporating interactive web tools to enable users to estimate fire emissions for simulated fire severity maps with usermodified climate and fuels management. For California's Fifth State Climate Assessment, the Pyregence Consortium is simulating individual wildfire events (location, date, size, and severity fraction) for future climate, land use and fuels management scenarios, mapping simulated severity at 30 m across the state, and using the framework presented here for estimating emissions from individual large fire simulations through 2100.

Data availability statement
Data used for this study can be obtained from our Git-Hub repository or downloaded following the R code. All the data generated from this study are available at the following URL: https://doi.org/10.6071/M3QX18 The data that support the findings of this study are openly available at the following URL/DOI: https://doi.org/10.6071/M3QX18.

Acknowledgments
This work was funded by the California Energy Commission EPC-18-026, National Oceanic and Atmospheric Administration Climate Program-NOAA NA170AR4310284, California Department of Insurance 18028CA-AM 1, Strategic Growth Council of California CCR20021, UC Lab Fees LFR-20-651032, and University of California Research Initiatives UCOP MRPI 20170261-01. M H acknowledges support from the California Department of Forestry and Fire Protection as part of the California Climate Investments Program, Grant #8GG14803. This research was supported in part by the USDA Forest Service, Rocky Mountain Research Station, Aldo Leopold Wilderness Research Institute. The findings and conclusions in this publication are those of the authors and should not be construed to represent any official USDA or U.S. Government determination or policy. Code and data is also hosted by Pyregence: GitHub.com/pyregence/ and https://data.pyregence.org/wg4/WBSE/.

Code availability
The code to implement our methods is available at the following GitHub repository: https://github.com/qxu6/WBSE.git.