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).

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).

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 1A simplified diagram of the four-dimensional variational data assimilation we used to estimate CO₂ surface fluxes over Australia.

Download

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 2Representation 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.

Table 1Physics parameterizations used in WRF model setup.

Download Print Version | Download XLSX

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:

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 2Reference of the TCCON stations used in this work for evaluation of our inverse model system.

Download Print Version | Download XLSX

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 3TCCON 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.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 3Summary of the configuration of the OCO-2 MIP (version 9) design.

Download Print Version | Download XLSX

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 4Bias 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.

Download

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 5Time 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.

Download

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 6Prior 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 7Posterior 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 8Annual 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 9Monthly 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 10Box 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.

Download

Figure 11CMAQ 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.

Download

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 12Box 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.

Download

Figure 13CMAQ 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.

Download