Spatiotemporal variations of the δ(O₂ ∕ N₂), CO₂ and δ(APO) in the troposphere over the western North Pacific

We analyzed air samples collected on board a C-130 cargo aircraft over the western North Pacific from May 2012 to March 2020 for atmospheric δ(O₂ $/$ N₂) and CO₂ amount fraction. Observations were corrected for significant artificial fractionation of O₂ and N₂ caused by thermal diffusion during the air sample collection using the simultaneously measured δ(Ar $/$ N₂). The observed seasonal cycles of the δ(O₂ $/$ N₂) and atmospheric potential oxygen (δ(APO)) varied nearly in opposite phase to the cycle of the CO₂ amount fraction at all latitudes and altitudes. Seasonal amplitudes of δ(APO) decreased with latitude from 34 to 25^∘ N, as well as with increasing altitude from the surface to 6 km by 50 %–70 %, while those of the CO₂ amount fraction decreased by less than 20 %. By comparing the observed values with the simulated δ(APO) and CO₂ amount fraction values generated by an atmospheric transport model, we found that the seasonal δ(APO) cycle in the middle troposphere was modified significantly by a combination of the northern and southern hemispheric seasonal cycles due to the interhemispheric mixing of air. The simulated δ(APO) underestimated the observed interannual variation in δ(APO) significantly, probably due to the interannual variation in the annual mean air–sea O₂ flux. Interannual variation in δ(APO) driven by the net marine biological activities, obtained by subtracting the assumed solubility-driven component of δ(APO) from the total variation, indicated a clear influence on annual net sea–air marine biological O₂ flux during El Niño and net air–sea flux during La Niña. By analyzing the observed secular trends of δ(O₂ $/$ N₂) and the CO₂ amount fraction, global average terrestrial biospheric and oceanic CO₂ uptakes for the period 2012–2019 were estimated to be (1.8±0.9) and (2.8±0.6) Pg a⁻¹ (C equivalents), respectively.

1 Introduction

Atmospheric ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratios have been observed since the early 1990s, for the primary application of constraining the marine and terrestrial exchange of CO₂ (Keeling and Shertz, 1992). For this purpose, observations of the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio have been carried out at many surface stations and on commercial cargo ships (e.g., Bender et al., 2005; Manning and Keeling, 2006; Tohjima et al., 2008, 2019; Goto et al., 2017). The ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio varies in opposite phase with the CO₂ amount fraction due to the terrestrial biosphere exchanges and fossil fuel combustion, of which respective O₂ : CO₂ exchange ratios (oxidative ratio (OR), −Δy(O₂)Δy(CO₂)⁻¹) are about 1.1 and 1.4, respectively (Keeling, 1988; Severinghaus, 1995), where y stands for the dry-amount fraction of gas, as recommended by Cohen et al. (2007). Using the OR value of 1.1 for terrestrial biospheric activities, atmospheric potential oxygen (APO) is defined by y(APO) =y(O₂) + 1.1y(CO₂) (Stephens et al., 1998). While APO is conserved for terrestrial biospheric activities, the air–sea exchange of O₂ is much faster than that of CO₂ since the air–sea CO₂ exchange is highly suppressed by the carbonate buffer system in seawater (e.g., Keeling et al., 1993). Therefore, APO can be used to evaluate air–sea O₂ fluxes associated with marine biological and physical processes (e.g., Nevison et al., 2012).

Aircraft observation serves as a useful platform for measuring altitude-dependent APO driven by spatially integrated air–sea O₂ fluxes at the surface. Aircraft observations of the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio have been conducted in the past (e.g., Sturm et al., 2005; Steinbach, 2010; Ishidoya et al., 2008a, 2012, 2014; van der Laan et al., 2014; Bent, 2014; Morgan et al., 2019; Birner et al., 2020; Stephens et al., 2018, 2021). Sturm et al. (2005) observed a vertical gradient and seasonal cycle in the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio in the altitude range of 0.8–3.1 km over Perthshire, United Kingdom, for the period 2003–2004. Longer-term observations of the tropospheric ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio have also been carried out by Ishidoya et al. (2012) and van der Laan et al. (2014). They conducted aircraft observations at altitudes of 2, 4 and above 8 km over Japan during 1999–2010 and at altitudes of 0.1 and 3 km over western Russia during 1998–2008, respectively, and provided additional evidence of seasonal cycles and secular changes in the tropospheric ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio. However, there were uncertainties associated with artificial fractionations of O₂ and N₂ in Ishidoya et al. (2012) and van der Laan et al. (2014). The $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratio and/or stable isotopic ratios can be used to evaluate natural and artificial molecular-diffusive fractionations of O₂ and N₂ (e.g., Kawamura et al., 2006; Ishidoya et al., 2013); however Ishidoya et al. (2012) and van der Laan et al. (2014) did not observe them.

Ishidoya et al. (2014) and Stephens et al. (2021) observed ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ simultaneously and were able to correct for thermally diffusive artificial fractionation of ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, using coefficients of 3.54 and 3.77 of $\mathrm{Ar}/{\mathrm{O}}_{\mathrm{2}}$, respectively. Using the corrected ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio, Ishidoya et al. (2014) were able to observe spatiotemporal variations in the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio from the surface to the middle troposphere over the western North Pacific around Japan, on monthly scheduled flights, for the period May 2012 to April 2013. Similarly, Stephens et al. (2021) were able to provide a better picture of much wider area distributions of the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio, from 0–14 km and 87^∘ N to 85^∘ S, using measurements from a series of aircraft campaigns such as five HIAPER Pole-to-Pole Observations (HIPPO) campaigns in 2009–2011 (Wofsy et al., 2011) and the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ Ratio and CO₂ Airborne Southern Ocean (ORCAS) study in 2016 (Stephens et al., 2018).

Stephens et al. (2021) also conducted continuous observations of O₂ mole fraction using a vacuum ultraviolet (VUV) absorption detector (Stephens et al., 2003). They adjusted the continuous O₂ data to the simultaneously observed flask-based ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio corrected for the artificial fractionation using the $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratio, since the artificial fractionations for the continuous observations were more significant than the fractionation for the flask sampling due to the lower flow rate. Based on the continuous O₂ data, Morgan et al. (2019) reported summertime vertical gradients of the atmospheric ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio and the CO₂ amount fraction through the atmospheric boundary layer over the Drake Passage region of the Southern Ocean, to evaluate the air–sea ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{CO}}_{\mathrm{2}}$ flux ratios in the region. Aircraft observations of $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ are also used to evaluate gravitational separation of the atmospheric components, which is an indicator of the Brewer–Dobson circulation (e.g., Ishidoya et al., 2013), in the lowermost stratosphere under the condition that the artificial fractionation is reduced sufficiently (Ishidoya et al., 2008; Birner et al., 2020).

In this study, as an update to Ishidoya et al. (2014), we present 9-year-long ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio variations observed in the troposphere over the western North Pacific. Measurements were carried out on monthly scheduled cargo aircraft flights with a fixed flight route, and the thermally diffusive artificial fractionations on the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio were corrected using the simultaneously measured $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratio. Using these corrected values, we made precise evaluation of the altitude–latitude distributions of seasonal cycle, vertical profile and year-to-year variation along the flight route. We mainly focus our discussions on the variations in the tropospheric APO and CO₂ amount fraction with the aid of a 3-D atmospheric chemistry-transport model. We also estimate average terrestrial biospheric and oceanic CO₂ uptakes for the period 2012–2019, using long-term trends of the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio and the CO₂ amount fraction.

2 Methods

The C-130 cargo aircraft flies once per month from Atsugi Base (35.45^∘ N, 139.45^∘ E), Kanagawa, Japan, to Minamitorishima, Japan, a small coral atoll (MNM; 24.29^∘ N, 153.98^∘ E). Two types of C-130, C-130H and C-130R, have been used for the observations. The cruising altitude is about 6 km, and 24 air samples were pressurized into 1.7 L silica-coated titanium flasks to an absolute pressure of 0.4 MPa during the flight. A set of 17–20 samples were collected during the level flight while others were collected during the descent portion at MNM. Flask air sampling for the C-130H flight is manually operated by two JMA personnel on board the aircraft. Therefore, a diaphragm pump is modified to be operated by hand, without the electric power supply. Air samples are collected from the air-conditioning system in the C-130H aircraft. Fresh air outside the aircraft is compressed by an engine and fed into the air-conditioning ducts by passing it through a pneumatic system and air cycle packs. A Teflon tube with a diameter of $\mathrm{1}/\mathrm{4}$ in. is inserted into the air-conditioning blowing nozzle upstream of the recirculation fan for the air sampling, so that the sample air is not contaminated with the cabin air. The open end of the air sampling intake is situated at wing root and faces the front side of the aircraft. Unfortunately, details of the air sampling line from the inlet to flask sampler have not been given to researchers from the Japan Ministry of Defense. It is noted that we only use a Teflon tube upstream of the diaphragm pump to avoid absorption and/or permeation of O₂ due to a pressurization of the Teflon tube. In this regard, we found measured values of ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio for the air samples corrected in January 2016 were anomalously low, when a Teflon tube was used downstream of the diaphragm pump. Therefore, we exclude the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ ratio values in January 2016 from the analyses in this study. During other flights, we use a flexible tube made of stainless steel downstream of the diaphragm pump. Details of the air sampling method have been described elsewhere (Tsuboi et al., 2013).

Figure 1 Locations of C-130 aircraft air sampling for the period May 2012–March 2020 (circles). Locations of Minamitorishima (MNM), Ioto Island and Atsugi Base are also shown by stars. Latitudinal and vertical distributions are taken during the level flight and descent portion at MNM, respectively. In total, 2206 of air samples were collected, and 1783 of them were analyzed for ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratios, as well as the stable isotopic ratios of N₂, O₂ and Ar.

The flask air samples were brought back to the Japan Meteorological Agency (JMA) and analyzed for CO₂, CH₄, CO and N₂O amount fractions. The observational results are reported by Niwa et al. (2014). The CO₂ amount fraction was measured using a non-dispersive infrared analyzer (Licor, LI-7000) with a precision of better than ±0.07 µmol mol⁻¹ (Tsuboi et al., 2013) against the World Meteorological Organization (WMO) mole fraction scales for CO₂ (Zhao and Tans, 2006). The dataset is posted on the WMO's World Data Centre for Greenhouse Gases (WMO/WDCGG; https://doi.org/10.50849/WDCGG_0001-8002-1001-05-02-9999; Saito, 2016, updated, 2022). After the JMA analyses, the flasks were sent to the National Institute of Advanced Industrial Science and Technology (AIST) to measure ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratios, as well as the stable isotopic ratios of N₂, O₂ and Ar (Ishidoya et al., 2014). In this study, we present the measured data obtained from the air samples collected for the period May 2012–March 2020. In Fig. 1, we show all the locations where the air samples were collected on board C-130 aircraft during the observation period and the locations of MNM, Atsugi Base and Ioto Island (24.76^∘ N, 141.29^∘ E), Japan.

The values of δ(O₂ $/$ N₂) and δ(Ar $/$ N₂) and stable isotopic ratios of N₂, O₂ and Ar (δ(¹⁵N), δ(¹⁸O) and δ(⁴⁰Ar)) are reported per meg (one per meg is equal to $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{6}}$):

$$\begin{array}{}\text{(1)}& {\displaystyle }\mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)={\displaystyle \frac{R_{\mathrm{sample}}\left({}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)-R_{\mathrm{standard}}\left({}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}{R_{\mathrm{standard}}\left({}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}},\text{(2)}& {\displaystyle }\mathit{\delta }\left(\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}\right)={\displaystyle \frac{R_{\mathrm{sample}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)-R_{\mathrm{standard}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}{R_{\mathrm{standard}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}},\text{(3)}& {\displaystyle }\mathit{\delta }\left({}^{\mathrm{15}}\mathrm{N}\right)={\displaystyle \frac{R_{\mathrm{sample}}\left({}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)-R_{\mathrm{standard}}\left({}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}{R_{\mathrm{standard}}\left({}^{\mathrm{15}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}{/}^{\mathrm{14}}{\mathrm{N}}^{\mathrm{14}}\mathrm{N}\right)}},\text{(4)}& {\displaystyle }\mathit{\delta }\left({}^{\mathrm{18}}\mathrm{O}\right)={\displaystyle \frac{R_{\mathrm{sample}}\left({}^{\mathrm{18}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}\right)-R_{\mathrm{standard}}\left({}^{\mathrm{18}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}\right)}{R_{\mathrm{standard}}\left({}^{\mathrm{18}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}{/}^{\mathrm{16}}{\mathrm{O}}^{\mathrm{16}}\mathrm{O}\right)}},\text{(5)}& {\displaystyle }\mathit{\delta }\left({}^{\mathrm{40}}\mathrm{Ar}\right)={\displaystyle \frac{R_{\mathrm{sample}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}\right)-R_{\mathrm{standard}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}\right)}{R_{\mathrm{standard}}\left({}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}\right)}},\end{array}$$

where the subscripts “sample” and “standard” refer to the values of the sample and standard air, respectively. The values of δ(O₂ $/$ N₂), δ(Ar $/$ N₂), δ(¹⁵N), δ(¹⁸O) and δ(⁴⁰Ar) of the air samples were determined against our primary standard air (cylinder no. CRC00045) using a mass spectrometer (Thermo Scientific Delta-V) (Ishidoya and Murayama, 2014) with a respective reproducibility of about 5, 8, 1.5, 3 and 13 per meg (1σ). The scale based on the primary standard air is our original scale and called the “EMRI/AIST scale” in Aoki et al. (2021).

Figure 2 Measured δ(Ar $/$ N₂), δ(¹⁸O) of O₂ and δ(⁴⁰Ar) plotted against δ(¹⁵N) of N₂ for all the collected air samples (solid dots). Least-squares regression lines fitted to the data are shown as solid lines, while the relationships expected from mass-dependent fractionation of air molecules are shown by dotted lines.

As already discussed in Ishidoya et al. (2014), the measured values of δ(O₂ $/$ N₂) are contaminated by significant artificial thermally diffusive fractionation of O₂ and N₂ during the air sample collection process on board the aircraft. Figure 2 shows the relationships between δ(Ar $/$ N₂), δ(¹⁸O) and δ(⁴⁰Ar) with δ(¹⁵N) for all the air samples analyzed in this study. It was found that δ(Ar $/$ N₂), δ(¹⁸O) and δ(⁴⁰Ar) change linearly in proportion to δ(¹⁵N), and the linear regression analyses gave respective slopes of (16.3±0.1), (1.58±0.01) and (2.69±0.05) per meg (per meg)⁻¹ for the δ(Ar $/$ N₂)/δ(¹⁵N), δ(¹⁸O) $/$ δ(¹⁵N) and δ(⁴⁰Ar) $/$ δ(¹⁵N) ratios. These ratios agree well with the ratios of (16.2±0.1), (1.55±0.02) and (2.75±0.05) for δ(Ar $/$ N₂) $/$ δ(¹⁵N), δ(¹⁸O) $/$ δ(¹⁵N) and δ(⁴⁰Ar) $/$ δ(¹⁵N), respectively, determined from the laboratory experiments on the effect of thermally diffusive fractionations on δ(Ar $/$ N₂), δ(¹⁵N), δ(¹⁸O) and δ(⁴⁰Ar) (Ishidoya et al., 2013). Therefore, we decided to correct for the thermally diffusive fractionation of O₂ and N₂ on the observed δ(O₂ $/$ N₂) using the following equation (Ishidoya et al., 2014):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{\mathit{\delta }}_{\mathrm{cor}.}\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)& ={\mathit{\delta }}_{\mathrm{meas}.}\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)-{\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{2}}}\\ & \cdot {\mathit{\alpha }}_{\mathrm{Ar}}^{-\mathrm{1}}\cdot \mathrm{\Delta }{\mathit{\delta }}_{\mathrm{meas}.}\left(\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}\right).\end{array}\end{array}$$

Here, δ_cor.(O₂ $/$ N₂) and δ_meas.(O₂ $/$ N₂) denote the corrected and measured δ(O₂ $/$ N₂), respectively. The coefficients ${\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{2}}}=(\mathrm{4.57}\pm \mathrm{0.02})$ and ${\mathit{\alpha }}_{\mathrm{Ar}}=(\mathrm{16.2}\pm \mathrm{0.1})$ are the δ(O₂ $/$ N₂)/δ(¹⁵N) and δ(Ar $/$ N₂) $/$ δ(¹⁵N) ratios respectively, determined from the laboratory experiments as described by Ishidoya et al. (2013). The value for ${\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{2}}}$ is not directly reported in Ishidoya et al. (2013), but it was obtained from the same laboratory experiment (Fig. 2 in their study). The ${\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{2}}}/{\mathit{\alpha }}_{\mathrm{Ar}}$ ratio of $\mathrm{4.57}/\mathrm{16.2}$ is close to the Keeling et al. (2004) diffusion factor for (Ar $/$ N₂) $/$ (O₂ $/$ N₂) and results in the same tracer δ(O₂ $/$ N₂)^∗ (Stephens et al., 2021). Δδ_meas.(Ar $/$ N₂) is the deviation of the measured δ(Ar $/$ N₂) from its reference point that is determined using the annual mean value of δ(Ar $/$ N₂) in 2013 observed at the surface in Tsukuba (36^∘ N, 140^∘ E), Japan (Ishidoya and Murayama, 2014). Therefore, the effect of the seasonal δ(Ar $/$ N₂) cycle on δ_cor.(O₂ $/$ N₂) was not excluded in this study. This could lead to an over- and under-correction of the surface δ_cor.(O₂ $/$ N₂) by about 2 per meg in the summertime and wintertime, respectively. Considering the measurement reproducibility of δ_meas.(O₂ $/$ N₂) of 5 per meg, uncertainty of the correction of 2 per meg (8 per meg $\times {\mathit{\alpha }}_{{\mathrm{O}}_{\mathrm{2}}}{\mathit{\alpha }}_{\mathrm{Ar}}^{-\mathrm{1}}$; 8 per meg is the measurement reproducibility of δ_meas.($\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$)) and the δ(Ar $/$ N₂) seasonality, the overall uncertainty of δ_cor.(O₂ $/$ N₂) was evaluated to be less than 6 per meg. It is noted that the annual mean value of δ(Ar $/$ N₂) also shows interannual variation with a peak-to-peak amplitude of 9 per meg a⁻¹ and secular trend of 0.75 per meg a⁻¹ (Ishidoya et al., 2021). They result in uncertainties of δ_cor.(O₂ $/$ N₂) by 2.5 and 0.1 per meg a⁻¹ for the interannual variation and secular trend, respectively, which are significantly smaller than those of the observed δ_cor.(O₂ $/$ N₂) discussed below. The details of the correction of artificial fractionations of O₂ and N₂ are given in Ishidoya et al. (2014). It is noted that the correction using δ(¹⁵N), which is stable in the atmosphere over long time periods, is also suitable to obtain δ_cor.(O₂ $/$ N₂). Since the uncertainty of the corrected δ_cor.(O₂ $/$ N₂) using δ(¹⁵N) is ±7 per meg ($\pm \mathrm{1.5}\times \mathrm{4.57}$), which is larger than that using δ(Ar $/$ N₂) (±2 per meg), we determined to use δ(Ar $/$ N₂) rather than δ(¹⁵N) in the present study.

We used δ(APO) for the detailed analyses of the air–sea exchange of O₂. δ(APO) was calculated from the observed δ_cor.(O₂ $/$ N₂) and CO₂ amount fraction:

$$\begin{array}{}\text{(7)}& \mathit{\delta }\left(\mathrm{APO}\right)={\mathit{\delta }}_{\mathrm{cor}.}\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)+{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{B}}}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}}y\left({\mathrm{CO}}_{\mathrm{2}}\right)-\mathrm{2000}\times {\mathrm{10}}^{-\mathrm{6}},\end{array}$$

where y(CO₂) is the dry-amount fraction of CO₂, α_B is the OR of 1.1 for terrestrial biospheric activities, X(O₂) of 0.2093–0.2094 is the amount fraction of atmospheric O₂ (Tohjima et al., 2005; Aoki et al., 2019) and 2000 is an arbitrary reference. We use 0.2094 for the calculation of δ(APO). From the definition, δ(APO) is conserved for terrestrial biospheric activities and driven not only by the air–sea exchange of O₂, N₂ and CO₂ but also by fossil fuel consumption, of which the OR value is generally larger than 1.1 (Keeling, 1988). In order to investigate the observed δ_cor.(O₂ $/$ N₂) variations, we used a three-dimensional atmospheric transport model NICAM-TM (Niwa et al., 2011) to simulate the CO₂ amount fraction and δ(APO) using surface O₂, N₂ and CO₂ fluxes. NICAM-TM is based on the Nonhydrostatic ICosahedral Atmospheric Model (NICAM; Satoh et al., 2008, 2014), and its tracer transport version has been used for atmospheric transport and flux inversion studies of greenhouse gases (e.g., Niwa et al., 2012). The horizontal model resolution used in this study had a mean grid interval of about 112 km. The model was driven by nudging the horizontal winds towards the Japanese 55-year Reanalysis data (JRA-55; Kobayashi et al., 2015).

The surface fluxes incorporated into NICAM-TM were the air–sea fluxes of O₂, N₂ and CO₂ and also CO₂ and O₂ fluxes from fossil fuel combustion and the terrestrial biosphere. The air–sea O₂ and N₂ fluxes were the climatological seasonal anomalies taken from the TransCom experimental protocol (Blaine, 2005; Garcia and Keeling, 2001). The fluxes were computed to give the seasonal component and the annual mean values at every grid point to be zero. The air–sea CO₂ flux was obtained from the monthly sea surface CO₂ flux climatology of Takahashi et al. (2009). The CDIAC fossil fuel database was used for the fossil fuel CO₂ flux (Andres et al., 2016; Gilfillan et al., 2019). The inversion flux reported by Niwa et al. (2012) was used for the terrestrial biosphere CO₂ flux. Model-based changes in δ(APO), CO₂ amount fraction and δ(O₂ $/$ N₂) (Δδ(APO), Δy(CO₂) and Δδ(O₂/N₂)) were calculated using the following equations (e.g., Nevison et al., 2008; Tohjima et al., 2012) in units per meg, micromoles per mole (µmol mol⁻¹) and per meg, respectively:

$$\begin{array}{}\text{(8)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }\mathit{\delta }\left(\mathrm{APO}\right)& =\left({\displaystyle \frac{\mathrm{\Delta }{y}^{\mathrm{SA}}\left({\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}}-{\displaystyle \frac{\mathrm{\Delta }{y}^{\mathrm{SA}}\left({\mathrm{N}}_{\mathrm{2}}\right)}{X\left({\mathrm{N}}_{\mathrm{2}}\right)}}\right)\\ & +\left({\displaystyle \frac{-{\mathit{\alpha }}_{\mathrm{F}}\cdot \mathrm{\Delta }{y}^{\mathrm{FF}}\left(\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\right)+{\mathit{\alpha }}_{\mathrm{B}}\cdot \mathrm{\Delta }{y}^{\mathrm{FF}}\left(\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}}\right)\\ & +{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{B}}\cdot \mathrm{\Delta }{y}^{\mathrm{OC}}\left(\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}},\end{array}\text{(9)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }y\left({\mathrm{CO}}_{\mathrm{2}}\right)& =\mathrm{\Delta }{y}^{\mathrm{FF}}\left({\mathrm{CO}}_{\mathrm{2}}\right)+\mathrm{\Delta }{y}^{\mathrm{OC}}\left({\mathrm{CO}}_{\mathrm{2}}\right)\\ & +\mathrm{\Delta }{y}^{\mathrm{TB}}\left({\mathrm{CO}}_{\mathrm{2}}\right),\end{array}\text{(10)}& {\displaystyle }\begin{array}{rl}\mathrm{\Delta }\mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)& =\left({\displaystyle \frac{\mathrm{\Delta }{y}^{\mathrm{SA}}\left({\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}}-{\displaystyle \frac{\mathrm{\Delta }{y}^{\mathrm{SA}}\left({\mathrm{N}}_{\mathrm{2}}\right)}{X\left({\mathrm{N}}_{\mathrm{2}}\right)}}\right)\\ & +{\displaystyle \frac{-{\mathit{\alpha }}_{\mathrm{F}}\cdot \mathrm{\Delta }{y}^{\mathrm{FF}}\left(\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}}\\ & +{\displaystyle \frac{-{\mathit{\alpha }}_{\mathrm{B}}\cdot \mathrm{\Delta }{y}^{\mathrm{TB}}\left(\mathrm{C}{\mathrm{O}}_{\mathrm{2}}\right)}{X\left({\mathrm{O}}_{\mathrm{2}}\right)}},\end{array}\end{array}$$

where Δy(O₂), Δy(N₂) and Δy(CO₂) are changes in dry-amount fractions of the respective gases calculated using NICAM-TM. The superscripts “SA”, “FF”, “OC” and “TB” denote the seasonal anomaly of the air–seas O₂ and N₂ flux, CO₂ flux from fossil fuel combustion, ocean and terrestrial biosphere, respectively. X(O₂) and α_B have the same meaning as in Eq. (7), while X(N₂) is the dry-amount fraction of N₂ in the atmosphere, and α_F is the global average OR for fossil fuel combustion. In this study, we adopted X(N₂)=0.7808 and α_F=1.37. The α_F was calculated from the fossil fuel and cement production emissions by fuel type for the period 2012–2019, reported by the Global Carbon Project (GCP; Friedlingstein et al., 2020), and the oxidative ratios for the different fuel type were taken from Keeling (1988). It should be noted that we assume initial amount fractions of y(O₂) and y(N₂) are equal to X(O₂) and X(N₂), respectively, in Eqs. (8) and (10). We can rewrite Eq. (8) as

$$\begin{array}{}\text{(11)}& \mathrm{\Delta }\mathit{\delta }\left(\mathrm{APO}\right)=\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{SA}}\left(\mathrm{APO}\right)+\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{FF}}\left(\mathrm{APO}\right)+\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{OC}}\left(\mathrm{APO}\right).\end{array}$$

Finally, δ(O₂ $/$ N₂), y(CO₂) and δ(APO) obtained from NICAM-TM are reported as respective deviations of Δδ(APO), Δy(CO₂) and Δδ(O₂ $/$ N₂) from arbitrary reference points, in other words, reported on different scales from those of observations. The δ(APO) simulation run incorporating the above-mentioned surface fluxes is referred to as the “control run”. In this calculation, δ(APO) driven by an annual mean air–sea O₂ and N₂ fluxes (hereafter referred to as the “δ^AM(APO)”), which were estimated by Gruber et al. (2001), was ignored. In other words, the δ(APO) obtained from the NICAM-TM control run ignores the contribution of not only interannual variations but also spatial distribution in annual mean air–sea O₂ and N₂ fluxes. If the control run represents other components other than δ^AM(APO) correctly, then the contribution of δ^AM(APO) was evaluated by subtracting the simulated δ(APO) from the observed δ(APO). In Sect. 3.2, we will discuss interannual variation in this context.

Figure 3 (a) Measured values of δ(O₂ $/$ N₂) and δ(Ar $/$ N₂) for all the air samples collected on board the C-130 aircraft. δ(Ar $/$ N₂) values observed at Tsukuba, Japan, are also shown as gray circles (Ishidoya et al., 2021). Shaded areas denote the periods when C-130R aircraft was used, and C-130H aircraft was used in other periods. The color bar indicates the altitude where the air sampling was carried out. (b) δ_cor.(O₂ $/$ N₂) corrected for artificial fractionation by applying Eq. (6) (see text). The color bar indicates the altitude where the air sampling was carried out. (c) Vertical profiles of detrended δ_cor.(O₂ $/$ N₂) separated by season (DJF: December to February, MAM: March to May, JJA: June to August, and SON: September to November), obtained by subtracting a linear secular trend, fitted to the δ_cor.(O₂ $/$ N₂) in (b), from each δ_cor.(O₂ $/$ N₂) value. The color bar indicates the latitude where the air sampling was carried out.

Figure 4 (a) δ_cor.(O₂ $/$ N₂), CO₂ amount fraction and δ(APO) observed at the altitude of (6.1±0.5) (±1σ) km at various latitudes over the western North Pacific. Best-fit curves to the data (solid lines) and secular trends (dashed lines) are also shown (Nakazawa et al., 1997a). (b) Same as in (a) but for calculated values obtained from the NICAM-TM control-run. The data from NICAM-TM are reported on arbitrary scales different from those of observations.

3 Results and discussion

3.1 Latitudinal and vertical distributions of δ_cor.(O₂ $/$ N₂), CO₂ amount fraction and δ(APO)

Figure 3a shows the measured δ(O₂ $/$ N₂) and δ(Ar $/$ N₂) for all the air samples observed in this study, and Fig. 3b shows the δ_cor.(O₂ $/$ N₂) corrected values by applying Eq. (6) to the measured values. As seen from the figures, significant artificial fractionations found in the measured δ(O₂ $/$ N₂) are reduced dramatically by the correction. It can also be seen from Fig. 3a that the fractionations have become smaller since 2018, especially at higher altitude. It should be noted that these noted changes across 2018 could be at least partly related to changes in the aircraft type from C-130H to C-130R used for flask sampling. The periods when the C-130R was used are indicated by the pale red shade. No systematic data gaps were found in the δ_cor.(O₂ $/$ N₂) time series across 2018, so that we successfully corrected the fractionation of O₂ and N₂ both for C-130H and C-130R. To examine the validity of the significant corrections applied to the δ_meas.(O₂ $/$ N₂), we also show the vertical profiles of detrended δ_cor.(O₂ $/$ N₂) separated by season in Fig. 3c, obtained by subtracting a linear secular trend from the data in Fig. 3b. As seen from the figure, seasonal δ_cor.(O₂ $/$ N₂) variations with summertime maxima are clearly seen at all altitudes, which supports validity of the correction. We also presented some examples that synoptic-scale variations in δ(O₂ $/$ N₂) can also be reproduced using δ_cor.(O₂ $/$ N₂) in our previous study (Fig. 4 in Ishidoya et al., 2014). Since details of the air sampling line from the inlet to flask sampler have not been informed to researchers, it is difficult to explain the cause(s) of the much lower variability in δ(Ar $/$ N₂) data for periods of C-130R. The visible difference is the location of the air-conditioning blowing nozzles of C-130H and C-130R. The former and the latter are attached at the roof in the cockpit and at the front of the assistant driver's seat, respectively, and air sampling tubes are inserted into the nozzles.

Figure 5 (a) δ_cor.(O₂ $/$ N₂), the CO₂ amount fraction and δ(APO) observed in the troposphere over MNM at various altitudes. Best-fit curves to the data (solid lines) and secular trends (dashed lines) are also shown (Nakazawa et al., 1997a). (b) Same as in (a) but for calculated values using NICAM-TM. The data from NICAM-TM are reported on arbitrary scales different from those of observations.

Figure 4a shows variations in the δ_cor.(O₂ $/$ N₂), CO₂ amount fraction and δ(APO) observed in the layer (6.1±0.5) km (±1σ) (hereafter referred to as “middle troposphere”) at five latitudes over the western North Pacific. Best-fit curves to the data and secular trends obtained using a digital filtering technique (Nakazawa et al., 1997a) are also shown. Using this filtering technique, the average seasonal cycles were approximated by the sum of the fundamental and its first harmonic with the respective periods of 12 and 6 months. The residuals obtained by subtracting the approximated average seasonal cycle from the data were interpolated linearly to calculate daily values of δ_cor.(O₂ $/$ N₂), CO₂ amount fraction and δ(APO). The daily values were smoothed by the 26th-order Butterworth filter with a cutoff period of 36 months to derive the long-term trend. The long-term trend thus obtained was subtracted from the data, and the average seasonal cycle was determined again from the residuals. These steps were repeated until an unchangeable long-term trend obtained. The observational data deviated from the best-fitted curves more than ±3σ are excluded from the analyses and not shown in the figures discussed in this study (2 % of the observational data are excluded in total). As can be seen in Fig. 4a, secular decreases in δ_cor.(O₂ $/$ N₂) and δ(APO) and increases in the CO₂ amount fraction, accompanied by prominent seasonal cycles, were observed at each latitude. The secular changes in δ_cor.(O₂ $/$ N₂) and the CO₂ amount fraction can be attributed mainly to O₂ consumption and CO₂ emission resulting from fossil fuel combustion. The seasonally dependent air–sea O₂ flux and the terrestrial biospheric activity contribute towards the observed seasonal δ(O₂ $/$ N₂) cycle, while the terrestrial biospheric activity is the main contributor to the seasonal CO₂ amount fraction cycle (e.g., Keeling et al., 1993; Keeling and Manning, 2014). The average rates of change in the observed δ_cor.(O₂ $/$ N₂), CO₂ amount fraction and δ(APO) at the four latitudes shown in Fig. 4a were ($-\mathrm{24.2}\pm \mathrm{0.4}$) per meg a⁻¹, (2.43±0.05) $\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$ and ($-\mathrm{11.3}\pm \mathrm{0.4}$) per meg a⁻¹, respectively, for the observational period. General features of the observed variations in δ(O₂ $/$ N₂), the CO₂ amount fraction and δ(APO) are well reproduced by the control run of NICAM-TM, as shown in Fig. 4b. Figure 5a shows variations in δ_cor.(O₂ $/$ N₂), the CO₂ amount fraction and δ(APO) observed over MNM. Clear secular trends and prominent seasonal cycles of δ_cor.(O₂ $/$ N₂), the CO₂ amount fraction and δ(APO) are distinguishable, similar to those in Fig. 4a. We also note in Fig. 5a that the seasonal amplitudes of δ_cor.(O₂ $/$ N₂), the CO₂ amount fraction and δ(APO) decrease with increasing altitudes. These features are also reproduced by the control run of NICAM-TM (Fig. 5b).

Figure 6 (a) Average seasonal cycles of δ(APO) and the CO₂ amount fraction observed in the troposphere at various latitudes over the western North Pacific. Dashed lines denote the average values throughout the observation period. (b) Same as in (a) but for calculated values obtained from the NICAM-TM control run and the corresponding results from the without-SH-flux run (see text). (c) Same average seasonal cycles of the observed and simulated δ(APO) in (a) and (b), along with the detrended observed values at each latitude. Error bands (shaded) indicate standard deviations of the observed δ(APO) from the average seasonal cycles (±1σ). (d) Same as in (c) but for the CO₂ amount fraction.

Figure 6a shows average seasonal cycles of δ(APO) and the CO₂ amount fraction observed at five latitudes over the western North Pacific. They show clear seasonal cycles with summertime maxima in δ(APO) and minima in CO₂. However, the amplitude of seasonal δ(APO) cycle decreases significantly toward the lower latitudes, with seasonal maxima and minima clearly occurring earlier than those of the CO₂ cycle. Figure 6b shows the corresponding average seasonal cycles of δ(APO) and the CO₂ amount fraction obtained from the control run of NICAM-TM. Earlier appearances of the seasonal maxima and minima in δ(APO) than those in the CO₂ amount fraction are reproduced by NICAM-TM. The seasonal amplitude of the simulated δ(APO) also decreases toward the lower latitudes, in agreement with the observation, but is underestimated. As for the CO₂ amount fraction, the simulated seasonal cycles agree well with the observations. In Fig. 6c and d, we also show the same average seasonal cycles of the observed and simulated δ(APO) in Fig. 6a and b, respectively, along with the detrended observed values at each latitude for the convenience of visual comparison between the observed and simulated results.

Figure 7 Same as in Fig. 6 but for average seasonal cycles of δ(APO) and the CO₂ amount fraction in the troposphere over MNM.

Figure 7a shows average seasonal cycles of δ(APO) and the CO₂ amount fraction observed over MNM. The δ(APO) seasonal cycle varies in opposite phase to the CO₂ amount fraction. However, amplitudes of the seasonal δ(APO) cycle decrease significantly with the higher altitude, with seasonal minima clearly occurring earlier than those of the CO₂ cycle. These salient characteristics are reproduced generally by the NICAM-TM control run (Fig. 7b). Similar to Fig. 6, we also show Fig. 7c and d for visual comparison between the observed and simulated results.

In order to identify and explore some of the cause(s) that gave rise to the observed differences in the latitudinal and altitudinal changes in the seasonal cycles between δ(APO) and the CO₂ amount fraction shown in Figs. 6a and 7a, we carried out additional NICAM-TM simulations. In the calculation, we used the same fluxes that were used in the control run but for the northern hemispheric flux only for the TransCom seasonal climatology (hereafter referred to as “without-SH-flux run”). It is well known that the opposing phase of the seasonal δ(APO) cycles between the Northern and Southern Hemisphere is due to seasonal changes in the air–sea O₂ (N₂) flux with summertime maxima (e.g., Keeling et al., 1998; Tohjima et al., 2012). Therefore, in comparing the control run with the without-SH-flux run, we decided to evaluate changes in the seasonal δ(APO) cycle by superimposing the anti-phase seasonal cycles through the interhemispheric mixing of air. On the other hand, the seasonal CO₂ amount fraction cycle is much smaller in the Southern Hemisphere than that in the Northern Hemisphere (e.g., Nakazawa et al., 1997b). Therefore, it is expected that the seasonal cycle in the CO₂ amount fraction in the Northern Hemisphere does not change significantly by the interhemispheric atmospheric mixing.

Figure 8 (a) Latitudinal distribution of average ratios of the observed seasonal δ(APO) and CO₂ amount fraction amplitudes calculated relative to the values observed at 25.5^∘ N in the troposphere over the western North Pacific throughout the observation period (black filled circles). Error bands (shaded) indicate year-to-year variations (±1σ). The corresponding results calculated using the NICAM-TM control run (blue triangles) and the without-SH-flux run (green triangles) are also shown. Amplitude of seasonal APO cycle was calculated as a difference between an average value during July to September and that during February to April considering its broad peak, while that of CO₂ amount fraction was evaluated as a difference between seasonal maximum and minimum values. (b) Same as in (a) but for vertical distribution over MNM relative to the corresponding values at 1.3 km, and amplitude of seasonal APO cycle was evaluated as a difference between an average value during July to October and that during February to April. The average fractions of surface seasonal cycles obtained from continuous observations of δ(O₂ $/$ N₂) and the CO₂ amount fraction at MNM since December 2015 and the corresponding results calculated using NICAM-TM control run are also plotted.

The seasonal δ(APO) cycles obtained from the without-SH-flux run are also shown in Figs. 6b and 7b. As seen from Fig. 6b, latitudinal differences in the seasonal δ(APO) amplitudes from the without-SH-flux run are clearly smaller than those from the control run. In addition, appearances of the maxima and minima of the seasonal δ(APO) cycle from the without-SH-flux run are later than those from the control run, pushing the timing closer to those of seasonal CO₂ amount fraction cycles. The seasonal δ(APO) cycles from the without-SH-flux run in Fig. 7b also show similar features. These results suggest that the interhemispheric mixing of air modifies the seasonal δ(APO) cycles significantly, especially in the higher altitude in the lower latitude region. To compare the observed latitudinal and altitudinal variation in the seasonal δ(APO) cycle amplitude with that calculated by NICAM-TM in more detail, Fig. 8a shows a latitudinal distribution of average fractions of the observed seasonal δ(APO) and CO₂ amount fraction amplitudes relative to the 25.5^∘ N values. The decrease in the seasonal amplitude of δ(APO) toward the lower latitude is about 50 %, while that of the CO₂ amount fraction is less than 10 %; these are well reproduced by the control run of NICAM-TM. On the other hand, the without-SH-flux run yields a decrease in the δ(APO) amplitude by 20 % toward the lower latitude, which is significantly smaller than that from the δ(APO) from control run but slightly larger than that of the CO₂ amount fraction. These results indicate that the SH makes a significant contribution to the amplitude and phase of the seasonal δ(APO) cycle at lower latitude.

Figure 8b shows average fractions of the seasonal amplitudes of δ(APO) and the CO₂ amount fraction with altitude, relative to the corresponding values at 1.3 km. The average fractions of surface seasonal cycles, obtained from continuous observations of δ(O₂ $/$ N₂) and the CO₂ amount fraction at MNM since December 2015 (updated from Ishidoya et al., 2017), are also plotted. The seasonal amplitude of δ(APO) decreases rapidly with altitude by about 70 %, while that of the CO₂ amount fraction decreases by less than 20 %. The altitudinal dependence of the seasonal δ(APO) amplitude is not consistent with that of CO₂. This could be attributed to the fact that the seasonal CO₂ cycle is driven mainly by terrestrial biospheric activities in the northern midlatitudes and high latitudes and the seasonal minimum appears in summer, which makes the seasonal CO₂ cycle uniform from the surface to the middle troposphere due to continental strong convection in summer. On the other hand, the seasonal APO cycle is driven by air–sea flux, so that altitudinal homogenization of the seasonal APO cycle due to convection is weaker than that of CO₂. These features are well reproduced by the control run of NICAM-TM. The altitudinal decrease in the seasonal δ(APO) amplitude is also reproduced by the without-SH-flux run, although a slight underestimation is found above 5 km. Therefore, the altitudinal decrease in the seasonal δ(APO) amplitude over MNM is mainly due to an attenuation of the seasonal air–sea O₂ and N₂ fluxes around MNM with altitude, with some influence from the interhemispheric atmospheric mixing. Consequently, over the western North Pacific region, the altitude–latitude distribution of seasonal δ(APO) cycles is more sensitive to the atmospheric transport processes associated with interhemispheric air mixing and vertical attenuation of surface signal, compared with those of seasonal CO₂ amount fraction cycles.

Figure 9 (a) Latitudinal distribution of average deviations of the annual mean values of δ(APO) and the CO₂ amount fraction relative to those at 25.5^∘ N in the troposphere over the western North Pacific throughout the observation period (black filled circles). Error bands (shaded) indicate year-to-year variations (±1σ). The corresponding results calculated using NICAM-TM for the control run (blue triangles) and the without-SH-flux run (green triangles) are also shown. (b) Same as in (a) but for vertical distribution over MNM relative to the corresponding values at 1.3 km.

We also compared the observed and simulated annual mean values of δ(APO) and the CO₂ amount fraction. Figure 9a shows average deviations of the middle tropospheric annual mean values of δ(APO) and the CO₂ amount fraction at each latitude, relative to the 25.5^∘ N values. The deviation values of δ(APO) and the CO₂ amount fraction are well reproduced by the control and without-SH-flux runs. Therefore, the surface fluxes of O₂, N₂ and CO₂ in the Northern Hemisphere are the main contributors to the observed latitudinal variations in Fig. 9a. As discussed in connection with Eq. (11), we ignored δ^AM(APO) in our δ(APO) simulation using NICAM-TM, which is a component of δ(APO) driven by annual mean air–sea O₂ and N₂ fluxes. Therefore, the results of our simulation suggest that δ^AM(APO) does not affect significantly the latitudinal variations in the annual mean values of the middle tropospheric δ(APO) at 25–34^∘ N. Figure 9b shows the altitude deviations of the annual mean values of δ(APO) and the CO₂ amount fraction, relative to their corresponding values at 1.3 km over MNM. The observed profile of the CO₂ amount fraction is well reproduced by NICAM-TM. On the other hand, the average vertical gradient of δ(APO) profiles, obtained from the control and without-SH-flux runs of NICAM-TM, seems to be slightly larger than the observation. This may be due to the ignored contribution of δ^AM(APO); in that case, it may be that the sea area around MNM emits O₂ to the atmosphere throughout the observation period. Moreover, it is clearly seen from the figure that the interannual variation in the observed δ(APO) profiles is much larger than that in the NICAM-TM simulations. This may also be due to the ignored contribution of δ^AM(APO), and we discuss interannual variations in the observed and simulated δ(APO) in the section below.

Figure 10 (a) Same secular trends in Fig. 4a and their annual change rates along with the deseasonalized values of δ(APO) at each latitude over the western North Pacific. Error bands (shaded) indicate standard deviations of the deseasonalized δ (APO) from the secular trends (±1σ). (b) Same annual change rates in panel (a). The corresponding change rates for the APO obtained from the control run of NICAM-TM (middle) and those of δ^AM(APO) obtained by subtracting the calculated δ(APO) from the observed δ(APO) (bottom) are also shown.

3.2 Interannual variations in δ(APO) and its implication to global air–sea O₂ flux and CO₂ budget

In this section, we discuss causes of the interannual variations found in the middle tropospheric δ(APO) observed over the western North Pacific. Figure 10a shows the same secular trends in Fig. 4a and their annual change rates along with the deseasonalized values of δ(APO) at each latitude. The change rates observed at 29–34^∘ N show interannual variation with maxima around early 2015 and mid-2019, with a minimum at all latitudes in early 2017. The corresponding change rates obtained from the control run of NICAM-TM are also shown in Fig. 10b. The change rates obtained from NICAM-TM at 29–34^∘ N show interannual variations in phase with the observed rates, although the amplitudes are smaller by about 80 %. The interannual variations observed at 24–28^∘ N are also larger than the simulated values. Therefore, it is likely that interannual variation in the δ^AM(APO), which was not incorporated into NICAM-TM, is a main contributor to the observed interannual variations at various latitudes. In this connection, it is possible that the interannual variation in the global air–sea CO₂ flux could also contribute to the larger interannual variation in the observed δ(APO) since the NICAM-TM model incorporated only the monthly sea surface CO₂ flux climatology to calculate δ^OC(APO). However, the global air–sea CO₂ flux reported by the GCP (Friedlingstein et al., 2020) showed an interannual variation of 0.07 Pg a⁻¹ during 2012–2019, corresponding to 0.2 per meg a⁻¹ of δ^OC(APO), which is much smaller than the interannual variation in the observed δ(APO) shown in Fig. 10.

By assuming that all other components besides δ^AM(APO) are well represented in the NICAM-TM control run, we subtracted the change rates simulated by NICAM-TM from the observed rates, to extract the interannual variations due only to the δ^AM(APO). The calculated change rates of δ^AM(APO) are shown at the bottom of Fig. 10b. The change rates show similar interannual variations to the observed rates, but the latitudinal differences are smaller. This suggests that the interannual variations driven by δ^AM(APO) do not differ significantly as a function of latitude. In the following discussion, we make a bold assumption that an average of the change rates of δ^AM(APO) shown in Fig. 10b is a global average.

Figure 11 Anomaly in the average annual change rate of δ^AM(APO) shown at the bottom of Fig. 10 (thick black line). Anomaly in the change rate of APO driven only by the solubility change, expected from the observed surface δ($\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$) at Tsukuba, Japan (δ_therm(APO); red line, see text), and that driven by the net marine biospheric activities (δ_netbio(APO); green line), obtained by subtracting the change rate of δ_therm(APO) from that of δ^AM(APO), are shown. Anomaly in the global air–sea O₂ flux corresponding to the change rate of δ_netbio(APO) is also shown (see text). The time series of the NINO.WEST index (black open circles) and the annual average values of the index (thin black line) are shown at the bottom of the figure.

An anomaly of the average interannual variation of δ^AM(APO) change rate is shown in Fig. 11 (black line). In this figure, we also plotted a similar anomaly of interannual variation of the δ(APO) change rate due to solubility change (red line, hereafter referred to as “δ_therm(APO)”). The δ_therm(APO) was calculated from the δ(Ar $/$ N₂) measurements observed at Tsukuba (36^∘ N, 140^∘ E), Japan (Ishidoya et al., 2021), by multiplying a coefficient of about 0.9 derived from differences in the solubility in O₂ and Ar (Weiss, 1970). It should be noted that this is a rough approximation of the δ_therm(APO), which is a combination of air–sea fluxes of O₂, N₂ and CO₂ caused by solubility changes. As discussed in Ishidoya et al. (2021), the interannual variation in the δ(Ar $/$ N₂) change rate is in phase with the global ocean heat content reported by ocean temperature measurements (e.g., Levitus et al., 2012). This suggests that δ_therm(APO) is also driven by changes in the solubility of the global seawater. By subtracting δ_therm(APO) from δ^AM(APO), we estimated interannual variation of the δ(APO) change rate due to marine biological activities (green line, hereafter referred to as “δ_netbio(APO)”). It is expected that δ_netbio(APO) is driven by marine biological activities, not only in the surface mixed layer but also through a ventilation of subsurface low-O₂ waters.

Both the δ_therm(APO) and δ_netbio(APO) show significant interannual variations, roughly in opposite phase with each other. Moreover, the change rate of δ_netbio(APO) varies in opposite phase with the NINO.WEST (Japan Meteorological Agency; https://www.data.jma.go.jp/gmd/cpd/db/elnino/index/ninowidx.html, last access: 13 May 2020), which is an index of the El Niño–Southern Oscillation (ENSO). The ENSO is in the El Niño and La Niña phase, respectively, during the period with the negative and positive NINO.WEST index. Therefore, the δ_netbio(APO) tends to increase and decrease during El Niño and La Niña, respectively. This is consistent with Eddebbar et al. (2017), who examined global and tropical air–sea O₂ flux responses to ENSO, based on the Community Earth System Model (CESM). They reported that the upper ocean loses O₂ to the atmosphere during El Niño and gains O₂ during La Niña, mainly due to changes in ventilation of low-O₂ waters in the tropical Pacific, the region that has a dominant influence over the interannual variation in global air–sea O₂ flux (McKinley et al., 2003). By assuming the interannual variation in the δ_netbio(APO) represents a global average, and assuming a one-box atmosphere with 5.124×10²¹ g for the total mass of dry air (Trenberth, 1981), 28.97 g mol⁻¹ for the mean molecular weight of dry air, and respective fractions of 0.2093 and 0.7808 for O₂ and N₂ in the atmosphere, we estimated an interannual variation in the global air–sea O₂ flux due to marine biological activities (right axis of the green line in Fig. 11). The peak-to-peak amplitude of the O₂ flux is found to be about 300 Tmol a⁻¹, which is almost consistent with that of global APO flux estimated using the δ(APO) data from Scripps stations (Keeling and Manning, 2014) and a global atmospheric transport inversion (Rödenbeck et al., 2008; Eddebbar et al., 2017). Therefore, the seasonal cycles and interannual variations obtained from δ_cor.(O₂ $/$ N₂) have a coherent signal we can explain, although we recognize the number of simplifying assumptions is extensive to derive δ_netbio(APO), and the aircraft observations are from a comparatively small region. This helps to prove that the δ_cor.(O₂ $/$ N₂) is of good quality, despite the significant correction for enormous sampling artifacts.

By assuming the average secular trends of the middle tropospheric δ_cor.(O₂ $/$ N₂) and the CO₂ amount fraction observed in this study to be global average values, we were able to estimate the global CO₂ budget. The equations for separating out the global net terrestrial biospheric and oceanic CO₂ uptake are given by Keeling and Shertz (1992) firstly as

$$\begin{array}{}\text{(12)}& B={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{F}}}{{\mathit{\alpha }}_{\mathrm{B}}}}F+{\displaystyle \frac{\mathrm{1}}{\mathrm{0.471}}}{\displaystyle \frac{X\left({\mathrm{O}}_{\mathrm{2}}\right)}{{\mathit{\alpha }}_{\mathrm{B}}}}{\displaystyle \frac{\mathrm{d}\mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)}{\mathrm{d}t}}-{\displaystyle \frac{Z_{\mathrm{eff}}}{{\mathit{\alpha }}_{\mathrm{B}}}},\end{array}$$

and

$$\begin{array}{}\text{(13)}& \begin{array}{rl}O& =-{\displaystyle \frac{\mathrm{1}}{\mathrm{0.471}}}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(y\left({\mathrm{CO}}_{\mathrm{2}}\right)+{\displaystyle \frac{X\left({\mathrm{O}}_{\mathrm{2}}\right)}{{\mathit{\alpha }}_{\mathrm{B}}}}\mathit{\delta }\left({\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\right)\right)\\ & +{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{B}}-{\mathit{\alpha }}_{\mathrm{F}}}{{\mathit{\alpha }}_{\mathrm{B}}}}F+{\displaystyle \frac{Z_{\mathrm{eff}}}{{\mathit{\alpha }}_{\mathrm{B}}}}.\end{array}\end{array}$$

Here, B, F and O (in Pg a⁻¹, C equivalents) are the global terrestrial biospheric CO₂ uptake, the anthropogenic CO₂ emitted from fossil fuel combustion and cement manufacturing, and the oceanic CO₂ exchange, respectively; dδ(O₂ $/$ N₂)dt⁻¹ (per meg a⁻¹) and dy(CO₂)dt⁻¹ ($\mathrm{\mu }\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}{\mathrm{a}}^{-\mathrm{1}}$) are the observed change rates in atmospheric δ_cor.(O₂ $/$ N₂) and the CO₂ amount fraction, respectively; 0.471 converts CO₂ emissions of 1 Pg (C equivalents) to micromoles per mole of atmospheric CO₂; X(O₂) is the standard mole fraction of O₂ in air (0.2094), α_F and α_B are the OR for global average fossil fuel combustion and net terrestrial biospheric activities, respectively; and Z_eff (Pg a⁻¹) represents the net effect of oceanic O₂ outgassing on the oceanic and terrestrial biospheric CO₂ uptakes, which has been considered since Bender et al. (2005) and Manning and Keeling (2006). As described in the “Methods” section, we use α_F of 1.37 calculated from the fossil fuel and cement production emissions by fuel type reported by the GCP.

Long-term change in Z_eff is caused mainly by stratification of the ocean and the decrease of O₂ solubility in seawater due to a secular increase in the global ocean heat content (e.g., Bopp et al., 2002). However, as discussed above for Fig. 11, the ocean O₂ outgassing shows significant interannual variation, which makes it difficult to estimate year-to-year variations in the global CO₂ budget from Eqs. (12) and (13). In this regard, Tohjima et al. (2019) estimated the terrestrial biospheric and the oceanic CO₂ uptakes using their δ(O₂ $/$ N₂) and CO₂ amount fraction data, by changing the time period to obtain average secular change rates. They reported that the CO₂ uptakes estimated using the change rates averaged over a longer period greater than 5 years were consistent with those reported by the GCP (Le Queéreé et al., 2018) within ±0.5 Pg a⁻¹, while those estimated using annual change rates scattered significantly.

It seems that the averaging period needed to reduce the interannual variations in δ_therm(APO) and δ_netbio(APO) in Fig. 11 is about 4–5 years, similar to Nevison et al. (2008) and Tohjima et al. (2019). Therefore, we estimated average terrestrial biospheric and oceanic CO₂ uptake throughout the observation period (2012–2019, 8 years) to reduce the interannual variation in Z_eff sufficiently. Using the global ocean (0–2000 m) heat content data from the National Oceanographic Data Center (NOAA)/National Centers for Environmental Information (NCEI) (updated from Levitus et al., 2012, https://www.nodc.noaa.gov/OC5/3M_HEAT_CONTENT/, last access: 14 January 2020) and the same ratio of air–sea O₂ (N₂) flux to air–sea heat flux used in Manning and Keeling (2006), we adopted (0.6±0.6) Pg a⁻¹ for Z_eff for the period 2012–2019. We assumed 100 % uncertainty for Z_eff following Manning and Keeling (2006). The value of F, (9.7±0.5) Pg a⁻¹, was obtained from emissions from fossil fuel combustion and industrial processes by the GCP (Friedlingstein et al., 2020). Using the average secular trends of δ_cor.(O₂ $/$ N₂) and the CO₂ amount fraction for the observational period of the present study, the respective terrestrial biospheric and oceanic CO₂ uptakes were estimated to be (1.8±0.9) and (2.8±0.6) Pg a⁻¹ for the period 2012–2019. These values agree well with the corresponding CO₂ uptake of (1.8±1.1) and (2.6±0.5) Pg a⁻¹ reported by the GCP (Friedlingstein et al., 2020). It is noted that the terrestrial biospheric CO₂ uptake by the GCP is calculated by subtracting the CO₂ emission due to land-use change ((1.6±0.7) Pg a⁻¹) from their estimated total land CO₂ uptake ((3.4±0.9) Pg a⁻¹).

4 Summary

Regular air samples were taken on C-130 cargo aircraft flights from Atsugi Base to MNM, and air samples have been collected during the level flight and during the descent portion at MNM. In this paper, we have presented the analytical results of the air samples for the CO₂ amount fraction, δ(O₂ $/$ N₂), δ(Ar $/$ N₂), δ(¹⁵N) of N₂, δ(¹⁸O) of O₂ and δ(⁴⁰Ar) for the period May 2012–March 2020. The relationships of δ(Ar $/$ N₂), δ(¹⁸O) and δ(⁴⁰Ar) with δ(¹⁵N) indicate a significant artificial fractionation due to thermal diffusion during the air sample collection.

The δ_cor.(O₂/N₂) values, corrected for the artificial fractionation using δ(Ar $/$ N₂), and the δ(APO) values derived from δ_cor.(O₂ $/$ N₂) were shown to have clear seasonal cycles nearly in opposite phase to the cycle of the CO₂ amount fraction from the surface to 6 km along the latitudinal path from 25.5 to 33.5^∘ N. We then used a three-dimensional atmospheric transport model NICAM-TM that was driven by the air–sea fluxes of O₂, N₂ and CO₂, along with fluxes of CO₂ and O₂ from fossil fuel combustion, to interpret some of the characteristic features we observed in the seasonal cycles and vertical profiles of APO and the CO₂ amount fraction.

Seasonal amplitudes of δ(APO) and the CO₂ amount fraction decreased toward the lower latitude from 34.5 to 24.5^∘ N by about 50 % and less than 10 %, respectively: these features were reproduced by the corresponding ratios from the control run of NICAM-TM. On the other hand, the without-SH-flux run underestimated the latitudinal change in the δ(APO) amplitude, which indicated that the seasonal cycle of the mid-tropospheric δ(APO) was modified significantly by a combination of the northern and southern hemispheric seasonal cycles through the interhemispheric atmospheric mixing. The decrease in the δ(APO) seasonal amplitude was about 70 % with altitude from the surface to 6 km, while that of the CO₂ amount fraction was less than 20 %. These features were also reproduced well by the control run of NICAM-TM.

The observed decrease in the annual mean values of the CO₂ amount fraction with altitude was reproduced by the control run of NICAM-TM. On the other hand, the average vertical gradient of the δ(APO) profiles was slightly overestimated by the NICAM-TM simulations, while the simulated interannual variation was underestimated. This may be due to the ignored contribution of δ^AM(APO), which is a component of δ(APO) driven by annual mean air–sea O₂ and N₂ fluxes.

The interannual variations in the middle tropospheric δ^AM(APO) were estimated by subtracting the simulated δ(APO) by the NICAM-TM control run from the observed δ(APO). We also estimated the solubility-driven component of δ(APO) (δ_therm(APO)) from the δ(Ar $/$ N₂) observed at Tsukuba, assuming its interannual variation was driven by changes in the globally averaged solubility of the seawater. The interannual variation in δ(APO) driven by marine biological activities (δ_netbio(APO)) was calculated by subtracting δ_therm(APO) from δ^AM(APO). The δ_netbio(APO) showed significant interannual variations in the opposite phase to that of δ_therm(APO), and the change rate varied in opposite phase with the NINO.WEST. Therefore, the δ_netbio(APO) values obtained in this study tended to increase and decrease with El Niño and La Niña, respectively, which is in agreement with Eddebbar et al. (2017), who examined responses of the global and tropical air–sea O₂ flux to ENSO based on the CESM.

By assuming the observed secular trends of the middle tropospheric δ_cor.(O₂ $/$ N₂) and CO₂ amount fraction to be representative of global average values, we estimated terrestrial biospheric and oceanic CO₂ uptakes to be (1.8±0.9) and (2.8±0.6) Pg a⁻¹, respectively, for the period 2012–2019. These values agree well with the corresponding CO₂ uptake values of (1.8±1.1) (land) and (2.6±0.5) Pg a⁻¹ (ocean) reported by the GCP.

Additionally, our study has shown that our aircraft observation and the method we used to correct artificial fractionation of O₂ and N₂ are useful in evaluating interhemispheric air mixing processes based on the seasonal δ(APO) cycle, as well as interannual variations in the global air–sea O₂ flux, and in calculating global CO₂ budgets based on the long-term trends of δ_cor.(O₂ $/$ N₂) and the CO₂ amount fraction.