Global carbon budgets estimated from atmospheric O₂∕N₂ and CO₂ observations in the western Pacific region over a 15-year period

Time series of the atmospheric O₂∕N₂ ratio and CO₂ mole fraction of flask samples obtained from the National Institute for Environmental Studies' (NIES's) flask sampling network are presented. The network includes two ground sites, Hateruma island (HAT; 24.05^∘ N, 123.81^∘ E) and Cape Ochiishi (COI; 43.17^∘ N, 145.50^∘ E), and cargo ships regularly sailing in the western Pacific. Based on temporal changes in fossil-fuel-derived CO₂ emissions, global atmospheric CO₂ burden and atmospheric potential oxygen (APO), which were calculated from the observed O₂∕N₂ ratio and CO₂ mole fraction according to APO = O₂$+\mathrm{1.1}\times {\mathrm{CO}}_{\mathrm{2}}$, we estimated the global carbon sinks of the ocean and land biosphere for a period of more than 15 years. In this carbon budget calculation, we adopted a correction for the time-varying ocean O₂ outgassing effect with an average of 0.54 PgC yr⁻¹ for 2000–2016. The outgassing effect, attributed mainly to global ocean warming, was evaluated under the assumption that the net ocean gas flux is proportional to the change in the ocean heat content for the 0–2000 m layer. The resulting oceanic and land biotic carbon sinks were 2.6±0.7 and 1.5±0.9 PgC yr⁻¹, respectively, for a 17-year period (2000–2016) and 2.4±0.7 and 1.9±0.9 PgC yr⁻¹, respectively, for a 14-year period (2003–2016). Despite the independent approaches, the sink values of this study agreed with those estimated by the Global Carbon Project (GCP) within a difference of about ±0.4 PgC yr⁻¹. We examined the carbon sinks for an interval of 5 years to assess the temporal trends. The pentad (5-year) ocean sinks showed an increasing trend at a rate of 0.08±0.02 PgC yr⁻² during 2001–2014, while the pentad land sinks showed an increasing trend at a rate of 0.23±0.04 PgC yr⁻² for 2001–2009 and a decreasing trend at a rate of $-\mathrm{0.22}\pm \mathrm{0.04}$ PgC yr⁻² during 2009–2014. Although there is good agreement in the trends of the pentad sinks between this study and that of GCP, the increasing rate of the pentad ocean sinks of this study was about 2 times larger than that of GCP.

1 Introduction

In spite of various international efforts to reduce anthropogenic carbon dioxide (CO₂) emissions, the atmospheric CO₂ levels observed around the world have shown a steady increase and exceeded the benchmark of 400 parts per million (ppm) mole fraction in past years (Betts et al., 2016). The Carbon Dioxide Information Analysis Center (CDIAC) reported that global fossil-fuel-derived CO₂ emissions in recent years still increased gradually and rose to 9.9 PgC yr⁻¹ by 2014 (Boden et al., 2017). Under these circumstances, the Paris Agreement adopted at the 21st Conference of Parties (COP) to the United Nations Framework Convention on Climate Change in 2015 aimed to reduce the anthropogenic greenhouse gas emissions to maintain the increase in global mean surface temperatures at well below 2^∘ by 2100 and, if possible, to limit the increase to 1.5^∘. To achieve this goal, it is crucially important to quantitatively understand the natural sink strengths or land biosphere and ocean sinks. A variety of approaches have so far been applied to the quantification of ocean or land sinks or both, including process-based land and ocean models, bottom-up emission estimates based on flux measurements, and top-down estimates based on atmospheric measurements. Developing process-based models to enhance the accuracy of the global carbon budget is crucially important because they are expected to predict the future global carbon cycle in a warmer world. However, carbon budget estimates based on observations are still important in validating and improving the process-based models.

The budget estimation based on atmospheric CO₂ and O₂ observations is a simple and straightforward approach, and it has historically settled the controversy of whether the land biosphere is a net carbon sink or source (Keeling and Shertz, 1992). Although several techniques based on the interferometer (Keeling, 1988), mass spectrometry (Bender et al., 1994), paramagnetic analyzer (Manning et al., 1999), fuel cell analyzer (Stephens et al., 2007), vacuum ultraviolet absorption photometer (Stephens et al., 2003) and so on have been developed to detect the parts-per-million-level changes in the atmospheric O₂ concentration, the accurate quantification of the O₂ change is still challenging. The carbon budget is evaluated by simultaneously solving the mass balance equations of the atmospheric CO₂ and O₂ as follows (Manning and Keeling, 2006):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{2}}=F-B-O,\text{(2)}& {\displaystyle }\mathrm{\Delta }{\mathrm{O}}_{\mathrm{2}}=-{\mathit{\alpha }}_{\mathrm{F}}F+{\mathit{\alpha }}_{\mathrm{B}}B+Z,\end{array}$$

where ΔCO₂ and ΔO₂ represent the changes in the atmospheric CO₂ and O₂ burdens based on atmospheric observations, respectively, F represents the fossil-fuel-derived emissions, and B and O represent the uptake by the land biosphere and the ocean, respectively. α_F and α_B are the −O₂∕C exchange ratio for the globally averaged fossil-fuel combustions and land biotic processes, respectively. The estimated value for α_B is about 1.1 (Severinghaus, 1995), and that for α_F is about 1.4 (Keeling, 1988). These equations mean that the CO₂ and O₂ fluxes associated with fossil-fuel combustion and land biotic processes are tightly coupled. In contrast, the ocean CO₂ uptake, O, and ocean O₂ emissions, denoted as Z, are decoupled because the ocean acts as a carbon sink by physicochemically dissolving the CO₂. Since the values of F and α_F can be evaluated from energy statistics (Keeling, 1988), we can evaluate the ocean and land uptake by solving the above equations if we could evaluate the value of Z.

The global carbon budget can also be related to tracer atmospheric potential oxygen (APO), which is defined by the equation of APO = O₂$+{\mathit{\alpha }}_{\mathrm{B}}\times$ CO₂ (Stephens et al., 1998). Since the APO is defined to be invariant with respect to the land biotic exchange, the secular trend in the APO is determined by fossil-fuel combustions which cause a gradually decreasing trend in APO and the air–sea gas exchange. Combining Eq. (1), multiplied by α_B, and Eq. (2), in accordance with the APO definition, results in the following equation for the APO budget (Manning and Keeling, 2006):

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }\mathrm{APO}=-\left({\mathit{\alpha }}_{\mathrm{F}}-{\mathit{\alpha }}_{\mathrm{B}}\right)F-{\mathit{\alpha }}_{\mathrm{B}}O+Z.\end{array}$$

Since observation sites for atmospheric O₂ are still limited compared with those for atmospheric CO₂, Manning and Keeling (2006) proposed an alternative approach that the global carbon budgets could be obtained by simultaneously solving Eqs. (1) and (3) and using globally averaged CO₂ data based on the measurements of the Global Monitoring Division of the National Oceanic and Atmospheric Administration's Earth System Research Laboratory (NOAA/ESRL/GMD). This approach, making maximum use of the available data, is expected to give the most reliable estimation. This APO approach has been adopted for the estimation of the global carbon budget based on atmospheric O₂ and CO₂ measurements (e.g., Manning and Keeling, 2006; Tohjima et al., 2008; Ishidoya et al., 2012; Goto et al., 2017).

To evaluate the carbon budgets based on O₂ and CO₂ measurements, we need to quantify the magnitude of Z and its temporal variation, if possible. It is considered that Z has a large interannual variability because observed trends of APO generally show large interannual variations which would result in unrealistic variations in the ocean uptake if the variability in Z were rather small (e.g., Bender et al., 2005). Probably, an imbalance of the air–sea seasonal O₂ exchanges, outgassing flux associated with primary production in spring and summer, and ingassing flux associated with ocean ventilation in autumn and winter cause the interannual variations in Z. The results of ocean model simulations also support this mechanism (e.g., McKinley et al., 2003; Nevison et al., 2008). Therefore, it is difficult to estimate the short-term carbon budgets unless the temporal variations in Z are accurately evaluated. Additionally, as for long timescales, it is considered that the present ocean acts as an O₂ source because of the global ocean warming (Keeling and Garcia, 2002). The increase in surface ocean temperature not only reduces the solubility of gases in seawater but also strengthens the ocean stratification, which reduces the ventilation of interior water masses. The reduction of ventilation reduces the ingassing flux of O₂. In contrast, the reduction of ventilation also causes a reduction of the nutrient supply from deep water, which might decrease the primary production and O₂ outgassing in summer. Therefore, the influence of the ocean warming on the net air–sea gas exchange is rather complicated.

Unfortunately, there is little observational evidence to quantify the magnitude of Z and its temporal variations. The long-term average values of Z are inferred under the assumption that Z is proportional to the air–sea heat flux (Keeling and Garcia, 2002). The change in the global ocean heat content has been evaluated based on the large dataset of ocean observations (Levitus et al., 2012). Keeling and Garcia (2002) estimated the ratio of the O₂ flux to heat flux from the relationship between the dissolved O₂ corrected for the mineralization effect and the potential temperature. This approach was basically adopted by most of the studies to evaluate the long-term global carbon budgets (e.g., Bender et al., 2005; Manning and Keeling, 2006; Tohjima et al., 2008). On the other hand, Ishidoya et al. (2012) evaluated the instantaneous variations in the land and ocean sinks based on the APO data at Ny-Ålesund, Svalbard and Syowa, Antarctica, for the period 2001–2009 by calculating the interannual variation in Z from the temporal variation in the ocean heat content. They concluded that the above-mentioned Z values adequately suppressed artifacts caused by the imbalance of the seasonal air–sea O₂ exchange.

We have been conducting air sampling into glass flasks for the measurement of the atmospheric O₂∕N₂ ratio and CO₂ mole fraction at two ground sites in Japan since the late 1990s (Tohjima et al., 2003) and evaluated the global carbon budgets for up to 7 years (1999–2005) based on the APO data from the flask observations (Tohjima et al., 2008). To extend the observation area, we started additional flask sampling aboard commercial cargo ships regularly sailing in the Pacific region in 2002 (Tohjima et al., 2005b, 2012). About a decade has passed since we previously reported the global carbon budgets, and now we have more than a 15-year-long record of atmospheric O₂∕N₂ and CO₂ of the flask samples. In this study, we estimated the ocean and land biotic carbon sinks for over a decade by using the temporal changes in the APO based on these flask data. In addition, we sequentially computed the ocean and land sinks for an interval of 5 years and examined the changing trends of both sinks. In these budget calculations, we estimated the values of Z for the corresponding period by using the temporal changes in the global ocean heat content. Finally, the estimated ocean and land carbon sinks of this study were compared with those of the Global Carbon Project (GCP).

2 Data and analysis

2.1 Flask sampling locations

We started air samplings for the measurement of the atmospheric O₂∕N₂ ratio and CO₂ mole fraction at two monitoring stations located on Hateruma island (HAT; 24.05^∘ N, 123.81^∘ E) in July 1997 and at Cape Ochiishi (COI; 43.17^∘ N, 145.50^∘ E) in December 1998 (Tohjima et al., 2003). In addition, we have been collecting air samples from the Pacific regions by using commercial cargo vessels equipped with automated flask sampling systems (Tohjima et al., 2005b; 2012). The shipboard flask samplings were started between Japan and North America in December 2001, between Japan and Australia and New Zealand in December 2001, and between Japan and southeastern Asia in September 2007. The flask sampling sites are depicted in Fig. 1. Unfortunately, the shipboard data in southeastern Asia, the northern North Pacific (north of 30^∘ N) and the eastern North Pacific are spatiotemporally rather sporadic. (See the inserted figure in Fig. 1 showing time–latitude plots of the shipboard flask samples.) Thus, in the following analysis, we only used the dataset obtained at HAT, COI and the western Pacific region between 40^∘ S and 30^∘ N and between 130 and 180^∘ E.

Air samples were collected in glass flasks hermetically sealed by two glass valves with Viton O-rings. The volumes of the flasks were 2 L for the samplings at HAT and COI and 2.5 L for the shipboard samplings. It should be noted that glass flasks with a volume of 1 L were also used in the early period from the start to March 2006, and only 1 L flasks were used from the start to January 1999 at HAT.

Figure 1 Map showing the air sampling locations in the Pacific region. The light blue and red squares represent the monitoring stations of COI and HAT, respectively. The orange, blue and green circles correspond to the positions where flask samplings were taken aboard cargo ships in southeastern Asia, western Pacific and eastern Pacific, respectively. The inserted figure shows the time–latitude distribution of onboard flask samples. The flask data from COI, HAT and the regions in the black rectangle (130–180^∘ E, 40^∘ S–30^∘ N) were used in the budget calculations.

2.2 O₂∕N₂ and CO₂ analytical methods

In this study, we used a gas chromatograph (GC) equipped with a thermal conductivity detector (TCD) for the measurements of atmospheric O₂ (Tohjima, 2000). In this GC–TCD method, O₂∕N₂ ratios of sample air and working reference air were alternately measured, and the atmospheric O₂ change was determined as the relative difference in the O₂∕N₂ ratio from an arbitrary reference. We used the delta notation according to Keeling and Shertz (1992) to express the relatively small difference in the O₂∕N₂ ratio as follows:

$$\begin{array}{}\text{(4)}& \mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)={\displaystyle \frac{{\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)}_{\mathrm{sam}}}{{\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)}_{\mathrm{ref}}}}-\mathrm{1},\end{array}$$

where subscripts “sam” and “ref” refer to sample and reference, respectively, and the δ(O₂∕N₂) value multiplied by 10⁶ is expressed in per meg units. The change in 1 µmol of O₂ per mole of dry air changed the O₂∕N₂ ratio by 4.77 per meg, which corresponds to a 1 ppm change in the atmospheric trace gas abundance. APO was calculated from the CO₂ mole fraction (X${}_{{\mathrm{CO}}_{\mathrm{2}}}$) in parts per million and δ(O₂∕N₂) in per meg according to

$$\begin{array}{}\text{(5)}& \mathit{\delta }\mathrm{APO}=\mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)+{\mathit{\alpha }}_{\mathrm{B}}{X}_{{\mathrm{CO}}_{\mathrm{2}}}/S_{{\mathrm{O}}_{\mathrm{2}}}-\mathrm{1850},\end{array}$$

where $S_{{\mathrm{O}}_{\mathrm{2}}}$ is the mole fraction of O₂ in the air ($S_{{\mathrm{O}}_{\mathrm{2}}}=\mathrm{0.2094}$, Tohjima et al., 2005a) and the value of 1850 is an arbitrary reference point of δAPO in per meg. The values of δ(O₂∕N₂) were determined against the National Institute for Environmental Studies (NIES) O₂∕N₂ scale (Tohjima et al., 2008). Its temporal stability is examined in the following section.

A nondispersive infrared (NDIR) analyzer (LI-COR, Lincoln, Nebraska, model LI-6252) was used for the CO₂ measurement of the flask samples. The CO₂ mole fractions were determined against the NIES 09 scale, which is based on a set of gravitationally prepared CO₂-in-air standard gases (Machida et al., 2011). The relationship between the NIES 09 scale and the NOAA scale were repeatedly compared through the WMO round-robin intercomparison program. The results showed that the differences of the NIES 09 scale from the NOAA scale were kept within ±0.15 ppm during the period from 1996 to 2014 (https://www.esrl.noaa.gov/gmd/ccgg/wmorr/wmorr_results.php?, last access: 16 July 2019).

2.3 O₂∕N₂ scale stability

As details of the NIES O₂∕N₂ scale are given elsewhere (Tohjima et al., 2008), here we describe only briefly the outline of the scale and add some new information below. The zero point of NIES O₂∕N₂ scale was related to an ambient dry air stored in a high-pressure cylinder (HDA-1). The O₂∕N₂ scale was maintained by three cylinders during the first 4 years (1997–2001) of our O₂∕N₂ measurement program. In 2001, the NIES O₂∕N₂ scale was transferred to 11 other high-pressure cylinders (five 10 L cylinders and six 48 L cylinders), in which the δ(O₂∕N₂) values were carefully determined against the original O₂∕N₂ scale. Another primary reference gas (48 L cylinder, CQB-07080) was added in 2002, and now 12 primary reference gases keep the NIES O₂∕N₂ scale. The air samples delivered from glass flasks or high-pressure cylinders were measured against the working reference airs stored in 48 L aluminum cylinders, from which the δ(O₂∕N₂) values were repeatedly determined against the individual primary gas cylinders at intervals of a few months. The working reference gas cylinders were replaced by new ones every 1–2 years.

Figure 2 Temporal changes in the O₂∕N₂ ratio of primary reference gases relative to the NIES O₂∕N₂ scale. The differences of the O₂∕N₂ ratio from the average are plotted for HDA-2 along with the differences of the O₂∕N₂ ratio from the initial values for the individual cylinders except HDA-2. Solid and dashed horizontal bars in the lower part of the figure indicate the periods when the working reference gases were used.

Figure 2 is the extended version of the previously reported figure (Fig. 1 in Tohjima et al., 2008), showing the temporal changes in the O₂∕N₂ ratio of primary reference gases relative to the NIES O₂∕N₂ scale. In the figure, the deviations of the O₂∕N₂ ratio from the average value for another high-pressure cylinder (HDA-2) and from the initially determined values for the second set of 12 cylinders are plotted. The averages and the standard deviations (1σ) of the differences for the individual 12 cylinders range from −4.2 to 3.3 per meg and from 3.1 to 5.2 per meg, respectively. The changing rates of the deviations for the 12 reference gases during 2001–2017, determined by least square linear regression, range from −0.34 to 0.2 per meg yr⁻¹. Solid and dashed horizontal bars in the bottom of the figure indicate the durations of use of the individual working reference gases.

To assess the stability of the NIES O₂∕N₂ scales, we have continued to measure the reference gases in two 48 L aluminum cylinders (CQB-15645 and CQB-15649) since 2003, which are independent of the reference gases for the NIES scale. The results are shown in Fig. 3. The average changing rates for the whole period, evaluated by linear regression analysis, are $-\mathrm{0.14}\pm \mathrm{0.06}$ per meg yr⁻¹ for CQB-15645 and $-\mathrm{0.05}\pm \mathrm{0.06}$ per meg yr⁻¹ for CQB-15649. Therefore, we conclude that the stability of the NIES O₂∕N₂ scale has been maintained within ±0.2 per meg yr⁻¹, at least during the period of 2003–2016. However, this stability test cannot exclude the possibility that the O₂∕N₂ ratios of the reference gases drift across all the cylinders rather uniformly. There are several mechanisms that affect the O₂∕N₂ ratios of the gases within the high-pressure cylinders, including corrosion of the inner surface, leakage, thermal diffusion and gravitational fractionation. Keeling et al. (2007) assessed carefully and comprehensively the influences of those potential mechanisms on the long-term stability of the O₂∕N₂ ratio of the reference gases and obtained an estimated uncertainty of ±0.4 per meg yr⁻¹. We also treated the reference cylinders, which were kept horizontally in a thermally insulated box, with the greatest care (Tohjima et al., 2008). Therefore, we adopted the value of ±0.4 per meg yr⁻¹ as the long-term drift of the reference gases caused by the above degradation effects. Consequently, we assumed that the total uncertainty of the long-term stability of the O₂∕N₂ reference scale was ±0.45 per meg yr⁻¹ (or $({\mathrm{0.2}}^{\mathrm{2}}+{\mathrm{0.4}}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$) in this study.

Figure 3 Temporal changes in the O₂∕N₂ ratio of reference airs in two aluminum cylinders which are independent of the NIES primary reference gases, relative to the NIES O₂∕N₂ scale. The dashed lines represent the linear regression lines.

2.4 Carbon budget calculation

The ocean uptake and land uptake, O and B, respectively, are given by the following equations (Manning and Keeling, 2006; Tohjima et al., 2008):

$$\begin{array}{}\text{(6)}& {\displaystyle }O=\left[-\left({\mathit{\alpha }}_{\mathrm{F}}-{\mathit{\alpha }}_{\mathrm{B}}\right)F-\left({\displaystyle \frac{S_{{\mathrm{O}}_{\mathrm{2}}}}{\mathit{\beta }}}\right)\times \mathrm{\Delta }\mathrm{APO}+Z_{\mathrm{eff}}\right]\times {\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{B}}}},\text{(7)}& {\displaystyle }B=\left[{\mathit{\alpha }}_{\mathrm{F}}F+\left({\displaystyle \frac{S_{{\mathrm{O}}_{\mathrm{2}}}}{\mathit{\beta }}}\right)\times \mathrm{\Delta }\mathrm{APO}-\left({\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{B}}}{\mathit{\beta }}}\right)\times \mathrm{\Delta }{X}_{{\mathrm{CO}}_{\mathrm{2}}}-Z_{\mathrm{eff}}\right]\times {\displaystyle \frac{\mathrm{1}}{{\mathit{\alpha }}_{\mathrm{B}}}},\end{array}$$

where β is the coefficient converting petagrams of carbon to parts per million CO₂ in the atmosphere (β=0.470 ppm PgC⁻¹; Tohjima et al., 2008), and Z_eff represents the net effect of the air–sea O₂ and N₂ exchange on the atmospheric O₂∕N₂ ratio. In these equations, O, B, F and Z_eff are given in units of petagrams of carbon per year; ΔAPO in units of per meg per year; and $\mathrm{\Delta }{X}_{{\mathrm{CO}}_{\mathrm{2}}}$ in units of parts per million per year. Note that F and α_F include the CO₂ emissions associated with cement manufacturing. The values of α_F were calculated from the CO₂ emission amounts and the −O₂∕CO₂ molar exchange ratios of the individual fuel types (Keeling, 1988). Since α_F slightly varies from year to year, the value of α_F for the relevant period was used for the carbon budget calculations. The values of Z_eff were calculated based on the effects of global ocean warming and anthropogenic nitrogen deposition in accordance with the approach of Keeling and Manning (2014). The details of the Z_eff calculation is discussed in the following section.

We used the same dataset of the fossil-fuel-derived CO₂ emissions and the global average of the atmospheric CO₂ mole fractions as those used by the GCP for the estimation of global carbon budget 2018 (Le Quéré et al., 2018). The fossil CO₂ emissions from fossil-fuel combustion and cement production were basically based on the dataset from the CDIAC and other energy statistics (Boden et al., 2017). The change in the atmospheric CO₂ burden was calculated based on the global observation by the US National Oceanic and Atmospheric Administration Earth System Research Laboratory (NOAA/ESRL; Dlugokencky and Tans, 2018). The temporal variations in the fossil CO₂ emissions and the atmospheric CO₂ accumulation rate are depicted in Fig. 4.

Annual means of APO centered on 1 January were computed by using the same procedure as Tohjima et al. (2008). First, smooth-curve fits to the data were computed in accordance with the methods of Thoning et al. (1989), with a cut-off frequency of 4.6 cycles yr⁻¹. Then the flask APO data were modified to represent the values at the center of the individual months by shifting them in parallel with the smooth-curve fits. This procedure aimed to reduce the influence from biases of the sampling timings within the individual months. The monthly averages were calculated from the modified APO data. When there were no flask data in the monthly time frame, the monthly average of the smooth curve was used. The annual means were calculated from the consecutive 12 monthly averages from July to June of the following year. The standard errors of the differences between the flask data and the smooth-curve fits for the corresponding annual periods were adopted as the uncertainties for the annual averages. The averages and ranges (minimum–maximum) of the errors for the annual means of APO were 0.8 per meg (0.6–1.3 per meg) for HAT and 1.1 per meg (0.9–1.4 per meg) for COI.

In this study, we adopted the value of 1.10 of Severinghaus (1995) for α_B, in accordance with a series of previous studies (e.g., Bender et al., 2005; Manning and Keeling, 2006; Ishidoya et al., 2012; Keeling and Manning, 2014; Goto et al., 2017). However, several studies (e.g., Randerson et al., 2006; Worrall et al., 2013), investigating the elemental compositions of organic matters in soil and plants, indicated that the value of 1.1 is rather large for the globally averaged net −O₂∕CO₂ exchange ratio for the terrestrial biosphere. These studies suggest that the value of 1.05 is much more appropriate for α_B. Although it is beyond the scope of this study to discuss which value is better for α_B, it is useful to mention that the use of 1.05 for α_B results in larger decreasing rates of APO by about 5 % and an increase in land sinks and a decrease in ocean sinks by about 0.06 PgC yr⁻¹ on average in the following results. Considering recent reports about the global net −O₂∕CO₂ exchange ratio, Keeling and Manning (2014) revised the uncertainty of α_B to be increased, from ±0.05 (Severinghaus, 1995) to ±0.10. Thus, we also adopted ±0.10 for the uncertainty of α_B in this study.

Figure 4 Temporal changes in fossil-fuel CO₂ emissions (red squares), atmospheric CO₂ accumulation rate (orange circles) and ocean outgassing effect Z_eff divided by land biotic −O₂∕CO₂ exchange rate (1.1; purple diamonds). The 5-year averages of Z_eff∕1.1 used for the pentad (5-year) carbon sink calculations are also depicted as purple lines.

2.5 Evaluation of outgassing effect (Z_eff)

As is discussed in the Introduction, today's ocean is considered to act as a net source of atmospheric O₂ because of global ocean warming, which also affects the air–sea N₂ exchange. Since the atmospheric O₂ change is measured as the change in the atmospheric O₂∕N₂ ratio, the outgassing effect caused by the global ocean warming, which is denoted as Z_gow, should include the influences from not only ocean O₂ outgassing but also ocean N₂ outgassing. Assuming that the relationship is proportional between the gas fluxes and heat fluxes across the air–sea interface, Manning and Keeling (2006) gave the equation for Z_gow as

$$\begin{array}{}\text{(8)}& Z_{\mathrm{gow}}=\left({\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}-{\displaystyle \frac{S_{{\mathrm{O}}_{\mathrm{2}}}}{S_{{\mathrm{N}}_{\mathrm{2}}}}}{\mathit{\gamma }}_{{\mathrm{N}}_{\mathrm{2}}}\right)\times Q\times m_{\mathrm{c}}\times {\mathrm{10}}^{-\mathrm{15}},\end{array}$$

where Q represents the changing rate of the global ocean heat storage (in units of J yr⁻¹), ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$ and ${\mathit{\gamma }}_{{\mathrm{N}}_{\mathrm{2}}}$ are ratios of the gas flux to heat flux between the air and sea (in units of mol J⁻¹), $S_{{\mathrm{N}}_{\mathrm{2}}}$ is the mole fraction of atmospheric nitrogen ($S_{{\mathrm{N}}_{\mathrm{2}}}=\mathrm{0.7809}$; Tohjima et al., 2005b), and m_c is the atomic mass of carbon (m_c=12.01). Z_gow is given in units of petagrams of carbon per year.

The primary mechanism that affects the air–sea gas exchange is a reduction of gas solubility caused by the increase in the ocean temperature. Therefore, the ratio of the gas flux to heat flux derived from the above thermal effect can be evaluated from the temperature dependence of gas solubility in the seawater and the specific heat of the seawater. Since the air–sea N₂ exchange is predominantly driven by the thermal effect, we adopted the estimated ${\mathit{\gamma }}_{{\mathrm{N}}_{\mathrm{2}}}$ of 2.2 nmol J⁻¹ in this study, in accordance with previous studies (Keeling and Garcia, 2002; Manning and Keeling, 2006).

In contrast to the air–sea N₂ exchange, the changes in the ocean circulation and ocean primary production also affect the air–sea O₂ exchange, as mentioned in the Introduction. Examining the ratio of the seasonal ocean outgassing of O₂ to the seasonal ocean heating and the negative linear relationship between the dissolved O₂ concentrations corrected for ocean biological processes and the potential temperature in the main thermocline based on archived global observation data, Keeling and Garcia (2002) obtained the estimate of 4.9 nmol J⁻¹ for ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$. The value of ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$ was also investigated by using ocean biogeochemical models to revise the global carbon budgets based on O₂ observations (e.g., Plattner et al., 2002; Bopp et al., 2002). Keeling et al. (2010) summarized the model-based values of ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$ ranging from 5.9 to 6.7 nmol J⁻¹. On the other hand, Stendardo and Gruber (2012) examined a huge archived dataset of observations in the northern Antarctic Ocean during the past 5 decades and obtained changing ratios of O₂ inventory to heat content of $-\mathrm{4.3}\pm \mathrm{2.4}$ nmol J⁻¹ in the upper 700 m and $-\mathrm{1.6}\pm \mathrm{1.9}$ nmol J⁻¹ between 700 and 2750 m. These basin-scale ocean O₂ ∕ heat changing ratios seem to suggest that the global ocean acts as a net O₂ source due to global ocean warming. Therefore, in this study, we used the value of 4.9 nmol J⁻¹ for ${\mathit{\gamma }}_{{\mathrm{O}}_{\mathrm{2}}}$ in accordance with previous studies.

To compute Q, we used the estimates of the world ocean heat content (OHC) based on a variety of oceanographic data (Levitus et al., 2000, 2012). Time series of OHC for the 0–700 and 0–2000 m layers are available from the NOAA's website (https://www.nodc.noaa.gov/OC5/3M_HEAT_CONTENT/, last access: 19 February 2019). In previous carbon budget estimations based on atmospheric O₂∕N₂ measurements, the values of Q were estimated from the OHC for the 0–700 m layer (e.g., Manning and Keeling, 2006; Tohjima et al., 2008; Ishidoya et al., 2012). Levitus et at. (2012) showed, however, that the ocean heat storage of the 700–2000 m layer contributes to about one-third of the total heat storage of the 0–2000 m layer. Keeling and Manning (2014) estimated the value of Q by considering not only Q for the depths above 700 m but also Q for the depths below 700, which contributed to 30 % of the total value of Q. Therefore, the time series of OHC for the 0–2000 m layer was used in this study. Note that since the annual average of OHC for the 0-2000 m layer is available only after 2005, we used the pentad (5-year) averages before 2005.

In addition to the ocean warming effect, Keeling and Manning (2014) introduced recently another ocean outgassing effect caused by atmospheric deposition of excess anthropogenic nitrogen to the open ocean. The excess nitrogen is considered to enhance the ocean biotic production of organic matter, which is associated with the O₂ production. Keeling and Manning (2014) evaluated the anthropogenic nitrogen-induced outgassing as being about 0.1×10¹⁴ mol O₂ yr⁻¹. Since the outgassing effect caused by the anthropogenic nitrogen deposition is small but rather significant, we adopted the effect as Z_anthN, with a magnitude of 0.12 PgC yr⁻¹ ($=\mathrm{0.1}\times {\mathrm{10}}^{\mathrm{14}}$ mol O₂ yr${}^{-\mathrm{1}}\times \mathrm{12.01}$ gC mol⁻¹). Eventually, the total outgassing effect, Z_eff, is expressed as the summation of Z_gow and Z_anthN:

$$\begin{array}{}\text{(9)}& Z_{\mathrm{eff}}=Z_{\mathrm{gow}}+Z_{\mathrm{anthN}}.\end{array}$$

The time series of the annual Z_eff divided by α_B are depicted in Fig. 4 as purple lines. The value of Z_eff/α_B ranges from 0.2 to 1.0 Pg yr⁻¹, and the 19-year average for 1998–2016 is 0.52 Pg yr⁻¹. There are not many differences in Z_eff between this study and previous studies (e.g., Manning and Keeling, 2006; Keeling and Manning, 2014; Tohjima et al., 2008; Ishidoya et al., 2012; Goto et al., 2017).

3 Results and discussion

3.1 Observed O₂∕N₂, CO₂ and APO

The time series of the atmospheric CO₂ mole fraction, O₂∕N₂ ratio and APO of the air samples collected at HAT and COI and aboard cargo ships sailing between 40^∘ S and 30^∘ N in the western Pacific are depicted in Fig. 5 together with the smooth-curve fits. The ship data were binned into 10^∘ latitudinal bands (40–30^∘ S, 30–20^∘ S, 20–10^∘ S, 10–0^∘ S, 0–10^∘ N, 10–20^∘ N and 20–30^∘ N). Note that there are no data gaps with more than 50 days in the time series at HAT and COI, while the time series of the cargo ships have data gaps during the 7-month period from October 2006 to April 2007. The ship data during the 7-month period were significantly contaminated by the inboard air due to the failure of diaphragm of the sampling pump.

Figure 5 Time series of the atmospheric CO₂ mole fraction (a), O₂∕N₂ ratio (b) and APO (c) of the flask samples obtained from the NIES flask sampling network shown in Fig. 1. Observed data from COI, HAT and cargo ships operating between 40^∘ S and 30^∘ N were used for the global carbon budget calculation. The time series of CO₂, O₂∕N₂ and APO are offset by 20 ppm, 150 per meg and 100 per meg, respectively, to allow them to be plotted on the same panels. The numbers on the y axis represent the values for the data at HAT.

The temporal variations in the annual means of the atmospheric CO₂, O₂∕N₂ and APO, centered on 1 January, are shown in Fig. 6a, b and c, respectively, where linear trends obtained from least square fitting to the data of HAT were subtracted from the individual time series to emphasize the interannual variations. The standard errors of the annual means for HAT and the 10–0^∘ S bin are depicted as vertical bars for typical examples. Note that the annual means of the atmospheric CO₂, O₂∕N₂ and APO and the corresponding standard errors for HAT, COI and the 10^∘ bins are summarized in Table S1–S9 in the Supplement. The annual means of CO₂ and O₂∕N₂ show latitudinal gradients of northward increase and southward increase, respectively, because the fossil-fuel-derived CO₂ emissions and O₂ consumptions occur predominantly in the northern midlatitudes. In contrast, the highest values of the annual mean APO were generally observed around the Equator, as previously reported (Battle et al., 2006; Tohjima et al., 2005b, 2012). This equatorial peak is mainly attributed to large-scale air–sea gas exchanges: ingassing in the midlatitudes and high latitudes and outgassing in the equatorial region. Conducting atmospheric inversion analyses based on the APO data from the Scripps observation network, Rödenbeck et al. (2008) suggested anomalous outgassing of APO from the equatorial region during El Niño periods, while Tohjima et al. (2015) found a suppressed equatorial peak during El Niño periods based on the western Pacific observations. Eddebbar et al. (2017) reconciled these conflicting results by predicting the existence of a zonal dipole-like El Niño–Southern Oscillation (ENSO) response in the equatorial Pacific based on several ocean process-based models and an atmospheric transport model. These results suggest that enhanced zonal coverage of the atmospheric observations in the equatorial Pacific is needed to constrain the full basin-scale ENSO response. We can see a considerable suppression of the equatorial peak during the strong 2015–2016 El Niño event in Fig. 6c, which was not reported in Tohjima et al. (2015). Any detailed discussion about the temporal variation of the equatorial peak during the 2015/2016 El Niño event is, however, beyond the scope of this study and will be given elsewhere.

Figure 6 Temporal variations of the (a) annual mean CO₂, (b) annual mean O₂∕N₂, (c) annual mean APO and (d) annual changing rate of APO based on the flask samples collected from HAT, COI and cargo ships in the western Pacific (40^∘ S–30^∘ N). The differences in the annual means from the linear trends fitted to the data at HAT are depicted in the figures to emphasize the interannual variations. Vertical bars for the plots of HAT and 10–0^∘ S bin correspond to the standard errors of the annual means.

Figure 6d shows the time series of the annual changes in the annual mean APO, which are the annual changing rates of APO for a 1-year interval (Δt=1 year). As you can see, there are considerable differences in the annual changing rates among the observation sites in the same years; the standard deviations range from 1.6 to 4.4 per meg and the average is 2.8 per meg. We also depict the averages of the annual changing rates of APO of HAT and COI and of all the shipboard data as a thick grey line in Fig. 6d. Note that these average annual changing rates of APO were used for calculation of the global carbon budget in the following sections. The average annual changing rates show also a large interannual variability with a standard deviation of 4.7 per meg yr⁻¹ for the entire observation period.

The differences in the changing rate of APO among the sites in the same years decrease with the increase in the interval for the calculation (Δt) as shown in Fig. 7, where the average (red circles) and the minimum and maximum (red dashed lines) standard deviations of the changing rates are plotted against Δt. The differences among the sites decrease almost inversely with Δt; the average standard deviation for the Δt of 5 years is 0.54 per meg yr⁻¹. The temporal variability in the changing rate also decreases inversely with Δt as depicted in Fig. 7 (blue circles); the standard deviation is reduced to 1.2 per meg yr⁻¹ for the Δt of 5 years. The above results seem to suggest that the temporal variability in the APO fluxes exceeds the spatial variability. As is indicated by Eq. (3), the temporal variability in the APO changing rate should be attributed mostly to those in O and Z_eff. Therefore, the above results also indicate that an interval of 5 years could suppress the temporal variability in Z_eff to the level of ±1.2 per meg yr⁻¹, which corresponds to a carbon budget of about ±0.5 PgC yr⁻¹.

Table 1 Comparison of changing rate of the atmospheric CO₂, O₂ and APO.

^* Average and standard deviation of the changing rates for the three sites (ALT, MLO and SPO) listed in Table 1 of Goto et al. (2017) are given in this table.

Figure 7 Relationship between the standard deviation (SD) of the APO changing rate and the time interval to calculate the changing rate. The red circles represent the averages of the standard deviations of the APO changing rates from the different sites for the same year. The dashed lines represent the minimum and maximum of the standard deviations. The blue circles represent the temporal variability in the average APO changing rate of the different sites.

The changing rates of the atmospheric CO₂, O₂∕N₂ and APO for several combinations of time periods and the observed data (HAT, COI and shipboard) are summarized in Table 1. Here, the uncertainties of the changing rates were computed from the uncertainties of the corresponding annual means at both ends of the periods, the estimated uncertainty of the O₂∕N₂ scale stability (±0.45 per meg yr⁻¹; Sect. 2.3) and the uncertainty in the O₂∕N₂ span sensitivity (±3 %; see below). The time periods of 2000–2010, 2001–2010 and 2001–2014 were selected to compare the observational results of this study with those of the Scripps Institution of Oceanography (SIO) and Tohoku University (TU; Keeling and Manning, 2014; Ishidoya et al., 2012; Goto et al., 2017). As discussed in the above section, the differences in the long-term changing rates of APO between HAT, COI and shipboard data are less than 0.3 per meg yr⁻¹, while the increasing rates of CO₂ and the decreasing rates of O₂∕N₂ for HAT are slightly larger than those for COI and other sites. The monitoring station of HAT is located at the marginal region of continental East Asia, and the anthropogenic CO₂ emissions from China often influence the observations at HAT during winter due to the East Asian monsoon (c.f. Minejima et al., 2012; Tohjima et al., 2010, 2014). Additionally, for the period of 2000–2014, the fossil-fuel-derived CO₂ emissions from China show a rapid increase associated with the unprecedented economic growth. These situations may explain the rather large increase in CO₂ and decrease in O₂∕N₂ at HAT. In contrast to CO₂ and O₂∕N₂, the emissions from fossil-fuel combustion and land biotic processes contribute less to the APO variations, resulting in relatively small differences in the long-term APO changing rates among the sites.

It should be noted that the decreasing rates of APO of our study are 0.5–1.1 per meg yr⁻¹ smaller than those of SIO and TU. Except for the differences of the observation sites, we can offer two explanations for the discrepancy. First, the calculation methods of the changing rate adopted by Goto et al. (2017) are different from those adopted in this study, which might partially explain the discrepancy. This is understandable when comparing the changing rates of CO₂, O₂∕N₂ and APO for the individual studies. In this study, the APO changing rates are almost consistent with those calculated from the CO₂ and O₂∕N₂ changing rates according to the APO definition. However, in the study of Goto et al. (2017), the CO₂ and O₂∕N₂ changing rates of Ny-Ålesund give APO decreasing rates of −9.4 per meg yr⁻¹, which is 0.7 per meg yr⁻¹ smaller than the originally reported values. Second, inter-laboratory comparison of flask samples and high-pressure cylinders suggests a possibility that the span sensitivity of the O₂∕N₂ measurements of NIES is about 3 % lower than that of SIO, which can almost explain the differences in the APO decreasing rates. However, to obtain an accurate conclusion, we need many more studies.

Table 2 Comparison of global carbon budgets based on APO with those from GCP^a,b.

^a Figures are given (in units of per meg yr⁻¹ for ΔAPO, mol mol⁻¹ for α_F and PgC yr⁻¹ for the others). ^b Figures in parentheses represent the uncertainties. ^c ΔAPO is based on the data from HAT, COI and cargo ships (40^∘ S–30^∘ N). The uncertainty in parentheses includes the uncertainty associated with the observations, stability in the O₂∕N₂ scale, uncertainty derived from limited sampling and uncertainty in the O₂∕N₂ span sensitivity (see text). ^d These figures are computed from the dataset summarized by the GCP. The uncertainties are ±0.2 PgC yr⁻¹ for the atmospheric CO₂ and ±5 % for the fossil-fuel emissions, ±0.5 PgC yr⁻¹ for the ocean sinks, and ±0.9 PgC yr⁻¹ for the land sinks (Le Quéré et al., 2018). ^e The uncertainties for α_F are ±0.04 mol mol⁻¹ (Tohjima et al., 2008). ^f The values of Z_eff include both global ocean warming and anthropogenic nitrogen deposition effects, and the uncertainties are assumed to be ±100 % (see text).

3.2 Calculation of global carbon budgets

The rates of the global carbon uptake by the ocean and land biosphere were calculated from the average changing rate of APO based on observations at COI, HAT and on cargo ships (40^∘ S–30^∘ N). The results for several time periods are summarized in Table 2 together with the average changing rates of APO, globally averaged atmospheric CO₂ accumulation rates, fossil-fuel-derived CO₂ emission rates and the ocean outgassing effect divided by α_B. The 19-year (1998–2016), 17-year (2000–2016) and 14-year (2003–2016) periods correspond to the individual maximum observation periods for HAT, COI and the western Pacific, respectively. For example, the estimated ocean and land sinks for 2000–2016 were found to be 2.6±0.7 and 1.5±0.9 PgC yr⁻¹, respectively.

The uncertainties in the parameters used for the carbon budget calculation (Eqs. 6 and 7), which are also listed in Table 2, are briefly discussed here. Note that in this study the estimated uncertainties are ±1σ. Since the ocean outgassing effect is rather speculative, we assumed that the values of Z_eff for the individual periods had ±100 % uncertainties, in accordance with previous studies (e.g., Manning and Keeling, 2006; Tohjima et al., 2008). We adopted uncertainties of ±5 % for the fossil-fuel-derived CO₂ emission rate and ±0.2 PgC yr⁻¹ for the atmospheric CO₂ increasing rate from Le Quéré et al. (2018). As for the uncertainties of the observed APO changing rates, we adopted the standard deviations among the sites shown in Fig. 7 (±0.37 per meg yr⁻¹ for longer than 10 years and ±0.54 per meg yr⁻¹ for 5 years). The estimated uncertainty of the O₂∕N₂ scale stability (±0.45 per meg yr⁻¹) discussed in Sect. 2.3, the uncertainty of the O₂∕N₂ span sensitivity (±3 %) and the uncertainty in the globally averaged APO associated with the limited atmospheric sampling (±0.2 PgC yr⁻¹) discussed in Nevison et al. (2008) were also included in the calculation of the uncertainties in ΔAPO. The uncertainties of α_B and α_F were ±0.10 (Keeling and Manning, 2014) and ±0.04 (Tohjma et al., 2008), respectively. Finally, these uncertainties were propagated to the ocean and land sink uncertainties, in accordance with Eqs. (6) and (7).

We compared our global carbon budget estimations with those of the GCP (Global Carbon Budget, 2018) updated by Le Quéré et al. (2018). In the GCP carbon budget assessment, the ocean and land sinks were estimated by combining multiple results from a variety of models, including global ocean biogeochemistry models (GOBMs) and dynamic global vegetation models (DGVMs). Since the sum of the model-based ocean and land sinks was not necessarily balanced with the difference between fossil-fuel emissions and atmospheric accumulation, Le Quéré et al. (2018) listed the discrepancies as budget imbalances. The ocean sinks, the land sinks which are net land sinks computed as the differences between land uptake and emissions associated with land-use change, and the budget imbalances for the corresponding periods are listed in Table 2. Note that the uncertainties of the sinks of the GCP are ±0.5 PgC yr⁻¹ for ocean and ±0.9 PgC yr⁻¹ for land. The carbon sinks of this study and the GCP for the three long periods are consistent with each other; the largest difference in sink strength is 0.35 PgC yr⁻¹, which is smaller than the uncertainties associated with the individual estimations. A 3 % higher span sensitivity of the O₂∕N₂ measurements, which corresponds to the difference in the span sensitivity between SIO and NIES (Sect. 3.1), would result in an increase and decrease of 0.27 PgC yr⁻¹ in the ocean and land sinks, respectively, for the three long periods. Although these changes would enlarge the differences in sinks between the GCP and this study, the differences are still not significant given the uncertainties of both this study and the GCP.

The carbon budgets for four pentad periods (2000–2004, 2004–2008, 2008–2012 and 2012–2016) are also listed in Table 2. Here, we consider that the 5-year interval effectively reduced the apparent errors caused by the imbalance of the seasonal ocean O₂ fluxes, as discussed in Sect. 3.1. Again, the discrepancies of the pentad ocean and land sinks between this study and the GCP are within ±0.5 PgC yr⁻¹, which is also less than the estimated uncertainties. The land sink during 2008–2012 and the ocean sink during 2012–2016 of this study are about 0.5 PgC yr⁻¹ larger than those of the GCP. These discrepancies in the carbon sinks, partly explained by the rather large values of the carbon budget imbalances of the GCP estimation, might give a clue about how to partition the imbalance values between the land and ocean sinks.

Examining the temporal variations in the pentad sink strengths of this study, we found a gradual increase in the ocean sinks for the latter three pentad periods and a rapid increase and decrease in the land sinks for the former and latter two pentad periods, respectively. The pentad averages of the GCP sinks seem to show similar temporal variations: a steady increase in the ocean sinks for the whole period and a rapid increase and decrease in the land sinks for the former and latter two pentad periods, respectively. These results suggest that the carbon sinks for the pentad periods can be used to evaluate the temporal changes. In the following section, we will examine the temporal change in the carbon sinks in more detail.

3.3 Temporal change in the carbon sinks

Figure 8 shows the temporal variations in the ocean and land sinks for the annual (dashed red lines) and pentad (red lines) intervals calculated from the average of the APO changing rates based on the observations from HAT, COI and cargo ships in the western Pacific. The uncertainties for the pentad sinks (±1σ), which were calculated as described in the previous section, are shown as grey shading. To clearly understand the effect of Z_eff correction, the pentad sinks without corrections are also depicted as purple lines in the figure. Additionally, the annual and pentad sinks of the GCP are also depicted in the figures for comparison. Although the annual sinks show considerable variability, especially for the first several years with peak-to-peak differences of more than 10 PgC yr⁻¹, the variability in the pentad sinks is effectively suppressed. Only the pentad budgets for 2000 show a rather large ocean uptake (3.21 PgC yr⁻¹) and a rather weak land emission (0.56 PgC yr⁻¹), which are depicted as dotted lines. These anomalous values may be explained by the fact that the influences from the considerable drawdown of APO in 2000–2001 cannot be compensated for in the pentad APO changing rate for 2000. Hamme and Keeling (2008) reported that the APO drawdown in 2000–2001, which was also observed in the SIO observations, may be attributed to deep ventilation associated with the unprecedented cooling of the western Pacific, and the variations in the ocean heat content exerted only secondary influence. Therefore, we do not use the anomalous pentad ocean and land sinks for 2000 in the following discussions.

Figure 8 Temporal variations in (a) ocean and (b) land biospheric sinks estimated from APO variations of this study (red) and process-based models of GCP (blue). The thin dashed lines represent the annual sinks, and the thick lines represent the pentad sinks. The purple lines represent the pentad sinks based on APO without ocean outgassing correction (Z_eff), and the light blue lines represent the sinks of GCP with the imbalance sinks added. The uncertainty associated with the pentad sinks with Z_eff corrections is shown as shaded area.

The pentad ocean sinks show an overall increasing trend, although there is a dip in the ocean sink centered on 2004–2005 by about 0.6 PgC yr⁻¹. Nevison et al. (2008) suggested that a decadal or longer period is needed to suppress the influence of the interannual variation in the ocean O₂ flux on the carbon sink estimation within ±0.1 PgC yr⁻¹ based on an ocean ecosystem model and an atmospheric transport model. In addition, the pentad APO changing rate still contains an uncertainty corresponding to ±0.5 PgC yr⁻¹ as discussed in Sect. 3.1. Therefore, the anomalous dip in the ocean sink for 2004–2005 might be an error caused by the anomalous ocean O₂ flux variations. The increasing rate of the ocean sink during 2001–2014, determined by linear regression, is 0.08±0.02 PgC yr⁻², which is larger than that of the GCP, which was 0.04±0.01 PgC yr⁻². Although the temporal variability in the ocean sink in the GCP study is rather suppressed, which is attributed to the rather coarse resolution of the GOBMs (Le Quéré et al., 2018), a much larger decadal and sub-decadal variability has been reported in the ocean sink estimations based on archived data of the observed surface partial pressure of CO₂ (pCO₂; Landschützer et al., 2015; DeVries et al., 2017). Results from the Surface Ocean pCO₂ Mapping (SOCOM) initiative show that the decadal linear trend of the global ocean sink enhancement over 2001–2011 based on pCO₂ data and selected mapping methods is about 0.8 PgC yr⁻¹ per decade (Rödenbeck et al., 2015). For a detailed comparison, the global ocean sinks based on pCO₂ observations and interpolation techniques (Landschützer et al., 2016; Rödenbeck et al., 2014) for the period of 1990–2017 are plotted in Fig. 9 together with the ocean sinks of both this study and the GCP. Note that the extended pCO₂-derived ocean sinks were given as supplementary data of Le Quéré et al. (2018), and those sinks were uniformly inflated by 0.78 PgC yr⁻¹ to compensate for the pre-industrial steady-state source of CO₂ derived from riverine input of carbon to the ocean (Resplandy et al., 2018). Both the GCP- and pCO₂-derived ocean sinks show changes in the trends before and after 2001, while the magnitude of the changes in the pCO₂-derived sinks is larger. The increasing rates determined by linear regression during 2001–2014 are 0.08±0.01 PgC yr⁻² for Landschüzer et al. (2016) and 0.07±0.02 PgC yr⁻² for Rödenbeck et al. (2014), which are more consistent with the rate found in this study. Therefore, our result seems to support a previous conclusion that the recent increase in the ocean sinks exceeds the increasing trend of ocean sink expected only from the atmospheric CO₂ increase (Landschützer et al., 2015; DeVries et al., 2017).

Figure 9 Comparison of the temporal variations in the ocean sinks based on the APO data of this study (red), global ocean biogeochemistry models (GOBMs) of GCP (blue), and pCO₂ data of Landschützer et al. (2016; light blue) and Rödenbeck et al. (2014; orange). The dashed lines represent the regression lines for the corresponding data during 2001–2014. Note that the pCO₂-based ocean sinks are adjusted for the pre-industrial ocean CO₂ emissions (±0.78 PgC yr⁻¹) caused by riverine CO₂ input to the ocean (Resplandy et al., 2018).

In contrast, the pentad land sinks of both this study and the GCP study show an increasing trend during 2001–2009 followed by a decreasing trend during 2009–2014, although the range of variations of this study is about 2 times larger than that of the GCP. The linear trends for the former period are 0.23±0.04 PgC yr⁻² for this study and 0.10±0.03 PgC yr⁻² for the GCP, and those for the latter period are $-\mathrm{0.22}\pm \mathrm{0.04}$ PgC yr⁻² for this study and $-\mathrm{0.12}\pm \mathrm{0.04}$ PgC yr⁻² for the GCP. An enhancement of the land uptake during the 2000s has been reported recently by several studies based on atmospheric inversions and biosphere models (Keenan et al., 2016; Kondo et al., 2018; Piao et al., 2018). Although there is an ongoing discussion about the detailed mechanisms of the enhanced net land uptake, the accelerated land uptake may partially explain the stagnation of the growth rate of the atmospheric CO₂ in the 2000s despite the increasing anthropogenic CO₂ emissions. Examining the atmospheric inversion studies and the previous version of the GCP (Le Quéré et al., 2015), in which the net land uptake was computed as residuals among the other carbon budget components, Piao et al. (2018) found that the linear increasing trend of the net land carbon sink during 1998–2012 was 0.17±0.05 PgC yr⁻². The linear trend of this study during 2001–2009 is close to the above value within the uncertainty. Although the corresponding linear trend of the latest GCP estimation is about half that of the present study, the sum of the net land sink and the budget imbalances of the GCP, plotted as light blue lines in Fig. 8, shows a much larger increasing trend, 0.20±0.03 PgC yr⁻², which is almost identical to our trend.

The land sinks of both this study and the GCP study exhibit decreasing trends for the period 2009–2014, which are partially compensated for by the steady increase in the ocean uptake. The atmospheric accumulation rate of CO₂ significantly increased in 2015 and 2016 (see Fig. 5), when one of the strongest El Niño events occurred. Studies based on atmospheric CO₂ observations from stations and by satellite indicated that a reduction in biospheric uptake and an increase in biomass burning contributed to the CO₂ increase during the El Niño event (Chatterjee et al., 2017; Patra et al., 2017). The decreasing trend of the pentad land uptake also reflects the change in the global carbon cycle associated with the El Niño event.

Shown as discrepancies between the pentad sinks of the GCP with and without budget imbalances (Fig. 8), the magnitude of the budget imbalances seems to increase after 2007. For the pentad sinks centered on 2007, 2008 and 2009, the ocean and land sinks of this study agree with those of the GCP without and with the budget imbalances, respectively. For the pentad sinks between 2010 and 2014, both the ocean and land sinks of this study are plotted between those of the GCP with and without the budget imbalances. Although we cannot show any definitive partitioning of the budget imbalance between ocean and land sinks because of a rather large uncertainty associated with the sink estimations, the above results seem to suggest that the budget imbalances are allocated to land sinks for the former period and to both sinks for the latter period.

From the above discussions, we feel that a 5-year duration effectively suppresses to some extent the anomalous variations in the carbon budget estimations based on APO, which are considered to be caused by the imbalance of the seasonal air–sea O₂ exchange. Probably, the 5-year average suppresses the variability in Z_eff to a level of ±0.5 PgC yr⁻¹, as discussed in Sect. 3.1. To reduce uncertainty in the carbon budget estimation, we need more effort to improve the quantification of the net O₂ outgassing associated with global ocean warming because the quantification of the Z_eff at this state is still very speculative. Applying the approach of Stendardo and Gruber (2012), who examined the long-term changes in dissolved O₂ and heat content by using archived oceanographic data of the Atlantic Ocean, to other ocean basins would improve our understanding of the long-term net ocean ratio of the O₂ flux to the heat flux.

4 Conclusions

We evaluated the global carbon budgets based on the APO data computed from the O₂∕N₂ and CO₂ of the flask samples collected in the Pacific region since 1997. In the carbon budget calculation, we corrected the ocean and land sinks with the ocean O₂ outgassing effect, Z_eff, based on the ocean heat increment for the 0–2000 m layer. Eventually, we obtained the following conclusions:

1. The long-term oceanic and land biotic carbon sinks were 2.6±0.7 PgC yr⁻¹ and 1.5±0.9 PgC yr⁻¹, respectively, for a 17-year period (2000–2016), and 2.4±0.7 and 1.9±0.9 PgC yr⁻¹, respectively, for a 14-year period (2003–2016). These long-term carbon sinks agreed well with those of the latest GCP estimation (Le Quéré et al., 2018); the differences of the individual estimations are less than ±0.35 PgC yr⁻¹.

2. The ocean and land sinks for the four pentad (5-year) periods (2000–2004, 2004–2008, 2008–2012 and 2012–2016) of this study also showed good agreement with those of the GCP within a difference of ±0.5 PgC yr⁻¹. The land and ocean sinks of this study showed larger values by about 0.5 PgC yr⁻¹ than those of the GCP for 2008–2012 and 2012–2016, respectively, when rather large carbon budget imbalances (> 0.5 PgC yr⁻¹) were found. Therefore, the discrepancies in the sinks between this study and the GCP might give a clue about how to partition the imbalance values between the land and ocean sinks.

3. Calculating the carbon budgets for the pentad periods consecutively, we examined the changing trend of the ocean and land sinks during a 14-year period (2001–2014). In general, the changing trends of both land and ocean sinks of this study agreed well with those of the GCP, although the range of variations of this study was about 2 times larger than that of the GCP study. The pentad ocean sinks showed an overall increasing trend for the entire period (2001–2014), with a linear increasing rate of 0.08±0.02 PgC yr⁻². This increasing rate was about 2 times larger than that for the GCP ocean sinks (0.04±0.01 PgC yr⁻²) but was consistent with that for the global ocean sinks based on pCO₂ observations and interpolation techniques (Landschützer et al., 2016; Rödenbeck et al., 2014). In contrast, the pentad land sinks showed an increasing trend for 2001–2009 and a decreasing trend for 2009–2014. The linear trends of the land sinks for this study and the GCP (in parentheses) were 0.23±0.04 PgC yr⁻² (0.10±0.03 PgC yr⁻²) for the former period and $-\mathrm{0.22}\pm \mathrm{0.04}$ PgC yr⁻² ($-\mathrm{0.12}\pm \mathrm{0.04}$ PgC yr⁻²) for the latter period. Enhancement of the land carbon uptake was reported also by previous studies (Keenan et al., 2016; Kondo et al., 2018; Piao et al., 2018). In addition, the recent decreasing trend of the land uptake was found to be partially related to the global carbon cycle variation associated with the strong El Niño event in 2015 and 2016.