Was Australia a sink or source of CO₂ in 2015? Data assimilation using OCO-2 satellite measurements

In this study, we present the assimilation of data from the Orbiting Carbon Observatory-2 (OCO-2) (land nadir and glint data, version 9) to estimate the Australian carbon surface fluxes for the year 2015. To perform this estimation, we used both a regional-scale atmospheric transport–dispersion model and a four-dimensional variational assimilation scheme. Our results suggest that Australia was a carbon sink of −0.41 ± 0.08 PgC yr⁻¹ compared to the prior estimate 0.09 ± 0.20 PgC yr⁻¹ (excluding fossil fuel emissions). Most of the carbon uptake occurred in northern Australia over the savanna ecotype and in the western region over areas with sparse vegetation. Analysis of the enhanced vegetation index (EVI) suggests that the majority of the carbon uptake over the savanna ecosystem was due to an increase of vegetation productivity (positive EVI anomalies) amplified by an anomalous increase of rainfall in summer. Further from this, a slight increase of carbon uptake in Western Australia over areas with sparse vegetation (the largest ecosystem in Australia) was noted due to increased land productivity in the area caused by positive rainfall anomalies. The stronger carbon uptake estimate in this ecosystem was partially due to the land surface model (CABLE-BIOS3) underestimating the gross primary productivity of the ecosystem. To evaluate the accuracy of our carbon flux estimates from OCO-2 retrievals, we compare our posterior concentration fields against the column-averaged carbon retrievals from the Total Carbon Column Observing Network (TCCON) and ground-based in situ monitoring sites located around our domain. The validation analysis against TCCON shows that our system is able to reduce bias mainly in the summer season. Comparison with surface in situ observations was less successful, particularly over oceanic monitoring sites that are strongly affected by oceanic fluxes and subject to less freedom by the inversion. For stations located far from the coast, the comparison with in situ data was more variable, suggesting difficulties matching the column-integrated and surface data by the inversion, most likely linked to model vertical transport. Comparison of our fluxes against the OCO-2 model intercomparison (MIP) was encouraging. The annual carbon uptake estimated by our inversion falls within the ensemble of the OCO-2 MIP global inversions and presents a similar seasonal pattern.

1 Introduction

Australia's carbon budget has been investigated by several large scale global ecosystem models (Sitch et al., 2015, Carbon cycle model intercomparison project (TRENDY)) and by the REgional Carbon Cycle Assessment and Processes project (RECCAP) (Haverd et al., 2013a). However, although they have contributed to a more refined knowledge of the Australian carbon cycle, the estimated fluxes still diverge significantly. In the latest RECCAP report (Haverd et al., 2015), the net biome production (NBP) estimate for the country was a net carbon source of 59 ± 35 TgC yr⁻¹ between 1990–2011. A large component of the uncertainty in this carbon budget was attributed to the estimate of net primary productivity (NPP) over grassland (Haverd et al., 2013b), with a large contribution to the land cover type they used to force their simulations (e.g. the Advanced Very High Resolution Radiometer (AVHRR) (1990–2006); Donohue et al., 2009) and the Moderate Resolution Imaging Spectroradiometer (MODIS) (2000–2011). Given this uncertainty, it is essential to bring any other observations we have to bear on the Australian carbon balance.

Data assimilation (also called atmospheric transport inversion), along with an increase of remotely sensed concentrations of carbon dioxide CO₂ data, has been revolutionary for quantifying land–ocean–atmosphere CO₂ flux exchange in the last decade. Satellite data from the Greenhouse Gases Observing Satellite (GOSAT) (Yokota et al., 2009), launched in 2009, and the Orbiting Carbon Observatory-2 (OCO-2) (Eldering et al., 2017), launched in 2014, have been used by several studies (Basu et al., 2013; Chevallier et al., 2014; Deng et al., 2014; Maksyutov et al., 2013; Crowell et al., 2019) to infer carbon CO₂ sources and sinks at continental scales. Few regional studies have been performed and “none in Australia”, while the global inversions show large differences for this region. For example, a study based on six satellite-based inversions using GOSAT (Chevallier et al., 2014, Fig. 1) shows that Australia was a carbon sink (−0.7 PgC yr⁻¹) for 2010. For the same year, Basu et al. (2013) inferred it to be a net carbon source (0.4 PgC yr⁻¹)¹.

The accuracy of flux inversions using global atmospheric transport models has been the subject of discussion due to errors related to modelled transport (Chevallier et al., 2014; Basu et al., 2018). Transport model error in global inversions often emerges because inversions are run at horizontal resolutions of 1–5^∘. Increasing the model resolution (Law et al., 2004) potentially reduces the representation errors found in global-scale models. Regional-scale inversions arose about a decade ago. They rely on mesoscale transport models (run at 1^∘ down to 10 km resolution); for example, Broquet et al. (2011) performed a regional-scale variational inversion of the European biogenic CO₂ fluxes at 50 km resolution. Another example of regional-scale inverse modelling is found in Villalobos et al. (2020), who performed an inversion at 81 km resolution over Australia. Finer-resolution models have the potential to be more successful, since they can offer a better representation of surface CO₂ fluxes and variability, as well as a better simulation of the processes driving high-frequency variability of transport (Schuh et al., 2010).

Australia has recently been subject to attention from the global carbon cycle community due to a large terrestrial carbon sink anomaly recorded in 2011 (Poulter et al., 2014). Poulter et al. (2014) found that Australia's flux anomaly was −0.66 for 2011 (relative to the 2003–2012 mean). Trudinger et al. (2016) also found a similar carbon sink anomaly for this period (ranging between −0.40 to −0.61 PgC yr⁻¹). These studies suggest that Australia’s ecosystems might act as strong sinks of CO₂ in the future during extreme wet periods. However, the efficiency and the spatial distribution of these carbon sinks remain largely uncertain (Ma et al., 2016). Some studies, i.e. Ma et al. (2016), found that the anomalous carbon uptake recorded in Australia in 2011 rapidly diminished in the following year (∼ 0.08 PgC yr⁻¹), suggesting that semi-arid ecosystem can act as carbon sink in the short term but not over longer periods compared to tropical forest ecotypes. An important unanswered question in carbon cycle research remains regarding the carbon sink strength of semi-arid ecosystems in non-wet years.

In this study, we present a regional inversion to infer CO₂ fluxes over Australia for 2015 based on the Community Multiscale Air Quality (CMAQ) model and OCO-2 satellite retrievals. In 2015, Australia was affected by the El Niño–Southern Oscillation (ENSO), and although some parts of the continent were impacted by rainfall deficiency, in other regions such as northern and southeastern Australia rainfall was above average (Annual climate statement, Bureau of Meteorology, 2015).

This paper is structured into five sections. Section 2 describes the flux inversions system and the datasets used. Section 3 presents the main results of the Australian carbon budget, as well as an analysis of the enhanced vegetation index (EVI) and rainfall anomalies, and a comparison between our posterior CO₂ concentration and the Total Carbon Column Observing Network (TCCON) and in situ measurements. In Sect. 4, we present a discussion of our results, as well as a comparison of our optimized fluxes against the ensemble mean of nine different global inversions that participate in the OCO-2 model intercomparison (MIP). In Sect. 5, we summarize our findings.

2 Methodology and data

To estimate the Australian CO₂ surface fluxes for 2015, we followed the same four-dimensional variational assimilation scheme described in Villalobos et al. (2020). In this section, we will present a brief description of the system and an update of all changes we made to the data used for our inversion.

2.1 Bayesian inverse system

Finding the optimal value (x^a) of the CO₂ flux estimates involves identification of the best fits between both observations (y) and a prior (or background) estimate (x^b) of these fluxes (Ciais et al., 2010; Rayner et al., 2019). Using Bayes' theorem and under the hypothesis of unbiased Gaussian-distributed errors of x^b and y, the best estimate of x^a (likelihood maximum a posteriori) is equivalent to finding the minimum of the cost function J(x) shown in Eq. (1). Notation in this study follows Rayner et al. (2019).

$$\begin{array}{}\text{(1)}& \begin{array}{rl}J\left(\mathit{x}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[(\mathit{x}-{\mathit{x}}^b{)}^T{\mathbf{B}}^{-\mathrm{1}}(\mathit{x}-{\mathit{x}}^b)\right]\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\left(\mathbf{H}\right(\mathit{x})-\mathit{y}{)}^T{\mathbf{R}}^{-\mathrm{1}}(\mathbf{H}\left(\mathit{x}\right)-\mathit{y})\right]\end{array}\end{array}$$

In Eq. (1), H, represents the application of the forward model and the “observation operator”, which allows us to map the model variables (e.g. fluxes) to observations. R represents the error covariance matrix of the observations y, including the transport model error. R is defined as a diagonal matrix (details Sect. 2.5). x represents the control vector of unknowns. x includes not only CO₂ surface fluxes but also initial and boundary conditions (details Sect. 2.2). B is the associated error covariance matrix of x^b, boundary and initial concentrations, and includes off-diagonal terms. In these off-diagonal values, we only include spatial and non-temporal correlations of the prior fluxes (details of the structure of the prior error covariance matrix can be found in Sect. 2.4 in Villalobos et al., 2020).

We calculate the minimum of J(x) by an iterative process and not by an analytical expression. This numerical problem requires the value of the cost function gradient ∇_xJ(x).

$$\begin{array}{}\text{(2)}& {\mathrm{\nabla }}_{x}J={\mathbf{B}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}^b\right)+{\mathbf{H}}^T\left({\mathbf{R}}^{-\mathrm{1}}\left[\mathbf{H}\left(\mathit{x}\right)-\left(\mathit{y}\right)\right]\right)\end{array}$$

We compute H^T using the adjoint of the CMAQ model (version 4.5.1; Hakami et al., 2007). We can see in Eq. (2) that H^T is applied to the vector ${\mathbf{R}}^{-\mathrm{1}}\left(\mathbf{H}\left(\mathit{x}\right)-\mathit{y}\right)$, which is often called the “adjoint forcing”, and represents the error-weighted differences between the forward model and the observed concentrations. Applying the adjoint model to the adjoint forcing, running backward in time from t_i−1 to t₀, allows us to construct the gradient of the cost function, ∇_xJ(x). The algorithm that our inverse system uses to optimize the J(x) is the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS-B), implemented in the SciPy Python module (Byrd et al., 1995). Figure 1 shows a simplified version of how our inversion system works to find the optimal values of CO₂ surface fluxes.

The error statistics of x^a are embodied in the posterior error covariance matrix (A). In this study, A was computed by a series of observing system simulation experiments (OSSEs) carried out by Villalobos et al. (2020, Sect. 2.4). However, here we adjusted the prior and observation uncertainties assumed in Villalobos et al. (2020) by a factor of 1.2. We made this adjustment to satisfy the theoretical assumption in the variational optimization, which indicates the value of the cost function in its minimum has to be approximately equal to half of the number of observations (for more details, see Sect. 3.1). In general, errors assumed in the inversion are not Gaussian and independent but rather have errors correlated in time and space (including flat biases) that render the statistical assumptions made in deriving the estimation method invalid and lead to a higher cost function than expected. A description of how the prior and observation uncertainties were assumed in our study is found in Sect. 2.3 and 2.5. Appendix D (Figs. D1 and D2) shows the spatial distribution of the prior and posterior that we reference in this study.

Figure 1 A simplified diagram of the four-dimensional variational data assimilation we used to estimate CO₂ surface fluxes over Australia.

2.2 Defining the control vector

In our data assimilation system, we solve for monthly average surface fluxes at 81 km grid-cell-scale resolution as the multipliers of the principal eigenvectors of the prior error covariance matrix B, computed as ${\mathbf{W}}^T{\mathbf{w}}^{-\mathrm{1}/\mathrm{2}}(\mathit{x}-{\mathit{x}}^b)$, where W was defined as the matrix of eigenvectors and w as a diagonal matrix of corresponding eigenvalues (Villalobos et al., 2020, Sect. 2.2). In order to avoid the impact of the initial conditions (ICs) and boundary conditions (BCs) on our assimilated fluxes, we also solved them within the control vector x. We did not optimize them in the same way as the fluxes in order to not increase the control vector size, so we treat the unknowns related to the BCs and ICs as scaling factors of the emissions added to the CMAQ model. Lateral BCs were solved as eight boundary regions divided by the upper and lower boundary areas within the CMAQ domain (south, east, north and west). In Fig. 2, we provide a representation of these boundaries. In this figure, we can see that our study domain not only covers the Australian continent (AUS) but additionally other countries such as Indonesia (IND), Papua New Guinea (PNG) and New Zealand (NZ). The extension of this domain was created as an extra precaution to minimize the influence of the boundaries over Australia.

Figure 2 Representation of the horizontal Weather Research and Forecast model (WRF) domain (black rectangle) and CMAQ model domain (dark blue rectangle). Boundary components (south, east, north and west) are represented as between the outer domain of CMAQ domain and the dotted dark blue lines. Land biosphere emissions incorporated over Australia are represented by the small dotted blue lines (CABLE model in BIOS-3 setup). Outside this area, land biosphere emissions come from the CABLE global product.

Lower boundary layers were defined to cover from σ=1 to σ=0.96, which correspond (on average) to a pressure of 972.5 hPa, while the upper boundary layer was defined to cover from 972.5 up to 50 hPa. As mentioned before, our inversion system solves for these lateral boundary components, while surface fluxes are also being optimized. Boundary conditions are provided to the CMAQ model as daily averages, but we optimize them as monthly averages. BC and IC datasets were taken from the CAMS global CO₂ atmospheric inversion product (version v19r1) (Frédéric Chevallier, personal communication, 2019). Uncertainties for the initial condition were set at 1 % (approximately 4 ppm), and the uncertainties in the lateral boundary conditions were assumed as the standard deviation (1σ uncertainty) of CAMS concentration data in the perimeter of the boundary.

2.3 Prior information and its uncertainties

We updated the prior CO₂ fluxes described in Villalobos et al. (2020, Sect. 2.4). Biosphere carbon fluxes were derived using a modified version of the Community Atmosphere-Biosphere Land Exchange model (CABLE) (Haverd et al., 2018), which was forced by Australian regional drivers and observations (BIOS3 setup). CABLE consists of a biophysical core: the Carnegie–Ames–Stanford approach, carbon, nitrogen, phosphorus (CASA-CNP) biogeochemical model (Wang et al., 2010), the POP module for woody demography and disturbance-mediated landscape heterogeneity (Haverd et al., 2013d), and a module for land use and land management (POPLUC; Haverd et al., 2018). For our case, Vanessa Elizabeth Haverd (personal communication, 2020) ran the CABLE model in the BIOS3 setup (hereafter CABLE-BIOS3) at a resolution of 0.25^∘. We calculated 3-hourly biosphere CO₂ fluxes by combining two datasets: daily net ecosystem exchange (NEE) fluxes with 3-hourly gross primary production (GPP). Given that the BIOS3 product did not cover our whole CMAQ model domain, we also incorporated monthly biosphere fluxes from CABLE-POP global simulations as shown in Fig. 2. These CABLE-POP simulations were used in the carbon cycle model intercomparison project (TRENDY version 8) for the 2019 global carbon budget (Friedlingstein et al., 2019). Biosphere flux uncertainties in our system were assumed to be equal to the NPP simulated by CABLE, with a ceiling of 3 gC m⁻² d⁻¹ following Chevallier et al. (2010).

Anthropogenic fluxes were created by the combination of two different inventory datasets: the Open-source Data Inventory for Anthropogenic CO₂ (ODIAC) (Oda et al., 2018) and the Emissions Database for Global Atmospheric Research (EDGAR) version 5 (Crippa et al., 2020). The combination of these two anthropogenic inventories (each used to cover different source sectors) was necessary because the version of the ODIAC selected did not contain emissions from aviation. The EDGAR emissions combined with ODIAC were aviation climbing and descent, aviation cruise, and aviation landing and take-off datasets. Aviation emissions were also distributed across the vertical layers of the CMAQ domain. EDGAR is a gridded product with spatial resolution of $\mathrm{0.1}{}^{\circ }\times \mathrm{0.1}{}^{\circ }$ with monthly temporal resolution. ODIAC (version 2019) is also a gridded product, which has a spatial resolution of 1×1 km. Monthly ODIAC fluxes were modified by incorporating a diurnal-scale factor estimated by Nassar et al. (2013). The ODIAC data product selected did not include bunker emissions. Fossil fuel carbon emission uncertainties were created by multiplying the emissions dataset by a factor of 0.44. This factor was calculated by a linear regression between the mean fluxes and the spread of an ensemble of 25 realizations of posterior emissions estimated by the Fossil Fuel Data Assimilation System (FFDAS) (Asefi-Najafabady et al., 2014).

Prior ocean fluxes were taken from the CAMS greenhouse gas flux inversion (version v19r1) (Frédéric Chevallier, personal communication, 2019). The prior fluxes that CAMS uses in its inversion also includes EDGAR emissions over the ocean; thus, we did not include this anthropogenic flux over the ocean to avoid double counting. We assumed that the ocean uncertainties were uniform and set up a value of 0.2 gC m⁻² d⁻¹ over ocean, as in Chevallier et al. (2010). We also used monthly fire emissions from the Global Fire Emission Database (GFED) v4.1 (van der Werf et al., 2017), which includes small fire emissions. Fire emission uncertainties were assumed as 20 % of the biomass burning carbon emissions. All datasets mentioned above (terrestrial biosphere exchange, fossil fuel, fires and ocean fluxes) were interpolated to the spatial resolution of the CMAQ model.

As described in Villalobos et al. (2020, Sect. 2.4), we included spatial correlations into our prior error covariance matrix B following Basu et al. (2013, Sect. 3.1.1). The correlation length between grid points over land was assumed to be 500 km and over ocean 1000 km. We assume that fossil fuel uncertainties were not correlated, so we only use the diagonal values of the matrix. In our eigen-decomposition of B, the eigen-spectrum (eigenvectors of the covariance matrix) retains 99 % of the explained variance (eigenvalues).

2.4 Atmospheric transport model

The inversion was based around the CMAQ modelling system (version 5.3) and its adjoint (version 4.5.1; Hakami et al., 2007). The CMAQ modelling system is an Eulerian (gridded) mesoscale chemical transport model (CTM). We added CO₂ into the CMAQ model as an inert chemical species, whose concentration is determined by atmospheric transport, fluxes, initial and boundary concentrations. The CMAQ model was driven by meteorological fields from the Weather Research and Forecast model (WRF) Advance Research Dynamical Core WRF-ARW (henceforth WRF) version 4.1.1 (Skamarock et al., 2008), the data of which were processed by the Meteorology-Chemistry Interface Processor (MCIP) version 4.2 (Otte and Pleim, 2010). WRF configuration details are shown in Table 1. Our WRF model was set up at a spatial resolution of 81 km with 32 vertical layers from the surface up to 50 hPa. The numerical simulation was carried out on a single domain (i.e. non-nested). WRF initial conditions were taken from the ERA-Interim global atmospheric reanalysis (Dee et al., 2011), which has a resolution of approximately 80 km on 60 vertical levels from the surface up to 0.1 hPa. Sea surface temperatures were obtained from the National Centers for Environmental Prediction/Marine Modeling and Analysis Branch (NCEP/MMAB). The WRF model was run with a spin-up period of 12 h.

(Morrison et al., 2009)(Iacono et al., 2008)(Iacono et al., 2008)(Monin and Obukhov, 1954)(Tewari et al., 2007)(Janjić, 1994)(Grell and Dévényi, 2002)

Table 1 Physics parameterizations used in WRF model setup.

2.5 OCO-2 satellite information and its uncertainties

We assimilated satellite data from OCO-2 level 2 (lite file version 9) for 2015, which is distributed by the National Aeronautics and Space Administration (NASA) (OCO-2 Science Team/Michael Gunson, Annmarie Eldering, 2018). OCO-2 was launched in 2014 and since then has provided nearly global coverage of column-averaged dry air mole fraction of CO₂ that has been used by several carbon cycle researchers to estimate surface carbon fluxes at global and regional scales. OCO-2 provides data in three modes: nadir, glint and target mode. In nadir mode, OCO-2 instrument points straight down at the surface of the Earth (surface solar zenith angle is less than 85^∘); in glint mode, OCO-2 instrument points to the bright glint spot on Earth where solar radiation is directly reflected off the Earth's surface (local solar zenith angle is less than 75^∘); and target mode, the instrument is configured to scan about a particular point on the ground as it passes overhead. In this study, we used the combination of both land nadir and land glint observations (LNLG), because there are no systematic offsets between the two datasets (O'Dell et al., 2018). We also performed an inversion using the combination of land (nadir and glint) and ocean glint observations (LNLGOG). However, this inversion was treated as an independent experiment (see Appendix F, Table F1), and the assimilated fluxes estimated by using LNLGOG were not included in our main results. We decided to not incorporate them because ocean glint retrievals still have undetermined biases (O'Dell et al., 2018) that may complicate or confound the Australia flux estimation. We discussed the impact of assimilation LNLGOG in the validation of our inversion with independent data (see Sect. 3.5 for more details). We only used OCO-2 retrievals with quality flag 0 and only soundings that were bias corrected, as described by Kiel et al. (2019). The spatial distributions of OCO-2 soundings (LNLG and LNLGOG) across the CMAQ domain for 2015 are shown in Appendix C, Figs. C1 and C2, respectively.

Given that multiple OCO-2 soundings cross one grid cell over the CMAQ domain, we had to average them before doing any comparison with the CMAQ model simulations. This averaging process was carried out in two steps. First, we averaged all OCO-2 soundings that fall within 1 s intervals, and then these 1 s averages were averaged again within the CMAQ vertical columns (approximately 11 s average) across 81 km × 81 km grid-cells. The 1 s weighted averaging process is described in detail in Villalobos et al. (2020, Sect. 2.3). In summary, to obtain the uncertainties of these 1 s averaging processes, we considered three different forms of uncertainty calculation, similar to Crowell et al. (2019). First, we averaged OCO-2 uncertainties assuming that these were correlated in a 1 s span (uncertainties defined as σ_s). Second, given that the average of OCO-2 uncertainties is sometimes low because they neglect systematic errors, we also used the spread of the OCO-2 retrievals in the 1 s average (uncertainties defined as σ_r). Third, we also defined baseline uncertainties (defined as σ_b) for cases where the number of soundings was not enough to compute a realistic spread. The values for our baseline uncertainties were assumed to be 0.8 ppm over land and 0.5 ppm over ocean. Finally, we selected the maximum value between these three uncertainties (σ_s, σ_r and σ_b). We also added (in quadrature) to this term 0.5 ppm as the contribution of the model uncertainty (defined as σ_m).

Solving the cost function shown in Eq. (1) requires convolving the vertical levels of the CMAQ model with the retrieval profile from OCO-2. For this, we used Eq. (3) derived by Connor et al. (2008) as follows:

$$\begin{array}{}\text{(3)}& x_{{\mathrm{CO}}_{\mathrm{2}}}^m=x_{{\mathrm{CO}}_{\mathrm{2}}}^a-\sum _j{\mathit{h}}_j{\mathit{a}}_{{\mathrm{CO}}_{\mathrm{2}},j}{\mathit{x}}_a+\sum _j{\mathit{h}}_j{\mathit{a}}_{{\mathrm{CO}}_{\mathrm{2}},j}{\mathit{x}}_j^m,\end{array}$$

where x^a is the OCO-2 a priori, h is a vector of pressure weight, h_j is the mass of dry air in layer j divided by the mass of dry air in the total column, ${\mathit{a}}_{{\mathrm{CO}}_{\mathrm{2}}}$ is the averaging kernel of OCO-2, x_a is the OCO-2 a priori profile, and x^m is the simulated profile from the CMAQ model. In our inversion system, the OCO-2 averaging kernel is defined on 20 pressure levels and we interpolate these to the CMAQ vertical levels.

2.5.1 TCCON measurements

To validate our posterior CO₂ CMAQ concentrations, we used ground-based remote-sensing data from the Total Carbon Column Observing Network (TCCON) (Wunch et al., 2011). There are three TCCON stations in our domain (see Table 2 for references and Fig. 3 for coordinate locations). A TCCON instrument is a Fourier transform spectrometer (FTS) developed to record direct solar spectra in the near-infrared spectral region. TCCON provides accurate and precise column-averaged concentrations of CO₂ and other greenhouse gases. This instrument represents the “gold standard” for surface-based remote-sensing estimates of the total-column concentration of these gases. Data from TCCON are widely used by carbon cycle researchers, in particular for global flux inversion (Byrne et al., 2020) and validation of satellite data products (such as from OCO-2). To perform a quantitative comparison against CMAQ simulations, we averaged all the TCCON retrievals to create hourly average XCO2 values, which were consistent with the CMAQ hourly simulations. After calculating the average of these retrievals, we interpolated the TCCON column averaging kernels and TCCON a priori CO₂ profile to the CMAQ vertical levels. After the interpolation, we followed Eq. (3) to compute the TCCON column-mixing ratios simulated by CMAQ. The statistical analysis of CMAQ model–TCCON differences was based on monthly mean concentration, which were calculated by taking local time averages (10:00–14:00 LT, Australian local time, locally referred to the monitoring site location), where the solar radiation intensity is most stable (Kawasaki et al., 2012).

Griffith et al. (2014a)Griffith et al. (2014b)Sherlock et al. (2014)

Table 2 Reference of the TCCON stations used in this work for evaluation of our inverse model system.

2.5.2 Ground-based in situ measurements

Additional datasets used to validate our posterior concentrations were taken from four ground-based in situ monitoring sites forming part of the Global Atmosphere Watch (GAW) Programme of the World Meteorological Organisation (WMO): Cape Grim, Gunn Point, Burncluith and Ironbark. Coordinates of these locations are shown in Fig. 3. All these datasets were supplied by Zoë M. Loh (personal communication, 2019) at hourly temporal resolution. For the comparison with our model simulation, we used hourly data from these monitoring sites, but the monthly mean averaged data shown in Sect. 3.5.2 were calculated using local time averages between midday and 17:00 LT (Australian local time, locally referred to the monitoring site location).

Measurements of atmospheric CO₂ concentration at the Gunn Point, Ironbark and Burncluith sites were made continuously at high frequency (∼ 0.3 Hz) using CSIRO Picarro cavity ring-down spectrometers (model G2301 at Gunn Point and Ironbark, and G2401 at Burncluith) all with inlets placed at the height of 10 m. Details of the Ironbark and Burncluith installation are given by Etheridge et al. (2016) and are broadly similar to the installations elsewhere, including Gunn Point. Cape Grim also operates a Picarro G2301 analyser, with the inlet positioned at a height of 70 m.

The instrumental precision for these analysers is better than ± 0.1 ppm for CO₂ (Etheridge et al., 2014), and all measurements are calibrated to the WMO X2007 CO₂ mole fraction scale (Zhao and Tans, 2006), ensuring comparability between all measurements used.

Cape Grim is a significant monitoring station in the GAW programme because it samples air with some of the least recent anthropogenic and terrestrial influence in the world, representing hemispheric background concentrations. These air masses, known as “baseline”, have blown straight off the Southern Ocean and have often been used in modelling studies. However, in this study, we used all Cape Grim data because our inversion assimilates only data that were collected over land and carry terrestrial signals.

Figure 3 TCCON and in situ location sites. Red dots indicate TCCON locations. TCCON Darwin and Wollongong are located over Australia, while TCCON Lauder is located in New Zealand. Blue dots represent in situ location around Australia (Gunn Point, Burncluith, Ironbark and Cape Grim).

2.6 Auxiliary data

In this study, we also use auxiliary data such as the EVI, rainfall and GPP from the CABLE-BIOS3 model to understand the difference between the prior and posterior fluxes over Australia in 2015.

2.6.1 The EVI

To understand if there was higher than expected growth of vegetation across Australia in 2015, we evaluated the monthly EVI anomalies relative to the long-term mean from 2000–2014. We used the EVI from the MODIS MOD13C1 version 6 data product from the NASA satellite Terra (Didan, 2014). This gridded EVI MODIS product has a temporal resolution of 16 d composite and 0.05^∘ spatial resolution. We constructed the EVI anomalies by subtracting the long-term mean (2000–2014) for each month of 2015. The spatial distribution of the EVI anomalies is shown in Fig. S1 in the Supplement. EVI measures the greenness of vegetation and can be used as a proxy for monitoring the density or productivity of the vegetation biomass. EVI indices range from −0.2 to 1, where values less than 0 indicate a lack of green vegetation or arid areas. These monthly EVI MODIS products were regridded to the CMAQ grid to calculate the spatial correlation between prior and posterior flux differences (see Sect. 3.3).

2.6.2 Australian Water Availability Project (AWAP)

Monthly rainfall data were taken from the Australian Water Availability Project (AWAP), Bureau of Meteorology (Jones et al., 2009). We used data for the period 2000–2015. AWAP data are obtained from a spline interpolation technique, which interpolates all available in situ rainfall observations onto grid-cells of 0.05^∘ (more details can be found in Jones et al., 2009). AWAP rainfall anomalies were calculated in the same way as EVI anomalies, by subtracting their long-term mean from 2000 to 2014 (see Fig. S2 in the Supplement).

Table 3 Summary of the configuration of the OCO-2 MIP (version 9) design.

2.6.3 MODIS GPP

We compared the MODIS Terra GPP MOD17A2H version 6 product for 2015 (Running et al., 2015) against CABLE-BIOS3 model GPP predictions (see Appendix E, Fig. E1). The MODIS GPP product has a spatial resolution of 500 m and a temporal resolution of eight days. The 8 d composite was averaged to monthly resolution and aggregated to the CMAQ grid for comparison with the CABLE-BIOS3 model GPP.

Figure 4 Bias and root mean square error (RMSE) between OCO-2 and the prior and posterior concentrations simulated by CMAQ model. Orange and purple circles represent prior and posterior concentration biases, and orange and blue bars represent the RMSE (Units: ppm). The top edge of the box represents the 75th percentile and the bottom edge represents the 25th percentile. The top and bottom whiskers represent the 95th and 5th percentiles.

2.6.4 MIP in situ and OCO-2 satellite-derived fluxes

For validation, we compared our posterior Australian biosphere CO₂ flux estimates (excluding fossil fuel) against the ensemble monthly mean of nine OCO-2 satellite-based and in situ global inversions (see Sect. 4 for details). In situ and OCO-2 satellite-derived fluxes were consolidated by the OCO-2 MIP (Crowell et al., 2019), which used the OCO-2 satellite version 7 product. In this study, we used the latest update of the OCO-2 MIP product (Peiro et al., 2021), which uses OCO-2 data lite file version 9, with an improved bias correction approach (Kiel et al., 2019) compared to the version 7 product (Crowell et al., 2019). Within the MIP design, in situ carbon flux estimates are derived by utilizing five collections in ObsPack observations (Masarie et al., 2014). A description of these data can be found in Peiro et al. (2021, Sect. 2.3).

Table 3 shows a summary of these global inversions. This table shows that MIP global inversions were performed using different prior flux estimates, and the transport models were run at different spatial scales. Within MIP, prior estimates also include fossil fuel data, which was fixed and derived from ODIAC. With regard to fires estimates, they use different versions of the GFED dataset. Some of them used version 4, while other modellers use version 3. The main difference between these two datasets is that GFED version 3 does not include small fire-burned areas.

3 Results

3.1 Inversion evaluation: analysis of the residual between CMAQ simulation and OCO-2

As described in Eq. (1), the main purpose of the inversion is to optimize fluxes by minimizing the mismatch between the model simulation and observations. In order to evaluate the performance of the inversion, we compared the CO₂ concentrations obtained when forcing the CMAQ model with the prior and posterior fluxes (for convenience, we will call these the prior and posterior CO₂ concentrations, respectively). Figure 4 shows the bias and root mean square error (RMSE) between the prior and posterior CMAQ simulations against the OCO-2 observations for 2015. This figure shows that the biases and RMSE in the posterior concentration were reduced by the inversion and indicate the inversion system leads to an overall improvement of the representation of OCO-2 observations. Our findings indicate that the prior concentrations overestimate OCO-2 from March to April and from July to September. Prior biases in these months were reduced by more than 90 %. In March, for example, the monthly mean bias was reduced from 0.56 to 0.05 ppm, with a decrease in the RMSE from 1.11 to 0.84 ppm. In April, we see similar results to those for March; the prior bias was reduced from 0.40 (RMSE = 1.05) to 0.03 (RMSE = 0.88). On the other hand, in January, February, May and December, prior biases were negligible, showing a good agreement with OCO-2. In a consistent system, we know that the theoretical value of the cost function at its minimum should be close to half the number of assimilated observations, assuming all error statistics are correctly specified (Tarantola, 1987, p. 211). In our inversion, after iteration 27, we obtained a cost function of 4392.15, which was close to half of the total number of OCO-2 assimilated observations for 2015 (N=9556). In general, we also see a modest reduction in the prior RMSE each month during 2015, and its variability is proportional to the number of assimilated observations. Thus, a slight prior RMSE decrease corresponds to a month with a reduced number of OCO-2 data available.

3.2 Australian carbon flux estimate

In this section, we will only discuss the results of carbon fluxes that were assimilated by the combination of LNLG OCO-2 retrievals, not the carbon fluxes estimated using LNLGOG observations. We decided not to discuss the results based on LNLGOG because ocean glint observations (version 9) still have undetermined biases (O'Dell et al., 2018) that might contaminate the Australian carbon flux estimate. However, we include these findings in the Appendix F, Table F1. Adding ocean OCO-2 glint observations to our inversion system does not significantly alter the terrestrial annual mean flux estimate for Australia (−0.36 PgC yr⁻¹) compared to estimates made by only using OCO-2 LNLG observations (−0.41 PgC yr⁻¹). Results based on LNLGOG will be further discussed in Sect. 4.

Figure 5a represents the terrestrial prior and posterior annual mean flux for 2015 (excluding fossil fuel). As mentioned previously, our assimilated carbon fluxes using LNLG indicate that the Australian annual terrestrial flux for 2015 was a slight carbon sink of −0.41 ±  0.08 PgC yr⁻¹ (1σ uncertainty) compared to the prior terrestrial estimate of 0.09 ± 0.17 PgC yr⁻¹. Our prior fossil fuel estimate from ODIAC and EDGAR, which is about 0.06 ± 0.01 PgC yr⁻¹ (mostly constant for each month in 2015) over Australia represents only 25 % of the annual posterior flux. We decided to exclude these emissions from our analysis because variations in land uptake cause most of the variation in our posterior fluxes.

Figure 5b shows the seasonal cycle of the prior and posterior fluxes along with its uncertainties. As mentioned in Sect. 2.3, the prior and posterior uncertainties included in Fig. 5a and b were calculated from an ensemble of five different OSSEs adjusted by a factor of 1.2 in this study. We also plotted the spatial distribution of the prior and posterior fluxes (Figs. 6 and 7), and the difference between them (Appendix A, Fig. A1).

Figure 5b shows that the posterior flux estimates generally refine the prior with the exception of March and the period July to September. In January and February, the posterior fluxes were not modified much by the inversion. In January, for example, the terrestrial posterior flux was −0.84 ± 0.18 PgC yr⁻¹ compared to the prior −0.89 ± 0.75 PgC yr⁻¹. The agreement follows from the small residual between prior simulated concentration and observation (Fig. 4). From April to May, we also see the posterior is shifted from the prior, although not significantly considering the prior uncertainty. In April, for instance, the prior flux (0.24 ± 0.69 PgC yr⁻¹) was slightly shifted to a posterior carbon sink (−0.28 ± 0.25 PgC yr⁻¹). However, these two estimates do not disagree because they fall within 1σ uncertainties.

Figure 5 Time series of monthly mean prior (orange dots) and posterior (blue dots) CO₂ fluxes and their uncertainties in PgC yr⁻¹ over Australia for 2015. The dashed orange and blue line represents a smooth line for the prior and posterior fluxes, respectively.

As mentioned in the previous paragraph, March, July, August and September were the exceptions to this general agreement. In March, we see a prior flux of 0.12 ± 0.74 PgC yr⁻¹ compared to the posterior carbon sink of −0.82 ± 0.17 PgC yr⁻¹. The difference between the posterior (Fig. 7b) and the prior flux (Fig. 6b) at grid-cell scale (see Appendix A, Fig. A1c) suggests that most of the posterior sink comes from the north and southeast corner of Australia. July represents the month where the posterior is most shifted from the prior. In this month, we see a posterior flux of −1.75 ± 0.34 PgC yr⁻¹ compared to the prior flux of 0.09 ± 0.62 PgC yr⁻¹. The spatial distribution of the posterior and prior flux difference at grid-cell scale for July (Fig. A1g) indicates that the shift largely comes from northern and southeastern Australia. The stronger posterior sink seen in July decreased in August (−0.93 ± 0.27 PgC yr⁻¹) and September (−0.78 ± 0.20 PgC yr⁻¹), and changed sign in October and November. In November, the posterior flux was 1.75 ± 0.31 PgC yr⁻¹ compared to the prior, which was 0.53 ±  0.58 PgC yr⁻¹. The largest difference in this month is found in the north and on the southeast coast of Australia (Appendix A, Fig. A1k). The carbon release from land in northern Australia is likely attributed to a combination of fire anomalies (Fig. S3k in the Supplement) and the lack of rainfall seen in Australia in 2015 (Fig. S2k in the Supplement). In December, we see that the posterior source seen in November changed to a posterior carbon neutral (0.003 ± 0.15 PgC yr⁻¹). A further analysis which explains the reasons for this shift is given in the following section.

Figure 6 Prior fluxes derived by the CABLE model in the BIOS3 setup in combination with fires emissions selected by GFED for 2015 (fossil fuel emissions are excluded).

Figure 7 Posterior fluxes assimilated using LNLG OCO-2 satellite observations for 2015 (fossil fuel emissions are excluded).

3.3 Spatial patterns of the EVI and rainfall anomalies in Australia

To investigate why our inversion led to a higher carbon uptake (relative to the prior flux) in some months in 2015, we studied the spatial pattern of monthly EVI anomalies and rainfall anomalies. EVI anomalies were calculated relative to 2000–2014 over Australia from the MOD13C1 version 6 data product and rainfall anomalies from AWAP data relative to 2000–2014 (Figs. S1–S2 in the Supplement). We indicated in Sect. 3.2 that the posterior sink recorded in January and February agrees with the prior estimates but disagrees with the prior flux estimates from March to September, where the most considerable difference is seen in March and from July to September. From the inversion viewpoint, the significant shift between prior and posterior fluxes occurs because the prior column average simulated by the CMAQ model overestimates the column-average retrieval by OCO-2 in these periods (see Fig. 4).

As indicated in Sect. 3.2, most of the posterior carbon uptake seen in March comes from the northern part of Australia (except coastal regions) and the southeastern area (Appendix A; Fig. A1c). We found that the higher posterior uptake relative to the prior in the northern part of the continent was not attributed to an increase of greenness in vegetation (negative EVI anomalies). These results suggest that most of the posterior carbon uptake observed in March is likely associated with positive EVI anomalies seen in January and February affected by the positive rainfall anomalies recorded in January. High anomalous rainfall in January is not unexpected because it is the wet season in the northern region of Australia (tropical monsoonal climate).

The spatial pattern of the difference between the posterior and prior flux estimate recorded in July indicates that the majority of the posterior carbon uptake estimated by the inversion comes from the southeastern and northern region of Australia (Appendix A; Fig. A1g). We found that the posterior sink estimated in southeastern Australia was likely driven by a higher-than-expected greenness of vegetation (Fig. S1g in the Supplement), probably induced by anomalously positive rainfall in that period (Fig. S2g in the Supplement). We cannot conclude the same results in the northern region that we found in the south. We see in Appendix A; Fig. A1g that positive EVI anomalies were slight compared to the one found in southeastern Australia. In the following section, we will show that the underestimation of the GPP by the CABLE-BIOS3 model might be the likely reason for the difference between prior and posterior in this region.

An increase in carbon uptake estimated by our inversion in August comes from the northern and southern regions of Australia (with the exception of coastal areas in the southeastern corner of Australia), which mainly shows a release of carbon (Appendix A; Fig. A1h). The release of carbon by the land in this coastal region is likely attributed to a decrease in land productivity (Fig. S1h in the Supplement). The subtle decrease of photosynthesis activity in the coastal area is likely associated with a decrease of rainfall seen in June and July (Fig. S2f and g in the Supplement).

In September, the posterior carbon uptake primarily comes from the southeast corner of Australia (with a slight exception seen in the southeast and east coast of Australia), which shows a release of carbon into the atmosphere (Appendix A; Fig. A1i). The carbon uptake seen in the southeast of Australia aligns with a higher-than-usual increase in land productivity, as reflected by the positive EVI anomalies in that region (Fig. S1h in the Supplement), likely benefited by the positive rainfall anomalies seen in August in that area. In September, we also see that positive EVI anomalies were not as strong as in July and August. These findings are probably associated with the fact that rainfall in September decreased considerably for most parts of the country, where rainfall was lower than average (negative rainfall anomalies) (Fig. S2i in the Supplement).

Anomalies in EVI in Australia are closely related to fluctuations in rainfall, which is one of the most important drivers of ecosystem dynamics and productivity. This is the case in (e.g. semi-arid) regions where rainfall is the limiting factor for plant growth, which is indeed the case in much of Australia. These results are consistent with findings of previous studies (e.g. Weltzin et al., 2003) that Australia's semi-arid ecosystems are water resilient and can respond to favourable rainfall conditions by capturing large amounts of carbon.

3.4 Australian carbon flux estimate classified by bioclimatic zones

To understand which Australian ecosystem contributed most to our posterior carbon sink estimate, we divided the continent into six bioclimatic classes: tropical, savanna, warm temperate, cool temperate, Mediterranean and sparsely vegetated (Fig. 8). We used the same six bioclimatic regions at a 0.05^∘ spatial resolution as in Haverd et al. (2013a). The six bioclimatic classes used in this study correspond to an aggregation of the 18 agroclimatic zones generated by Hutchinson et al. (2005). The climatic classification in Hutchinson et al. (2005) was adapted from an existing global agroclimate classification (Hutchinson, 1992), which was refined and closely aligned with natural vegetation formations and common land uses across Australia using 182 weather climate stations and the Interim Biogeographic Regionalisation for Australia (IBRA). In Fig. 8, we can see that Australian tropical land only covers the northern coastal part of Australia. Savanna extends across the northern tropics to the southeastern subtropical zone. Warm temperate land covers the southeast Australian coast, while cool temperate land covers the southeastern corner of Australia. The Mediterranean region is confined to the southwestern corner of Australia and the gulf region of South Australia. The sparsely vegetated ecosystem represents the biggest ecosystem over Australia, which extends from the northern subtropical zone to southern Australia.

Figure 8 also includes the prior and posterior annual flux aggregated into these bioclimatic regions. It is evident that savanna and sparsely vegetated ecosystems were the regions across Australia that most contribute to posterior carbon sink estimated for 2015. The annual posterior carbon flux for savanna was −0.17 ± 0.03 PgC yr⁻¹ compared to the prior annual flux (0.09 ± 0.11 PgC yr⁻¹), and the annual posterior carbon sink over sparsely vegetated was even higher (−0.25 ± 0.07 PgC yr⁻¹) compared to prior annual flux (−0.01 ± 0.11 PgC yr⁻¹). These results are not unexpected because the sparsely vegetated ecosystem represents the largest bioclimatic region in Australia, and a slight shift of carbon fluxes across this area causes a significant impact on the total annual flux for this ecoregion and for total annual flux estimated for Australia.

Figure 8 Annual prior and posterior flux estimates aggregated into six bioclimatic classes (tropics, savanna, warm temperate, cool temperate, Mediterranean and sparsely vegetated) over the Australian region. Fossil fuel emissions are excluded.

Figure 9 shows the monthly time series of the prior and posterior terrestrial flux aggregated into these bioclimatic regions. Over the savanna ecosystem (Fig. 9b), our inversion indicates that from January to June, this ecosystem acted as carbon sink. In February, in this ecosystem, we see that the prior sink (−0.48 ±  0.40 PgC yr⁻¹) strengthens to a posterior of −1.07 ± 0.10 PgC yr⁻¹. The stronger carbon sink (relative to the prior) from January to March coincides with an increase of greenness in vegetation (positive EVI anomalies) in this ecosystem (see Fig. S1a and b in the Supplement), benefited by anomalous rainy conditions in January (see Fig. S2a in the Supplement). Thus, it seems that the anomalous increase of rainfall in northern Australia in January benefits the increase in vegetation growth and carbon uptake recorded in February. However, it is difficult to draw conclusions about the posterior carbon uptake seen in months subsequent to March because of unfavourable raining conditions and negative EVI anomalies in these periods.

Another noticeable difference between prior and posterior flux estimates over the savanna is seen in July and August. In July, we cannot conclude if the prior was a sink or source of carbon (0.19 ± 0.28 PgC yr⁻¹). However, our inversion indicates that savanna was acting as a carbon sink of −0.35 ± 0.11 PgC yr⁻¹. In August, the prior source (0.25 ±  0.20 PgC yr⁻¹) becomes a posterior carbon sink of −0.22 ±  0.08 PgC yr⁻¹. To understand the difference between the prior and posterior estimate in this period, we calculated the GPP estimated by the CABLE-BIOS3 model and the GPP generated by MODIS (see Appendix E; Fig. E1b). The temporal correlation between CABLE-BIOS3 and MODIS GPP was moderate (R = 0.69). According to MODIS estimates, the CABLE-BIOS3 GPP is overestimated from January to March and underestimated from May to October. The underestimation of the GPP flux by the CABLE-BIOS3 model might explain why we find a stronger posterior sink estimated by our inversion in this category.

Over the warm temperate region, from February to April, our posterior estimate suggests a carbon source (Fig. 9c). For this period, we cannot determine if the prior flux estimate was a carbon sink or source due to its uncertainty range. In February, the prior flux (−0.05 ±  0.08 PgC yr⁻¹) becomes a posterior carbon source of 0.17 ± 0.06 PgC yr⁻¹. In March, the prior estimate was nearly neutral (0.04 ± 0.05 PgC yr⁻¹) compared to the posterior carbon source estimate (0.17 ± 0.05 PgC yr⁻¹). The reduced carbon uptake estimated by the inversion in this period does not agree with the positive EVI anomalies seen in this region; however, it is likely that this extra carbon release to the atmosphere is related to an increase of leaf respiration in response to high temperatures recorded in 2015 for the majority of Australia (Annual climate statement, Bureau of Meteorology, 2015). Another possible reason for the relatively small shift from the prior in this period was most likely because the CABLE-BIOS3 GPP overestimates MODIS GPP (Appendix E1, Fig. E1c). For the warm temperate category, the correlation of CABLE-BIOS3 and MODIS GPP flux is high (R = 0.86).

We also see a subtle disagreement between prior and posterior estimates over the cool temperate ecosystem in April and May (Fig. 9d). In this period, our posterior estimate indicates that this category was a stronger carbon source than the prior flux estimate. In April, the inversion strengthened the prior source (0.12 ± 0.1 PgC yr⁻¹) to a posterior of 0.47 ± 0.05 PgC yr⁻¹. In May, we cannot define if the prior was a sink or a source (0.06 ± 0.09 PgC yr⁻¹); however, our assimilated fluxes indicate this category was acting as a posterior carbon source (0.36 ±  0.04 PgC yr⁻¹). The most likely reason for a larger carbon release in this period is related to negative EVI anomalies seen across this ecosystem. While it is true that April and May see positive EVI anomalies (Fig. S1d and e in the Supplement), in April, we notice predominantly negative EVI anomalies in the southern corner of the Australian (mainland) and Tasmania. The analysis of GPP between the CABLE-BIOS3 model and MODIS also shows some discrepancies (see Appendix E; Fig. E1d). For this category, in general, the CABLE-BIOS3 GPP is overestimating the productivity of the land for the whole year. For example, the absolute difference between both GPP datasets in April and June is about 0.2 PgC yr⁻¹. For the cool temperate category, the correlation between CABLE-BIOS3 and MODIS GPP is moderate (R = 0.73).

Another disagreement between the prior and posterior terrestrial flux estimate is seen over the Mediterranean ecoregion in August (Fig. 9e). Our posterior estimate is a flux of −0.35 ± 0.08 PgC yr⁻¹ compared to the prior of −0.12 ± 0.12 PgC yr⁻¹. An increase in vegetation productivity may also be the reason for the increase in carbon uptake by this category (positive EVI anomalies; Fig. S1h in the Supplement). This larger carbon uptake was likely a consequence of an increase of rainfall in this category (greater than 60 % on average (relative to the mean for 2000–2014) for some areas of this ecosystem; Fig. S3h in the Supplement). We also found that CABLE-BIOS3 underestimates MODIS GPP by 0.2 PgC yr⁻¹.

We found a noteworthy discrepancy between the prior and posterior flux estimate over sparsely vegetated ecosystem from March to September (Fig. 9f). In this period, in general, the absolute difference between the prior and posterior mean was around 0.4 PgC yr⁻¹. The largest difference found in July and September was about 0.6 PgC yr⁻¹. This highly unexpected and counterintuitive difference is not because we see a significant increase of positive EVI anomalies across this ecoregion (Fig. S1 in the Supplement). On the contrary, it is because a “small shift” in the carbon fluxes over this large ecosystem causes an important impact on the total carbon net flux calculated for the whole country. We clearly demonstrate this fact in Appendix B, Fig. B1. This figure in the Appendix shows the fluxes divided by area. In the western region of this category, we see evident positive EVI anomalies which start from April and last all the way through to September, which again line up with positive rainfall anomalies in that period.

Analysis of the GPP also shows the stronger posterior sink estimated by our inversion might be associated with a underestimation of the GPP by CABLE-BIOS3 in this category. The absolute difference between the CABLE-BIOS3 GPP and MODIS GPP was almost the same between May and September, with a range of 0.8–1.1 PgC yr⁻¹. This underestimation in GPP also suggests an underestimation of the land productivity. In this same category the posterior sink estimated in September disappears in October. For this period, our posterior source estimate (0.14 ± 0.08 PgC yr⁻¹) did not change much from the prior (0.18 ± 0.13 PgC yr⁻¹). In November and December, our posterior source was strengthened by the inversion. In November, we estimated a carbon source of 0.21 ± 0.08 PgC yr⁻¹ in comparison with the prior, which was 0.12 ± 0.13 PgC yr⁻¹. The extra carbon release estimated by the inversion in November might likely be associated with the combination of fires (Fig. S4k in the Supplement) located in the west and central northwestern region of Australia (Fig. S4k in the Supplement) and due to high temperatures recorded across Australia in summer (Annual climate statement, Bureau of Meteorology, 2015). These conditions certainly intensified the wildfires seen in that period.

In summary, our results showed that OCO-2 produced a shift in the carbon flux (relative to the prior) over the savanna and sparsely vegetated region. We found strong negative correlations (R>0.8) at grid-cell scale between the EVI anomalies and the posterior and prior difference in northern Australia (savanna ecosystem) and in the western region of the sparsely vegetated ecosystem, which align with the spatial pattern of rainfall in that area. These results suggest that our OCO-2 inversion might likely be better at capturing the anomalies in comparison with the biosphere land model.

Figure 9 Monthly time series of the Australian land biosphere prior and posterior CO₂ flux and their uncertainties in PgC yr⁻¹ aggregated over six bioclimatic regions. The prior and posterior estimates do not include fossil fuel emissions.

3.5 Evaluation of the inversion with independent data

In this section, we evaluate the accuracy of our posterior fluxes by comparing the residual between the prior and posterior concentrations simulated by CMAQ against TCCON and in situ observations. In this comparison, we simulate the posterior concentration with fluxes that were not only assimilated by nadir and glint satellite observations (LNLG) but also by the combination of both land nadir and land and ocean glint observations (LNLGOG). We decided to examine whether biases in our posterior concentration could improve when incorporating glint ocean observations into the inversion.

3.5.1 Comparison with TCCON observations

As mentioned in Sect. 2.5.1, we selected TCCON observations from three different sites (Darwin, Wollongong and Lauder; see Fig. 3). The comparison between the monthly mean column average from TCCON sites and the prior and posterior column-averaged concentration simulated by CMAQ, including both bias and RMSE are shown in Figs. 10 and 11, respectively.

In Fig. 10a, we see that in late spring, summer and early autumn in Australia (November to March) the posterior column-average simulated by CMAQ model (LNLG) is in better agreement with TCCON Darwin estimates than the prior. In this period, prior mean biases were reduced by approximately 30 %–80 %. For example, in November and December, the prior concentration biases were reduced from −0.25 (RMSE = 0.51) to 0.11 (RMSE = 0.42) and from −0.34 (RMSE = 0.48) to −0.02 (RMSE = 0.31), respectively. Large seasonal differences (approximately 1 ppm) are seen between the TCCON observation, the prior and posterior column-average concentrations (LNLG) from June to September. Despite the fact that we see an improvement of the prior biases in this period, assimilating OCO-2 data does not significantly reduce them. The remaining posterior concentration biases of about +1 ppm might be explained by spurious OCO-2 soundings affected by biomass burning aerosols seen in that period. In northern Australia, winter occurs in dry season, and it is highly impacted by wildfires (see Fig. S4 in the Supplement). OCO-2 spectrometers measure reflected sunlight from the Earth's surface, and regions heavily affected by fires can lead to a modification of the light path length because the instrument struggles to distinguish between photons reflected by intermediate scatterers and photons reflected from the Earth's surface (O'Dell et al., 2018). In terms of the posterior bias improvement when fluxes have been assimilated by OCO-2 LNLGOG data, we can see that the improvements are negligible, and in some periods such as January or May, the posterior biases get worse. This result suggests that the uncharacterized OCO-2 glint ocean bias degrades the performance of the inversion. We also found that our posterior column-average concentrations were better correlated with TCCON in comparison with the prior concentration (Appendix G, Table G1).

At Wollongong site, in general, we see a consistent overestimation of the prior and the posterior column average (LNLG) simulated by CMAQ (Fig. 10b). We observe a relatively slight reduction of the prior biases in February, March and November (spring, summer and early autumn in Australia, Fig. 11b). In November, for example, prior negative biases of about −0.74 ppm (RMSE = 1.22) were reduced to −0.40 ppm (RMSE = 1.13). The small reduction of the biases in this period is likely associated with strong winds coming from the ocean to the TCCON station (Fig. S7 in the Supplement). Wollongong TCCON site is strongly affected by ocean fluxes, which are less restricted by our inversion when we only use LNLG observations. In late autumn and winter at Wollongong, we see high significant positive posterior biases (range between 1.1 and 1.61 ppm). Biases in winter season are likely related to OCO-2 than TCCON biases. It has been found that passive satellite instruments have difficulties measuring at high and middle latitudes in winter because the Sun stays low in the sky (Wunch et al., 2017). A low solar altitude angle corresponds to a high solar zenith angle and high air mass, which means it takes longer for the sunlight to reach the satellite instrument. Biases related to high air masses (“long path length”) can be obtained because the absorption spectra tend to saturate at the line centre, causing the column line shape of the absorption line to be more sensitive (Jacobs et al., 2020). TCCON retrievals contain an air-mass-dependent bias, which is corrected using the method described by Wunch et al. (2011) and Deutscher et al. (2010). To evaluate whether any residual air-mass-dependent bias is present in the TCCON retrievals, which would cause the seasonal posterior biases seen here, we filtered the TCCON dataset to contain only selected solar zenith angles ≥ 40 and ≤ 50 (see Appendix G, Table G3). We found only a slight improvement of the posterior biases, which means TCCON retrieval bias is not likely the reason for the biases seen in winter. Similar to the Darwin site, we did not find an improvement of the prior biases by adding ocean glint data to the inversion. Besides, ocean glint data (as they are shown in Appendix C, Fig. C2d–h) are quite sparse around the Wollongong site, providing little constraint on carbon fluxes around this location.

Similar to the Darwin and Wollongong sites, we also found a systematic overestimation of our posterior column-average (LNLG) concentration at the TCCON Lauder site (Fig. 10c). Prior biases are lower than 0.8 ppm. A slight improvement was seen in June, July, September and November (Fig. 11c). In June and July (winter season) the reduction of the biases was only about 7 %. Improvement of the biases in November and September was better (10 % and 18 %), and we did not find improvement in the correlation for these months (see Appendix G, Table G4). The small or negligible improvement of the prior biases at this site is likely due to a combination of New Zealand's size and shape, the prevailing wind direction, and the fact that we do not allow much freedom for ocean fluxes by specifying a small prior uncertainty. Adding ocean glint observation to the inversion did not improve the accuracy of the biases at this site. The sparseness of OCO-2 soundings over the ocean around New Zealand in the period from May to September might explain the lack of improvement in this bias. Higher-resolution models and smaller correlation lengths (allowing more flexibility in spatial fluxes) would be required for good performance over New Zealand.

Figure 10 Box plot diagrams show the monthly mean average of CO₂ concentration at the Darwin (a), Lauder (b) and Wollongong (c) TCCON sites for 2015. The top edge of the box represents the 75th percentile and the bottom edge represents the 25th percentile. The top and bottom whiskers represent the 95th and 5th percentiles. The horizontal black line shows the median and the circle indicates the mean. Mean values are indicated by blue circles and median values by the black line.

Figure 11 CMAQ prior and posterior concentration bias and root mean square error at the (a) Darwin, (b) Wollongong and (c) Lauder TCCON sites for 2015. Yellow-green and green boxes represent prior posterior concentration biases, and coral and purple bars represent the RMSE. Box plots represent the 75th percentile and the bottom edge represents the 25th percentile. The top and bottom whiskers represent the 95th and 5th percentiles. The horizontal black line shows the median and the circle indicates the mean. Mean values are indicated by black circles and median values by the black line.

3.5.2 Comparison with in situ measurements

Figures 12 and 13 show the comparison between ground-based in situ measurements (Gunn Point, Burncluith, Ironbark and Cape Grim) and our prior and posterior concentrations simulated by CMAQ at the surface.

As illustrated in Fig. 12a, the inversion using only LNLG OCO-2 observations does not match Gunn Point observed concentrations well except in September. Most biases are negative, indicating that the posterior simulation at the surface of the CMAQ model underestimates the observations. The prior concentration indicates a better agreement, but biases are still significant. One possible explanation for the large negative posterior biases in January, February, March and December might be related to strong westerly winds that blow from the ocean to this site location (see Fig. S9 in the Supplement). Using only LNLG observations restricts the ability of our inversion system to optimize ocean fluxes, primarily because the ocean uncertainties set up in the inversion were relatively low compared to the uncertainties assigned over land. However, our results suggest that adding glint observations to the system improves the posterior biases at this site. These results are not unexpected because Gunn Point is a coastal site largely affected by ocean carbon fluxes. In February, and when the wind comes from the ocean, the posterior bias using LNLGOG shows a significant improvement compared to prior bias concentration. Here, we see a reduction of the bias from 1.93 (RMSE = 4.21) to 0.65 (RMSE = 3.79).

In winter (June to August), we see that the posterior biases using either OCO-2 LNLG or LNLGOG show no improvement, and the prior biases are in better agreement with the observations. One possible explanation might be related to the fact the column-integrated CO₂ measurements are less sensitive to near-surface dynamics compared to in situ measurements (Lauvaux and Davis, 2014), or to remaining bias in the OCO-2 data. Despite the fact that version 9 has an improvement in the bias correction, in the recent study performed by OCO-2 MIP (Peiro et al., 2021) shows that LNLG data still have large negative latitudinal biases in the Southern Hemisphere. Another potential explanation could be associated with an inaccurate representation of vertical transport within the planetary boundary layer in winter by the CMAQ model. Incorrect vertical transport might lead to erroneous horizontal distributions of air masses (Lauvaux and Davis, 2014). Therefore, correcting the prior column-average simulated by CMAQ to match OCO-2 might not improve near-surface simulations.

Improvements of the bias using LNLG observations at Ironbark are only seen in January, February, May, September and November (Fig. 12b). We found high negative posterior biases in June and July. The negative posterior bias in June, −2.79 ppm (RMSE = 3.53), might be associated with the small number of OCO-2 soundings located around Ironbark (see Appendix C, Fig. C1f) and the wind direction in that region. We can see in Fig. S7 (see the Supplement) that prevailing winds blow from the southeast, an area with no OCO-2 soundings to constrain fluxes. In July, posterior biases are larger than prior ones (−0.35 to −2.33 ppm). Again, biases in winter might be associated with error in the transport model or remaining biases in LNLG OCO-2 observations. At this site, we did not see an improvement of the posterior bias when we added glint ocean data to the inversion.

Results for the Burncluith station (Fig. 3c) are similar to those for Ironbark. This is not surprising given the stations' proximity. The posterior simulation performs better in July, October and December at Burncluith than Ironbark.

The posterior LNLG simulation at Cape Grim, shown in Fig. 12d, is in better agreement with the observations than the prior concentrations for the austral autumn and early winter of 2015. By contrast, high posterior negative biases are seen from September to December (> 2 ppm). This seasonality of bias is likely related to the seasonality of wind direction. The predominantly northerly flow in winter brings air from mainland Australia where fluxes have been constrained by OCO-2 observations. The southerly flow later in the year brings air from the Southern Ocean, unconstrained by observations. Similar to the Gunn Point site, we found that adding ocean glint observations to the inversion improved the prior mean concentration bias considerably. In January, for example, we see a reduction of the bias from −2.31 (RMSE = 3.38) to −0.54 (RMSE = 2.46) ppm using LNLGOG in the inversion. Again, these findings are not unexpected because Cape Grim is an oceanic station strongly influenced by oceanic carbon fluxes.

Figure 12 Box plot diagrams show the monthly mean average of CO₂ concentration at (a) Gunn Point, (b) Ironbark, (c) Burncluith and (d) Cape Grim for 2015. For details of what the different components of the box plot represent, see the caption of Fig. 10.

Figure 13 CMAQ prior and posterior concentration bias and root mean square error at (a) Gunn Point, (b) Ironbark, (c) Burncluith and (d) Cape Grim for 2015. For details of what the different components of the box plot represent, see the caption of Fig. 11. Note that bias and RMSE in Gunn Point in July exclude the highest concentration value for that period, which was 558.408 ppm; we did this to better represent the figure. Prior and posterior concentration biases for this period were −4.9 (RMSE = 16.8) and −5.1 ppm (RMSE = 17.1), respectively.

4 Discussion

We saw that semi-arid ecosystems in Australia, such as savanna and areas with sparse vegetation, are responsible to some extent for the stronger carbon sink (relative to the prior flux) recorded in 2015. We associated this carbon uptake with an increase in vegetation productivity (positive EVI anomalies) and an underestimation of the GPP flux by the CABLE land surface model. We speculate that the leaf area index (LAI) estimated by the land surface model CABLE-BIOS3 fails to capture the abrupt response of the terrestrial biosphere to rainfall over areas with sparse vegetation. This hypothesis could be tested by comparing CABLE-BIOS3 LAI with satellite LAI. However, this is beyond the scope of this study and will be taken up in a forthcoming article.

Figure 14 Panel (a) shows the ensemble annual mean of carbon fluxes derived by the in situ (green dot) and OCO-2 (LNLG) MIP (grey dot) global inversions, and the annual mean of the posterior fluxes estimated in this study (blue dot). The green and grey error bar represents the annual ensemble spread of the global models. In contrast, the blue error bar represents the uncertainties in the posterior flux calculated by different OSSEs estimated by Villalobos et al. (2020). Panel (b) shows the ensemble mean seasonal cycle of MIP in situ (green dots) and OCO-2 (LNLG) (grey dots) and the seasonal cycle of the posterior fluxes estimated in this study (blue dots) (all fluxes in this comparison are without fossil fuel emissions).

We compared our findings against the ensemble mean of nine in situ and OCO-2 (LNLG) MIP global inversions (AMES, PCTM, CAMS, CMS-Flux, CSU, CT, OU, TM5−4DVAR, UT) for 2015 (Fig. 14). We can see in Fig. 14a that the OCO-2 MIP and our posterior inversion suggest that Australia was a carbon sink for 2015. Our flux annual mean carbon flux estimate for Australia (−0.41 ± 0.08 PgC yr⁻¹) falls within the annual ensemble mean estimate of the OCO-2 MIP (LNLG) flux inversion (−0.23 ± 0.12 PgC yr⁻¹). Contrary to these findings, the ensemble mean of the in situ MIP inversions suggest that Australia was a carbon source of 0.20 ± 0.22 PgC yr⁻¹ for 2015. However, this fact cannot be concluded with high confidence because of the spread ensemble mean of in situ MIP global inversions.

In terms of seasonality, we can see in Fig. 14b that our inversion produces a similar seasonal pattern to the ensemble monthly mean of OCO-2 MIP (LNLG) (except for July) and produces an almost identical flux estimate for several months in 2015 (Fig. 14b). For example, in February and March, the monthly ensemble OCO-2 MIP (LNLG) estimate was −0.77 ± 0.48 and −0.96 ± 0.21 PgC yr⁻¹ compared to our posterior flux estimate, which was −0.77 ± 0.14 and −0.82 ±  0.17 PgC yr⁻¹, respectively. These findings give confidence that the posterior carbon fluxes estimated in this study are reliable.

Figure 15 As in Fig. 14, but in this case, the monthly mean posterior flux (blue dot) for July has been interpolated between June and August.

As we saw in Sect. 3.5, the validation of our inversion against Ironbark and Burncluith sites suggests that the anomalous sink seen in July might likely be related to errors in the transport model. If we interpolate July between June and August, we reproduce a monthly mean (−0.62 ± 0.34 PgC yr⁻¹), which is closer to the ensemble mean of the OCO-2 MIP (−0.47 ± 0.58 PgC yr⁻¹), considering the range of the spread the model and the uncertainties of our posterior fluxes (see Fig. 15b). By doing this interpolation, we shift our posterior annual flux from −0.41 to −0.32 PgC yr⁻¹, which is also closer to the annual ensemble mean of OCO-2 MIP (−0.23 PgC yr⁻¹).

The individual analysis of seasonal variations for the nine global carbon flux estimates, either derived by in situ or OCO-2 (LNLG) observations, show a large disagreement between them (Appendix I; Figs. I1 and I2). However, the variation of the seasonal cycle between in situ global flux inversions is more evident. We can see in Fig. I1 that the seasonal cycle derived by in situ global inversions over Australia is highly uncertain. One reason for the large disagreement between the in situ global inversions in Australia is the sparsity of observations. There around six existing Australian monitoring stations, and not all are operational (Ziehn et al., 2016). Besides, these global in situ inversions rely on measurements that come from monitoring stations such as Cape Grim, a station designed to sample background maritime air masses much of the time, thus providing minimal constraint on Australian fluxes (Haverd et al., 2013c). The OCO-2 MIP disagreement is likely driven by the choice of the prior flux, transport and data assimilation methodology used in the inversion (Crowell et al., 2019). For example, prior flux estimates used in global inversions rely on biosphere models such as CASA (van der Werf et al., 2017) or ORCHIDEE (Krinner et al., 2005). These models do not simulate well the NPP for grasslands (Wang et al., 2016) and hence underestimate the seasonality of the net ecosystem exchange for important ecosystems such as savanna and sparsely vegetated. This last point is critical for flux estimates over Australia because most of the land ecosystem is grassland and shrubs.

The analysis of the peak-to-peak seasonal variability of the nine OCO-2 global inversions individually shown in Table 4 indicates that, in general, November was the month in Australia with the highest carbon release to the atmosphere, similar to our posterior estimates. However, there is no unanimous agreement between them about the month with the largest carbon uptake. Analysing the ensemble mean of OCO-2 MIP, we can see in Fig. 14b that February and March were the months with the largest carbon uptake (−0.77 and −0.96 PgC yr⁻¹, respectively), a close estimate to our posterior fluxes, whose values were −0.77, −0.82 PgC yr⁻¹.

To further analyse our results, we also assess agreement between our posterior monthly spatial maps and these 10 global inversions individually. We plotted monthly maps for each global inversion (see Sect. 7 in the Supplement). We found that our posterior flux distribution across Australia agree well with at least four global inversions (TM5, CAMS, PCTM, AMES; Figs. S14, S15, S16 and S18 in the Supplement). We believe that this intercomparison is valuable for Australia because it shows that our results are reliable with a better spatial resolution than global inversion.

The comparison with OCO-2 MIP strengthens our confidence that our inversion is capturing fluxes across Australia. It supports the surprising result that Australia was a carbon sink in 2015 despite the significant El Niño (ECMWF, 2020). El Niño is only one of several large-scale drivers of the Australian climate, and we have already noted positive rainfall anomalies associated with some strong sinks.

Table 4 Summary of the peak-to-peak amplitude of our posterior terrestrial fluxes, prior fluxes and terrestrial fluxes from nine different in situ and OCO-2 (LNLG) MIP global inversions. Units are PgC yr⁻¹

¹ Figure I1, Appendix I, shows the seasonal cycle of the Australian carbon fluxes derived by the models contributing to in situ MIP inversions. ² Figure I2, Appendix I, shows the seasonal cycle of the Australian carbon fluxes derived by the models contributing to the OCO-2 (LNLG) MIP global inversions.

Evaluating the carbon fluxes assimilated in this study through the validation of our assimilated posterior field is difficult. We saw in Sect. 3.5.1 that most of the improvement of the prior concentration biases was seen in TCCON site (mainly the Darwin site in northern Australia) in the summer season compared to in situ observations. It is difficult to validate our posterior concentration field against monitoring stations located in coastal areas, such as Gunn Point or Cape Grim. These sites are strongly affected by oceanic fluxes, which are less restricted by the inversion when they are assimilated by only using LNLG observations. We demonstrated that adding OCO-2 glint observations to the inversion improves the biases considerably at these sites but not for TCCON sites or sites located far away from the ocean such as Burncluith and Ironbark.

We also found that adding OCO-2 ocean glint data to the inversion does not significantly alter the annual carbon sink estimated for the continent (0.36 PgC yr⁻¹) compared to the estimate made by only using LNLG OCO-2 observations (0.41 PgC yr⁻¹), suggesting that adding ocean glint observations does not strongly drive the continental carbon budget.

To assess the impact of biases in the lateral boundaries of the CMAQ domain, we performed two sensitivity experiments. In both experiments, we add a constant offset of 0.25 ppm to each grid cell of the BCs. In the first experiment we solve for the BCs and use LNLG. This induces a bias in our posterior annual flux of −0.8 PgC yr⁻¹. Solving for the BCs and using LNLGOG observations reduces the bias further to −0.4 PgC yr⁻¹. Adding 0.25 ppm everywhere is an extreme test since the global assimilation fields we use are unlikely to have such systematic errors against data they assimilate. The results do highlight the importance of solving for the BCs in a regional inverse system and also the importance of large domains with enough observations in a buffer region around our area of interest. These results are also reinforced by the good agreement that our assimilated fluxes have with OCO-2 MIP.

There are still several methodological choices that are somewhat arbitrary in this study. The most important is the implied spatial resolution. This is determined by the correlation length used in the prior uncertainty as much as the resolution of CMAQ. Villalobos et al. (2020) showed the impact of this correlation length on posterior uncertainty, and our choice makes a compromise between the information available from observations and avoidance of aggregation errors (Kaminski et al., 2001). A more important limitation is the restriction to 1 year. This will be addressed in a forthcoming study extending over the OCO-2 dataset.

Higher-resolution flux inversions assimilating satellite retrievals of greenhouse gas concentrations, as illustrated by this study, will be increasingly important in a world seeking climate solutions and a better understanding of the global carbon cycle. They will likely play a role in not just addressing questions of scientific interest but also in ongoing monitoring and assessment of emission targets. Australia, as a large and geographically isolated land mass, with a terrestrial biosphere highly responsive to climate drivers, offers an ideal testing ground for such flux inversions. The overall success of this study suggests great promise, especially in regions with sparse in situ networks.

5 Conclusions

We performed a four-dimensional variational data assimilation inversion to estimate Australian CO₂ fluxes for 2015. The inversion was based around the Community Multiscale Air Quality (CMAQ) transport–dispersion model and satellite data from the Orbiting Carbon Observatory-2 (OCO-2) (land nadir and glint data, version 9). Our regional inversion suggests that Australia was a carbon sink of −0.41 ± 0.03 PgC yr⁻¹ compared to the prior estimate of 0.09 ± 0.20 PgC yr⁻¹. We found a higher-than-average increase of land productivity (relative to 2000–2014) over the savanna ecosystem (northern Australia) during summer, leading to most of the carbon uptake in this ecosystem. The sparsely vegetated ecosystem is the most extensive ecosystem over Australia and also showed a slight increase of land productivity in the autumn and winter seasons in the western region of Australia, which was also driven by an increase of vegetation productivity in response to positive rainfall anomalies in this period. We also found that the higher carbon uptake by our inversion (relative to the prior) was due to an underestimation of GPP simulated by the CABLE-BIOS3 model.

Evaluation with the TCCON Darwin site shows that our inversion is able to reduce biases mainly in the summer period compared to the winter season. Reduction of the biases at TCCON Lauder and Wollongong showed a very slight systematic decrease, mostly because both sites are strongly affected by ocean winds and there is a reduced number of OCO-2 soundings passing over these sites in some periods in 2015. Posterior column-integrated simulations at coastal monitoring sites are challenging to validate because they are strongly influenced by ocean fluxes, which were assigned small uncertainties in our inversion. Comparison with in situ data was also a challenge mainly over oceanic monitoring stations such as Cape Grim and Gunn Point sites, which are also strongly impacted by ocean fluxes. Comparison with monitoring stations over land such as Ironbark and Burncluith also shows difficulties in simultaneously matching column-integrated and surface data, most likely linked to model vertical transport. The scarcity of in situ observations across the Australian continent, mainly over the savanna and sparsely vegetated ecosystem, restricts our ability to conclude with confidence whether the stronger carbon sink (relative to the prior) found in those ecosystems is real or not. However, the comparison with the annual and monthly ensemble means of the OCO-2 MIP is encouraging and supports our results.

Appendix A: Spatial pattern of the differences between posterior and prior fluxes

Figure A1 Spatial pattern of the differences between posterior and prior fluxes for 2015.

Appendix B: Monthly time series of the Australian land biosphere prior and posterior CO₂ flux over six bioclimatic regions

Figure B1 Monthly time series of the Australian land biosphere prior and posterior CO₂ flux and their uncertainties in gC m⁻² yr⁻¹ over six bioclimatic regions. The prior and posterior estimates do not include fossil fuel emissions.

Appendix C: Spatial distribution of OCO-2 data across Australia

Figure C1 Spatial distribution of OCO-2 soundings (LNLG) over the CMAQ domain for 2015.

Figure C2 Spatial distribution of OCO-2 soundings (LNLGOG) over the CMAQ domain for 2015.

Appendix D: Spatial distribution of the prior and posterior uncertainties across Australia

Figure D1 Prior uncertainties accounting for the major terms in the CO₂ budget (anthropogenic fluxes, fires, land and ocean exchange).

Figure D2 Posterior uncertainties calculated by OSSEs in Villalobos et al. (2020).

Appendix E: Time series of BIOS GPP and MODIS GPP

Figure E1 Time series of monthly mean of CABLE-BIOS3 GPP and MODIS GPP.

Appendix F: Posterior fluxes assimilated by using LNLGOG satellite observations

Table F1 Australia terrestrial carbon fluxes estimated using LNLG and LNLGOG for 2015 (units: PgC yr⁻¹).

Appendix G: TCCON comparison

Table G1 Analysis of the residual between CMAQ prior and posterior simulations and the TCCON Darwin site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Table G2 Analysis of the residual between CMAQ prior and posterior simulations and the TCCON Wollongong site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Table G3 Analysis of the residual between CMAQ prior and posterior simulations and the TCCON Wollongong^* site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

^* Wollongong TCCON data are filtered by solar zenith angles ≤ 40 and ≥ 50^∘.

Table G4 Analysis of the residual between CMAQ prior and posterior simulations and the TCCON Lauder site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Appendix H: In situ comparison

Table H1 Analysis of the residual between CMAQ prior and posterior simulations and the Gunn Point site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Table H2 Analysis of the residual between CMAQ prior and posterior simulations and the Ironbark site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Table H3 Analysis of the residual between CMAQ prior and posterior simulations and the Burncluith site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Table H4 Analysis of the residual between CMAQ prior and posterior simulations and the Cape Grim site for 2015: averaged bias, RMSE and Pearson's coefficient (R).

Appendix I: Australian fluxes derived by MIP in situ and OCO-2 (LNLG) global inversions for 2015

Figure I1 Time series of monthly mean carbon fluxes derived by MIP in situ (IS) global inversion for 2015.

Figure I2 Time series of monthly mean carbon fluxes derived by OCO-2 MIP (LNLG) global inversion for 2015.