The Carbon Cycle of Southeast Australia During 2019–2020: Drought, Fires, and Subsequent Recovery

.

BYRNE ET AL.

10.1029/2021AV000469
2 of 20 Burrows, 2002). However, the region is experiencing more frequent, extensive and severe fires (Pitman et al., 2007;Stephens et al., 2013), a trend that is expected to continue with climate change (Abatzoglou et al., 2019;Di Virgilio et al., 2019;Dowdy et al., 2019;Perkins-Kirkpatrick & Gibson, 2017). Despite the adaptations of Australian ecosystems to fire, these changing fire regimes have been shown to impact tree mortality D. M. J. S. Bowman et al., 2014) and threaten the persistence of some forest biomes in Australia (D. M. J. S. Bowman, Williamson, Price, et al., 2020;Fairman et al., 2016;D. M. J. S. Bowman et al., 2014), including their carbon stores. Thus, monitoring the response of ecosystems in southeast Australia to extreme drought, heat and fire is critical for understanding how the carbon balance of this region will evolve under climate change.
Since 2017, southeast Australia has been in drought, with the 2017-2019 period having the largest 3 year rainfall deficit since 1900 (King et al., 2020). These conditions have been most extreme during 2019, which was the hottest and driest year recorded in southeast Australia (Abram et al., 2020;Bureau of Meteorology, 2020), provoking one of the worst bushfires seasons in recorded history Collins, Bradstock, et al., 2021;D. M. J. S. Bowman et al., 2021;Deb et al., 2020;King et al., 2020;Nolan et al., 2020;Ward et al., 2020). These extreme conditions subsided in early February 2020 with heavy rainfall and cooler conditions that persisted throughout the austral autumn. This combination of drought and fire followed by heavy rainfall imparts a large and complex perturbation on the carbon cycle of the region and impacted forested regions that cover much of the southeast coast and mountainous regions, and more arid savanna, grassland and cropland ecosystems further inland ( Figure 1).
The impact of extreme drought and heat events on ecosystems are complex and challenging to monitor. Ecosystem responses are sensitive to the specific characteristics of the event, such as the intensity and timing (Bastos et al., 2014;De Boeck et al., 2011;Denton et al., 2017;Frank et al., 2015), legacy effects from previous disturbances Longo et al., 2020) and vegetation type (Turner et al., 2020;Zhang et al., 2016). The recent expansions of space-based observing systems of carbon-cycle-relevant quantities are now providing the opportunity for finer scale quantification of carbon cycle perturbations and more detailed understanding of the response of ecosystems to extreme drought, heat, and fire (Byrne, Liu, Lee, et al., 2020;Byrne, Liu, Bloom, et al., 2020;Byrne et al., 2019;Turner et al., 2020;Yin, Bloom, et al., 2020). In this study, we utilize space-based 3 of 20 observations to provide a comprehensive analysis of the carbon cycle perturbations due to extreme drought, heat and fire during the 2019-2020 growing season in southeast Australia.
We combine observations from multiple satellites to quantify the carbon cycle anomalies within southeast Australia. We employ TROPOspheric Monitoring Instrument (TROPOMI) CO column abundance measurements (Borsdorff et al., 2018;Landgraf et al., 2016) to quantify biomass burning emissions. Anomalies in net ecosystem exchange (NEE, which is defined as the residual between TER and GPP) are obtained by combining top-down constraints on net surface-atmosphere CO 2 fluxes from column-averaged dry-air mole fractions of CO2 (XCO 2 ) measurements from the Orbiting Carbon Observatory 2 (OCO-2) Eldering et al., 2017) with estimates of GPP anomalies from FluxSat (Joiner & Yoshida, 2020), which produces GPP from MODerate-resolution Imaging Spectroradiometer (MODIS) reflectances trained against FLUXNET sites.
The combination of these newly available observations offers a unique opportunity to monitor individual components of the carbon cycle anomalies across southeast Australia during 2019-2020. Specifically, we aim to answer: How much CO 2 was released to the atmosphere due to drought and biomass burning, respectively? How did this event impact forest and non-forest ecosystems differently? What were the differences in carbon cycle perturbations between burned and unburned ecosystems? And how does 2019-2020 compare with previous years? To that end, we first quantify biomass burning emissions of CO using the TROPOMI observations, which are converted to CO 2 emissions (Section 3.1). Next, an anomaly in atmospheric CO 2 (ΔCO 2 ) is derived using the OCO-2 measurements (Section 3.2). This top-down constraint is then combined with estimates of GPP anomalies using FluxSat to derive NEE anomalies over the 2019-2020 growing season (Section 3.3). We then synthesize these estimates and present the evolution of carbon cycle anomalies over the 2019-2020 growing season (Section 4), and compare this event with the regional carbon budget over the 2010-2019 period (Section 5). This is followed by a discussion of our biomass burning emission estimates in the context of previous bottom-up and top-down estimates (Section 6.1), the implications of this extreme event for the carbon cycle of southeast Australia (Section 6.2), and the uncertainties and remaining challenges in estimating carbon fluxes from extreme events (Section 6.3). Finally, we provide our conclusions in Section 7.

Environmental and Geographical Data
Environmental and geographical data are used to help interpret the carbon cycle anomalies. We examine the covariations of carbon cycle anomalies with variations in soil temperature and soil moisture from ERA5-Land reanalysis (Munoz Sabater, 2019), generated using Copernicus Climate Change Service Information 2020. For this analysis, we calculate the area-weighted soil moisture and temperature over the top 1 m of soil. Vegetation land cover is obtained from the MODIS land cover dataset (MCD12C1; Friedl & Sulla-Menashe, 2015) and elevation data are obtained from ETOPO1 (Amante & Eakins, 2009). Figure 2 shows a schematic diagram of the methods used to estimate biomass burning emissions and anomalies in NEE (ΔNEE). Biomass burning CO 2 emissions are estimated from TROPOMI X CO measurements (Section 3.1). First, emissions of CO are estimated through flux inversion analyses that assimilate TROPOMI X CO measurements. Then CO emissions are converted to CO 2 emissions using emission factors.

CO 2 Flux Estimates
Estimates of ΔNEE are obtained through combining several different data sources. First, we infer a top-down CO 2 anomaly signal (ΔXCO 2 ) due to anomalies in biosphere-atmosphere CO 2 fluxes (Section 3.2). Then we subtract the ΔXCO 2 signal due to biomass burning emissions, giving ΔXCO 2 due to ΔNEE. This provides a constraint on the magnitude of ΔNEE. Finally, we estimate the spatiotemporal structure of ΔNEE by combining the atmospheric CO 2 constraints with FluxSat GPP (Section 3.3). Note that the CO 2 flux and atmospheric XCO 2 are related to fluxes using a chemical transport model (Section 3.1.1).
Atmospheric chemical transport simulations and flux inversions are performed with the Greenhouse Gas Framework-Flux (GHGF-Flux) inversion system. GHGF-Flux is a flux inversion system developed under the NASA Carbon Monitoring System Flux (CMS-Flux) project (https://cmsflux.jpl.nasa.gov), and inherits the chemistry transport model from GEOS-Chem and its adjoint (Henze et al., 2007;Liu et al., 2014). Chemical transport is driven by the Modern Era Retrospective Analysis for Research and Applications, Version 2 meteorology BYRNE ET AL.

10.1029/2021AV000469
4 of 20 produced with version 5.12.4 of the Goddard Earth Observing System (GEOS) atmospheric data assimilation system (Gelaro et al., 2017). To perform tracer transport, these fields are regridded to the desired horizontal resolution and archived with a temporal resolution of three hours except for surface quantities and mixing depths, which have a temporal resolution of one hour. Flux inversions are performed using 4-D variational assimilation (4D-Var), with the details provided in the subsections.

Biomass Burning Emissions
Atmospheric CO inversions have been shown to be an effective top-down approach for estimating fire carbon emissions (Langenfelds et al., 2002;Liu et al., 2017;Yin, Bloom, et al., 2020;Yin et al., 2015Yin et al., , 2016Zheng et al., 2019). Here, we perform atmospheric CO inversions to estimate biomass burning emissions by assimilating TROPOMI retrievals of (X CO ). TROPOMI is a grating spectrometer aboard ESA's Sentinel-5 Precursor (S-5P) satellite that measures Earth reflected radiances (Veefkind et al., 2012). CO total column densities are retrieved in the shortwave infrared (around 2.3 μm) using the Shortwave Infrared CO Retrieval algorithm (Landgraf et al., 2016). Retrieved CO total column densities are then converted to dry-air mole fractions of CO (X CO ) using the dry-air surface pressure and hypsometric equation. The column averaging kernel is similarly converted to mole-fraction space.
Biomass burning CO emissions are estimated using one-way nested flux inversions over Australia (100°−177.5° E and 0°−60° S) at 0.5° × 0.625° spatial resolution. Nested flux inversions are performed from 5 November 2019 through 14 January 2020 (to cover the period with the majority of fires) and assimilate TROPOMI X CO super-obs (aggregated observations) to optimize scaling factors for each gridcell over the entire period. Details Combining the mean NEE seasonal cycle over this period with a chemical transport model, we simulate the expected 2019-2020 baseline atmospheric CO 2 fields given climatological fluxes. Then, the difference between the actual 2019-2020 measurements and the expected XCO 2 gives the anomaly in atmospheric XCO 2 (shown in blue shaded area). ΔNEE is then estimated from combining all of the constraints. The spatiotemporal structure of ΔNEE is based on FluxSat gross primary productivity (shown in green), while the magnitude is derived from combining the top-down and biomass-burning-derived ΔCO 2 estimates (shown in purple). BYRNE ET AL. Eight nested flux inversions are performed, which vary in prior biomass burning emissions, quantities optimized, and boundary conditions (Table 1). Differences in flux inversion configuration are employed to test the sensitivity of posterior fluxes to the inversion set-up. We employ two different biomass burning emissions data sets as prior CO fluxes, namely the Global Fire Emissions Database version 4 (GFED4.1s; van der Werf et al., 2017) and Global Fire Assimilation System (GFAS; Kaiser et al., 2012). GFED4.1s provides estimates of biomass burning using MODIS 500 m burned area (Giglio et al., 2013), 1 km thermal anomalies, and 500 m surface reflectance observations to statistically estimate burned area associated with small fires (Randerson et al., 2012). GFAS v1.2 provides estimates of daily biomass burning emissions by assimilating MODIS fire radiative power observations (Di Giuseppe et al., 2018;Kaiser et al., 2012). For both data sets, we incorporate the impact of the diurnal cycle based on Mu et al. (2011). The inversions also differ by either prescribing or optimizing diurnal variations on biomass burning emissions. Finally, inversions are either run using boundary conditions from a global TROPOMI flux inversion or with these boundary conditions adjusted by adding 10 ppb (roughly equivalent to the mean data-model difference) at all levels and times to test the sensitivity of the nested CO inversion to lateral boundary conditions.
Video 1 shows the spatial distribution of the mean posterior fluxes and X CO measurements across southeast Australia. Biomass burning emissions were most concentrated in forest ecosystems along the coast and further inland along the border between New South Wales and Victoria. Posterior CO emissions are increased for all inversion configurations, with a posterior mean CO emission estimate of 15.6 TgC (range: 9.7-24.3 TgC), relative to prior emission estimates of 11.4 TgC for GFED and 5.8 TgC for GFAS over the growing season. The largest source of spread among posterior fluxes is due to the prior biomass burning flux employed, with GFED-based inversions giving larger posterior emissions than GFAS-based inversions (see Figure S2 in Supporting Information S1).
The performance of the nested CO flux inversions are evaluated by comparing the posterior CO fields with the TROPOMI X CO measurements, independent X CO measurements from the nearby Wollongong  and Lauder (Pollard et al., 2017(Pollard et al., , 2019 Total Column Carbon Observing Network (TCCON; Wunch et al., 2011) sites, the Cross-track Infrared Sounder (CrIS) and surface-based flask and in situ measurements at the nearby Cape Grim (CGO), Baring Head (BHD), and Lauder (LAU) sites. CrIS is a Fourier Transform Spectrometer aboard the satellite Suomi-NPP and has a spectral resolution of 0.625 cm −1 and a ground pixel diameter of 14 km at nadir. CrIS and TROPOMI make collocated measurements because Suomi-NPP and Sentinel 5p are in a tandem orbit with a roughly 10 min separation. However, CrIS takes measurements in both day and night. The retrieval of CO uses the MUlti-SpEctra, MUlti-SpEcies, Multi-SEnsors (MUSES) algorithm (Fu et al., 2016) that is based on the optimal estimation method with heritage from the tropospheric emission spectrometer (K. W. Bowman et al., 2006). We generate X CO measurements from version 1.8 of the L2 tropospheric CO profile product, and compare posterior CO fields against daytime and nighttime X CO measurements separately.
As trace gas emissions from fires are impacted by pyroconvective motions (that are not well represented in chemical transport models), we evaluate the posterior fluxes with two sets of model runs that release the CO emissions at different model levels. In one set of runs, we release the emissions at the surface (as was done in the inversion), while in the second set we release CO emissions at the injection height (mean altitude of maximum injection) simulated by a plume rise model (IS4FIRES; Rémy et al., 2017), which was provided with the GFAS emission data. Here we provide a brief summary of the evaluation, while a detailed evaluation of the flux inversions is presented in Text S1 in Supporting Information S1. Posterior fluxes generally show better agreement with the TROPOMI, TCCON, CrIS, and the in situ/flask measurements. This is true for all measurements and a subset of measurements that are biomass burning sensitive. However, posterior CO fluxes tend to underestimate X CO for biomass-burning-sensitive measurements (but less so than the prior). This residual mismatch is likely related to transport model errors, as the modeled observations often show differences in plume structure (Video 1). Furthermore, the transport model underestimates vertical motions around the bushfires, which were impacted by pyroconvection. The impact of weak modeled vertical motions can be seen in Video 1c,d. The column averaging kernel for TROPOMI shows greater sensitivity to CO between 400 hPa and the surface, while CrIS shows greater sensitivity to CO in the upper troposphere. Both TROPOMI and CrIS show mean X CO mole fractions greater than 200 ppb in southeast Australia for the duration of the biomass burning over November-January. However, 6 of 20 posterior data-model mismatches are much less positive for TROPOMI than for CrIS, implying that vertical motions are underestimated and the CO emissions do not reach the upper troposphere to the levels observed.
Finally, to estimate CO 2 biomass burning emissions we apply the ratio of CO 2 to CO emission factors (that are constant in time). We apply the emission factors from the biomass burning database used as the prior (e.g., either GFAS or GFED). The emission ratios are variable by vegetation type, but aggregating for fires across Australia gives effective CO 2 /CO emission ratios of 12.01 gC gC −1 for GFED and 11.30 gC gC −1 for GFAS. Differences are primarily driven by differences in emission factors for forest emissions, but are within the natural variation of emission factors reported by Akagi et al. (2011) (see Text S2 and Figure S4 in Supporting Information S1) and reported for Australian forests (Table S5 in Supporting Information S1) (Guérette et al., 2018;Paton-Walsh et al., 2014). The impact of emission factor uncertainty is further discussed in Section 6.3.2.

Atmospheric ΔCO 2 Signal Simulation
We simulate the biomass burning XCO 2 anomaly signal (ΔXCO 2 BB) by running the nested chemical transport model. The ΔXCO 2 BB signal is calculated by performing simulations with climatological fluxes and with the climatological fluxes plus the biomass burning estimates, then taking the difference between these two simulations at the OCO-2 and TCCON measurements locations to isolate the signal due to biomass burning. We simulate OCO-2 good-quality land (land glint and land nadir) and ocean glint super-obs (aggregated to 0.5° × 0.5° resolution grids following Liu et al. (2017), with the additional requirement that there must be a minimum of three OCO-2 observations within each 0.5° × 0.5° grid box per track). For TCCON measurements, we use all good quality data.

Top-Down ΔCO 2 Signal
The top-down estimate of ΔXCO 2 (ΔXCO 2 top−down) is calculated based on the data-model difference between OCO-2 and TCCON measurements and simulated CO 2 fields based on climatological NEE and ocean fluxes and ODIAC fossil fuel emissions (Oda & Maksyutov, 2011;Oda et al., 2018).
Climatological NEE and ocean fluxes are generated through CO 2 flux inversion analyses as the average over the period 2010-2018. Fluxes over 2010-2014 are taken as the mean GOSAT + surface + TCCON inversion of Byrne, Liu, Lee, et al. (2020). To generate climatological fluxes over 2015-2019, we perform a flux inversion at 4° × 5° assimilating OCO-2 measurements and surface-based CO 2 measurements concurrently and use the identical inversion set-up to Byrne, Liu, Lee, et al. (2020). For surface measurements, we use version 6.0 of the GLOBALVIEW plus package (Masarie et al., 2014;Schuldt et al., 2020). For OCO-2 measurements, we use ACOS b10 land (land glint and land nadir) and ocean glint retrievals aggregated into super-obs at 0.5° × 0.5° resolution grids following Liu et al. (2017), with the additional requirement that there must be a minimum of three OCO-2 observations within each 0.5° × 0.5° grid box per track. We use all data that pass the quality flag filter. This 2015-2019 flux inversion is referred to as the "IS + LNLGOG" inversion.
Calculations of the ΔXCO 2 top−down signal are performed with the one-way nested grid over Australia. First, we generate boundary conditions by performing a simulation at 2° × 2.5° with regrided optimized NEE and ocean fluxes and prescribed fluxes from the 4° × 5° flux inversion for 2019-2020. Then we run the nested model with the climatological NEE and ocean fluxes and ODIAC fossil fuel emissions (we impose 2018 ODIAC emissions for 2019 and 2020 based on the availability of data). Simulated CO 2 fields are sampled for OCO-2 and TCCON observations from 1 October 2019 through 31 January 2020. Finally, we calculate the ΔXCO 2 top−down anomaly signal as the data-model mismatch for these simulated observations.

NEE Anomaly Estimate
NEE anomalies (ΔNEE) over the 2019-2020 growing season are estimated by combining the constraints on GPP from FluxSat Version 2 (Joiner & Yoshida, 2020) with the constraints on the net CO 2 flux from the top-down ΔXCO 2 top−down signal and biomass-burning-ΔXCO 2 BB . The spatial and temporal structure of ΔNEE is assumed to be directly proportional to ΔGPP from FluxSat, while the magnitude of the ΔNEE is inferred from the atmospheric ΔXCO 2 signal.

of 20
FluxSat estimates GPP based on Nadir BRDF-Adjusted Reflectances (NBAR) from the MODIS MYD43D product (Schaaf et al., 2002). The GPP estimates are calibrated with the FLUXNET 2015 GPP derived from eddy covariance flux measurements at Tier 1 sites (Joiner & Yoshida, 2020). The native spatiotemporal resolution of FluxSat GPP is daily on a 0.05° × 0.05° grid. For our analysis, we regrid spatially to 0.1° × 0.1° while retaining daily temporal resolution. We calculate ΔGPP from FluxSat as the difference between fluxes for 2019-2020 relative to a 2010-2018 mean.
To estimate ΔNEE, we assume ΔNEE ∝ − ΔGPP. Empirical evidence from the OzFlux eddy covariance network supports this assumption. Li et al. (2017) found that ΔNEE can be expressed linearly as a function of ΔGPP with reasonable accuracy, and obtain the relationships of ΔNEE = −0.24 ΔGPP for non-forest ecosystems, where anomalies in GPP and respiration are correlated, but ΔNEE = − 0.8 ΔGPP for forest ecosystems, where GPP and respiration do not co-vary. To estimate the magnitude of ΔNEE, we simulate the OCO-2 observed XCO 2 anomaly signal due to ΔGPP (ΔXCO 2 GPP) using the same approach as was used for biomass burning (See Section 3.1.1). We invert a magnitude of ΔNEE through regressions of ΔXCO 2 NEE against an observationally constrained anomaly in XCO 2 : Note that β is included to account for possible small residual biases from the observations or model. Initially, we attempted a multivariate regression to solve this for forest and non-forest ΔXCO 2 NEE individually but recovered unrealistic negative coefficients for forests. The ΔXCO 2 NEE for forests is relatively small and may be impacted by errors in biomass burning emissions and transport, potentially limiting our ability to differentiate forest and non-forest ΔNEE. To avoid these unphysical values, we prescribe the ratio of ΔNEE between forest and non-forest ecosystems. Following from Li et al. (2017), we perform one set of regressions using, However, due to the large CO 2 biomass burning emissions over this event, it is possible that ΔNEE and ΔGPP may diverge from this relationship. Therefore, we also perform a second set of regressions using the relationship: We perform a series of linear regressions using Equation 1 to estimate "α," the parameter that relates ΔNEE and ΔGPP. This regression is performed a total of 32 times by varying the emission height of biomass burning emissions between the surface and injection height, the posterior biomass burning emissions estimated by the eight TROPOMI flux inversions, and the parameterization relating forest and non-forest ΔNEE using Equations 2 and 3. Table 2 shows the statistics of α for these 32 regressions. The best estimate of α is then calculated as the mean of the truncated distribution of the 32 α values, with the largest and smallest two values removed, and the range of the truncated distribution is taken as the uncertainty. This gives an α of 0.41 (0.23-0.66) for forest ecosystems, which is half the value of Li et al. (2017), and 0.23 (0.13-0.35) for non-forest ecosystem, which is almost identical to the value of Li et al. (2017). These α values are applied to estimate ΔNEE over the entire growing season.
A comparison of ΔXCO 2 top−down and the simulated ΔXCO 2 GPP + ΔXCO 2 BB signal for TCCON and OCO-2 measurements is shown in the Figures S6 and S7 in Supporting Information S1. The flux estimates found here are generally consistent with these top-down data sets, although there is considerable scatter between different TC-CON sites and OCO-2 viewing modes. Data-model mismatches for individual retrievals are strongly impacted by retrieval errors and model errors in simulating the observations, however, aggregating to 0.2 ppm intervals in the flux signal reveals strong positive correlations (R 2 > 0.9, Figure S6 in Supporting Information S1). Similarly, simulated boundary layer CO 2 at CGO and LAU shows improved agreement with the measurements when the flux anomalies are included, while results at BHD are mixed (Table S6 in Supporting Information S1).

Carbon Cycle Anomalies Over the 2019-2020 Growing Season
The climate anomalies over the 2019-2020 growing season can be partitioned into two phases. Warm-dry conditions dominated the region during the austral spring and early summer (October through January), when there were a number of biomass burning events, primarily in the forested regions. This was followed by a cooler-wetter 8 of 20 period during February through May (Figures 1a and 1b). Video 2 shows the evolution of ΔNEE and biomass burning over the 2019-2020 growing season. During the warm-dry phase, GPP was suppressed across the region, falling below the range of observed GPP over the 2010-2018 period (2.0 gC m −2 day −1 for October-January 2019-2020 vs. 3.0-4.3 gC m −2 day −1 over 2010-2018). Suppression of productivity occurred uniformly across southeast Australia during October-January (Figure 3), impacting both forest and non-forest ecosystems. This is followed by a large-scale recovery in GPP to above average values during February-May, when cooler-wetter conditions dominate. This recovery was relatively uniform across the region with the exception of burned areas (indicated by hatching in Figure 3), which show suppression of GPP during February-May that is similar to October-January. Figure 4 shows the timeseries of ΔGPP for forest and combined non-forest ecosystems (includes cropland, grassland, shrubland, and savanna ecosystems) over southeast Australia (145.5-154.5 E and 28.5-38.5 S). We divide forests into burned and unburned regions using a threshold of 50 gC m −2 of biomass burning emissions over the 2019-2020 growing season for each 0.1° × 0.1° grid cell. For non-forested regions, GPP was suppressed during October-January (54% below mean), but rapidly recovered to above average when cooler-wetter conditions dominate (33% above mean for February-May). In the unburned forested regions, GPP was suppressed during October-January (23% below mean), with a partial recovery during February-May (8% below mean). In contrast, the burned forests showed a larger reduction in GPP during October-January (37% below mean) that persisted throughout February-May (31% below mean). Similar reductions are found for MODIS near-infrared reflectance of terrestrial vegetation (NIR V ) and solar induced fluorescence (SIF) measurements from TROPOMI and OCO-2 for these vegetation types (see Text S3 and Figure S5 in Supporting Information S1). The similar reduction in NIR V suggest that structural changes in vegetation are partially responsible for the reductions in GPP (He et al., 2020;Y. Sun et al., 2015;Yoshida et al., 2015), and are consistent with site level observations of foliar death in eucalypt forests during 2019-2020 . In total, 166 TgC (range: 113-236 TgC) of CO 2 was released through biomass burning and 33 TgC (range: 19-52 TgC) was released due to anomalies in NEE over October-May (Table 3).

Impact of 2019-2020 Anomalies on the Regional Carbon Budget
To contextualize the carbon loss over the 2019-2020 growing season, it is useful to compare this period to the long-term mean. Here, we compare the estimated net biosphere exchange (NBE, sum of NEE and biomass burning) for 2019-2020 relative to CO 2 flux inversions spanning the 2010-2019 growing seasons. For simplicity, we will refer to growing seasons as YY 1 /YY 2 (e.g., "19/20"), which encompass July of YY 1 through June of YY 2 .  Figure 5 shows NBE estimates for southeast Australia over the period 10/11 through 19/20. We find that the magnitude of the 19/20 NBE anomaly (mean: 186 TgC year −1 , range: 118 to 275 TgC year −1 ) significantly exceeds NBE variability over the 10/11 to 18/19 period, confirming the extreme magnitude of this event. The mean annual net NBE sink over 10/11-18/19 is found to be −9.5 TgC year −1 (range: −16.1 to −3.4 TgC year −1 ). However, this mean value is a small residual of considerable inter-annual variations (standard deviation of 40 TgC year −1 ), ranging from sink of −73 TgC year −1 (range: −114 to −41 TgC year −1 ) in 10/11, driven by a strong La Niña (Poulter et al., 2014), to a source of 57 TgC year −1 (range: 28 to 99 TgC year −1 ) during the 18/19 drought. The magnitude of interannual variations in NBE are independently confirmed by FluxSat-based ΔNEE (calculated using the regressions from Section 3.3).
Interannual variations in NBE are found to be closely associated with climate variability. Strong correlations are obtained with soil temperature and moisture for the flux inversions (R 2 = 0.69/0.69 for temp/moist) and FluxSat-based ΔNBE (R 2 = 0.60/0.87 for temp/moist). We further examine the relationship between climate variability and uptake over the 18 year period covering 01/02-18/19 for both forests and non-forests with FluxSat GPP ( Figure S8 in Supporting Information S1). Over this longer period, we find that soil moisture variability is strongly correlated with variability in FluxSat GPP for forests (R 2 = 0.77) and non-forests (R 2 = 0.88). However, soil temperature is not correlated with GPP for forests (R 2 = 0.06) and moderately correlated with non-forests (R 2 = 0.57), suggesting that moisture availability is the primary driver of interannual variations in productivity.
This relationship between carbon uptake and climate variability has significant implications for the recovery of these ecosystems. The rate at which the 19/20 carbon loss will be re-absorbed may depend strongly on climate variability and change. Based on the mean NBE estimate of −9.5 TgC year −1 (range: −16.1 to −3.4 TgC year −1 ) over 10/11 to 18/19, we would expect an anomalous carbon release of 199 TgC year −1 (range: 131-288 TgC year −1 ) to be recovered in ∼21 years (range: 8-85 years). However, having ideal cool-wet conditions like 10/11 could shorten this to six years, while hot-dry conditions could prevent a full recovery indefinitely. This highlights the importance of both quantifying the sensitivity of these ecosystems to climate variability, and accurately projecting regional climate changes for predicting the recovery of these ecosystems. The emission estimates calculated in this study are "top-down," in that they are based on observations of the emitted trace gases in the atmosphere. Thus, the different approaches are complementary, and consistency between top-down and bottom-up estimates provides increased confidence in emission estimates. Our estimate of CO 2 emissions overlaps with all existing burned-area-based estimates of biomass burning CO 2 over southeast Australia (Table 4), providing increased confidence in these estimates. However, our estimated range suggests larger emissions than provided by the GFAS radiative-power-based method, suggesting the GFAS underestimates biomass burning over southeast Australia during 2019-2020. A similar top-down estimate of biomass burning CO 2 emissions was reported by van der Velde et al. (2021), and agrees well with the estimates reported here.

Implications for Southeast Australia
The 2019-2020 carbon loss of 199 TgC (range: 131-288 TgC) significantly exceeds annual CO 2 flux anomalies of any year since 2010 (∼5σ anomaly) and exceeds total annual Australian CO 2 emissions from fossil fuel . This demonstrates the impact that extreme events can have on the regional carbon budget, and suggests changes in the frequency of extreme heat, fire weather and drought could have a strong impact on the regional carbon balance.
During the study period, there was a robust recovery for unburned ecosystems, suggesting that drought-and heat-induced carbon losses over southeast Australia will be rapidly re-absorbed. This is consistent with recent modeling work suggesting resilience to drought in southeast Australian forests . However, this rapid recovery may be unusual, due to heavy rainfall and below average temperatures during the autumn, which strongly modulate productivity in dryland ecosystems (Haverd et al., 2017;Huxman et al., 2004). Furthermore, there may be drought-induced damages to these ecosystems that are not captured in this analysis, such as drought-induced tree mortality, which has the potential to impact species and biomass composition (Batllori et al., 2020;Burton et al., 2021;Fensham et al., 2019). Recovery in burned ecosystems was much more muted, consistent with major structural damage, preventing a rapid recovery when favorable conditions return.
For the years ahead, the speed and extent of carbon uptake will depend strongly on the climate conditions. The top-down 2010-2019 mean annual sink of −9.5 TgC year −1 (range: −16.1 to −3.4 TgC year −1 ) suggests that a full recovery of carbon pools will take 21 years (range: 8-85 years). However, the regional net annual flux showed large interannual variations closely linked with variability in temperature and moisture, with cooler-wetter years being associated with increased uptake. This is consistent with site-level observations showing that rainfall is   In addition to the carbon sequestration, it is important to consider ecosystem recovery. Severe fires can have legacy impacts on ecosystem function even after carbon stocks have been largely regenerated. Severe fires in Eucalyptus forests have been shown to include persistent changes to canopy structure (Karna et al., 2019), increase tree mortality Etchells et al., 2020) and cause changes in understory composition and structure (D. M. J. S. Bowman et al., 2014;Fairman et al., 2016;Fletcher et al., 2014;Pellegrini et al., 2021). In particular, fire-induced tree mortalities are generally higher in smaller-younger cohorts Bowd et al., 2021). For some species, such as Mountain Ash, ∼20 years are required to reach maturity and produce Note. The median and range of α are given for regressions using the eight posterior biomass burning estimates for simulations that vary in the emission height and forest/non-forest parameterization. The bottom row gives the mean and range for the truncated distribution of all simulations, wherein we remove largest and smallest two outliers from the 32 simulations performed by varying biomass burning emissions, emission height, and the forest/non-forest parameterization. 13 of 20 viable seeds , potentially threatening the survival of these species if subject to recurrent fire for several decades.
Finally, the impact of disturbance and ecosystem recovery should be considered within the context of ongoing climate change. This region is experiencing more frequent extreme heat and fire events (Abram et al., 2020;Bradstock et al., 2014;Sharples et al., 2016), a trend that is expected to continue with climate change (Abatzoglou et al., 2019;Di Virgilio et al., 2019;Dowdy et al., 2019;Herold et al., 2021;Perkins-Kirkpatrick & Gibson, 2017). Projections of regional trends in drought are less certain. However, several studies suggest drought intensity and frequency may increase in the coming years over much of southeast Australia (Herold et al., 2021;J. Wang et al., 2021;Kirono et al., 2020;Ukkola et al., 2020), with Herold et al. (2021) finding 1-in-20 year droughts may become 1-in-5 year events by 2060-2079.
Trends in climate variability will likely have a number of impacts on the carbon cycle of the region. A major risk for forest ecosystems is recurrent fires during recovery from the previous fires, with studies finding significant negative impacts on ecosystem function for both obligate seeder (D. M. J. S. Bowman et al., 2014) and resprouter-dominated communities (Collins, Hunter, et al., 2021;Fairman et al., 2017). This risk may be compounded by longer recovery periods after fire due to frequent extreme heat and drought events. In particular, if the inter-fire interval decreases below the recovery time, permanent carbon loss will be experienced by these ecosystems, potentially leading to major changes in the ecosystem structure and fire regimes of the region (Boer et al., 2016).
Detecting permanent changes in the regional carbon budget will require sustained monitoring of the regional carbon budget through a combination of expanding top-down constraints (Crisp et al., 2018), as presented in this work, in addition to continued and improved site-level monitoring (D. M. J. S. Bowman, Williamson, Yebra, et al., 2020).

Uncertainties in Estimating Carbon Flux
In this analysis, we have calculated drought-induced NEE anomalies and biomass burning CO 2 anomalies over southeast Australia during 2019-2020 that are consistent with observed X CO , XCO 2 and FluxSat GPP. Still, there are remaining challenges in quantifying carbon cycle perturbations, leading to large uncertainties in the estimates presented here.

Model Transport
Accurate representation of atmospheric transport of CO and CO 2 from biomass burning remains a major challenge (Eastham & Jacob, 2017). Rapid pyroconvective motions are not well represented in our model simulations. This leads to errors in simulated X CO fields relative to the observations and systematic errors in flux inversions.
In our analysis, we performed sensitivity analyses by evaluating the posterior CO fields for emissions released at the surface and at an estimated plume injection height (emitted at up to 6 km in altitude, Text S1 and Figure S1 in Supporting Information S1), and found that the posterior emissions better matched independent CO observations in both cases. Still, Modeled CrIS X CO , which are most sensitive to the upper troposphere, showed weak sensitivity to biomass burning emissions despite the fact that biomass burning species were observed in the stratosphere (Hirsch & Koren, 2021;Khaykin et al., 2020;Schwartz et al., 2020). This suggests that modeled vertical motions are too weak and do not fully capture the vertical structure of biomass burning species produced by strong pyroconvective motions. Such systematic errors are challenging to address, but one possible avenue of future study would be to utilize weak constraint 4D-Var (Stanevich et al., 2019), which would allow for optimizing both surface fluxes and the atmospheric state. Accounting for the total CO change throughout the column would provide a quantitative assessment of the impact of systematic transport errors on CO emission estimates. Another avenue of future work could be to improve the   representation of pyroconvective motions in transport models. As these motions are sub-grid scale for typical chemical transport models, this would most likely require prescribing vertical mass fluxes calculated by a high-resolution cloud resolving model.

CO 2 /CO Emission Ratio
To estimate biomass burning CO 2 emissions from estimated CO emissions, the CO 2 /CO emission ratio needs to be precisely and accurately known. However, there is considerable uncertainty in this value, with recent reported values for Australian forests ranging from 8.59 ± 1.16 gC gC −1  to 12.65 ± 2.34 gC gC −1 (Guérette et al., 2018) (Table S7 in Supporting Information S1). We incorporated some of this uncertainty by applying different emission ratios for the GFAS (9.44 gC gC −1 for forests) and GFED (11.91 gC gC −1 for forests) based biomass burning estimates. Comparison of simulated and measured XCO 2 and X CO retrievals at Wollongong and Lauder supports the emission ratios employed here ( Figure S9 in Supporting Information S1). For Wollongong, we found an observed XCO 2 /X CO ratio of 0.014 ppm ppb −1 (range: 0.011-0.036 ppm ppb −1 ) and a simulated XCO 2 /X CO ratio of 0.017 ppm ppb −1 (range: 0.009-0.021 ppm ppb −1 ), while the dynamic range of biomass-burning-impacted measurements at Lauder was not sufficient to provide a strong constraint on the emission ratio (note that XCO 2 /X CO ratios are not directly comparable to emission ratios due to chemical loss of CO). Still, we acknowledge that uncertainty in the CO 2 /CO emission ratio remains a major challenge in estimating CO 2 biomass burning emissions from CO flux inversion analyses.

Data Gaps
It is also notable that the largest biomass burning enhancements of XCO 2 were not observable by OCO-2 or TC-CON sites due to the presence of co-emitted aerosols (J. Wang et al., 2020). Rapid deployment of aircraft campaigns that observe the chemical composition of the biomass burning plumes would help mitigate these sampling biases. The serendipitous occurrence of the Atmospheric Carbon and Transport-America (ACT-America) flight campaign during the 2019 Midwest floods provided supporting evidence of the flood-induced CO 2 flux anomalies estimated by Yin, Byrne, et al. (2020), resulting in increased confidence in those estimates.

Estimating ΔNEE
Due to atmospheric mixing and the relatively sparse sampling of XCO 2 by OCO-2, it is not possible to fully resolve the spatial and temporal structure in ΔNEE. Thus, we utilized the spatiotemporal structure of ΔGPP to predict the spatiotemporal structure of ΔNEE for forests and non-forest to regularize the problem. Although, this linear relationship is generally supported by eddy-covariance measurements within Australia (Li et al., 2017), there are likely many cases where this linearity breaks down. We do not account for this source of systematic error in our analysis, suggesting that the uncertainties may be larger and more systematic than estimated here.

Unaccounted for Carbon Fluxes
Finally, we note that, we only quantify land-atmosphere CO 2 fluxes in this study, and that a full accounting of the carbon stock changes due to this event would need to incorporate lateral carbon fluxes. Intense rainfall following immediately after fire likely increased runoff of ash and debris to waterways, leading to a number of record fish kills in estuarine sites located downstream of burned areas (Silva et al., 2020). Thus, there may have been considerable export of carbon to waterways and the ocean, but this has not been quantified to our knowledge. In addition, we only include estimates of biomass burning emissions of CO 2 . We estimate an additional 15-29 TgC emitted as CO from biomass burning.

Conclusions
Extreme events play a major role in the carbon cycling of ecosystems, but quantifying the impact of these events on the carbon budget remains challenging. Incorporating a variety of space-based observations, we have provided a comprehensive accounting of biosphere-atmosphere CO 2 flux anomalies due to drought, heat, and fire over southeast .5 E and 28.5-38.5 S) during the 2019-2020 austral growing season. In total, biomass burning released 113-236 TgC of CO 2 and anomalies in October-May NEE reduced carbon uptake by 19-52 TgC. Carbon losses were found to be most severe in forested regions and were dominated by biomass burning emissions. Unburned forests and non-forest ecosystems recovered to mean or greater productivity when 15 of 20 cooler-wetter conditions dominated during the late austral summer and autumn, however, primary productivity remained suppressed in burned regions.
The carbon loss over 2019-2020 is found to significantly exceed interannual variability in the regional carbon uptake over 2010-2019 from a set of top-down estimates (∼ 5σ anomaly), highlighting the extreme nature of this event. In the years to come, these ecosystems are expected to largely recover lost carbon stocks. However, the speed of recovery may be strongly regulated by climate variability and change, with reduced uptake during hot and dry conditions. This has important implications for the future carbon budget of the region. Climate-change-driven increases in the frequency of extreme heat and drought events will increase the recovery period after fires and decrease the inter-fire interval. If the recovery period becomes longer than the inter-fire interval then permanent carbon losses are likely.
This analysis finds that space-based remote sensing of trace gases and MODIS reflectances provide strong constraints on carbon cycle anomalies produced by extreme events. Still, there are remaining challenges that result in significant uncertainties in inferred fluxes. For inferring biomass burning estimates from X CO measurements, resolving pyroconvective tracer transport remains a major challenge and source of uncertainty. Aerosols co-emitted with biomass burning CO and CO 2 prevent total-column trace gas retrievals within much of the biomass burning plume. In addition, estimating CO 2 emissions from CO has considerable uncertainty, as does estimating the spatiotemporal structure of ΔNEE estimates. Addressing these sources of uncertainty, in addition to expanding spacebased trace gas observations will provide increasingly precise estimates of carbon release from extreme events.

Appendix A: Flux Inversion Configuration
The nested CO flux inversions are performed over a one-way nested domain of (100°−177.5° E and 0°−60° S) at 0.5° × 0.625° spatial resolution. Assimilated TROPOspheric Monitoring Instrument (TROPOMI) X CO super-obs are generated by aggregating measurements for each orbit with the quality flag ≥0.5 to the 0.5° × 0.625° spatial grid. The flux inversions optimize scaling factors to each model gridcell for prior biomass burning emissions from 5 November 2019 through 14 January 2020. Prior biomass burning emissions vary between flux inversions and are listed in Table 1. For the anthropogenic emissions, we combine off-line emission inventories from the Emission Database for Global Atmospheric Research (EDGAR 4.2) global model (Olivier & Berdowski, 2001) and several regional models including the US Environmental Protection Agency National Emission Inventory for 2008 in North America, the Criteria Air Contaminants inventory for Canada, the Big Bend Regional Aerosol and Visibility Observational Study Emissions Inventory for Mexico (Kuhns et al., 2003), the Cooperative Program for Monitoring and Evaluation of the Long-range Transmission of Air Pollutants in Europe inventory for Europe in 2000 (Vestreng, 2002) and the Streets Asia emissions inventory for 2000 (Streets et al., 2006). Monthly BioFuel emissions are from the EDGAR (Crippa et al., 2016), monthly shipping emissions from the International Comprehensive Ocean-Atmosphere Data Set (C. Wang et al., 2008), and hourly Biogenic emissions from Model of Emissions of Gases and Aerosols from Nature (Guenther et al., 2012).
Boundary conditions for the nested flux inversions are generated by performing a global inversion with GH-GF-Flux at 4° × 5° spatial resolution over the 3-month period from November 2019 through January 2020. The global inversion assimilates TROPOMI X CO super-obs (aggregated to 4° × 5° for measurements with quality flag equal to one) to optimize 14-day scaling factors for prior Global Fire Emissions Database version biomass burning emissions at each grid cell. Other prescribed emissions are identical to the nested flux inversion. Initial conditions for the global flux inversion are obtained from a global MOPITT X CO flux inversion. To test the sensitivity of inferred fluxes to the boundary conditions of the nested flux inversions, we generate a second set of boundary conditions that are identical to those from the global TROPOMI flux inversion but have CO increased by 10 ppb at all times and locations.

Conflict of Interest
The authors declare no conflicts of interest relevant to this study.

Data Availability Statement
The gridded daily estimates of ΔNEE and biomass burning emissions from this study can be downloaded from (https://doi.org/10.48577/jpl.rrk6-7453). GFED data were downloaded from https://www.globalfiredata.org/. GFAS data were downloaded from https://apps.ecmwf.int/datasets/. GFAS is generated using Copernicus Atmosphere Monitoring Service Information 2020, neither the European Commission nor ECMWF is responsible for any use that may be made of the information it contains. TCCON data were obtained from the TCCON Data Archive, hosted by CaltechDATA (https://tccondata.org). We downloaded version 10 of the ACOS OCO-2 lite files from the CO 2 Virtual Science Data Environment (https://co2.jpl.nasa.gov/). OCO-2 data were produced by the OCO-2 project at the Jet Propulsion Laboratory, California Institute of Technology, and obtained from the OCO-2 data archive maintained at the NASA Goddard Earth Science Data and Information Services Center. FluxSat data were downloaded from https://avdc.gsfc.nasa.gov/pub/tmp/FluxSat_GPP/. MODIS land cover data was downloaded from the EOSDIS Land Processes DAAC. ETOPO1 elevation data was downloaded from https://www. ngdc.noaa.gov. ERA5-Land data are obtained from the Climate Data Store (https://cds.climate.copernicus.eu). TROPOMI CO data were downloaded from http://www.tropomi.eu/data-products/carbon-monoxide. CrIS CO is provided by the NASA TRoposperhic Ozone and its Precursors from Earth System Sounding (TROPESS) and available from https://tes.jpl.nasa.gov. MODIS NIRv was calculated from MODIS NBAR measurements (MCD43A4), which were downloaded from the LP DAAC. TROPOMI L2 SIF data were downloaded from ftp:// fluo.gps.caltech.edu/data/tropomi/ungridded/SIF740nm. OCO-2 L2 SIF data are available from the GES DISC (https://disc.gsfc.nasa.gov). BHD CO/CO 2 and LAU CO data were downloaded from the World Data Centre for Greenhouse Gases (https://gaw.kishou.go.jp/). LAU CO 2 data were obtained from the site PIs. Paul Krummel, Ray Langenfelds, and Zoe Loh are thanked for supplying the CGO CO/CO 2 data. Australia fossil fuel combustion statistics were downloaded IEA's CO 2 Emissions from Fuel Combustion Highlights: https://www.iea.org/ data-and-statistics/data-product/co2-emissions-from-fuel-combustion-highlights.