Highly anomalous fire emissions from the 2019–2020 Australian bushfires

While it is widely recognized that extreme fires have been increasing under warming and drying climate, knowledge regarding the magnitude and intensity of extreme fires is very limited. Moreover, fire emissions reported by existing emissions inventories show large discrepancies due to different approaches and parameters. In this study, we analyzed the fire intensity and emissions magnitude of the 2019–2020 Australian bushfires using fire observations from multiple satellites. The results show that the bushfires were extreme in both their number and intensity, which were higher by a factor of 25 and 19, respectively, compared to the past two-decade seasonal mean. The 2019–2020 bushfires burned a total of 112.3 Tg biomass and released 178.6 ± 13.6 Tg CO2 (carbon dioxide), 1.71 ± 1.28 Tg PM2.5 (particulate matter with a diameter <2.5 μm), and 0.061 ± 0.04 Tg BC (black carbon) across eastern and southern Australia. The CO2 emissions are 35% of Australia’s greenhouse emissions from all sectors combined in 2020. Furthermore, the extreme fires in the most severe day and hour released 10% and 1.4% of the entire seasonal emissions, respectively. Our findings provide quantitative information for investigating the impacts of smoke emissions on air quality, ecosystem, and climate.


Introduction
Australia's 2019-2020 bushfire season greatly impacted communities, human health, the economy, and ecosystems (Filkov et al 2020, Bowman et al 2020a, Johnston et al 2021. The bushfires destroyed more than 9,000 buildings, resulting in at least 33 direct deaths (Filkov et al 2020). Meanwhile, smoke from fires blanketed eastern and southern states and substantially degraded air quality, causing exposure of more than 10 million people to hazardous air for months, thousands of hospitalizations, and ∼400 deaths indirectly (AIHW 2020, Johnston et al 2021). The bushfires disturbed ecosystems by burning flora, killing fauna, and devastating habitat .
The 2019-2020 bushfires exhibited record-breaking fire characteristics , Nolan et al 2020, Hirsch and Koren 2021. The bushfires scorched approximately 10 million ha across eastern and southern Australia (Bowman et al 2020a), where two eucalyptus forest dominated ecoregions (the temperate broadleaf and mixed forest ecoregion (TBMF) and the mediterranean forest, woodland and scrub ecoregion (MFWS)) are located. The burned area was larger by a factor of eight than the average for the past two decades, and likely topped the 170-year fire records since 1850 (Bowman et al 2020a). Meanwhile, the bushfires consumed a vast amount of fuels (Nolan et al 2020) and emitted massive amounts of greenhouse gases and aerosol emissions (Bowman et al 2020a, Hirsch andKoren 2021). Further, the bushfires generated ∼33 fire thunderstorms or pyrocumulonimbus clouds (pyroCbs) during the fire season (Kablick et al 2020, Bowman et al 2020b, which is very rare (Peterson et al 2018). Additionally, the bushfires injected smoke plumes into the stratosphere at an they show a discrepancy by a factor of 4-7 due to different assumptions and input data (Carter et al 2020, Pan et al 2020. For FRP-based emissions estimates, their accuracy depends largely on accurate estimation of fire radiative energy (FRE)-the temporal integration of FRP over a period, which requires high temporal resolution FRP to construct fire diurnal cycles (Freeborn et al 2009, Zhang et al 2012. All the currently available FRP-based emissions inventories (e.g., Quick Fire Emissions Dataset (QFED2.4), Fire Energetics and Emissions Research (FEER1.0), Global Fire Assimilation System (GFAS1.2), and Blended Global Biomass Burning Emissions (GBBEPx3.0)), however, use FRP, instead of FRE, from polar-orbiting sensors (e.g. MODIS) by simply assuming that FRP retrievals at a few satellite overpass times could represent diurnal fire activity (Darmenov and Silva 2015, Ichoku and Ellison 2014, Kaiser et al 2012. High spatiotemporal resolution FRP, which is not available from any single orbiting satellite, can be calculated by fusing high temporal resolution FRP from geostationary sensors and fine spatial resolution FRP from polar-orbiting sensors (Freeborn et al 2009, Zhang et al 2012. This motives this study to calculate fire emissions for the 2019-2020 bushfires by fusing the finest spatial and highest temporal resolution FRP observations that are available from the new generation satellite sensors.

Methods
We first measured the intensity of the 2019-2020 Australian bushfires by examining anomalies in fire intensity using two-decade 1 km FRP observations from the Moderate Resolution Imaging Spectroradiometer (MODIS) because MODIS onboard both Terra and Aqua satellites provides the longest consistent FRP records of scientific quality since 2002 (Giglio et al 2016). Then, we computed biomass consumption and fire emissions for the 2019-2020 bushfires by fusing the finest spatial resolution (375 m) FRP from the Visible Infrared Imaging Radiometer Suite (VIIRS) instrument and high temporal resolution (10 min) FRP from the Advanced Himawari Imager (AHI) at a spatial resolution of 2 km. Fire products and FRP from these sensors are introduced in detains in section 2.3.

Australia's ecoregions and forest map
Australia's terrestrial ecoregions were used to group fire activity. The terrestrial ecoregion map (figure 1(A), https://www.environment.gov.au/land/nrs/science/ibra/australias-ecoregions, last accessed on 2/4/2021) is developed by the Australian government department of Sustainability, Environment, Water, Populations, and Communities based on climate and vegetation. There are seven main ecoregions: (1) tropical broadleaf forest (TBLF), (2) tropical and subtropical grasslands, savannas and shrublands (TGSS), (3) deserts and xeric shrublands (DXES), (4) temperate broadleaf and mixed forest (TBMF), (5) mediterranean forest, woodland and scrub (MFWS), (6) temperate grasslands, savannas and shrublands (TEGS), and (7) montane grasslands and shrublands (MOGS). Because the MOGS ecoregion has a very small area and is located inside the TBMF ecoregion, for the sake of simplicity we merged MOGS ecoregion into the TBMF ecoregion. We adopted the latest Australian national forest map (a spatial resolution of 100 m) from Australia's State of the Forests Report 2018 (ABARES 2018). This map (figure 1(B)) was used to stratify biomass consumption and emissions from bushfires in forest and non-forest areas.

Three biomass-burning emissions inventories
Global fire emissions from GFASv1.2, GFED4s, and Fire INventory from NCAR (FINN v1.5 and v2.4) were collected to compare with our emissions estimates in this study. The GFASv1.2 estimates daily global fire emissions using MODIS FRP, biomass combustion coefficient, and emissions factors (Kaiser et al 2012). The GFED4s computes fire emissions using MODIS burned area data, modeled fuel loadings, combustion completeness, and emission factors (van der Werf et al 2017). After 2017, the GFED4s emissions are estimated using MODIS fire detections and empirical relationship between MODIS fire detections and historical GFED4s emissions estimates before 2017. Similarly, the FINN also estimates emissions using the burned area and fuel loading, where burned area in FINNv1.5 is derived from MODIS active fire detections and fuel loading is based on land cover types . According to the latest FINN Readme file (https://www.acom. ucar.edu/Data/fire/data/finn2/README_FINNv22_November2020.pdf), the FINN has been updated to version 2.4 recently by considering both 1 km MODIS and 375 VIIRS fire detections for estimating burned area. In this study, we obtained GFASv12 (https://apps.ecmwf.int/datasets/data/cams-gfas/; last access on 24 August 2021), GFED4s (https://www.globalfiredata.org/, last access on 24 August 2021), and both FINNv1.5 and v2.4 (http://bai.acom.ucar.edu/Data/fire/; last access on 24 August 2021). We also included CO 2 emissions reported recently by Shiraishi and Hirata (2021) for the 2019-2020 bushfires, which was calculated using burned area estimated from 1km MODIS fire detections and fuel loading from global above ground biomass (AGB) maps. This estimation was simply called Shiraishi and Hirata 2021 hereafter. The emissions inventories FEER and QFED, which are based on coefficients derived from or tuned by aerosol optical depth (AOD), were not included in this study because they are usually several times larger than other inventories (Carter et al 2020, Pan et al 2020.

Satellite active fire data
We used a two-decade record of MODIS active fire data, which is the longest satellite-derived fire record of scientific quality, to examine anomalies in the fire intensity of Australia's 2019-2020 bushfires. The MODIS sensor on Terra (launched in December 1999) and Aqua (launched in May 2002) satellites detects actively burning fires using the 4-and 11-μm bands approximately four times a day at 01:30, 10:30, 13:30, and 22:30 local time (Giglio et al 2016). The latest collection 6 MODIS level-2 active fire products (abbreviated MOD14 and MYD14 for Terra and Aqua satellites, respectively) at 1-km spatial resolution were obtained from NASA's Level-1 and Atmosphere Archive and Distribution System (LAADAS, https://ladsweb.modaps.eosdis.nasa. gov/) for the period 2002-2020. It provides for each fire detection the detection time, geolocation, detection confidence, and FRP (Megawatt, MW) that is calculated using the radiance of the fire pixel and its ambient nonfire pixels in the 4-μm band (Giglio et al 2016). The 4-μm band has a saturation temperature of ∼500 K that is reached only for extremely powerful fires (Giglio et al 2003).
We fused VIIRS and AHI active fire data to estimate emissions in the 2019-2020 fire season. VIIRS onboard both Suomi National Polar-Orbiting Partnership (Suomi NPP, launched in October 2011) and first satellite of the Joint Polar Satellite System (JPSS-01, subsequently named NOAA-20, launched in November 2017) crosses the equator at approximately 01:30 and 13:30 local time but NOAA-20 VIIRS observes fires over the same area ∼50 min earlier. The VIIRS active fire products at 375 m pixels from both Suomi NPP and NOAA-20 were obtained from NOAA's Center for Satellite Applications and Research (ftp://ftp.star.nesdis.noaa.gov/) for the period from November 2019-January 2020. These fire products provide for each fire detection the detection time, geolocation, confidence, persistent anomaly flag, and FRP (MW). The 375m VIIRS FRP is retrieved using the co-located radiance of the fire pixel and its background non-fire pixels in the 4-μm 750 m band that has a very high saturation temperature of 600 K, which leads to far fewer saturated fire detections than MODIS . Moreover, the 375 m VIIRS FRP has the finest spatial resolution among all existing satellite-based fire products. This resolution allows for the detection of many more very small and cool fires than the MODIS fire data at 1 km resolution , Li et al 2020a. For example, in the 2019-2020 fire season, a comparison of daily sum FRP from three sensors (Aqua MODIS (1 km), NPP VIIRS (375 m), and NOAA-20 VIIRS (375 m)) with close overpass times shows that the VIIRS sum FRP is on average ∼50% higher than MODIS sum FRP (figure S1 (available online at stacks.iop.org/ERC/3/105005/ mmedia)). This is the reason why 375 m VIIRS FRP was chosen to estimate fire emissions. The 375m VIIRS FRP would also be a better choice for examining anomalies of fire intensity if the data has a temporal coverage as long as the two-decade MODIS FRP. However, as a follow-on sensor of MODIS, VIIRS provides FRP for less than 10 years, which is insufficient for the purpose of anomaly examination. As the MODIS fire detection algorithm performs consistently for the whole two-decade observation period, the comparison of seasonal mean FRP during past two decades could reliably reveal interannual FRP anomalies (see section 2.3).
The Himawari-8 AHI observes fires at 2 km pixels across East and Southeast Asia and Australia at a temporal resolution of 10 min since 2015 (Bessho et al 2016). The official level-2 active fire product from Japan Aerospace Exploration Agency (JAXA 2020) (available at JAXA's P-tree system FTP, ftp://ftp.ptree.jaxa.jp) detects fires using the 4-and 11-μm AHI bands and was updated from the Beta version to version 1.0 in October 2020 (JAXA 2020). We used the version 1.0 to screen the cloud edge and coastline related false alarms in the Beta version, and calculated FRP using the radiance of a fire pixel and at its ambient non-fire background in the 4-μm band in the Beta version (see details in subsection 2.3 'Examination of the Anomalies in Fire Number and Intensity'). We did not use the FRP provided in version1.0, which was calculated using a bi-spectral method with radiance in the 4-and 2.3-μm bands (JAXA 2020). This is because it is very difficult to characterize non-fire background using the 2.3-μm band due to its high sensitivity to non-fire hot surfaces (e.g., newly burned areas) . Inaccurate characterization of non-fire background results in unreliable fire temperature and fire size, which are the sub-pixel parameters required for calculating FRP in the bi-spectral method .

Examination of the Anomalies in fire number and intensity
The total number of MODIS fire detections and summed FRP separately represent the extent of actively burning fires and fire-emitted energy at observation times. Thus, these two parameters during the season of November 2019-January 2020 were first compared with those in the past 17 seasons from 2002-2019. Because MODIS pixel size enlarges as scan angle increases from nadir to scan edge, which is known as the bow-tie effect that causes sizeable overlapping between adjacent scans in off-nadir view angles (Wolfe et al 1998), adjacent scans can observe the same fire twice (Freeborn et al 2014, Li et al 2018b. Therefore, we first removed the duplicate fire detections between consecutive MODIS granules following the method developed by Li et al (2018b). Then, for a given ecoregion, the number of fire detections and FRP were summed up for each season from November of one year to January of the next year during 2002-2020. Subsequently, the ecoregion-specific anomalies of the total number of fire detections and summed FRP for each season were calculated by subtracting the corresponding 2002-2019 means separately. Finally, the anomalies of the 2019-2020 fire season were compared with the past seasons.
Considering FRP as a proxy of fire intensity, we further investigated the extreme characteristics of fire intensity during the 2019-2020 fire season by comparing with the past seasons. MODIS fire detections with an FRP>1,500 MW are very rare globally, mostly observed in forest wildfires across North American boreal regions, the western United States, and eastern and southwestern Australia (Ichoku et al 2008). In eastern and southwestern Australia, fire regimes are typically characterized by intense fires (Archibald et al 2013). To define extreme FRP across Australia, we first examined the frequency distribution of the nearly two-decade FRP retrievals of fire pixels with a 4-μm brightness temperature (BT)500 K that is the approximate saturation temperature of the MODIS main fire detection band at 4 μm (channel 21). Very strong thermal signals are generally emitted from very powerful fires burning in a saturated fire pixel (Giglio et al 2003, Ichoku et al 2008. Because the FRP frequency distribution with a 4-μm BT500 K shows that the majority of the fire detections have an FRP1,600 MW (FRP value at the peak of the frequency distribution, see figure S2), we defined extreme fires using an FRP threshold of 1,600 MW. Then, the extreme FRP was further grouped into three categories: (1) 1,600-2,500 MW, (2) 2,500-5,000 MW, and (3)>5,000 MW. Finally, the number of fire detections with extreme FRP for the season from November-January was summed up for three categories in each ecoregion during 2002-2020. Additionally, the number of saturated fire detections (BT500 K) without valid FRP retrievals was also counted for each season and ecoregion. This is because these detections should theoretically have a valid FRP of at least 1,600 MW according to the frequency distribution.
Because the fire activity during November 2019-January 2020 mainly occurred in the eucalyptus forest dominated TBMF and MFWS ecoregions (figures 1 and 3), we only focused on these two ecoregions in comparing extreme FRP and estimating biomass consumption and fire emissions as described in the following subsection.

Calculation of biomass consumption and fire emissions
Biomass consumption and fire emissions were estimated using diurnal FRP that was fused from high temporal AHI and high spatial VIIRS observations.

Calibrating AHI FRP against and Fusing with VIIRS FRP
The AHI FRP was calculated first using parameters provided in the Beta version of the AHI fire product. This is because the FRP values in the version 1.0 AHI fire product are not reliable as described in section 'Satellite active fire data'. Specifically, the Beta version AHI fire data were first refined by screening false alarms using version 1.0 AHI fire data and BT-based tests. By comparing AHI fire data with MODIS and VIIRS fire detections across Australia during the study period, we found that false alarms in the Beta version fire data were mainly related to cloud edges and water edges (i.e., coastlines), which were removed in the version 1.0. Therefore, we excluded fire detections that appeared only in the Beta version by referencing version 1.0. After removing the cloud-and water-related false alarms, we further removed fire detections if they failed to meet two BT conditions (BT 4 280 K and BT 4 -BT 11 1 K) by following the BT-based tests in the Fire Thermal Anomaly (FTA) algorithm that was designed for the Meteosat SEVIRI (Spinning Enhanced Visible and InfraRed Imager) sensor (Roberts and Wooster 2008) and also applied successfully to AHI (Xu et al 2017). Note that BT 4 and BT 11 denote brightness temperature in the 4-and 11-μm bands. As a result, a total of 24% fire detections were considered as false alarms and excluded.
FRP at a 10-min interval was then calculated for the refined Beta version AHI fire detections using the Mid-Infrared (MIR)-radiance method (Wooster et al 2005): are the radiance of an AHI fire pixel and the background non-fire pixels in the 4-μm band, respectively; (Note that for a given fire pixel, the associated background pixels were determined by searching neighboring non-fire pixels in a window, where the window size is determined usually by balancing searching sufficient background pixels with algorithm run time). A is AHI pixel area; s is the Stefan-Boltzmann constant (5.6704×10 −8 Wm −2 K −4 ); and a is a sensor-specific constant (a =3.0×10 −9 W m −2 Sr −1 μm −1 K −4 ). The 4-μm band radiance was calculated from BT of the same band using the AHI calibration table (JMA 2015). The JAXA fire detection algorithm calculates background BT by averaging BT values of the non-fire background pixels in a window centered at each fire pixel (JAXA 2020). The background BT varies between 280 K and 334 K for AHI fire detections across the study area. Note that large difference in land cover type between a fire pixel and its non-fire background pixels could cause bias in FRP estimates, which could happen in all active fire products (e.g., MODIS and VIIRS active fire products) that calculate FRP using the MIR method (Schroeder et al 2010).
The AHI FRP values at 2-km pixels were further calibrated against the 375-m VIIRS FRP to mitigate the underestimation of 2-km AHI FRP due to missing observation of small and/or cool fires. To do this, VIIRS and AHI FRP retrievals were first aggregated to 0.25°grids separately. A grid size of 0.25°was chosen to match the grid size of commonly used large-scale global models, such as U.S. NWS (National weather Service) global forecasting system, and air quality forecast models (Eastham and Jacob 2017) (Note that emissions can also be generated at higher resolution for other applications). The gridded AHI FRP was then calibrated using the gridded and contemporaneous VIIRS FRP, where the contemporaneous FRP was referred to as the FRP observed by both VIIRS and AHI within 5 min. Thus, the calculation equation is: where FRP t A and FRP t A are the gridded AHI FRP before and after adjustment at observation time t, respectively; and r t is a calibration factor (unitless).
For each of grids where contemporaneous AHI and VIIRS FRP were available, the calibration factor r t was calculated for each AHI observation time t (with FRP>0) by temporally interpolating ratios of the gridded FRP differential (VIIRS FRP -AHI FRP) to AHI FRP linearly at two consecutive VIIRS observation times in a day (one before t and the other after t). The calibration assumes that the underestimation of AHI FRP is proportional to the magnitude of AHI FRP. The two consecutive VIIRS observation times were determined based on the times of valid VIIRS observations, two for Suomi NPP and two for NOAA-20. If only one VIIRS observation time was found for the periods before the earliest and after the latest VIIRS observation time, the ratio at the earliest or the latest observation time was used. Thus, the calibration factor r t varies spatially and temporally. Comparisons between contemporaneous gridded VIIRS FRP and AHI FRP (before and after applying calibration) are illustrated in figure S3.
VIIRS FRP and the calibrated AHI FRP were fused as: here FRP t f is fused AHI-VIIRS FRP; FRP t A is the calibrated AHI FRP as in equation (2); FRP t V is the gridded VIIRS FRP at observation time t; and w 1 and w 2 are fusing weights. When VIIRS FRP is available, w 0 1 = and w 1; 2 = and if only AHI FRP is available, w 1 1 = and w 0. 2 = 2.5.2. Calculating biomass consumption and emissions FRP diurnal cycle was reconstructed subsequently for each grid by filling temporal gaps between valid fused FRP values. Although AHI observes fires every 10 min, fire-observing gaps appear occasionally due to obscuration by clouds, smoke plumes, and forest canopies. FRP values in temporal gaps were filled by interpolating the adjacent valid fused FRP. For grids where only VIIRS observed fires, which were mostly too small/cool to be detected by AHI, at a single observation time, we assumed conservatively that fires had burned for one hour before and after the observation time separately. As a result, hourly or daily fire radiative energy (FRE, Megajoule (MJ)) was calculated by integrating the reconstructed FRP diurnal cycle: where FRP t f is the fused FRP from equation (3) for a given grid; and t1 and t2 are the first and last observation times of the reconstructed diurnal cycle, respectively.
Finally, biomass consumption and CO 2 , PM2.5, and BC emissions were calculated as: here E is the hourly or daily total mass of an emission species (kg); DM is dry biomass consumed (kg); F is the emission factor of the emission species (kg kg −1 ); FRE is from equation (4); and B is the FRE biomass combustion coefficient (B=0.368 kg MJ −1 ) (Wooster et al 2005). Because emissions factor (especially for aerosols) reported for Australia's temperate forest is very limited, this study used emissions factors (CO 2 , PM2.5, and BC) that were compiled by Andreae (Andreae 2019) based on recent studies of global fire emissions factors and the widely used emissions factors compiled by Akagi et al (2011). Uncertainties in emission factors, which are measured typically by standard deviation, are propagated to the associated emission estimates (equation (6)). Thus, standard deviations of emission factors were applied to estimate uncertainties in emission estimates. To quantify emissions from forest fires, biomass consumption and fire emissions were further classified into forest and nonforest based on Australia's national forest map (ABARES 2018).
The emissions estimates were compared with four different emissions estimates described in the section 2.2. To minimize the impact from different emissions factors, we used the CO 2 emission factors close to three global inventories (e.g.,<5% for temperate forest). Red and blue bars show positive and negative anomalies, respectively. Each fire season covers three months from November to January (next year). Temperate eucalyptus forests are mainly distributed in the TBMF and MFWS ecoregions. Note that fire detection number and summed FRP were not shown for the tropical broadleaf forest ecoregion and the temperate grasslands, savannas and shrublands ecoregion because fire activities were relatively limited during the study period.

Anomalies in Fire Number and Intensity
The 2019-2020 bushfires across the TBMF (temperate broadleaf and mixed forest) and MFWS (Mediterranean forest, woodland, and scrub) ecoregions burned in a way that has never been observed from the satellite record of the past 17 fire seasons from 2002-2019 (simply called two decades hereafter), in terms of the number of fires and fire intensity (figure 2). In the TBMF ecoregion, fire activity was very limited in the past two decades, with a mean MODIS fire detection number of 5,200 and a mean summed FRP of 0.4 Terawatt (TW). However, the 2019-2020 fire season had a total of 89,000 fire detections and a summed FRP of 7.9 TW, which were higher by a factor of 16 and 19 than the 2002-2019 means, respectively. The only other extreme fire season in Australia in the past two decades occurred during the 2002-2003 fire season. The 2019-2020 fires were two times more intense than fires in the 2002-2003 season with respect to both the number of fires and FRP. Moreover, the day of 4 January 2020 had the highest number of daily fire detections (6800) and summed FRP (1.03 TW) (figure 3(A)), which are 30% and 160% higher than the means for the past two decades. This most severe day accounted for ∼8% and ∼13% of fire observations and summed FRP in the 2019-2020 season, respectively. Additionally, fire observations on 4 January doubled relative to the prior day, and summed FRP increased by 300% correspondingly ( figure 3(A)). Similarly, in the MFWS ecoregion, the 2019-2020 fire season had a total of 14,300 fire detections and summed FRP of 2.6 TW, which are higher by a factor of 4.5 and 6.4 than the 2002-2019 means (2,600 and 0.35 TW), respectively. The highest number of daily fire detections (1060) and summed FRP (0.3 TW) were observed on 2 and 8 January 2020, respectively, and were approximately 40% and 86% of the means for the past two decades ( figure 3(B)).
The FRP during the 2019-2020 fire season could be much higher if including the FRP from MODIS saturated fire observations. MODIS fire observations can be saturated in very intense burnings, resulting in no valid FRP measurements. There were 14% of MODIS saturated fire detections during the three days (30 and 31 December 2019, and 4 January 2020) in the TBMF ecoregion, which led to the underestimation of summed FRP. In the MFWS ecoregion, likewise, FRP was underestimated on 19 December 2019 and 3 January 2020 because 23% of MODIS fire detections were saturated on these days.
In the TGSS (tropical and subtropical grasslands, savannas, and shrublands) and DXES (deserts and xeric shrublands) ecoregions, however, fire activity in the 2019-2020 season was lower relative to the past seasons (figure 2). The total number of fire detections for the 2019-2020 season was ∼5% and ∼53% smaller than the 2002-2019 means in these two ecoregions, respectively. The summed FRP was also ∼10% and 53% lower than the 2002-2019 means correspondingly.
Taking into account all ecoregions across Australia, the 2019-2020 season had 130% and 180% more fire detections and summed FRP than the means of the past two decades, respectively (figure 2). Compared to past individual fire seasons, the anomaly of the 2019-2020 season was slightly higher than the 2002-2003 season but The number of extreme fire observations in the 2019-2020 bushfires was also record breaking across the TBMF and MFWS ecoregions (figure 4). Extreme fire observations refer to fire detections with an FRP value greater than 1,600 MW, where the FRP threshold was determined based on the FRP frequency distribution of historical fire records across Australia (see details in Methods). In the TBMF ecoregion, the number of extreme fire observations during the 2019-2020 season was 50% more than the total number of those for the past two decades and larger by a factor of 19 than the seasonal average ( figure 4(A)). Their difference of extreme fire observations was as high as a factor of 25 if including the saturated fire detections that theoretically had an FRP1,600 MW. Compared to the second severest 2002-2003 fire season, the number of extreme fire observations for the 2019-2020 season was higher by a factor of 2.4. Similarly, in the MFWS ecoregion, the number of extreme fire observations during the 2019-2020 season was higher by a factor of 10 than the seasonal average for the past two decades and a factor of 2.5 than the second severest 2002-2003 fire season (figure 4(B)).  usually very limited), and eucalyptus woodland and mallee forest in MFWS. The bushfires in eucalyptus forest across these two ecoregions consumed a total of 103.8 Tg dry biomass that accounted for 92% of all biomass consumption and released 92% of CO 2 , 96% of PM2.5, and 93% of BC emissions accordingly.

Biomass consumption and fire emissions
Fire emissions during the 2019-2020 season also revealed dramatic temporal patterns in monthly, weekly, daily, and hourly variations (figure 6). The bushfires in November 2019, December 2019, and January 2020 separately contributed to 26%, 45%, and 29% of seasonal biomass consumption and emissions. During the week from 29 December 2019 to 4 January 2020, fires consumed a total of ∼35 Tg dry biomass and released ∼55±4.3 Tg CO 2 , ∼0.54±0.4 Tg PM2.5, and 0.019±0.012 Tg BC, which accounted for ∼31% of the total biomass consumption and emissions released during the 2019-2020 season. The three highest daily amounts of biomass consumed and fire emissions during the season occurred on 30 and 31 December 2019 and 4 January 2020, which accounted for 22% of seasonal fire emissions. The most severe fire was on 4 January with a biomass consumption of ∼11 Tg and emissions of ∼17±1.4 Tg CO 2 , ∼0.2±0.15 Tg PM2.5, and ∼0.006±0.004 Tg BC, which accounted for ∼10% of seasonal biomass consumption and fire emissions. The diurnal variation of fire emissions revealed that most of fire smoke was emitted during daytime and fire emissions peaked generally in the early afternoon (1:00 pm-3:00 pm) during the week (figure 7). For instance, fires during the peaking hour (2:00 pm) on 4 January 2020 accounted for 14.5% of the daily total biomass consumption and emissions or 1.4% of the total emissions of the entire season. Nevertheless, nighttime fires were equally remarkable; half of daily emissions on two severe days in December (30 and 31) were emitted during night.   Figure 8 shows the comparison of monthly total CO 2 emissions across the TBMF and MFWS ecoregions. In the three months from November 2019 to January 2020, the AHI-VIIRS FRP based monthly CO 2 (from this study) was 44.4±3.4 Tg (November), 82.2±6.3 Tg (December), and 52±3.9 Tg (January), with a total of 178.6±13.6 Tg for the whole period. In comparison with other inventories, the total mass of the AHI-VIIRS FRP based CO 2 was approximately 20% higher than FINNv1.5 (150 Tg) but differed by a factor of 2-3 from    FRP, the best and longest satellite-observed fire data, quantified the extreme anomaly of the 2019-2020 bushfires. The number of fire observations and summed FRP across the TBMF differed from the means for the past two decades by one order of magnitude, and were higher by a factor of two than the second-largest fire season (2002)(2003) on the record, which indicates the extremely abnormal extent of actively burning fires and emitted radiative energy. The extreme fire intensity in a single day occurred on 4 January 2020, where the number of fire observations and summed FRP were 30% and 160% higher than the means for the past two decades, respectively. Furthermore, the large number of extremely intense fire observations also demonstrated the extreme-intensity characteristics of the 2019-2020 bushfires ( figure 4). The number in the TBMF ecoregion was one order of magnitude larger in the 2019-2020 season than in the past two decades. Moreover, the saturated fire observations that indicate the much extremer intensity of fires were rarely seen in the past seasons but accounted for nearly 30% of all extreme FRP observations in the 2019-2020 season across the TBMF ecoregion. These extreme fires strongly imply substantially powerful disturbances to ecosystems .

Comparison with other fire emissions inventories
In the TGSS and DXES ecoregions, however, strong positive anomalies of fire number and intensity were not observed in the 2019-2020 fire season despite the widespread sever fire weather. The fire activity in these two ecoregions, especially DXES, was lower than the mean of past seasons. The weakened fire activity is likely explained by decreased fuel availability. Across these two regions, fuels are mainly savanna tall grasses (in TGSS) and sparse shrubs and grasses (in DXES) that usually accumulate during wet season before burning in dry season (Pausas and Ribeiro 2013). The year-long drought across the Australia continent likely limited vegetation growth before the 2019-2020 fire season (ABM 2019), which substantially reduced available fuels to burn.
The biomass consumption as well as greenhouse gases and aerosol emissions from the 2019-2020 bushfires are extremely high as observed by VIIRS and AHI. Because biomass consumption and fire emissions are a linear function of fire radiative energy (Freeborn et al 2011), their anomalies of extremities should be very similar as the summed FRP in the 2019-2020 fire season that was larger by an order of magnitude than the 2002-2019 mean (figure 2). The magnitude of fire emissions released from the 2019-2020 bushfires is further demonstrated by comparing it with Australia's industrial greenhouse gas emissions and other important fire emissions sources globally. First, the total amount of CO 2 emissions from the 2019-2020 bushfires in TBMF and MFWS is as high as 35±3% of Australia's annual CO 2 -equivalent emissions from all sectors for the year from July 2019 to June 2020 (513.4 Tg), and even higher than that of any industrial sectors (e.g., energy-related electricity and transport emissions) (ADISER 2020). Carbon emissions emitted from natural wildfires are not considered as a net source of carbon usually as post-fire vegetation regrowth sequesters carbon back from the atmosphere over time (Landry and Matthews 2017). Nevertheless, the amount of fire-released CO 2 emissions is approximately a third of Australia's national CO 2 -equivalent emissions in 2005, which is considered as the baseline greenhouse emissions for Australia in the Paris Agreement with the goal of reducing greenhouse emissions (ADISER 2020). Clearly, fire emissions are one important source of atmospheric greenhouse gas emissions in short term and will likely increase because of the projected increasing fire activity in warming climate scenarios ( Landry andMatthews 2017, Abatzoglou et al 2018). Second, Australia's 2019-2020 bushfire-released CO 2 emissions are very close to the CO 2 emissions (180 Tg) from wildfires across the western United States in the worst 2020 fire season (van der Werf et al 2017). Moreover, the 2019-2020 CO 2 emissions are ∼37% of the CO 2 emissions from deforestation fires across the Amazon rainforest in 2019, which has been considered as the worst fire season in the past decade (Escobar 2019).
Moreover, biomass burning was much severe on hourly, daily, and weekly time scales during the fire season. Generally, the temporal pattern of the daily fire emissions in the 2019-2020 fire season matches that of the satellite-based absorbing aerosol index (AAI) that measures the concentration of absorbing aerosols in fire smoke (Khaykin et al 2020), with the highest peaks during the week from 30 December 2019 to 4 January 2020 (figure 6). During this week, the bushfires accounted for more than 30% of the total biomass consumption and emissions during the entire season, which likely explains the fast elevated ground-measured PM2.5 (50-100 μg m −3 averaged by population, Borchers Arriagada et al 2020) and the rapid increase of fire-smoke-related hospitalizations across southeastern Australia in the following week (Johnston et al 2021) when emissions were transported to populated areas. The largest daily emission, which occurred on 4 January 2020, was ∼17 Tg CO 2 , ∼0.2 Tg PM2.5, and ∼0.006 Tg BC that explained ∼10% of the total emissions for the entire season, with the most extremely intense fire observations. On this day, moreover, the hourly fire emissions during the most severe hour accounted for 1.4% of the seasonal total emissions. The substantial amounts of energy and emissions released from these extreme fires during the week with the largest daily and hourly fire emissions are likely the major forcing agents that caused the occurrence of 18 observed pyroCbs (Kablick et al 2020) as well as the record-breaking injection height and spread distance of smoke plumes during the same period (

Implications of fire emissions for smoke monitoring
The extremely large amount of fire emissions has strong implications for monitoring smoke aerosols. The amount of aerosols is inferred typically from satellite-based and/or ground-based aerosol optical depth (AOD) that measures the degree to which aerosols cause attenuation of sunlight and implies air pollution, with a normal Figure 9. Aerosol optical depth (AOD) in southeastern Australia on 4 January 2020. The right panel shows spatial details of AOD over the region highlighted in pink in the left panel. In TBMF (temperate broadleaf and mixed forest) ecoregion, AOD is on average 12.8 (much greater than 5.0 -the maximum AOD that most satellites could retrieve) in 18% of actively burning fire areas in TBMF on 4 January 2020. Note that AOD was calculated from PM2.5 emission flux at a spatial resolution of 375m by following the method developed by Hoff and Christopher (Hoff and Christopher 2009), with a typical extinction coefficient of 3.14 m 2 g −1 and an integration time of 300 seconds. PM2.5 emission flux was computed using the 375m VIIRS FRP. value range of 0-1.0 (up to 5.0 for satellites and 7.0 for AERONET) (Levy et al 2013). Satellite sensors that retrieve AOD using visible wavelengths (e.g., MODIS, VIIRS, and AERONET) usually fail to detect AOD in extreme conditions when highly concentrated aerosols (e.g., thick fire smoke plumes and pyroCbs) absorb nearly all the incoming sunlight at visible wavelengths and scatter little back to sensors (Levy et al 2013). Until the 2019-2020 extreme bushfires, one could not imagine that emissions could be of this scale resulting in such high AODs. The PM2.5 emissions from the 2019-2020 bushfires enhanced our knowledge of extremely high AOD. Based on the PM2.5 emissions in this study and the AOD calculation method employed by many researchers (i.e., equation 10 in Hoff and Christopher (Hoff and Christopher 2009), with a typical extinction coefficient of 3.14 m 2 g −1 for smoke aerosol (Dobbins et al 1994) and an integration time of 300 s), we revealed that AOD was on average 12.8 (much greater than 5.0) in 18% of actively burning fire areas in TBMF on the most severe day 4 January 2020 ( figure 9). This is far beyond the observation capabilities of visible wavelength-based instruments. Meanwhile, the ultra-violet (UV)-based instrument (i.e., Ozone Monitoring Instrument or OMI on Aura) reported a maximum AOD of 13.0 (personal communication with Dr Hiren T. Jethva). Therefore, this study suggests that instruments retrieving AOD in different wavelength regions are needed to understand the physical and optical properties of smoke plumes in extreme fires like the 2019-2020 bushfires.

Uncertainty in AHI-VIIRS FRP-based emissions estimates
The uncertainty of the FRP-based emissions estimates in this study is associated with potential uncertainties of FRE estimate, FBCC value, and emissions factors. The FRE estimate can be affected by the omission errors of AHI fire detections due to two factors: (1) missing observations of small and/or cool fires relative to VIIRS, and (2) obscuration of clouds and thick smoke. The effect of the first factor was mitigated by calibrating AHI FRP against VIIRS (equation (2), illustrated in figure S3). The FRP observation gaps caused by clouds and smokes were interpolated using the closest valid FRP observations in reconstructing FRP diurnal cycle. The great agreement between the fused FRP diurnal cycles and VIIRS FRP suggests that our FRE estimates are reliable with high accuracy ( figure 10(a)). This is further discussed in comparison with GFASv1.2 emissions (see section 4.4).
Uncertainty in emissions estimation could also come from FRE biomass combustion coefficient (FBCC) but it is likely limited. FBCC was initially derived as a value of 0.368 kg MJ −1 in a controlled field fire experiment that was conducted to explore the relationship between rate of biomass consumption and fire-emitted instantaneous radiative energy (Wooster et al 2005). Several recent studies further investigated FBCC in landscape wildfires (Konovalov et  . Therefore, we believe that it is reasonable to use the experiment-based FBCC (0.368 kg MJ −1 ) for calculating emissions in Australia temperate forest fires and that FBCC-caused uncertainty in emissions estimation should be relatively limited.
The main source of uncertainty in our emissions estimates is likely from uncertainties in emission factors. This study used emission factors compiled recently by Andreae (Andreae 2019) based on emissions factors for global fires, including Australia fires (Guérette et al 2018). Given that the emission factor (e.g., g of PM2.5 per kg of biomass consumption) of a specific emissions species is derived usually as an average of reported emission factors from different studies, standard deviation (σ) is used as a measure of the uncertainty. In temperate forest where the 2019-2020 bushfires burned mostly, 1σ explains 8%, 78%, and 65% emission factors of CO 2 (1570±130 g kg −1 ), PM2.5 (18.5±14.4 g kg −1 ) and BC (0.55±0.36 g kg −1 ), respectively. Although Andreae's emission factors are average values, the CO 2 emission factor is very close to the airborne-based and ground-based CO 2 emission factor for temperate forest across southeastern Australia (Guérette et al 2018). As a result, the total emission mass varies from 165.4-192.6 Tg for CO 2 , 0.43-3.03 Tg for PM2.5, and 0.021-0.101 Tg for BC in the 2019-2020 fire season. It is clear that PM2.5 and BC have much larger uncertainties than CO 2 , which is attributed partly to the fact that measuring the emissions factors of aerosol emissions (e.g., PM2.5 and BC) is much more complex than gaseous emissions (e.g., CO 2 ) (Andreae 2019).

Sources of discrepancies among various emissions estimation
Large discrepancies in CO 2 emissions estimates for the 2019-2020 bushfires among existing emissions inventories indicate challenges in emissions estimation. First, the GFASv1.2 calculates emissions in the same way (equations (5) and (6)) as this study but uses 1 km MODIS FRP (Kaiser et al 2012). Because similar CO 2 emissions factors were used (difference <5%), the discrepancy between GFASv1.2 and this study is attributed to the FRE estimates and FBCC values. During the study period, the MODIS FRP based FRE (372 PJ) in GFASv1.2 was 22% higher than our FRE (305 PJ) calculated from the fused diurnal AHI-VIIRS FRP even though MODIS FRP ( Figure S1) is ∼50% lower than VIIRS FRP. Further, GFASv1.2 uses an FBCC of 0.49 kg MJ −1 for temperate forests, which is 33% higher than the FBCC (0.368 kg MJ −1 ) used in this study. These two sources can result in a CO 2 emission difference by a factor of 1.62, which is very close to the difference presented in figure 8. It should be noted that accurate FRE estimation requires FRP diurnal cycles. However, the GFASv1.2 calculates FRE using daily mean FRP averaged from the daily up-to-four MODIS FRP, which is overestimated generally (figure 10). Specifically, GFASv1.2 overestimated FRP in the morning and evening but underestimated in the afternoon relative to the AHI-VIIRS FRP ( figure 10(A)). Although some compensation from over and under estimations, GFASv1.2 highly overestimated FRP in the morning and evening, which lead to the overestimation of daily FRE in most days ( figure 10(b)). Moreover, FBCC values used in GFASv1.2 were derived by relating MODIS FRP to the biomass consumption from an earlier version of GFED (GFEDv3), which is subject to the accuracy of GFEDv3.
Second, the CO 2 estimates from the conventional method (burned area, fuel loading, combustion completeness and emission factors) (i.e., FINNv2.4 and Shiraishi and Hirata 2021) are overall larger than the FRP-based CO 2 estimates (i.e., this study and GFASv1.2), but show large discrepancies among themselves ( figure 8). Diagnosis of the difference between the FRP-based and the conventional method estimates is difficult because these two methods use different input parameters. The two methods are connected when FRE is considered as a function of burned area, fuel loadings, and combustion completeness, which highlights the advantage of the FRP-based method as it bypasses the uncertainties from three parameters in the conventional method. Indeed, the input parameters in the conventional method are the main sources of uncertainties in emissions estimates , as demonstrated in Liu et al (2020). For example, although the 1 km MODIS collection 6 active fire detection was used, Shiraishi and Hirata 2021 reports a total burned area of 86,000 km 2 , approximately two times of FINNv1.5 (40,300 km 2 ), across the study area. This large difference is because FINNv1.5 assumes a burned area of 0.75 km 2 for each 1 km MODIS fire detection and further scales it using vegetation cover fraction , while Shiraishi and Hirata 2021 assumes a burned area of 1 km 2 for each MODIS fire detection (Shiraishi and Hirata 2021). These assumptions ignore the physical principle of satellite fire detection algorithm, which allows MODIS to detect flaming fires (e.g., 1000 K) with a size of 0.0001 km 2 (0.01% fraction of a 1 km pixel) (Giglio et al 2003). When the 375m VIIRS active fire detections are also included to estimate burned area by FINN2.4, the total burned area from FINNv2.4 increases by 280% (relative to FINNv1.5) to 153,750 km 2 ,∼80% higher than that from Shiraishi and Hirata 2021. Nevertheless, the FINNv2.4 CO 2 mass is still 15% lower than the Shiraishi and Hirata 2021 estimation.
The fuel loading is certainly another major source of difference CO 2 estimates between FINN and Shiraisih & Hirata 2021. The FINN uses a static fuel loading lookup table summarized for global fires from published results , while Shiraishi and Hirata 2021 calculated fuel loadings from two global AGB maps (Shiraishi and Hirata 2021). Statistic fuel loadings do not consider spatial heterogeneity of fuels and also fail to account for fuel changes due to disturbances and vegetation regrowth (Gale et al 2021). The global AGB maps from satellite observations likely have limited capability of mapping surface dead and live fuel, which played a key role in the 2019-2020 bushfires (Nolan et al 2020), and their accuracy varies largely with availability of reference data and remains to be validated in many regions (Quegan et al 2017). Moreover, combustion completeness is difficult to measure for large wildfires especially as it can vary largely within the same fire with different fire behaviors (i.e., fire intensity and spreading speed) and fuel moisture (Veraverbeke and Hook 2013).
Third, it is difficult to pinpoint the potential factors that lead to overall larger GFED4s emissions than other inventories due to the lack of published or documented information on how GFED4s emissions of 2017 onward are computed exactly. According to the readme of GFED4 (https://www.geo.vu.nl/∼gwerf/GFED/GFED4/ Readme.pdf, last accessed on 30 August 2021), it merely mentions that GFED4s emissions after 2016 are generated using MODIS fire detections and empirical relationship between fire detections and historical GFED4s emissions. It is unknown if number of fire detections or FRP are used. Besides, fuel loading is another factor affecting emissions estimates. Nevertheless, it is interesting that GFED4s CO 2 is the lowest, differing from others by a factor of 3-10, in January 2020, although it is higher than the others for the whole study period.
The large discrepancies among various methods suggest that improving parameters in emissions estimation is most critical. For the FRP-based approach, as proposed in this study, the fusion of high spatial resolution fire observations from polar-orbiting satellites with high temporal resolution data from geostationary satellites is able to reconstruct diurnal FRP cycles and enhance FRE calculations, which can greatly reduce uncertainties in fire emissions estimates. For the conventional method, burned area can be accurately mapped using high resolution satellite data (e.g., Landsat-8 and Sentinel-2) and fuel loading can be reliably quantified by combining field-based data with the latest high resolution Lidar observations of forest structures (e.g., Global Ecosystem Dynamics Investigation (GEDI), Dubayah et al 2020).
To fully understand the accuracy of fire emissions estimation, it is urgent to develop robust validation approaches, which is very challenging. The common practice of validating fire emissions estimates is to simulate aerosols (e.g., AOD and PM2.5) using chemical transport models (CTMs) with fire emissions as the input and then compare the simulated values with the satellite-based and/or ground-based aerosol observations. However, such a validation practice is subject to uncertainties in model simulation of aerosols (Curci et al 2015). Alternatively, a direct comparison of fire emission estimates with high-resolution satellite observations of firereleased trace gases is of large potential if individual wildfires with fresh smoke plumes can be delineated (Li et al 2020b).
In summary, this study revealed the magnitude and intensity of an extreme fire event by (1) examining the fire intensity characteristics using the long-term MODIS FRP observations, and (2) quantifying the consumed biomass and released greenhouse gases and aerosol emissions by fusing the highest-temporal-resolution AHI FRP and finest-spatial-resolution VIIRS FRP. The fire intensity, biomass consumption, and emissions presented in this study are expected to complement other characteristics (e.g., fire area and smoke plumes) of the 2019-2020 bushfires and provide invaluable information for making sound policy for forest fuel management and emissions budget and understanding the impacts of fire emissions on air quality and climate. Because the frequency and intensity of extreme climate events have increased rapidly and will continue to increase under changing climatic conditions (Beniston et al 2007, Rummukainen 2012), our results for Australia's 2019-2020 bushfires provide new insight into the behavior of fire intensity and biomass-burning emissions, which is more critical than that of regular fire emissions.