Interannual variability in the Australian carbon cycle over 2015-2019, based on assimilation of OCO-2 satellite data

In this study, we employ a regional inverse modelling approach to estimate monthly carbon fluxes over the Australian continent for 2015–2019 using the assimilation of the total column-averaged mole fractions of carbon dioxide from the Orbiting Carbon Observatory-2 (OCO-2, version 9). Subsequently, we study the carbon cycle variations and relate their fluctuations to anomalies in vegetation productivity and climate drivers. Our five-year regional carbon flux inversion suggests that Australia 5 was a carbon sink averaging -0.46 ± 0.08 PgC yr−1 (excluding fossil fuel emissions), largely influenced by a strong carbon uptake (-1.04 PgC yr−1) recorded in 2016. Australia semi-arid ecosystems, such as sparsely vegetated regions (in central Australia) and savanna (in northern Australia), were the main contributors to the carbon uptake in 2016. These regions showed relatively high vegetation productivity, high rainfall and low temperature in 2016. In contrast to the large carbon sink found in 2016, the large carbon outgassing recorded in 2019 coincides with an unprecedented deficit of rainfall and higher than 10 average temperature across Australia. Comparison of the posterior column average CO2 concentration against the Total Carbon Column Observing Networks (TCCON) and in situ measurements offers limited insight into the fluxes assimilated with OCO2. However, the lack of these monitoring stations across Australia, mainly over ecosystems such as the savanna and areas with sparse vegetation, impedes us from providing strong conclusions. Comparison of our flux inversion to the ensemble mean carbon flux of the OCO-2 Multi-model Intercomparison Project (MIP) (2015-2018) agrees with our findings, and their results 15 also suggest that Australia was a strong carbon sink in 2016 (-0.73 ± 0.41 PgC yr−1). The analysis of the variability of the nine models that participate in the OCO-2 MIP also aligns with our findings, and it gives us the confidence to say that changes in rainfall and temperature drive most of the carbon flux variability across Australia. 1 https://doi.org/10.5194/acp-2022-15 Preprint. Discussion started: 11 February 2022 c © Author(s) 2022. CC BY 4.0 License.

inversion. Further details can be found in Villalobos et al. (2020Villalobos et al. ( , 2021.

Inversion set-up
Our regional inversion system optimizes monthly-mean gridded-based surface carbon emissions x using a four-dimensional variational data assimilation method, which was configured to use the Community Multi-scale Air Quality Model (CMAQ) (version, v5.3) and its adjoint (version 4.5.1; Hakami et al., 2007). Each year of the five-year period was run independently, 70 with a spin up of one month for each year. Our system optimizes CO 2 surface fluxes by finding the minimum of the cost function J(x) shown in Eq. 1. Notation in this study follows Rayner et al. (2019).
3-hourly NBP estimates as input for CMAQ (further details of how we constructed NBP can be found in Sect.2. 3 Villalobos et al., 2021). The prior error covariance matrix of the terrestrial biosphere flux from CABLE was assumed to be an approximation of the net primary productivity (NPP) following the approach of Chevallier et al. (2010) with a ceiling of 3 gC m −2 d −1 .
We assumed that these uncertainties were spatially correlated with length-scale 500 km over land following Basu et al. (2013). Within our inversion system, no temporal correlations were considered.

120
Fossil fuel emissions used here were based on two different inventory data sets: the Open-source Data Inventory for Anthropogenic CO 2 (ODIAC) (version 2019) (Oda et al., 2018) and the Emissions Database for Global Atmospheric Research (EDGAR) (Crippa et al., 2020). We added some missing sectors from the EDGAR inventory to ODIAC (such as aviation climbing and descent, aviation cruise, and aviation landing and take-off datasets). ODIAC is a global gridded product distributed at 0.1 • × 0.1 • spatial resolution over land, which uses power plant profiles (emissions intensity and geographical location) and 125 satellite-observed nighttime lights. We used ODIAC monthly fluxes and incorporated a diurnal scale factor to estimate diurnal CO 2 emission variability (Nassar et al., 2013). Given that the ODIAC product only covers the period from 2015 to 2018, we repeated the data from 2018 in 2019 but increased the value in each grid cell by 1.7%, which represents the mean annual growth rate of these emissions from 1970 to 2018. EDGAR is also gridded at 0.1 • × 0.1 • with monthly temporal resolution. There is no EDGAR gridded product for 2016-2019, so we repeated the 2015 product to cover the other years. We increased EDGAR 130 aviation emissions by 2.5%, which represents the mean growth rate in this emissions sector from 2016 to 2019. Fossil fuel prior uncertainties were assigned to be 0.44 times the value of the monthly fossil fuel estimates described above (see details in Sect. 2.3. in Villalobos et al., 2021). Errors in fossil fuel emissions were assumed to be uncorrelated.
Ocean flux estimates were selected from CAMS global data (version v19r1) (Chevallier, 2019). Ocean prior uncertainties were assumed to be 0.2gC m −2 d −1 and uniform across the ocean, as in Chevallier et al. (2010). Similar to correlations for 135 biosphere prior uncertainties, uncertainties of ocean fluxes were assumed to be correlated in space with length-scale 1000 km.
Fire prior emissions were selected from the Global Fire Emission Database (GFED), version 4.1s, which includes emissions from small fires. Fire emissions uncertainties were assumed to be 20% of the GFED emissions and correlated in space with length-scale 500 km, but not in time. The combination of all prior fluxes was regridded to the spatial resolution of the CMAQ model.
mode the 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 • ). In Target mode, the instrument points towards a specific location on the ground.

150
Target mode is use for validation, where the performance of the instrument is validated against ground-based observations from the Total Carbon Column Observation Network (TCCON) (Wunch et al., 2011).
In this study, the regional inversion was performed using the combination of both land (nadir and glint) (LNLG) OCO-2 observations (version 9). We used the combination of both datasets because it has been demonstrated by Miller and Michalak (2020) that combining both modes provides a stronger and better constraint of CO 2 fluxes at regional scales. Also, both datasets provides a very good coverage over the Australian region. Such spatial coverage offers good potential to help constrain regional biosphere CO 2 fluxes.
Given that the OCO-2 spatial resolution (1.29 km × 2.25 km) is higher than the CMAQ model grid cell (81 × 81 km), the OCO-2 data were averaged to the CMAQ model grid-level following a two-step process described in (Sect 2.3, Villalobos et al., 2021). The first step involves averaging all OCO-2 soundings across 1-second intervals, while the second step involves 165 averaging these 1-second averages into the CMAQ vertical column (approximately 11-seconds averages). The algorithm to estimate the uncertainties across 1-second averages follows Crowell et al. (2019). Here, we considered three different forms of uncertainty calculation. First, we assumed that uncertainties that fall within 1-second span were perfectly correlated in time and space (uncertainties defined as σ s ). Second, given that the average of OCO-2 uncertainties (σ s ) is relatively lower than the real OCO-2 uncertainties (mainly because they only consider the errors from measurement noise, and not systematic errors), we 170 also used the spread (standard deviation) of the OCO-2 retrievals in the 1-second average (uncertainties defined as σ r ). Third, we also considered a baseline uncertainty (defined as σ b ) for cases where the number of OCO-2 soundings was not enough to compute a realistic spread. 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 ). For each grid-cell, we also added (in quadrature) to this term 0.5 ppm as the contribution of the CMAQ model uncertainty (defined as σ m ). We also increase the 175 final observation uncertainty by a factor of 1.2 to satisfy the theoretical assumptions of the inversion (Villalobos et al., 2021).
We interpolated the retrieval OCO-2 profile to the CMAQ model vertical profile as described in (Sect. 2.6., Villalobos et al., 2020). Note that we only selected OCO-2 retrievals with quality flag "0" and bias-corrected data, as described by Kiel et al. (2019).
operates a Picarro G2301 analyser; however, the inlet is positioned at 70 m. The instrument precision for these spectrometers is better than 0.1 ppm (Etheridge et al., 2014) and all measurements are calibrated to the World Meteorological Organization (WMO) X2007 CO 2 mole fraction scale (Zhao and Tans, 2006), ensuring comparability between all measurements used. We 215 note that we used "baseline" and "non-baseline" data from Cape Grim. Baseline data is selected when winds blow straight off the Southern Ocean and have not been in recent contact with land. In this study, we used both datasets because our inversion only uses OCO-2 soundings located over land. in-situ sites (blue points). This map also shows a classification of six bioclimatic regions for Australia.

Australian bioclimatic classification
To understand which ecosystems contributed the most to the Australian inter-annual carbon flux variability between 2015-220 2019, we divided the continent into six bioclimatic classes: tropical, savanna, warm temperate, cool temperate, Mediterranean and sparsely vegetated (Fig. 1). We used the same six bioclimatic regions at a 0.05 • spatial resolution as in Haverd et al. (2013a). The classes were regridded over our CMAQ grid (81 × 81 km) resolution.

Climate data
In order to analyse the impact of climatic drivers on Australian terrestrial carbon cycle variability, we investigated the anomalies 225 of rainfall and temperature across Australia for the period 2015-2019. Rainfall data were selected from the Australian Water  (Jones et al., 2009) for the period 2015 to 2019. AWAP is a gridded product at 0.05 • resolution. It is generated by spline interpolation of in situ rainfall observations. We also used air temperature data at 2 m above the land surface from ERA5, the fifth generation of European Centre for Medium-Range Weather Forecasts (ECMWF) atmospheric reanalyses. The dataset selected from ERA5 was monthly and it was gridded at 230 0.25 degrees spatial resolution. We constructed 3-month running means of rainfall anomalies and air temperature anomalies relative to a mean across 2015-2019. These anomalies were calculated by subtracting their long-term mean (2015-2019) for each month from the raw time series and constructing the 3-month running mean on the resultant time series. We regridded the rainfall anomalies onto the grid of the CMAQ model in order to simplify the comparison with the estimated terrestrial carbon uptake from the flux inversions. Imaging Spectroradiometer (MODIS) MOD13C1 version 6 data product, which flies on board Terra, a NASA earth-observing satellite (Didan, 2014). The MODIS EVI is a gridded product, which has a temporal resolution of 16 days composite and 0.05degree spatial resolution. The EVI ranges from -0.2 to +1, where values less than 0 indicate a lack of green vegetation or arid areas. We calculate the 3-month running mean of EVI anomalies in Australia relative to the long-term mean from 2015-2019 245 and subtract the mean seasonal cycle. These monthly EVI MODIS products were also regridded into the CMAQ domain to calculate the temporal correlation between prior and posterior flux anomalies.

Gross Primary Productivity (GPP)
To understand the difference between posterior and prior fluxes, we compared the climatological seasonal cycle of the gross primary productivity (GPP) from the CABLE BIOS3 model against the remote-sensing based DIFFUSE model (Donohue et al.,250 2014), and the latest MODIS terra GPP product (MOD17A2H version 6) (Running et al., 2015) for the period 2015-2019. We also calculated 3-month running mean GPP anomalies for these three datasets.
The DIFFUSE GPP estimates are taken to be the product of the fraction of photosynthetically active radiation (PAR) absorbed by vegetation and the light-use efficiency. These datasets have a temporal resolution of 16-days at 250 m resolution.
Similar to the DIFFUSE estimates, the MODIS GPP product is based on a light-use efficiency approach and provides a cumu-255 lative 8-day composite product gridded at 500 m. For comparison, the CABLE BIOS3, DIFFUSE and MODIS GPP products were averaged to a monthly resolution and regridded over the CMAQ domain.

Global Atmospheric Inversions
We compared our Australian assimilated fluxes against nine independent global atmospheric inversions: AMES, PCTM, CAMS, CMS-Flux, CSU, CT, OU, TM5−4DVAR, UT (Sect. 4). These global inversions are part of the OCO-2 Model In-260 tercomparison Project (MIP) Peiro et al., 2021). In this study, we used the OCO-2 MIP flux version found in Peiro et al. (2021). In Peiro et al. (2021), the global inversions were performed using the assimilation OCO-2 (version 9, bias-corrected) from 2015-2018. A summary of these nine global inversions is given in Table 1, and a complete description of them and their input fields can be found in Peiro et al. (2021, Appendix A: model information). We can see in Table 1 that all global inversions were run using different inverse systems and were configured at different spatial resolutions with different 265 atmospheric transport models and prior fluxes. Some global inversion methods use a four-dimensional variational (4D-Var) approach, while others utilize the technique known as ensemble Kalman filter (EnKF) or Bayesian synthesis.

Inversion performance
In our inversion system, the L-BFGS algorithm iteratively adjusts the control vector until the cost function reaches an optimal 270 solution. In Bayesian inverse problems we require that the observational residuals (simulated − observed) and increments (posterior − background) are consistent with the assumed probability density functions (PDFs). This implies that the cost function should be approximately half the number of observations (Tarantola, 1987, p.211). Table 2 shows the analysis of convergence for our five-year inversion. In this table, we can see that for each year, the ratio between final cost function J(x) and observation was about 0.5, indicating that our system is self-consistent. In 2019, most of the prior biases were positive, except for January where biases were negative. For this period, prior concentration biases were -0.18 ppm, with a RMSE of 1.01 ppm, which were reduce by the inversion to -0.01 ppm and RMSE 0.87. Considerable positive biases are seen in August and September, which values were around 0.48 ppm, with a RMSE of 0.86 ppm. December was the only month that prior concentrations were in a good agreement with OCO-2 observations. Prior biases in this period were negligible (0.08 ppm, RMSE 0.98).

Australian posterior and prior monthly carbon flux estimates
Our five year inversion suggests that Australia was a carbon sink of -0.46 ± 0.09 PgC yr −1 compared to the prior flux estimate, which was 0.11 ± 0.17 PgC yr −1 . Within these estimates we only included the terrestrial part of the Australian carbon cycle, including fires but not fossil fuel emissions. We subtracted the fossil fuel emissions from the prior and posterior estimate because it only represents 25% of the total flux, and is subject to much less interannual variability than the biogenic fluxes. 305 Figure 3a shows the posterior and prior monthly mean terrestrial CO 2 fluxes for the period 2015-2019 along with their uncertainties. Posterior flux uncertainties from 2016 to 2019 were assumed to be the same as those calculated for 2015, which were estimated by five different observing system simulation OSSE experiments (see more details in Villalobos et al., 2021).
During the period 2015-2019, the posterior flux estimates show a stronger seasonal cycle compared to the prior flux estimate (Fig. 3). In terms of inversion results, we see that OCO-2 sees a strong seasonal biospheric carbon uptake each year between 310 June and September (winter and early spring in Australia), and a stronger carbon source from November to December (late spring and early summer in Australia). As is shown in Fig. 2, the stronger carbon uptake seen in winter and early spring occurs because the prior concentration simulated by CMAQ model overestimate OCO-2 observations in this period.
In August 2016, we see that the carbon uptake by Australia was about -2.92 ± 0.31 PgC yr −1 compared to the prior (-0.25 ± 0.53 PgC yr −1 ). This uptake rate was the largest seen in the period 2015-2019. However, this uptake was rapidly reduced in 315 September (-1.09 ± 0.23 PgC yr −1 ) and October (-1.27 ± 0.39 PgC yr −1 ), being reversed in November (1.12 ± 0.31 PgC yr −1 ) and December (1.10 ± 0.17 PgC yr −1 ). In July 2017, we again see a strong carbon uptake (-1.96 ± 0.39 PgC yr −1 ), which is rapidly reduced in September 2017 (-0.15 ± 0.23 PgC yr −1 ), and eliminated in October and November 2017. In 2018, the largest difference with the prior flux is seen from June to September. For these four months, we obtain an uptake rate of (-1.11 ± 0.3 PgC yr −1 ). This carbon uptake is rapidly released back to the atmosphere in November and December. In 2019, we see 320 that the largest carbon uptake observed in August (-1.77 ± 0.31 PgC yr −1 ) is returned to the atmosphere in December 1.50 ± 0.17 PgC yr −1 . This carbon release is the largest registered for the period 2015-2019. Later, in the discussion section, we will show that the stronger carbon uptake is likely related to an underestimated GPP simulated by the CABLE model over the largest ecosystem across Australia such as savanna and sparsely vegetated regions.   largely affected by the drought that covered large parts of the country in 2019. These conditions were also intensified by the wildfires, which were the worst on record. anomalies were negatively correlated with EVI anomalies (Fig. 6d), which indicate that an increase in vegetation productivity is associated with an increase of the carbon uptake by the terrestrial ecosystem. which coincide with negative posterior carbon anomaly in this period (∼ -0.35 PgC yr −1 ). In 2017, we see that the positive EVI anomalies recorded from January to April do not align with the larger than average posterior carbon release in this period. 370 We believe that the few OCO-2 soundings found in this period limit the potential of our inversion to constraint the surface fluxes in the savanna ecosystem (see Appendix A, Fig

390
As mentioned at the beginning of this section, we did not find correlations between posterior anomalies and EVI anomalies over the cool temperate category. But, we did find a moderate correlation between prior anomalies and EVI anomalies (R =-0.36) (Fig. 6g). We believe that this discrepancy is due to difficulty in estimating the posterior fluxes along coastal areas.
We should note that this category is extended in the south-eastern corner of Australia, which is strongly affected by strong winds that come from the southern ocean, and cloudy conditions that prevent to retrieve satellite observations in this area.   415 We saw in the previous section that the variability of the prior and posterior flux anomalies is associated with the variability of vegetation greenness (EVI as proxy of the vegetation productivity). A higher land productivity (i.e., positive EVI anomalies) leads to a higher uptake of CO 2 , and negative EVI anomalies cause the opposite effect. Such variability is more evident in some ecosystems than others. We observed that savanna and sparsely vegetated ecosystems account for the highest variability of our posterior fluxes. We also see some variability of the posterior fluxes over warm temperate, but not as strong as the categories mentioned above. In this section, we will examine how changes in rainfall and temperature drive most of the variability in these ecosystems, with a focus over the savanna, warm temperate and sparsely vegetated regions.  Fig. 8 shows the flux anomalies along with 3-month running mean of temperature anomalies. In general, we found that rainfall and temperature anomalies are negatively correlated, and Australia's semi-arid ecosystems are highly responsive to these climate drivers. It is evident that above-average rainfall and low temperatures lead to unexpectedly rapid vegetation growth over semi-arid ecosystems shown in An increase in the amount of rainfall might explain the larger than average negative posterior anomalies in this period. We see that the posterior negative flux anomalies remain in January 2016, but its intensity decreases as the amount of rainfall also decrease. Negative rainfall anomalies in 2016 are only seen from February to April. In the period May 2016 to December 2017, we again notice an increase of the amount of rainfall, which coincides with the anomalous posterior sink in that period.

435
For this category, the temporal correlation between the posterior flux anomalies and rainfall anomalies was weak (R = -0.15).
However, we plotted the temporal correlation calculated at each grid-point, and we discovered that in some areas over savanna, correlations are strongly negative, with values of about -0.4 to -0.8 (Fig. 9b).
In regards to temperature anomalies over savanna, we can see in Fig. 8b that they were below average by -0.5 • C in 2015.
Despite having estimated negative temperature anomalies recorded in 2015, we do not see a carbon sink anomaly in this 440 ecosystem as we expected. We believe that the savanna ecosystem responded strongly to the deficit of water seen in this period compared to the relatively small changes in the air temperature. Positive temperature anomalies recorded from March to June in 2016 (2 • C above average) lead to an increase in the amount of carbon released to the atmosphere in that period.
Temperature in 2017 was predominantly below average from January to June. In 2018, from February to March, we mostly found temperatures that were below average by -1 • C. In this period, negative temperature anomalies in combination with Positive air temperature anomalies (4 • C above average) coincide with the higher than average release of carbon recorded by the ecosystem in this period. 455 We also see climate relationships in 2016 and 2017 over sparsely vegetated areas (Fig. 7f). We see that the amount of rainfall in this category is above average from August 2016 to April 2017. These positive rainfall anomalies (25-30 mm) coincides with the negative posterior flux anomalies estimated in that period and with negative air temperature anomalies of about -1 • C to -2 • C. From May 2017 on-wards the positive precipitation anomalies are cancelled out and shift to being negative until the end of 2019. This deficit of rainfall coincides with a change in air temperature anomalies in this period.

460
In 2018, no significant changes in temperature were observed. In 2019, positive temperature anomalies led us to expect an increase in the ecosystem respiration. In this ecoregion, the time series temporal correlation between posterior anomalies and temperature anomalies was positive (R = 0.34). However, if we look at temporal correlations calculated at each grid-point in  Australia experiences a wet season (November to April), which is highly impacted by monsoonal rains and storms. Winter, the dry season in this region, is affected by fires. Some studies (i.e. Taylor et al., 2016) suggest that some OCO-2 retrievals can 480 be biased by clouds during the wet season and smoke aerosol plumes during the dry season, mainly because the OCO-2 cloud screening algorithms present some difficulty in identifying clouds near the surface. With regard to correlation analysis, we also found that the relationship between observations and the posterior simulations is improved in some periods (see Appendix A, Table B1).
Evaluation at the Wollongong site ( Fig. 10b) also shows systematic differences with our posterior concentrations. From in bias is negligible when the wind blows from the ocean to this site or not many OCO-2 soundings were found around the monitoring location. Improvement of correlation at Wollongong site are shown in Appendix A, (Table B2).
Unlike Wollongong site, we see a persistent overestimation of both the prior and posterior estimates at TCCON Lauder ( Fig. 10c). Posterior biases are, however, less than 1 ppm. Prior biases at the Lauder site were mainly reduced in winter and early spring. The reduction of the biases at this site were modest (about 10 to 25%). New Zealand is (relative to the 495 Australian mainland) much smaller and narrower along the south-west to north-east direction, and thus strongly affected by oceanic airflow. The smaller size means relatively few OCO-2 soundings are retrieved over this area. Ocean fluxes that affect New Zealand have less freedom to be modified due to the small prior uncertainties assumed in the inversion. Analysis of the correlation is shown in Appendix A, Table B3.  of the biases in some periods. For example, we note that prior biases are reduced in summer (e.g., January and February) and spring (e.g., September and October). In these seasons, improvements of the biases are greater than 50%. No improvement of 515 the biases are seen in June, July and August (winter season in Australia) 2016. In this period, the prior concentrations show a better agreement with the observations, with biases ranges between 0.54 to -0.46 ppm compared to the posterior biases range between -2.79 to -3.57 ppm). One possible explanation for these results could be associated with errors in the transport of the CMAQ model, such as a misrepresentation of the wind flows or errors related to the vertical mixing scheme. Transport errors related to vertical mixing are often associated with atmospheric turbulence near the surface where most CO 2 observations are 520 made. Therefore, correcting the prior column average simulated by CMAQ to match OCO-2, the correction imposed by the inversion near the surface might be wrong. For 2018, prior biases were not reduced but they were small (range -0.16 to 0.46 ppm), except February (-2.31 ppm) compared to posterior biases. Again, we suggest that high posterior biases might be related to the vertical transport of the CMAQ model. We also did not find much improvement in the correlations at this site (see Appendix A, Table C4).

525
Similar to Burncluith, Ironbark also has gaps in the data (Fig. 12c). The data for the Ironbark site spans January 2015 to April 18th to 20th , we also see a similar event that was not captured by the inversion, causing a posterior concentration bias of -6.44 ppm (RMSE = 9.82). During these 3 days the concentrations registered at Burnlcluith was greater than 450 ppm.
Cape Grim is the only site with a complete time-series of observations during this period (Fig. 12d). Like Gunn Point, Cape Grim is a coastal site, and it is affected by strong westerly winds that blow from the ocean into Tasmania. We can see in Fig associated with the the strong westerly and north-westerly winds that come from the ocean to Tasmania. As mentioned before, Cape Grim is a coastal station, whose aim is to record clean air that blows from the southern ocean and it is not representative of Tasmania's air mass. To further examine the difference between the prior and posterior terrestrial fluxes shown in Fig. 3, we calculated the climatological seasonal cycle of the prior, posterior (Fig. 14) and the gross primary productivity (GPP) fluxes derived from CABLE BIOS3, MODIS and the DIFFUSE model (Fig. 15) for the period 2015 and 2019 aggregated by the same six bio-climate ecosystems described in Section 3.3.1.
It is evident that savanna (Fig. 14b) and sparsely vegetated (Fig. 14d) ecosystems are the regions that show the largest 555 difference between prior and posterior flux estimates. Over the savanna region, the most evident difference is seen from June to September. The absolute difference between prior and posterior fluxes in this period is about 0.4 to 0.5 PgC yr −1 . According to MODIS and DIFFUSE GPP estimates, the stronger posterior sink observed in this period may be due to an underestimation of GPP simulated by the CABLE BIOS3 model (Fig. 15). The GPP estimated by MODIS from June to September was about 0.90 PgC yr −1 compared to CABLE BIOS3, which was 0.59 PgC yr −1 . DIFFUSE GPP estimates were around 0.68 PgC yr −1 .

560
Over the sparsely vegetated region, the seasonal discrepancy between the prior and posterior flux is more evident than for savanna. The seasonality of the posterior flux is stronger (relative to the prior estimate) from April to September. In this ecosystem, the largest absolute difference between the prior and posterior fluxes is seen from June to August (0.5 PgC yr −1 ). In July, for example, the inversion shifts the prior flux from -0.05 ± 0.09 PgC yr −1 to -0.56 ± 0.06 PgC yr −1 . Again, we can see in Fig. 14d that a possible reason of this shift may be associated with an underestimation of the GPP by CABLE BIOS3.

565
It is evident that MODIS and DIFFUSE GPP have a stronger seasonality compare to CABLE BIOS3 GPP. For example, from June to August, the CABLE BIOS-3 GPP was about 0.4 PgC yr −1 compared to DIFFUSE and MODIS, which were 0.9 and 1.3 PgC yr −1 respectively. We did not find a seasonal correlation between the prior fluxes and MODIS and DIFFUSE GPP fluxes (see Appendix D, Table D1), but we did find a positive correlation between the posterior fluxes and the GPP estimated by MODIS and DIFFUSE (R = 0.44 and R = 0.45 respectively).

570
Regarding the Tropics, warm temperate, cool temperate, and Mediterranean ecosystems, the seasonal correlation MODIS or DIFFUSE GPP estimates were stronger for the prior than for the posterior fluxes (Appendix D, Table D1). A stronger correlation between the prior flux and MODIS and DIFFUSE might be attributable to the fact that the assimilated coastal fluxes might be somehow less constraint by the inversion in these ecosystems, mainly because they are mostly influenced by ocean fluxes where the uncertainties have less freedom to be modified by the inversion.

Comparison between CABLE BIOS3 GPP and the GPP anomalies derived from MODIS and DIFFUSE model
We saw in the previous section that over the sparsely vegetated and savanna ecosystems, the seasonal cycle of our posterior flux estimate shows a larger amplitude compared to the prior estimate. We suggested that the amplitude difference between these fluxes may be due to an underestimate of GPP by the CABLE BIOS-3 model. In this section, we investigate whether the CABLE BIOS3 model is also underestimating the GPP flux anomalies. We did this analysis because even though the prior 580 and posterior flux anomalies, in general, follow similar patterns in each bioclimatic ecosystem, the amplitudes of these two flux anomalies estimates show discrepancies. Fig. 16 shows the time series of GPP flux anomalies for the period 2015 to 2019 classified by the same six bio-climate ecosystems described in Section 3.3.1 against the 3-month running mean of prior and posterior flux anomalies.
Over the savanna ecosystem (Fig. 16b), we see that the negative posterior flux anomaly shows a larger amplitude than the Over the sparsely vegetated ecosystem (Fig. 16f), we found that the GPP anomalies from CABLE BIOS3, MODIS and 600 DIFFUSE show similar patterns. However, we found that the correlation between the posterior anomalies and CABLE BIOS3, DIFFUSE and MODIS GPP anomalies were stronger than the one correlated with prior. For example, the correlations between the prior and CABLE BIOS3 GPP anomalies were -0.46 compared to the posterior (R = -0.61). Correlations between the posterior anomalies and DIFFUSE and MODIS GPP was also stronger than the prior, which value was -0.5 compared to the prior (R = -0.3) (more details in Appendix D, Table D2). These findings are significant for Australia because they suggest that 605 our OCO-2 inversion might likely be better at capturing the anomalies of this ecosystem (the largest ecosystem in Australia) compared to the biosphere land model ecosystem.
Overall, we found that the savanna and sparsely vegetated ecosystems dominate the inter-annual variability for the period 2015-2019. These results are consistent with Poulter et al. (2014) and Ahlström et al. (2017), who found that semi-arid ecosystems such as Australia are highly responsive to climate drivers such as rainfall and temperature. The result supports the broader view of carbon-cycle/climate interactions than the original focus on the wet tropics (e.g. Rayner and Law, 1999;Bousquet et al., 2000).

Comparison with the OCO-2 MIP global inversion
In Sect 3.1, we saw that validating the posterior concentrations against the current Australian GHG monitoring system is challenging, mainly because Australia has a small number of monitoring sites, especially given its large landmass. In addition, We also studied the carbon flux anomalies derived by the OCO-2 MIP and compared them with the prior and posterior flux anomalies (3-month running mean) that we have discussed throughout this study (Fig. 18). We see in Fig In terms of carbon release, we observe that almost all OCO-2 MIP inversions agree that the largest outgassing occurred in 635 November 2015. Our inversion places the maximum outgassing during the MIP study period in October 2017.
In summary, In summary, analysis of the inter-annual (peak-to-peak) variability shows that PCTM (3.05 PgC yr −1 ) and CSU (2.46 PgC yr −1 ) produce the largest amplitude of variability compared to our prior (0.70 PgC yr −1 ) and posterior anomalies (1.89 PgC yr −1 ). We also note that UO and AMES exhibit the lowest carbon amplitude of the variability, which values are 1.06 and 1.86 respectively. The disagreement between the global inversion and our study might also be driven by transport 640 differences Schuh et al., 2019). As a final discussion, we could say that previous inversion studies over Australia have been limited by the lack of in situ data.
OCO-2 certainly allows a quantum leap in resolution but this is still reasonably coarse, especially when one remembers that the prior covariance structures we use impose smooth variations up to the correlation length of 500 km. Instruments with scanning geometries which allow higher resolution observations such as OCO-3  may improve significantly the 645 available resolution of fluxes. This is particularly important when assessing the roles of drivers such as rainfall which may vary Table 3. Summary of the peak-to-peak amplitude of 3-month running mean posterior (black line) and prior (grey line) terrestrial carbon flux anomalies and 3-month running mean anomalies of nine different global transport models (Units PgC yr −1 ). on smaller scales. We also note continuing improvement in the OCO retrievals themselves which should allow joint assimilation of land and ocean measurements, hopefully improving the visibility of coastal fluxes and improving comparison with coastal in situ measurements such as Cape grim and Gunn Point as shown by Villalobos et al. (2021).

650
We estimated monthly carbon fluxes over Australia for 2015 to 2019, based on the assimilation of the Orbiting Carbon Observatory-2 (OCO-2) satellite data (land nadir and glint data, version 9). We investigated the effect of vegetation productivity (EVI anomalies as proxy) and climate driver variations such as rainfall and air temperature on the Australian terrestrial carbon flux variability. The mean of our five-year inversion suggests that Australia was a carbon sink of -0.46 ± 0.08 PgC yr −1 driven partly by large carbon uptake (-1.04 PgC yr −1 ) recorded in 2016 over the savanna and sparsely vegetated ecosystems. 655 We found that negative carbon flux anomalies recorded in this period over these ecosystems coincide with an increase in the vegetation greenness (positive EVI anomalies) driven by higher than average rainfall anomalies and lower than average air temperature anomalies. The 2017 sink over Australia also contributed to the 2015-2019 long term mean, but its contribution was not as significant as 2015 and 2016. Negative carbon flux anomalies recorded in 2017 also coincided with positive rainfall anomalies and temperatures below average in that period over areas with sparse vegetation. In 2018 we did not find significant 660 terrestrial flux anomalies across Australia, and 2019 mainly was affected by positive carbon flux anomalies, which also were in line with a deficit of rainfall and positive temperature anomalies.
Regarding validation of our inversion with independent data, we found it challenging to validate our posterior columnaveraged concentration with the current Australian monitoring sites. Despite the fact that for several periods between 2015-