Global atmospheric CO₂ inverse models converging on neutral tropical land exchange, but disagreeing on fossil fuel and atmospheric growth rate

2.1 Participating models

The list of participating inverse models is shown in Table 1 and more details are available in the Supplement. These include 10 different inverse modeling systems or system variants. The inversion systems differ in many aspects such as transport models, wind fields, analysis procedures, and subset of assimilated observations. The atmospheric chemistry-transport model (ACTM) system performed two inversions with different prescribed fossil fuel emissions (Saeki and Patra, 2017), one based on totals from the Carbon Dioxide Information Analysis Center (CDIAC; Boden et al., 2016) and another based on the International Energy Agency (IEA/OECD, 2016), which allows us to assess sensitivity to the FF prior only. This is also the case for the two CarbonTracker Europe versions, CTE2016-FT (Fast Track) and CTE2017-FT, where only the subset of observations and the FF prior are different (van der Laan-Luijkx et al., 2017).

It is worth noting that some inverse models are constructed in a similar framework. Some share the same transport model, such as TM5 that is used in four inversions, and some use the same meteorological fields. Five inverse systems nudge their forecast field to the ERA-Interim reanalysis (Dee et al., 2011). The two longest flux estimates, from CAMS (v16r1) and Jena (s85_v4.1), are used to reproduce the comparison with observations as in Stephens et al. (2007) over the T3L2 period (1992 to 1996). The Jena s85_v4.1 and s04_v4.1 inversions differ in their calculation periods and station sets used: Jena s85_v4.1 starts in 1985 using only 23 stations that cover this entire period, while s04_v4.1 uses many more sites (59) and starts in 2004. This also allows us to separate the impact of the number of sites assimilated over the most recent period.

Table 1List of the inverse modeling systems used in this study and general characteristics.

^a With updates documented at https://atmosphere.copernicus.eu/ (last access: 7 January 2019). ^b With updates documented at https://carbontracker.noaa.gov (last access: 7 January 2019).

Download Print Version | Download XLSX

2.2 The Global Carbon Budget 2016

The Global Carbon Project (GCP) gathers observational and model-based flux estimates from multiple organizations and research groups around the world to yearly report a global budget of atmospheric CO₂ (Le Quéré et al., 2016). GCP2016 is the most recent version with flux estimates forced to balance globally. The most recent version (GCP2017; Le Quéré et al., 2018) separated an explicit unknown ocean or land flux term, which prevents simple comparisons of the type presented here. Specifically, the land–ocean partitioning in GCP2016 is based on multiple observational constraints on the ocean flux for the 1990s, extrapolated forward with a suite of seven global ocean models. As pointed out in Le Quéré et al. (2018), there are considerable uncertainties in this extrapolation, with the estimated ocean–land partitioning for later decades dependent on the models. The GCP2016 atmospheric growth rate is derived from atmospheric CO₂ measurements at marine boundary layer (MBL) sites made by the US National Oceanic and Atmospheric Administration (NOAA) Earth System Research Laboratory (ESRL; Masarie and Tans, 1995; Dlugokencky and Tans, 2018). CO₂ emissions from land-use change (ELUC) are the net sum of all anthropogenic activities: deforestation, afforestation, logging, and shifting cultivation. Total emissions are estimated, following the bookkeeping method (Houghton, 2003; Houghton et al., 2012), with complementary interannual variability calculated from satellite data when available (Giglio et al., 2013; van der Werf et al., 2010). The average ELUC for the year 2009 to 2011 included here is estimated to be 0.85 PgC yr⁻¹ with an uncertainty of 0.5 PgC yr⁻¹. These emissions are added to the GCP2016 land sink for comparison to atmospheric inversion estimates.

Finally, the land sink is estimated in GCP2016 as a residual from all other components of the carbon budget. The GCP2016 method treats the riverine flux of carbon from land to ocean to atmosphere as separate components of the total air–land and air–sea fluxes and subtracts an estimate of this flux (0.45 PgC yr⁻¹; Jacobson et al., 2007) from the pCO2-based sea-to-air flux estimates to match estimates of the anthropogenic ocean sink alone. Because the land sink is a residual, this increase in the magnitude of the ocean sink results in a corresponding reduction by 0.45 PgC yr⁻¹ in the magnitude of the land sink in GCP2016. To compare to atmospheric inverse flux estimates, which represent the total air–sea and air–land fluxes, we have adjusted the GCP2016 ocean and land flux estimates by this same 0.45 PgC yr⁻¹, decreasing the ocean sink and increasing the land sink.

Note that we do not show GCP2016 estimates here as a truth metric against which to evaluate the models, but rather as one estimate of an internally consistent global budget that provides a useful reference for exploring axes of variability in our models and comparing to other community estimates.

Figure 1Reconstructed annual cycle in northern extratropical vertical CO₂ gradients, obtained from fits using two harmonics of the HIPPO data and correspondingly sampled model outputs, averaged over 20 to 90^∘ N (1000 to 800 hPa minus 800 to 400 hPa). The CO₂ average curtain observations for each of nine atmospheric transects have been added on the graph to illustrate the data uncertainties and temporal coverage, the y-axis error bar is derived from the range of disagreement among the three in situ instruments on board (QCLS, OMS, and AO2; see Supplement), and the line average is derived from the CO2.X merged dataset. The horizontal whiskers represent the time span of the flights contributing to each average. The observed line shown here is not a direct fit to the observation points, but rather comes from an average of fits to individual 100 hPa by 5^∘ latitude bins as described in the text.

Download

2.3 HIPPO observations and fitting procedures

The HIPPO project (Wofsy, 2011) used the NSF/NCAR Gulfstream V aircraft (GV) to conduct 5-month-long campaigns in different seasons over 3 years (2009–2011; see Supplement) that consisted of vertical profiling along North–South-Pacific transects between 87^∘ N and 67^∘ S. The five campaigns included nine transects of the NET Pacific. We exclude observations over North America conducted between Colorado and Alaska (Fig. S1 in the Supplement). HIPPO flew three different in situ CO₂ instruments and two whole air samplers with laboratory CO₂ measurements. We use the recommended CO2.X variable which comes primarily from the Harvard quantum cascade laser spectrometer (QCLS), gap filled during calibration sequences, and compare to the other systems to constrain potential systematic biases (see Supplement). We calculate the NET vertical gradient as the difference between the average from 1000 to 800 hPa for the lower troposphere (LT) and the average from 800 to 400 hPa for the upper troposphere (UT), spanning the latitude range from 20 to 90^∘ N. To do this, we first detrend the observations and model sampled along the flight-track output by subtracting a deseasonalized and smoothed long-term trend record from the fit of the Mauna Loa Observatory in situ measurement time series to provide a common reference for both observations and models, and we bin the observations by 100 hPa in pressure and 5^∘ in latitude bins. We then fit each bin with a curve using two harmonics and constant offset, and we average the resulting fits across boxes and pressure levels, with latitude weighting (see Supplement). Figure 1 shows the resulting daily fit of the annual cycle for the HIPPO observations and model simulations of the NET vertical gradient. Qualitatively, it shows that most models reproduce the CO₂ cycle well, with positive gradients in winter over a broad peak and negative gradients in summer over a narrower trough. The three CarbonTracker inversions (CT2016, CTE2016-FT, and CTE2017-FT) have somewhat lower seasonal gradient amplitude, while the two ACTM inversions (ACTM-IEA and ACTM-CDIAC) show larger amplitude. More quantitative details are given in Sect. 3.1. To illustrate the temporal coverage of the observations, we plot the measurements of the nine HIPPO transects in Fig. 1 as simple differences of the latitude-weighted average concentrations within the LT and UT boxes for each transect, while an example of a fit to an individual bin is shown in Fig. S1.

The QCLS instrument has a 1σ precision of 20 ppb (Santoni et al., 2014), and for all five CO₂ systems on the GV the instrumental precision is negligible for the large-scale average metrics we present here. More relevant sources of uncertainty are associated with the potential for altitude-dependent biases that might result from inlet or cabin-pressure effects, as well as misrepresentation of synoptic transport in the models. We estimate (i) uncertainty in the annual-mean NET vertical gradient metric by comparison of the five independent instruments and whole air samplers to be ± 0.15 ppm (see Supplement) and (ii) uncertainty on the individual HIPPO transect values to range from 0.02 to 0.48 ppm as shown by the vertical bars in Fig. 1. These values are derived from the maximum absolute differences between the sensors, which we conservatively treat as best-guess 1σ uncertainty estimates. These uncertainty estimates correspond to the vertical gradient as observed by the HIPPO flight tracks and calculated with the fitting procedure used here. Because we use model output along the flight tracks and treat model output and observations identically in our calculations, we do not include an estimate of potential spatial sampling bias, but we do use model output to assess the spatial representativeness of our calculated metrics with respect to full 150^∘ W transect and full zonal means in Sect. 4 of the Supplement (Figs. S5, S6). Also, because the models are driven by reanalysis winds, they should capture the position of synoptic systems and associated transport. However, the wind fields and model transport may be biased, which could result in different vertical gradients for reasons unrelated to the fluxes of interest. We have estimated synoptic variability in the vertical gradient metric and find a worst-case potential model synoptic sampling bias of ±0.06 ppm for the annual mean, ±0.14 ppm for JFM, and ±0.15 ppm for JAS (1σ; see Supplement).

3.1 Fluxes and posterior CO₂ comparisons with HIPPO

Each individual inversion system adjusts fluxes to fit the concentration fields with its given transport scheme and a priori source and sink information. Biases can appear in the retrieved posterior CO₂ resulting from errors in the estimated fluxes or from specific biases in transport to the location of the independent data (here in particular vertical transport to the upper atmosphere). We first evaluate if the spread of retrieved land fluxes over different zonal bands is correlated with NET vertical CO₂ gradients and if the modeled gradients match observations, as was previously done for the T3L2 models by Stephens et al. (2007).

Figure 2a presents the results for the HIPPO and model vertical gradients and model fluxes, broken into NET and T+SET regions for the years 2009–2011. The mean and relative spread of 10 simulations for the posterior annual mean NET land flux is −2.24 PgC yr⁻¹ ±0.29 PgC yr⁻¹ (13 %, 1σ). Aside from the ACTM-IEA simulation, all models are within the uncertainty range of 0.15 ppm or 50 % of the measured vertical gradient. This contrasts to the TransCom 3 Level 2 simulations which had an annual mean of −2.42 PgC yr${}^{-\mathrm{1}}\pm \mathrm{1.05}$ (43 %) PgC yr⁻¹ for NET land flux and disagreed with the observed vertical gradient by ∼0.5 ppm on average and as much as 1.3 ppm (186 %). As listed in Table 1, the inversions have significant differences in transport model, resolution, and driving meteorology and are converging despite these differences. In addition, there are no apparent relationships between vertical gradients and NET nor T+SET land fluxes. The standard deviation across 10 simulations on the difference between NET land and T+SET is 0.4 PgC yr⁻¹ while it was 2.1 PgC yr⁻¹ in T3L2 (Gurney and Denning, 2013; Gurney et al., 2004) and 1.28 PgC yr⁻¹ in RECCAP (Peylin et al., 2013), representing a steady and dramatic convergence of model estimates over the past 15 years. We reproduce the Stephens et al. (2007) annual mean figure in Fig. 2b, with the exception of showing T+SET instead of T, to highlight those differences. It is important to note that these results correspond to a different period and different models, with a smaller network of assimilated in situ network measurements and assimilation of monthly mean rather than discrete measurements. We took advantage of the two models that span the 1992–1996 period, CAMS (v16r1) and Jena (s85_v4.1), to further investigate differences from the T3L2 period. Those two models are quite close to the 2009–2011 vertical gradient observations (Fig. 2a), but they both overestimate the 1992–1996 vertical gradients (Fig. 2b). Notably, they fall along the lines fit to the T3L2 models in Fig. 2b, which could be a coincidence, but might also suggest that despite agreeing with the other models on the latitudinal flux distribution for 2009–2011 these models overestimate tropical sources and northern sinks during 1992–1996. This would require that these models be more dependent on vertical mixing biases in the earlier period. The different number of assimilated sites is one potential factor that might explain different biases in retrieved fluxes for these two periods, but this is not seen for the comparison of the two versions of the Jena model assimilating different numbers of sites during 2009–2011. It is worth noting that reanalyses of meteorological observations have noticeably improved thanks to a better representation of unresolved processes in global models, improved data assimilation methods, and the increasing availability of satellite data, which makes the reanalyses perform better in the 2000s than for the 1990s and earlier (Bauer et al., 2015; Gelaro et al., 2017). As an example, the assimilation of new observations from the constellation of COSMIC global positioning system radio occultation (GPSRO) satellites has led to a significant improvement in meteorological analyses and forecasts (Healy, 2008).

One concern is the spatial representativeness of the HIPPO measurements which were made over the Pacific Ocean while the light-aircraft observations used by Stephens et al. (2007) were mostly measuring profiles over land. We discuss this issue in the Supplement and show that across models HIPPO vertical gradients are significantly representative of the zonal mean for the 3-year mean and every year individually (Fig. S5). Seasonally (Fig. S6), it appears that the vertical gradients are representative of the parallel 150W for winter (JFM), spring (AMJ), and fall (OND) seasons, representative of the zonal mean for winter (JFM) and fall (OND), and representative of the zonal average over land only in boreal summer (JAS). We did find a significant correlation between vertical gradients defined by the HIPPO flight tracks and land zonal means during summer (JAS), when vertical gradients are weak.

Figure 2c and d show the vertical gradients and fluxes for 2009–2011 winter (JFM) and summer (JAS). The agreement between the models and HIPPO observations is not as strong as for annual means. The vertical gradient in the NET winter is reasonably well reproduced by nine models with differences lower than 0.36 ppm. The ACTM-IEA inversion is an outlier and overestimates by 0.94 ppm the winter season average vertical gradient. For ACTM, the global annual IEA emissions are less than CDIAC (Fig. 4c and d), which results in a weaker northern extratropical sink (Figs. 2s and 3a) that corresponds with a more positive LT–UT northern extratropical gradient (Figs. 2a and S2) and a more positive N–S gradient (Fig. S2), comparing just the two ACTM versions. Differences across inversion systems in Fig. S2 also depend on the transport and inversion scheme and the resulting spatial distribution of sources and sinks.

There are generally larger differences between observed and modeled vertical gradients in Northern Hemisphere summer (JAS), with only two models (ACTM-IEA and CAMS) within observation error bars, but the whole range of values is only 0.75 ppm. In this case a linear relationship (r²=0.4) is found between the modeled vertical gradient and the retrieved T+SET fluxes, but not for the NET flux. There is a significant relationship between HIPPO and the land-only zonal average vertical gradient and both are correlated with the T+SET fluxes (Fig. S7), but with a slope of 2.16 ppm yr PgC¹ for HIPPO while it is 0.93 ppm yr PgC¹ over land where the vertical gradients are bigger. This suggests that transport errors may be more critical in the summer season or that other factors compensate to obscure the relationship for these relatively coarse time averages in other seasons and for the annual means. While additional insights into model behavior could be gained from more detailed comparisons to individual models or in more controlled inversion ensembles, the varied nature of these inversion systems makes detailed analyses more challenging and beyond the scope of our current study.

For the annual means and winter there are no statistical relationships between the vertical gradients and the retrieved fluxes. This suggests that Northern Hemisphere vertical mixing errors do not play a major role in biasing the flux estimation across these models. However, the retrieved fluxes can still be biased because of the transport errors.

One potential limitation in our analysis could be the use of similar meteorological fields from the ECMWF base analysis and forecast cycle, which is the case for 5 out of 10 simulations. A careful comparison of model transport suggests that nudging to a particular reanalysis product does not imply identical tracer transport between the models (Locatelli et al., 2015; Orbe et al., 2017; Prather et al., 2008). The transport errors arise from resolved advection and heavily parameterized transport schemes such as convection and boundary layer mixing (Krol et al., 2018; Locatelli et al., 2015; Orbe et al., 2017). Qualitatively, we cannot distinguish the CO₂ vertical gradient from models using ERA-Interim winds from the five other models.

Figure 2Retrieved fluxes versus NET vertical gradients. (a) Annual mean NET land and T+SET land fluxes versus posterior NET vertical gradients (lower minus upper troposphere) from model output along HIPPO flight tracks and HIPPO observations (pink line) for the period 2009–2011. The shaded area represents an estimate of measurement uncertainty of ±0.15 ppm for the annual mean, as estimated in the Sect. S2 in the Supplement. Inverse model posterior concentration gradients and fluxes are shown as points (squares represent NET; triangles represent T+SET). The vertical axis represents the integrated annual mean land fluxes (PgC yr⁻¹). (b) Same as (a) but for 1992–1996 and showing TransCom 3 Level 2 models and our two current models that span this time period, showing dependence of posterior fluxes on transport and a large range of transport biases. Annual mean NET (red squares) and T+SET (blue triangles) land carbon fluxes for 1992–1996 estimated by the 12 T3L2 models plotted as a function of the models' post-inversion predicted mean vertical CO₂ gradients at 10 light-aircraft profiling sites (adapted from Stephens et al., 2007) with fluxes partitioned by TransCom region. The Jena (s85_v4.1) and the CAMS (v16r1) simulations have also been sampled at the same light-aircraft locations but their fluxes are partitioned at 20^∘ N and 20^∘ S. The crosses show our new best estimate of the fluxes estimated by the regression of all T3L2 models. The error bars on these points are estimated using the standard error of the regressions. (c) Same as panel (a) for January–February–March (JFM), and (d) same as panel (a) for July–August–September (JAS). For the seasonal plots, the width of the pink bar is 0.07 ppm for JFM and 0.17 for JAS. In panel (d), the black line represents the regression line, shown because the relationship is statistically significant at a 95 % confidence interval.

Download

3.2 The latitudinal distribution of retrieved land fluxes

In this section, we present the retrieved land flux partitioning between the NET and the T+SET, as shown in Fig. 3 and Table 2. Because the total sink is the sum of T+SET and NET, these lines have a slope of −1 and any deviation perpendicular to the lines indicates disagreement on the total land sink from the GCP2016 estimate. As noted in the previous section, inverse modeling results for the HIPPO period (2009–2011) are remarkably close to one another (Fig. 3a). These results converge on a NET land sink value slightly larger than 2 PgC yr⁻¹ ($-\mathrm{2.24}\pm \mathrm{0.29}$ PgC yr⁻¹) and a T+SET land sink of$-\mathrm{0.38}\pm \mathrm{0.31}$ PgC yr⁻¹. In Fig. 3, multi-model means are represented by blue diamonds and associated error bars are estimated by the standard deviation across models. The 2009–2011 period is marked by a large tropical land sink because of the strong La Niña event of 2011 (Bastos et al., 2013; Poulter et al., 2014). For these 3 years, the models clearly indicate a negative flux over the tropics and SET land. There are also increasing lines of evidence that the rate of deforestation and climate stress over tropics have been moderated in recent decades (e.g., 2000s), compared to the 1990s (Kondo et al., 2018), with a reduced change in tropical forest cover because the decrease in the South American deforestation has been compensated for by an increased Southeast Asian deforestation (Hansen et al., 2013).

In order to place these recent flux estimates in the context of previous studies, we show the flux estimates by the new models that also estimate fluxes for the earlier periods; two models have available outputs for the T3L2 period (1992–1996) and four for the RECCAP period (2001–2004), as shown in Fig. 3b and c. For Jena, one inversion (s85_v4.1) starts in 1985 and is constrained by only 23 atmospheric sites while the other (s04_v4.1) starts in 2004 and uses 59 sites. Interestingly, the difference between s85_v4.1 and s04_v4.1 for 2009–2011 is rather small (Fig. 3a), less than 0.15 PgC yr⁻¹.

According to GCP2016, the total land sink in 2009–2011 was around twice as large (around 3 PgC yr⁻¹) compared to that for 1992–1996 (around 1.7 PgC yr⁻¹) and 2001–2004 (around 1.3 PgC yr⁻¹). This is due to the combined effect of natural interannual variability as well as a long-term trend (Ballantyne et al., 2012). The retrieved total land fluxes for all study periods appear to be close to the corresponding GCP estimates with most models falling within the GCP2016 1σ uncertainty range. For the 2001–2004 period, the newer simulations move fluxes parallel to the GCP line in the direction of a weaker tropical source and a weaker NET sink relative to the original RECCAP estimates. For the 1992–1996 period, one of the two newer simulations shifts fluxes in that same direction, but not as far as suggested by Stephens et al. (2007).

However, we have revisited the Stephens et al. (2007) estimates, by considering the intercept of the regression lines with the aircraft observations rather than the mean of the three models nearest the annual mean observations and evaluating the error using the standard error of the linear regressions. The selection of three models by Stephens et al. (2007) was somewhat arbitrary as they did not directly overlap the observations and did not agree as well as other models seasonally. This new approach relying on the correlated signal from all models leads to a NET flux of $-\mathrm{1.7}\pm \mathrm{0.59}$ PgC yr⁻¹ and a T+SET flux of 0.15±0.66 PgC yr⁻¹, a similar shift in NET fluxes but only two-thirds of the shift in T+SET fluxes using the Stephens et al. (2007) subset of models, as shown in Fig. 3b.

For the RECCAP period, we used their Group 1 simulations (JENA, LSCE, MACC-II, CT2011_oi, CTE2013) identified in Peylin et al. (2013), four of which assimilated the observations at the sample time as opposed to using monthly means and all of which solved for fluxes at the resolution of the transport model or for small ecoregions over land. The T+SET flux estimate averaged over the RECCAP Group 1 models is 0.34±0.27 PgC yr⁻¹. This is nearly identical to the average of the new models from this study (0.34±0.27 PgC yr⁻¹; using CTE2016-FT, CTE2017-FT, CT2016, CAMS v16r1, and Jena s85_v4.1). Both estimate slightly positive T+SET fluxes that are only half of the RECCAP all-model average (0.93±0.90 PgC yr⁻¹). Our NET land sink estimates using newer models are less than the previous estimates in the original T3L2 and RECCAP studies for the 1992–1996 and 2001–2004 periods. Conversely, our new estimates suggest a change in the T+SET flux towards greater uptake and/or less emission for these periods; we found a decrease in the T+SET land flux by 0.71 PgC yr⁻¹ from 0.56±0.32 PgC yr⁻¹ for the 1994–2004 period compared to $-\mathrm{0.15}\pm \mathrm{0.43}$ PgC yr⁻¹ for the 2004–2014 period (Fig. S9). Then, to obtain a flux estimate less sensitive to year-to-year variability we calculate the fluxes for the full 11-year 2004–2014 period (Fig. 3d), for which we have five model estimates. For this longer period, the model spread is largely reduced, in particular for the NET land fluxes, and again we find near-neutral T+SET land fluxes. Taking all four of the estimation periods together (Table 2) all of our central estimates for T+SET are within 0.4 PgC yr⁻¹ of zero. The tropical land fluxes are $-\mathrm{0.2}\pm \mathrm{0.3}$ PgC yr⁻¹ for 2009–2011 and 0.0±0.12 PgC yr⁻¹ for 2004–2014. This implies a consistent uptake of carbon by intact tropical forests over several decades.

Figure 3Tropical and southern extratropical (T+SET) versus northern extratropical (NET) land fluxes for the periods (a) 2009–2011, (b) 1992–1996, (c) 2001–2004, and (d) 2004–2014. The new models used in this study are represented by squares and the average of the available or selected simulations is shown in blue with 1 standard deviation error bars. The pink line and shaded area represents the GCP2016 (river adjusted) estimates of the total land sink for the given period. (a) Results for the HIPPO period 2009–2011; (b) results for the T3L2 period 1992–1996. The TransCom 3 Level 2 outputs (Gurney et al., 2004) are shown in red, with the vertical gradient selected models from Stephens et al. (2007) as circles outlined in green and the rest as red squares outlined in black. The intercept of the regression line with the observed vertical gradient (Fig. 2) is used to define our best flux estimate with error bars estimated by the standard error of the linear regression. (c) Results for the RECCAP period 2001–2004. Also, from Peylin et al. (2013), model means and standard deviations are shown in pink for the subgroup 1 (Jena, LSCEa, MACC-II, CTE2013, CT2011_oi) and in gray for the subgroup 2 (MATCH, CCAM, TrC, NICAM). Panel (d) shows the results from our new set of models for the period 2004–2014.

Download

Table 2Previous and our new best estimates (in bold) of the latitudinal partitioning of land fluxes over four time periods. All values are in PgC yr⁻¹. Values are indicated by the model mean ± 1 standard deviation or 1σ error uncertainties. Regarding the T3L2 period (Gurney et al., 2004), our new estimate for the 1992–1996 period comes from the intercept of the fit lines with the observations in Fig. 2b, and the uncertainties on these values come from the standard error on these metrics from the fits. Regarding the RECCAP period (Peylin et al., 2013), our new estimate for the 2001–2004 period is the average of the five new models from this study.

Download Print Version | Download XLSX

3.3 Variation in retrieved global carbon budgets

Figure 4Synthesis of globally integrated fluxes for the 2009–2011 period, in PgC yr⁻¹. Each inversion is represented by a square and the model mean by a × symbol. The GCP2016 estimates are a pink diamond, which is sometimes hard to see because it is superimposed in each panel by the gray CAMS point. We have adjusted the GCP2016 ocean and land flux estimates by the riverine flux of carbon from land to ocean to atmosphere (0.45 PgC yr⁻¹; Jacobson et al., 2007; Le Quéré et al., 2018), decreasing the ocean sink and increasing the land sink. The magenta line and light-pink shaded area show the corresponding mass balance estimates from GCP2016. In each panel the line and equation shown represent the sum of the x and y variables, and thus the line has a slope of −1, and any deviation perpendicular to the line indicates disagreement on the sum. Here we use the total flux which by mass balance is the whole-atmosphere growth rate (see text), and for panels (a) and (d), the total flux – FF line also equals O + L, while for panels (b) and (c), the total flux line equals O + L + FF. Ellipses denote the variability around the model means of 1σ (darker gray) and 2σ (lighter gray). (a) Ocean versus land; (b) ocean versus land + FF; (c) ocean + land versus FF; (d) total flux versus $-\mathrm{1}\times$ FF.

Download

The global carbon budget partitioning for 2009–2011 is shown for our suite of models and for GCP2016 (river adjusted) in Fig. 4 with the model mean and GCP2016 reported in Table . In every panel of Fig. 4, the light-pink error band shows the constraint imposed by fixing the values to those of GCP2016, and the associated equation is shown on the graph. The pink diamond represents the GCP2016 estimate while the cross and the gray shaded area show the model mean and 1 standard deviation in darker and 2 standard deviations in lighter gray. For the models, the total flux is calculated as the subtraction of the ocean and land sink from the FF emissions. Note that by mass conservation the total flux equals the whole-atmosphere growth rate (WAGR), but that WAGR may differ from the MBL atmosphere growth rate (AGR) defined by surface stations, because of sampling biases or interannual variability in tropospheric mixing or stratosphere–troposphere exchange. GCP uses the MBL AGR (Dlugokencky and Tans, 2018) as an estimate of total flux and assigns uncertainty of ±0.19 PgC yr⁻¹ (Le Quéré et al., 2016) for recent decades, with speculation that the relative uncertainty should decrease when averaging multiple years. Note that, even though the CAMS results systematically align with the GCP2016 estimates in Fig. 4, the two are independent, except for the FF and for the atmospheric data that serve to estimate the total flux in GCP2016. By mass balance, the total annual flux must equal the total growth rate integrated over the entire atmosphere, and this is what we refer to as the total flux.

The integrated ocean versus land fluxes are presented in Fig. 4a. The equation for the range of ocean and land fluxes that would match FF and the total flux estimates from GCP2016 is also shown in Fig. 4a. The models and GCP2016 agree well on the ocean flux with a mean of $-\mathrm{2.04}\pm \mathrm{0.51}$ PgC yr⁻¹ over the 3 years of 2009–2011. The multi-model mean of the land flux is $-\mathrm{2.61}\pm \mathrm{0.42}$ PgC yr⁻¹. The GCP2016 land flux is $-\mathrm{3.04}\pm \mathrm{0.5}$ PgC yr⁻¹ and thus overestimates the model mean. The cloud of model ocean versus land flux estimates are rather scattered around the model mean with a correlation coefficient of only 0.51.

Table 3Global Carbon budget for 2009 to 2011, estimated by the Global Carbon Project 2016 (first row, with river adjustment) and by the suite of models from this study (second row); all values are in PgC yr⁻¹. Values are indicated by the model mean ±1σ error uncertainties, provided by GCP2016 or by the model standard deviation.

Download Print Version | Download XLSX

To better understand the reasons for these discrepancies, and specifically to investigate how much of the land spread in Fig. 4a is a result of differences in fossil fuel priors, we plotted the ocean flux versus the sum of land and FF emissions in Fig. 4b. This figure shows a tight correlation across models for these two parameters (r²=0.93). Given that prior uncertainties specified in the inversions for ocean fluxes are typically smaller than those for land and fossil emissions are fixed, this implies, for a given ocean and FF flux combination, the models are adjusting the land fluxes while matching CO₂ observations. While combining land and fossil fluxes together reduces the random scatter, it does not reduce the range of the continental fluxes, illustrating the fact that models do not simply compensate for biases in fossil priors with land fluxes, but rather that ocean fluxes are affected too (Saeki and Patra, 2017). Conversely, we plot the sum of ocean and land fluxes against FF emissions in Fig. 4c. This figure shows that the ocean + land total sink is largely controlled by the prescribed FF emissions. In general, the models use smaller fossil fuel sources than reported in GCP2016. Figure 4d compares the opposite of FF emissions versus the total flux, again defined by subtraction of the land and ocean fluxes from FF. The spread in models is not parallel to the line defined by the GCP2016 budget closure. We hypothesize that models that overestimate fossil emissions prioritize matching the spatial distribution of CO₂ and thus estimate overcompensating sinks. The spatial patterns of the different FF priors must also play a role, as well as the strength of the atmospheric constraint on annual timescales imposed by the inversion systems.

Taking the two extreme models the ACTM-CDIAC and TM5-4DVar estimates provide very different distributions of fluxes. ACTM-CDIAC suggests stronger land sinks, both over the NET and the T+SET regions, and a lower ocean sink while TM5-4DVar suggests the opposite. This leads to a range of around 2 PgC yr⁻¹ on the model ocean sink. Because of an intentionally different FF source, but with the same inversion system, the ACTM-CDIAC and ACTM-IEA retrieved land fluxes differ by slightly less than 1 PgC yr⁻¹ and ocean fluxes differ by 0.5 PgC yr⁻¹. Overall, this analysis suggests that errors in FF priors are larger than the uncertainty prescribed to them or, more specifically, the range of FF estimates used by leading inversions exceeds the uncertainty that GCP2016 places on the CDIAC estimates. This implies that uncertainties in FF emissions do not adequately consider potential regional biases (Peylin et al., 2011; Saeki and Patra, 2017; Thompson et al., 2016). The large spread of model results along the mass balance line in 4C highlights the need (i) to reduce uncertainty in estimates of FF emissions and (ii) to develop modeling systems that relax rigid FF prior constraints and observational systems that can support optimizing FF emission estimates. For the period 1980–2015, the total flux estimates from GCP2016 are estimated by the MBL AGR of Dlugokencky and Tans (2018). Only background sites that are located in the MBL are used in this calculation. Ballantyne et al. (2012) calculated a sampling error of 0.38 PgC yr⁻¹ (2σ) among the 40 sites and a GCP2017 estimate uncertainty of ±0.19 PgC yr⁻¹ (1σ) for the period 1980–2015 with respect to the total flux. We show the model-retrieved WAGR (equal to total flux) for each individual year in Fig. 5 along with the GCP2016 estimate and error bars. The total spread in the total flux from the inverse models over the 3 years of 2009–2011 equates to 1.38 PgC as shown in Fig. 5b. This is well outside of the uncertainty range estimated for the extrapolation of MBL measurements, implying several inversions are not rigidly constrained to match observed MBL AGR, even over periods of 3 years. Because CO₂ is variably mixed in different years and by different models in the troposphere and between the troposphere and the stratosphere, some inconsistency between the MBL-defined AGR and the total flux of CO₂ in the models might be expected. However, using CT2017 as a test case, the annual difference between the model total surface flux and the observed MBL growth rate over 2000–2016 has a standard deviation of 0.29 PgC yr⁻¹ and for 3-year averages within this period a standard deviation of only 0.10 PgC yr⁻¹, which is much smaller than the discrepancies shown in Fig. 5. Buchwitz et al. (2018) made a similar AGR comparison using CAMS output of total column and surface data and also found good agreement with differences of only ±0.2 PgC yr⁻¹ (1σ) on an annual basis. Another potential challenge to inversions having a consistent total flux during this time may be due to large interannual variability in natural fluxes, with rapid changes resulting from different climatic conditions from the moderate El Niño of 2009 to the strong La Niña of 2011 (Bastos et al., 2013; Poulter et al., 2014). This period has also been marked by rapid changes in emissions, related to lower emissions in 2009 during the financial crisis and a rapid increase in 2010 (Peters et al., 2011). However, Fig. 5 does not indicate that the model total flux estimates for the years 2009–2011 are more divergent than other years. Further work investigating these differences is needed but is beyond the scope of this study. In particular, the length of the assimilation window may have an impact. It may also be possible to force the inverse systems to agree, at least within the MBL, with the observationally defined AGR, and this may help to reduce model spread elsewhere.

Figure 5Modeled total flux (lines), equal to whole-atmosphere growth rate, that is the difference between the global FF emissions and the land and ocean fluxes. Atmospheric growth rate from GCP2016, derived from atmospheric CO₂ measurements made in the marine boundary layer by the NOAA ESRL flask network (Dlugokencky and Tans, 2018; Masarie and Tans, 1995) and GCP2016 assigned uncertainty (pink bands). (b) Shows the sum of the total flux for the 3 years (2009 to 2011).

Download