Using a speed-dependent Voigt line shape to retrieve O₂ from Total Carbon Column Observing Network solar spectra to improve measurements of XCO₂

High-resolution, laboratory, absorption spectra of the $a^{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{g}}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ oxygen (O₂) band measured using cavity ring-down spectroscopy were fitted using the Voigt and speed-dependent Voigt line shapes. We found that the speed-dependent Voigt line shape was better able to model the measured absorption coefficients than the Voigt line shape. We used these line shape models to calculate absorption coefficients to retrieve atmospheric total columns abundances of O₂ from ground-based spectra from four Fourier transform spectrometers that are a part of the Total Carbon Column Observing Network (TCCON). Lower O₂ total columns were retrieved with the speed-dependent Voigt line shape, and the difference between the total columns retrieved using the Voigt and speed-dependent Voigt line shapes increased as a function of solar zenith angle. Previous work has shown that carbon dioxide (CO₂) total columns are better retrieved using a speed-dependent Voigt line shape with line mixing. The column-averaged dry-air mole fraction of CO₂ (XCO₂) was calculated using the ratio between the columns of CO₂ and O₂ retrieved (from the same spectra) with both line shapes from measurements taken over a 1-year period at the four sites. The inclusion of speed dependence in the O₂ retrievals significantly reduces the air mass dependence of XCO₂, and the bias between the TCCON measurements and calibrated integrated aircraft profile measurements was reduced from 1 % to 0.4 %. These results suggest that speed dependence should be included in the forward model when fitting near-infrared CO₂ and O₂ spectra to improve the accuracy of XCO₂ measurements.

1 Introduction

Accurate remote sensing of greenhouse gases (GHGs) such as CO₂, in the Earth's atmosphere is important for studying the carbon cycle to better understand and predict climate change. The absorption of solar radiation by O₂ in the Earth's atmosphere is important because it can be used to study the properties of clouds and aerosols and to determine vertical profiles of temperature and surface pressure. Wallace and Livingston (1990) were the first to retrieve total columns of O₂ from some of the discrete lines of the $a^{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{g}}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ band of O₂ centred at 1.27 µm (which will be referred to bellow as the 1.27 µm band) using atmospheric solar absorption spectra from Kitt Peak National Observatory. Mlawer et al. (1998) recorded solar absorption spectra in the near-infrared (NIR) region to study collision-induced absorption (CIA) in the $a^{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{g}}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ band as well as two other O₂ bands. The spectra were compared to a line-by-line radiative transfer model and the differences between the measured and calculated spectra showed the need for better absorption coefficients in order to accurately model the 1.27 µm band (Mlawer et al., 1998). Subsequently, spectroscopic parameters needed to calculate the absorption coefficients from discrete transitions of the 1.27 µm band were measured in multiple studies (Cheah et al., 2000; Newman et al., 1999, 2000; Smith and Newnham, 2000), as was collision-induced absorption (CIA) (Maté et al., 1999; Smith and Newnham, 2000), while Smith et al. (2001) validated the work done in Smith and Newnham (2000) using solar absorption spectra.

The 1.27 µm band is of particular importance to the Total Carbon Column Observing Network (TCCON) (Wunch et al., 2011). TCCON is a ground-based remote sensing network that makes accurate and precise measurements of GHGs for satellite validation and carbon cycle studies. Using the O₂ column retrieved from solar absorption spectra, the column-averaged dry-air mole fraction of CO₂ (XCO₂) has been shown to provide better precision than when using the surface pressure to calculate XCO₂ (Yang et al., 2002). The O₂ column is retrieved from the 1.27 µm band because of its close proximity to the spectral lines used to retrieve CO₂, thereby reducing the impact of solar tracker mispointing and an imperfect instrument line shape (ILS) (Washenfelder et al., 2006a). To improve the retrievals of O₂ from the 1.27 µm band, Washenfelder et al. (2006a) found that adjusting the spectroscopic parameters in HITRAN 2004 (Rothman et al., 2005) decreased the air mass and temperature dependence of the O₂ column. These revised spectroscopic parameters were included in HITRAN 2008 (Rothman et al., 2009). Atmospheric solar absorption measurements from this band taken at the Park Falls TCCON site by Washenfelder et al. (2006a) were the first measurements to observe the electric-quadrupole transitions (Gordon et al., 2010). Leshchishina et al. (2011, 2010) subsequently used cavity ring-down spectra to retrieve spectroscopic parameters for the 1.27 µm band using a Voigt spectral line shape and these parameters were included in HITRAN 2012 (Rothman et al., 2013a). Spectroscopic parameters for the discrete spectral lines of the O₂ 1.27 µm band from HITRAN 2016 (Gordon et al., 2017) are very similar to HITRAN 2012 except that HITRAN2016 includes improved line positions reported by Yu et al. (2014).

Extensive spectral line shape studies have been performed for the O₂ A-band, which is centred at 762 nm and used by the Greenhouse Gases Observing Satellite (GOSAT) (Yokota et al., 2009) and the Orbiting Carbon Observatory-2 (OCO-2) satellite (Crisp et al., 2004) to determine surface pressure. Studies have shown that the Voigt line shape is inadequate to describe the spectral line shape of the discrete O₂ lines in the A-band. Dicke narrowing occurs when the motion of the molecule is diffusive due to collisions changing the velocity and direction of the molecule during the time that it is excited. This diffusive motion is taken into account by averaging over many different Doppler states, resulting in a line width that is narrower than the Doppler width (Dicke, 1953). Long et al. (2010) and Predoi-Cross et al. (2008) found it necessary to use a spectral line shape model that accounted for Dicke narrowing when fitting the discrete lines of the O₂ A-band. Line mixing, which occurs when collisions transfer intensity from one part of the spectral band to another (Lévy et al., 1992), was shown to be prevalent in multiple studies (Predoi-Cross et al., 2008; Tran et al., 2006; Tran and Hartmann, 2008). Tran and Hartmann (2008) showed that including line mixing when calculating the O₂ A-band absorption coefficients reduced the air mass dependence of the O₂ column retrieved from TCCON spectra. When fitting cavity ring-down spectra of the O₂ A-band, Drouin et al. (2017) found it necessary to use a speed-dependence Voigt line shape, which takes into account different speeds at the time of collision (Shannon et al., 1986), with line mixing to properly fit the discrete spectral lines of the O₂ A-band.

The need to include non-Voigt effects when calculating absorption coefficients for the O₂ 1.27 µm band was first shown in Hartmann et al. (2013) and Lamouroux et al. (2014). In Hartmann et al. (2013) and Lamouroux et al. (2014), Lorentzian widths were calculated using the requantised classical molecular dynamics simulations (rCMDSs) and used to fit cavity ring-down spectra with a Voigt line shape for some isolated transitions in the O₂ 1.27 µm band. The studies concluded that a Voigt line shape is insufficient for modelling the spectral lines of the O₂ 1.27 µm band and that effects such as speed dependence and Dicke narrowing should be included in the line shape calculation.

In this study, air-broadened laboratory cavity ring-down spectra of the O₂ 1.27 µm band were fitted using a spectral line shape that takes into account speed dependence. The derived spectroscopic parameters for the speed-dependent Voigt line shape were used to calculate absorption coefficients when fitting high-resolution solar absorption spectra. XCO₂ was calculated from O₂ total columns retrieved using the new absorption coefficients and CO₂ total columns retrieved with the line shape model described in Mendonca et al. (2016). These new XCO₂ values were compared to the XCO₂ retrieved using the Voigt line shape. Section 2 details the formulas used to calculate absorption coefficients using different spectral line shapes. In Sect. 3, we describe the retrieval of spectroscopic parameters from three air-broadened cavity ring-down spectra fitted with a speed-dependent Voigt line shape. For Sect. 4, the speed-dependent line shape along with the retrieved spectroscopic parameters is used to fit solar absorption spectra from four TCCON sites and retrieve total columns of O₂, which is compared to O₂ retrieved using a Voigt line shape. In Sect. 5, we investigate the change in the air mass dependence of XCO₂ with the new O₂ retrievals. In Sect. 6, we discuss our results and their implications for remote sensing of greenhouse gases.

2 Absorption coefficient calculations

2.1 Voigt line shape

The Voigt line shape is the convolution of the Lorentz and the Gaussian profiles, which model pressure and Doppler broadening of the spectral line respectively. The corresponding absorption coefficient, k, at a given wave number v becomes

$$\begin{array}{}\text{(1)}& k\left(v\right)=N{\sum }_j{S}_j\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\gamma }}_{{\mathrm{D}}_j}}}\right){\left({\displaystyle \frac{\ln \left(\mathrm{2}\right)}{\mathit{\pi }}}\right)}^{\mathrm{1}/\mathrm{2}}\left(\mathit{\text{Re}}\left[c\left(v,x_j,y_j\right)\right]\right),\end{array}$$

where N is the number density, S_j is the line intensity of spectral line j, ${\mathit{\gamma }}_{{\mathrm{D}}_j}$ is the Doppler half-width (HWHM), c is the complex error function, and

$$\begin{array}{}\text{(2)}& x_j={\displaystyle \frac{\left(v-v_j^o-P{\mathit{\delta }}_j^o\right)}{{\mathit{\gamma }}_{{\mathrm{D}}_j}}}{\left(\ln \left(\mathrm{2}\right)\right)}^{\mathrm{1}/\mathrm{2}},y_j={\displaystyle \frac{{\mathit{\gamma }}_{{\mathrm{L}}_j}}{{\mathit{\gamma }}_{{\mathrm{D}}_j}}}{\left(\ln \left(\mathrm{2}\right)\right)}^{\mathrm{1}/\mathrm{2}}.\end{array}$$

Here, $v_j^o$ is the position of the spectral line j, P is the pressure, and ${\mathit{\delta }}_j^o$ is the pressure-shift coefficient. The Lorentz half-width, ${\mathit{\gamma }}_{{\mathrm{L}}_j}$, is calculated using the following:

$$\begin{array}{}\text{(3)}& {\mathit{\gamma }}_{{\mathrm{L}}_j}\left(T\right)=P{\mathit{\gamma }}_{{\mathrm{L}}_j}^o{\left({\displaystyle \frac{\mathrm{296}}{T}}\right)}^n,\end{array}$$

where ${\mathit{\gamma }}_{{\mathrm{L}}_j}^o$ is the air-broadened Lorentz half-width coefficient (at reference temperature 296 K) and n is the exponent of temperature dependence. The Voigt line shape assumes that pressure broadening is accurately represented by a Lorentz profile calculated for the statistical average velocity at the time of collision.

2.2 Speed-dependent Voigt line shape

The speed-dependent Voigt line shape refines the pressure-broadening component of the Voigt by calculating multiple Lorentz profiles for different speeds at the time of collision. The final contribution from pressure broadening to the speed-dependent Voigt is the weighted sum of Lorentz profiles (weighted by the Maxwell–Boltzmann speed distribution) calculated for different speeds at the time of collision. The speed-dependent Voigt line shape (Ciuryło, 1998) with the quadratic representation of the Lorentz width and pressure shift (Rohart et al., 1994) is as follows:

$$\begin{array}{ll}\text{(4)}& {\displaystyle }& {\displaystyle }k\left(v\right)={\displaystyle }& {\displaystyle }N\left({\displaystyle \frac{\mathrm{2}}{{\mathit{\pi }}^{\frac{\mathrm{3}}{\mathrm{2}}}}}\right){\sum }_j{S}_j\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}{e}^{-V^{\mathrm{2}}}V\left({\tan }^{-\mathrm{1}}\left[{\displaystyle \frac{x_j-B{a}_{{\mathit{\delta }}_j}\left(V^{\mathrm{2}}-\mathrm{1.5}\right)+V}{y_j\left(\mathrm{1}+a_{{\mathit{\gamma }}_{{\mathrm{L}}_j}}\left(V^{\mathrm{2}}-\mathrm{1.5}\right)\right)}}\right]\right)\mathrm{d}V,\end{array}$$

where $a_{{\mathit{\gamma }}_{{\mathrm{L}}_j}}$ is the speed-dependent Lorentz width parameter (unitless) for line j, $a_{{\mathit{\delta }}_j}$ is the speed-dependent pressure-shift parameter (unitless), B is $\frac{{\left(\ln \left(\mathrm{2}\right)\right)}^{\mathrm{1}/\mathrm{2}}}{{\mathit{\gamma }}_{{\mathrm{D}}_j}}$, V is the ratio of the absorbing molecule's speed to the most probable speed of the absorbing molecule, and all other variables were defined before.

Figure 1 (a) Cavity ring-down absorption coefficients for O₂ measured at the three pressures indicated in the legend at approximately room temperature and a volume mixing ratio of 0.20720 (43). The difference between measured absorption coefficients and those calculated using (b) a Voigt line shape and (c) the speed-dependent Voigt line shape. Note the difference in scale between panels (b) and (c).

3 Fitting laboratory spectra

O₂, unlike CO₂ and CH₄, cannot produce an electric dipole moment and therefore should not be infrared active. However, O₂ has two unpaired electrons in the ground state that produce a magnetic dipole moment. Due to the unpaired electrons in the ground state ($X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$) the rotational state (N) is split into three components which are given by $J=N-\mathrm{1}$, J=N and $J=N+$1, while in the upper state (a¹Δ_g), J=N. When labelling a transition, the following nomenclature is used: $\mathrm{\Delta }N\left(N^{\prime \prime }\right)\mathrm{\Delta }J\left(J^{\prime \prime }\right)$ (Leshchishina et al., 2010), where ΔNis the difference between N^′ in the upper state and $N^{\prime \prime }$ in the lower state; ΔJ is the difference between J^′ in the upper state and $J^{\prime \prime }$ in the lower state. The magnetic transitions of $a^{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{g}}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ allow for ΔJ=0, ±1. This leads to nine observed branches: $P\left(N^{\prime \prime }\right)Q\left(J^{\prime \prime }\right)$, $R\left(N^{\prime \prime }\right)Q\left(J^{\prime \prime }\right)$ and $Q\left(N^{\prime \prime }\right)Q\left(J^{\prime \prime }\right)$ for ΔJ=0; $O\left(N^{\prime \prime }\right)P\left(J^{\prime \prime }\right)$, $P\left(N^{\prime \prime }\right)P\left(J^{\prime \prime }\right)$, and $Q\left(N^{\prime \prime }\right)P\left(J^{\prime \prime }\right)$ for $\mathrm{\Delta }J=-\mathrm{1}$; and $S\left(N^{\prime \prime }\right)R\left(J^{\prime \prime }\right)$, $R\left(N^{\prime \prime }\right)R\left(J^{\prime \prime }\right)$ and $Q\left(N^{\prime \prime }\right)R\left(J^{\prime \prime }\right)$ for ΔJ=1.

Absorption coefficients for three room temperature air-broadened (NIST Standard reference materal^® 2659a containing 79.28 % N₂, 20.720 (43) % O₂, 0.0029 % Ar, 0.00015 % H₂O and 0.001 % other compounds) spectra were measured at the National Institute of Standards and Technology (NIST) using the frequency-stabilised cavity ring-down spectroscopy (FS-CRDS) technique (Hodges et al., 2004; Hodges, 2005). The absorption spectra were acquired at pressures of 131, 99.3 and 66.9 kPa, at temperatures of 296.28, 296.34 and 296.30 K respectively. Figure 1a shows the three measured absorption spectra. A more detailed discussion of the present FS-CRDS spectrometer can be found in Lin et al. (2015).

Figure 2 The same as Fig. 1 but expanded to show four spectral lines in the P branch of the O₂ 1.27 µm band.

The spectra were fitted individually using a Voigt line shape (Eq. 1), with S_j, ${\mathit{\gamma }}_{{\mathrm{L}}_j}^o$ and ${\mathit{\delta }}_j^o$ for the main isotope of the magnetic dipole lines of the O₂ 1.27 µm band for lines with an intensity greater than $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{28}}$ cm⁻¹ (molecule cm⁻²)⁻¹. The spectroscopic parameters measured in Leshchishina et al. (2011) for the spectral lines of interest were used as the a priori for the retrieved spectroscopic parameters. The line positions were left fixed to the values measured in Leshchishina et al. (2011), and all other O₂ spectral lines (intensity less $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{28}}$ cm⁻¹ (molecule cm⁻²)⁻¹) were calculated using a Voigt line shape with spectroscopic parameters from HITRAN 2012 (Rothman et al., 2013a). Spectral fits were done using the lsqnonlin function in Matlab, with a user-defined Jacobian matrix. The Jacobian was constructed by taking the derivative of the absorption coefficients with respect to the parameters of interest. Using an analytical Jacobian instead of the finite difference method is both computationally faster and more accurate. The Voigt line shape was calculated using the Matlab code created by Abrarov and Quine (2011) to calculate the complex error function and its derivatives. To take collision-induced absorption (CIA) into account, a set of 50 Legendre polynomials were added together by retrieving the weighting coefficients needed to add the polynomials to fit the CIA for each spectrum. Figure 1b shows the residual (measured minus calculated absorption coefficients) when using a Voigt line shape with the retrieved spectroscopic parameters. The plot shows that residual structure still remains for all three spectra. The root mean square (rms) residual values for the spectra are given by the legend at the side of the plot.

Figure 3 The averaged measured (a) intensity, (b) Lorentz line width, (c) pressure shift and (d) speed-dependent pressure shift retrieved from the three cavity ring-down spectra of the 1.27 µm band of O₂. All data are plotted as a function of m, which is $m=-J$ for the P-branch lines, m=J for the Q branch and $m=J+\mathrm{1}$ for the R branch (where J is the lower state rotational quantum number) and the uncertainties shown are 2σ.

Figure 2 is the same plot as Fig. 1 but for the P(11)P(11), P(11)Q(10), P(9)P(9) and P(9)Q(8) spectral lines only. Figure 2b shows that for all four spectral lines there is a W-shaped residual at the line centre. The P(11)P(11) line was also measured by Hartmann et al. (2013) at pressures ranging from 6.7 to 107 kPa. Figure 5 of Hartmann et al. (2013) shows the P(11)P(11) line at a pressure of 66.7 kPa, which is approximately the pressure of the 66.9 kPa spectrum (blue spectrum in Figs. 1 and 2). When one compares the blue residual of the P(11)P(11) line in Fig. 2b to that of the residual of the left panel of Fig. 5 of Hartmann et al. (2013), one can see that the residuals are the same. Figure 6 of Hartmann et al. (2013) show that the amplitude of the residual increases with decreasing pressure, which is also seen in Fig. 2b. Figure 3 of Lamouroux et al. (2014) shows the same W residual for the P(9)P(9) lines and that the amplitude of the residual increases with decreasing pressure (although for lower pressures), consistent with the results for the P(9)P(9) line in Fig. 2b.

Figure 1c shows the residual when using the speed-dependent Voigt (Eq. 4) to fit each spectrum individually. Using Eq. (4) requires integration over all possible speeds, which is not computationally practical, so we employ the simple numerical integration scheme as was done by Wehr (2005). When fitting the spectra, parameters S_j, ${\mathit{\gamma }}_{{\mathrm{L}}_j}^o$, ${\mathit{\delta }}_j^o$, $a_{{\mathit{\gamma }}_{{\mathrm{L}}_i}}$ and $a_{{\mathit{\delta }}_j}$ were retrieved for lines of intensity greater than $\mathrm{7.0}\times {\mathrm{10}}^{-\mathrm{28}}$ cm⁻¹ (molecule cm⁻²)⁻¹, while all other O₂ lines were calculated using a Voigt line shape and spectroscopic parameters from HITRAN 2012 (Rothman et al., 2013b). The retrieved spectroscopic parameters are available in the Supplement. The Jacobian matrix was created by taking the derivative with respect to each parameter of interest, as was done with the Voigt fits. By taking speed-dependent effects into account, the residuals were reduced to 25 times smaller than those for the Voigt fit and the rms residuals (given in the legend of Fig. 1c) are 10 times smaller. However, some residual structure still remains, which is more evident in the in the Q and R branches than the P branch. Figure 2c shows the four lines in the P branch, as discussed when analysing the Voigt fits. A small residual W remains at the line centre as well as residuals from weak O₂ lines.

Figure 3 shows the averaged intensity Lorentz width coefficient, pressure shift coefficient and speed-dependent shift coefficient of the 1.27 µm O₂ band, retrieved from the three spectra, plotted as a function of quantum number m, which is $m=-J$ (where J is the lower state rotational quantum number) for the P-branch lines, m=J for the Q-branch lines and $m=J+\mathrm{1}$ for the R-branch lines. The intensity, Lorentz widths and pressure shifts show a m dependence for these parameters for the P and R sub-branches. The measured Lorentz widths and pressure shifts for the Q sub-branches show a m dependence but are not as strong as the P and R sub-branches. This is because the Q branch lines are broadened enough to blend with each other, since they are spaced closer together than the P or R branch lines. Figure 1c shows that some of the residual structure in the Q branch increases with pressure and is partly due to the blending of these transitions as the pressure increases. The weak O₂ absorption lines also blend in with the Q branch, contributing to the residual structure in Fig. 1c. We tried retrieving the spectroscopic parameters for the weak O₂ absorption lines, but since they were overlapping with the strong O₂ lines, it was not possible. Figure 4a shows the retrieved speed-dependent width parameter averaged over the three spectra, plotted as a function of m, showing that it increases with m. Error bars correspond to the 2σ standard deviation and are large regardless of sub-branch. Figure 4b shows the retrieved speed-dependent width for the PQ sub-branch for the different pressures. The speed-dependent width shows the same m dependence regardless of pressure, but also increases with decreasing pressure as is the case for sub-branches. It should be noted that the speed-dependent width parameter should be independent of pressure.

Figure 4 (a) The averaged measured speed-dependent width parameter of the 1.27 µm band of O₂ plotted as a function of m. (b) The measured speed-dependent width parameter for spectral lines that belong to the PQ sub-branch plotted as a function of m.

4 Fitting solar spectra

High-resolution solar absorption spectra were measured at four TCCON sites using a Bruker IFS 125HR FTIR spectrometer with a room temperature InGaAs detector at a spectral resolution of 0.02 cm⁻¹ (45 cm maximum optical path difference). The raw interferograms recorded by the instrument were processed into spectra using the I2S software package (Wunch et al., 2015) that corrects solar intensity variations (Keppel-Aleks et al., 2007), phase errors (Mertz, 1967) and laser sampling errors (Wunch et al., 2015), and then performs a fast Fourier transform to convert the interferograms into spectra (Bergland, 1969). The GGG software package (Wunch et al., 2015) is used to retrieve total columns of atmospheric trace gases. GFIT is the main code that contains the forward model, which calculates a solar absorption spectrum using a line-by-line radiative transfer model and an iterative non-linear least square fitting algorithm that scales an a priori gas profile to obtain the best fit to the measured spectrum. A priori profiles for GHGs are created by an empirical model in GGG that is based on measurements from the balloon-borne JPL MkIV Fourier transform spectrometer (FTS) (Toon, 1991), the Atmospheric Chemistry Experiment (ACE) FTS instrument aboard the SCISAT satellite (Bernath et al., 2005) and in situ GLOBALVIEW data (Wunch et al., 2011). Temperature and pressure profiles, as well as H₂O a priori profiles are generated from the National Centers for Environmental Prediction (NCEP) data. The calculations are performed for 71 atmospheric layers (0 to 70 km), so all a priori profiles are generated on a vertical grid of 1 km.

Figure 5 (a) The residuals (measured minus calculated) for a spectrum measured at Eureka on 27 March 2015 at a SZA of 81.32^∘. The red residual is the result of using the Voigt line shape and the blue is from using the qSDV. (b) The measured (red dots) and calculated values (black dots), with the qSDV spectrum, along with the gases included in the fit (refer to the legend to the right) in the spectral window.

In the current GGG software package (Wunch et al., 2015), the forward model of GFIT calculates absorption coefficients for the discrete lines of the O₂ 1.27 µm band using a Voigt line shape and spectroscopic parameters from Washenfelder et al. (2006a) and Gordon et al. (2010). To take CIA into account, absorption coefficients are calculated using a Voigt line shape and spectroscopic parameters from the foreign-collision-induced absorption (FCIA) and self-collision-induced absorption (SCIA) spectral line lists provided with the GGG software package (Wunch et al., 2015). Spectroscopic parameters in the FCIA and SCIA line lists were retrieved by Geoff C. Toon by fitting the laboratory spectra of Smith and Newnham (2000). This was done by retrieving the integrated absorption at every 1 cm⁻¹ of the spectrum and using a Voigt line shape, with fixed Lorentz width and no pressure shift. In GFIT, a volume scale factor is retrieved for the CIA and discrete lines separately so that the O₂ column is derived from the discrete lines of the 1.27 µm band only. Airglow is not considered when fitting the 1.27 µm band, since the spectrometer views the sun directly and airglow is overwhelmed by this a bright source. The continuum level and tilt of the 100 % transmission level is fitted using a weighted combination of the first two Legendre polynomials. Absorption coefficient for all other trace gases are calculated using a Voigt line shape and spectroscopic parameters from the atm.101 line list (Toon, 2014a) and solar lines are fitted using the solar line list (Toon, 2014b).

Figure 5 shows the spectral fit to a solar absorption spectrum recorded at Eureka on 27 March 2015, at a solar zenith angle (SZA) of 81.32^∘ (air mass of 6.3). This spectrum is an average of five Eureka scans. The TCCON standard is single scan but five scans were averaged to decrease the noise. The measured spectrum (red circles), calculated spectrum (black circles) and transitions from all gases in the window (coloured lines, refer to the legend for different gases) are shown in Fig. 5b. The residual obtained using a Voigt line shape to calculate the discrete lines of the O₂ 1.27 µm band is shown in red in Fig. 5a. The blue residual is the result of using a speed-dependent Voigt line shape with the spectroscopic parameters retrieved from fitting the absorption coefficients in Section 3. To decrease the amount of time it takes to calculate the absorption coefficients, the quadratic speed-dependent Voigt (qSDV) computational approach of Ngo et al. (2013) and Tran et al. (2013) was used instead of Eq. (4) since it requires the Voigt calculation only twice, while Eq. (4) requires numerical integration scheme with 33 iterations. The temperature-dependent parameter of the Lorentz width of the discrete lines of the O₂ 1.27 µm band reported in HITRAN 2012 was used to take temperature dependence into account for ${\mathit{\gamma }}_{{\mathrm{L}}_j}\left(T\right)$. There was only a slight improvement in the fit residuals with the new absorption coefficients (using the qSDV), as seen in Fig. 5a. Absorption coefficients calculated with the qSDV were used to retrieve total columns of O₂ from solar spectra recorded over a 1-year period at TCCON sites in Eureka (eu) (Nunavut, Canada) (Batchelor et al., 2009; Strong et al., 2017), Park Falls (pa) (Wisconsin, USA) (Washenfelder et al., 2006a; Wennberg et al., 2017a) , Lamont (oc) (Oklahoma, USA) (Wennberg et al., 2017b) and Darwin (db) (Australia) (Deutscher et al., 2010; Griffith et al., 2017). In total 131 124 spectra were fitted using the qSDV and the average root mean square (rms) residual of the fit only decreased by 0.5 %.

5 Impact of O₂ columns on XCO₂ measurements

The O₂ column retrieved from the 1.27 µm band with a Voigt line shape and spectroscopic parameters from the atm.101 line list (Toon, 2014a) has an air mass dependence such that the O₂ column retrieved increases as a function of the solar zenith angle (or air mass). Using spectra recorded from Eureka, Park Falls, Lamont and Darwin over 1-year periods, total columns of O₂ were retrieved using (1) a Voigt spectral line shape with spectroscopic parameters from the atm.101 line list and (2) the qSDV with the spectroscopic parameters determined in Sect. 3. Figure 6 shows the percent difference calculated as the column from the qSDV retrieval minus the column from the Voigt retrieval, which was then divided by the latter and multiplied by 100, plotted as a function of solar zenith angle (SZA). At the smallest SZA, the qSDV retrieves 0.75 % less O₂ than the Voigt, with the difference increasing to approximately 1.8 % as the SZA approaches 90^∘. Figure 7 shows XAIR for the entire data set plotted as a function of SZA. XAIR is the column of air (determined using surface pressure recorded at the site) divided by the column of O₂ retrieved from the spectra and multiplied by 0.2095, which is the dry-air mole fraction of O₂ in the Earth's atmosphere. Ideally XAIR should be 1 but when using O₂ retrieved with a Voigt line shape (Fig. 7a) to calculate XAIR the average XAIR value for the entire data set is 0.977. Using O₂ retrieved with the qSDV to calculate XAIR, the average value is 0.986, which is closer to the expected value of 1. However, XAIR has a dependence on SZA regardless of line shape used. Figure 7a shows that XAIR decreases as a function of SZA (evident at SZA >75^∘), which means that the retrieved column of O₂ increases as a function of SZA. Figure 7b shows that XAIR increases as a function of SZA (evident at SZA >70^∘), which means that the retrieved column of O₂ decreases as a function of SZA. Therefore retrieving total columns of O₂ with the qSDV changes the air mass dependence of the O₂ column, which in turn will impact the air mass dependence of XCO₂.

Figure 6 The percent difference between the O₂ column retrieved with the Voigt and qSDV line shapes for a year of measurements from Eureka (eu), Park Falls (pa), Lamont (oc) and Darwin (db).

Figure 7 (a) XAIR as a function of SZA calculated using the total column of O₂ retrieved using the Voigt line shape. Panel (b) is the same as panel (a) except the total column of O₂ was retrieved with the qSDV.

Figure 8 XCO₂ calculated from the CO₂ and O₂ columns retrieved from Park Falls spectra recorded on 18 June 2013. The CO₂ columns were retrieved using either the Voigt line shape or the qSDV with line mixing, while the O₂ columns were retrieved using either the Voigt or qSDV line shapes. XCO₂ was calculated in four ways: (1) both CO₂ and O₂ columns retrieved using the Voigt line shape (red), (2) CO₂ columns retrieved with the Voigt and O₂ columns retrieved with the qSDV (green), (3) CO₂ columns retrieved with the qSDV and line mixing and O₂ columns retrieved with the Voigt (cyan) and (4) CO₂ columns retrieved with the qSDV and line mixing and O₂ columns retrieved with the qSDV (blue). The top x axis is the SZA that corresponds to the hour on the bottom x axis.

5.1 Air mass dependence of XCO₂

Since the standard TCCON XCO₂ (and all other XGas) is calculated using the column of O₂ instead of the surface pressure, errors associated with the retrieval of O₂, such as the air mass dependence of the O₂ column, will affect XCO₂. Figure 8 is XCO₂ calculated for four different combinations pertaining to the two CO₂ column retrievals and the O₂ column retrievals. The CO₂ columns were retrieved with either a Voigt line shape (the standard GGG2014 approach) or the qSDV with line mixing as done in Mendonca et al. (2016), while the O₂ columns were retrieved with either a Voigt (the standard GGG2014 approach) or the new qSDV approach developed here. Figure 8 shows a spurious symmetric component to XCO₂ when the total column of O₂ is retrieved with the Voigt line shape, regardless of line shape used to retrieve CO₂. When the qSDV is used to retrieve total columns of O₂, the symmetric component of XCO₂ is dismissed regardless of line shape used to retrieve CO₂. This is because the air mass dependence of the column of O₂ retrieved using the qSDV is more consistent with the air mass dependence of the column of CO₂ (for both line shapes used to retrieve CO₂). Mendonca et al. (2016) showed that using the qSDV with line mixing results in better fits to the CO₂ windows and impacts the air mass dependence of the retrieved column of CO₂. When using a Voigt line shape the retrieved column amount of CO₂ decreases as air mass increases until the air mass is large (SZA of about 82^∘), at which point the retrieved column of CO₂ increases as the air mass increases, changing the shape of the air mass dependence of the CO₂ column. When the qSDV with line mixing is used, the retrieved column of CO₂ decreases as a function of air mass (up until the sun is above the horizon).

To correct for this, an empirical correction is applied to all TCCON XCO₂ (and XGas). The empirical correction determines the antisymmetrical component of the day's XCO₂, which is assumed to be the true variation of XCO₂ throughout the day, as well as the symmetrical component, which is caused by the air mass dependence of the retrieved column of the gas of interest and O₂. We can, therefore, represent a measurement as (Wunch et al., 2011)

$$\begin{array}{}\text{(5)}& y_i=\widehat{y}\left[\mathrm{1}+\mathit{\alpha }S\left({\mathit{\theta }}_i\right)+\mathit{\beta }A\left(t_i\right)\right],\end{array}$$

where $\widehat{y}$ is the mean value of XCO₂ measured that day, β is the fitted coefficient of the antisymmetric function A(t_i) and α is the fitted coefficient of the symmetric function S(θ_i). The antisymmetric function is calculated by (Wunch et al., 2011)

$$\begin{array}{}\text{(6)}& A\left(t_i\right)=\sin \left(\mathrm{2}\mathit{\pi }\left(t_i-t_{\mathrm{noon}}\right)\right),\end{array}$$

where t_i is the time of the measurement and t_noon is the time at solar noon, both in units of days. The symmetric function is calculated by (Wunch et al., 2011)

$$\begin{array}{}\text{(7)}& S\left({\mathit{\theta }}_i\right)={\left({\displaystyle \frac{{\mathit{\theta }}_i+\mathrm{13}{}^{\circ }}{\mathrm{90}{}^{\circ }+\mathrm{13}{}^{\circ }}}\right)}^{\mathrm{3}}-{\left({\displaystyle \frac{\mathrm{45}{}^{\circ }+\mathrm{13}{}^{\circ }}{\mathrm{90}{}^{\circ }+\mathrm{13}{}^{\circ }}}\right)}^{\mathrm{3}},\end{array}$$

where θ_i is the SZA in degrees. To determine α for the different line shapes, total columns of CO₂ were retrieved using the Voigt line shape (Wunch et al., 2015) and the qSDV with line mixing (Mendonca et al., 2016). Henceforth, we will refer to XCO₂ calculated from O₂ and CO₂ using the Voigt line shape as XCO₂ Voigt and the qSDV line shape as XCO₂ qSDV.

Figure 9 shows the average α calculated for each season at Darwin, Lamont and Park Falls. Eureka XCO₂ cannot be used to determine α because Eureka measurements do not go through the same range of SZAs as the other three sites due to its geolocation. The average α values derived from XCO₂ Voigt are represented by stars in Fig. 9, while the squares indicate XCO₂ qSDV. At all three sites, α is closer to 0 when the qSDV line shape is used in the retrieval compared to the Voigt retrieval, regardless of the season. The average α for XCO₂ Voigt calculated from a year of measurements from Darwin, Park Falls and Lamont is $-\mathrm{0.0071}\pm \mathrm{0.0057}$ and that for XCO₂ qSDV is $-\mathrm{0.0012}\pm \mathrm{0.0054}$.

For all four sites, $\mathit{\alpha }=-\mathrm{0.0071}$ is used to correct XCO₂ Voigt measurements. Figure 10a shows the XCO₂ Voigt anomalies plotted as a function of SZA. The data are expressed as the daily XCO₂ anomaly, which is the difference between the XCO₂ value and the daily median value, in order to remove the seasonal cycle. When XCO₂ is left uncorrected for air mass dependencies, XCO₂ decreases as a function of SZA up to approximately 82^∘ and increases as a function of SZA at angles greater than 82^∘. Figure 10b shows XCO₂ Voigt corrected for the air mass dependence. This air mass correction works well only up to a SZA of approximately 82^∘. Figure 10c is the same as Fig. 10a but for the uncorrected XCO₂ qSDV measurements, while Fig. 10d is the same as Fig. 10b but for the corrected XCO₂ qSDV measurements. When the air mass correction is applied to XCO₂ qSDV, there is a small difference between the corrected and uncorrected XCO₂ qSDV measurements, with the difference only noticeable for the Darwin measurements recorded at SZA >60^∘. For XCO₂ qSDV measurements taken at SZA >82^∘, XCO₂ does not increase with SZA as it does with the Voigt.

Figure 9 The average air-mass-dependent correction factor for XCO₂ derived from a year of spectra measured at Darwin, Lamont and Park Falls for different seasons. The dashed lines with stars are the α for XCO₂ retrieved using a Voigt line shape for both CO₂ and O₂ columns. The solid lines with squares are from XCO₂ retrieved using the qSDV for both CO₂ and O₂ columns.

Figure 10 (a) XCO₂ Voigt anomaly for a year of measurements from the four TCCON sites. The XCO₂ anomaly is the difference between each XCO₂ value and the daily median XCO₂. (b) The XCO₂ Voigt anomaly after the air mass dependence correction is applied to the XCO₂ Voigt data. (c) XCO₂ qSDV anomaly. (d) XCO₂ qSDV anomaly after correction for the air mass dependence.

Figure 11 (a) Correlation between TCCON and aircraft XCO₂ Voigt measurements for 13 TCCON sites. Each aircraft type is indicated by a different symbol given by the legend in the top-left corner. Each site is represented by a different colour given by the legend in the bottom right corner. The grey line indicates the one-to-one line and the dashed line is the line of best fit for the data. The slope of the line of best fit as well as the error on the slope are given in the plot. (b) The same as panel (a) but for XCO₂ qSDV.

Figure 12 (a–d) XCO₂ plotted as a function of day of the year for Eureka (2014), Park Falls (2013), Lamont (2010) and Darwin (2006). The mostly hidden red boxes are XCO₂ calculated from using a Voigt line shape in the retrieval and the blue boxes are from using the qSDV. (e–h) The difference between XCO₂ Voigt and XCO₂ qSDV.

5.2 Accuracy of XCO₂

To assess the accuracy of TCCON XCO₂ measurements, they are compared to aircraft XCO₂ profile measurements using the method described in Wunch et al. (2010). Figure 11a shows the comparison between the aircraft XCO₂ (Deutscher et al., 2010; Lin et al., 2006; Messerschmidt et al., 2010; Singh et al., 2006; Wofsy, 2011) measurements (legend at the top details the different aircraft) and TCCON XCO₂ Voigt measurements for 13 TCCON sites (given by the colour-coded legend at the bottom right). The grey line indicates the one-to-one line and the dashed line is the line of best fit. There is a bias of 0.9897±0.0005 given by the slope of the line of best fit in Fig. 11a for the XCO₂ Voigt measurements. Figure 11b is the same as Fig. 11a but for the XCO₂ qSDV measurements. The bias between the aircraft XCO₂ measurements and the XCO₂ qSDV measurements is 1.0041±0.0005 as given by the slope of the line of best fit in Fig. 11b. This increase in the slope can be explained by an increase in the retrieved column of CO₂ when using the qSDV with line mixing as shown in Mendonca et al. (2016) as well as combined with a decrease in the retrieved O₂ column due to using the qSDV. As discussed previously (Sect. 5) the decrease in the retrieved O₂ column is an improvement but the expected column of O₂ is still approximately 1.2 % higher (at the smallest SZA) than it should be. Therefore, the retrieved column of CO₂ is higher than it should be, and the slope would be greater if the retrieved column of O₂ was 1.2 % lower. Nevertheless, using the qSDV to retrieve total columns of CO₂ and O₂ reduces the difference between TCCON XCO₂ and aircraft XCO₂ measurements by 0.62 %.

TCCON XCO₂ measurements are divided by the scale factors (or bias determined in Fig. 11) to calibrate to the WMO scale. For all TCCON XCO₂ measurements retrieved with a Voigt line shape, the air mass correction is first applied to the data and the result is divided by the determined bias factor, 0.9897. Figure 12a to d shows XCO₂ Voigt (for Eureka, Park Falls, Lamont and Darwin) indicated by red square boxes in the plots. XCO₂ Voigt measurements taken at SZA >82^∘ have been filtered out because they cannot be corrected for the air mass dependence. The blue boxes are XCO₂ qSDV corrected for air mass dependence and scaled by 1.0041. No filter was applied to the XCO₂ qSDV measurements for SZA since the air mass dependence correction works at all SZA. Figure 12e–h shows the difference between XCO₂ Voigt and XCO₂ qSDV for Eureka, Park Falls, Lamont and Darwin. The mean differences for the data shown in Fig. 12e to h are 0.113±0.082, $-\mathrm{0.102}\pm \mathrm{0.223}$, $-\mathrm{0.132}\pm \mathrm{0.241}$ and $-\mathrm{0.059}\pm \mathrm{0.231}$ µmol mol⁻¹ (ppm) for Eureka, Park Falls, Lamont and Darwin respectively. The difference throughout the day at Park Falls, Lamont and Darwin varies between −0.6 and 0.2 µmol mol⁻¹ and is SZA dependent.

Figure 13 (a) XCO₂ from Park Falls retrieved from spectra recorded on 18 June 2013. Retrieved XCO₂ is plotted (1) with a Voigt line shape and corrected for the air mass dependence (red squares), (2) with the qSDV (cyan circles) and (3) with the qSDV and corrected for the air mass dependence (blue squares). (b) The difference between the corrected Voigt XCO₂ and the qSDV XCO₂ (cyan circles) and the difference between the Voigt corrected XCO₂ and the corrected qSDV XCO₂ (blue squares). The top x axis is the SZA that corresponds to the hour on the bottom x axis.

Figure 13a shows XCO₂ Voigt corrected for the air mass dependence, as well as XCO₂ qSDV, uncorrected and corrected for the air mass dependence. These XCO₂ measurements were retrieved from Park Falls spectra recorded on 18 June 2013. For all three XCO₂ measurements, the amount of XCO₂ decreases throughout the day. Figure 13b shows the difference between the corrected Voigt XCO₂ and the uncorrected qSDV XCO₂, as well as the difference between the corrected Voigt XCO₂ and the corrected qSDV XCO₂. The difference between the Voigt and the qSDV (corrected and uncorrected) shows that at the start and end of the day, more XCO₂ is retrieved with the qSDV, while at midday less is retrieved with the qSDV. The range in the differences seen in Fig. 12e to h varies with SZA throughout the day as shown in Fig. 13b.

6 Discussion and conclusions

Using cavity ring-down spectra measured in the lab, we have shown that the Voigt line shape is insufficient to model the line shape of O₂ for the 1.27 µm band, consistent with the results of Hartmann et al. (2013) and Lamouroux et al. (2014). By using the speed-dependent Voigt line shape when calculating the absorption coefficients, we were better able to reproduce the measured absorption coefficients than using the Voigt line shape. However, some residual structure remains as seen Figs. 1 and 2. This is partly due to the blending of spectral lines (i.e. line mixing) and the inability to retrieve the spectroscopic parameters for weak O₂ transitions. Fitting low-pressure spectra would help with isolating spectral lines and decreasing the uncertainty on the retrieved spectroscopic parameters for the Q branch lines.

Accurate measurements of the pressure shifts in the 1.27 µm band have been hard to obtain as shown in Newman et al. (1999) and Hill et al. (2003). While the retrieved pressure shifts show a dependence on quantum number m (Fig. 3c) as one would expect, this dependence is not as strong as the m dependence of the Lorentz widths (Fig. 3b). This can be explained by the fact that line mixing, which is shown to be important for the O₂ A-band, was not considered when fitting the cavity ring-down spectra. Neglecting line mixing usually produces an asymmetric residual in the discrete lines as well as a broad residual feature associated with the fact that collisions are transferring intensity from one part of the spectrum to another. By fitting a set of Legendre polynomials for CIA, we could be simultaneously fitting the broader-band feature associated with line mixing, while the retrieved pressure shifts and speed-dependent pressure shifts could be compensating for the asymmetric structure one would see in the discrete lines when neglecting line mixing. The remaining structure, as seen in Fig. 1c, could be due to neglecting line mixing, especially in the Q branch, where the spacing between spectral lines is small (in comparison to the P and R branches) and line mixing is most likely prevalent. The large error bars for the measured pressure shifts and speed-dependent pressure shifts as well as a deviation from a smooth m dependence of these parameters could be due to neglecting line mixing when fitting the lab spectra. Figure 3c and d show that the spectral lines that have large error bars and deviate from an expected m dependence belong mainly to the Q-branch spectral lines (which are mostly likely impacted by line mixing). To achieve the results obtained in this study it is best to use the parameters as they are instead of trying to apply an interpolation that depends on m or even omitting them, unless one tests these changes on atmospheric spectra that cover different range of conditions (i.e. seasons, dry/wet, SZA, geographical locations). It is evident that the parameters might be compensating for effects (such as line mixing) that were not included when fitting the lab spectra, and changing these parameters (or omitting them) could lead to degradation in the quality of the spectral fits of solar spectra and a change in the air mass dependence of the retrieved column of O₂, which would impact the air mass dependence of XCO₂.

The pressure dependence of the retrieved speed-dependent width parameter is an indication that Dicke narrowing needs to be taken into account, as shown by Bui et al. (2014) for CO₂. When both speed dependence and Dicke narrowing are present, a multispectrum fit needs to be used due to the correlation between the parameters (Bui et al., 2014). Domysławska et al. (2016) recommend using the qSDV to model the line shape of O₂ based on multiple line shape studies of the O₂ B-band. In these studies, a multispectrum fit to low-pressure (0.27–5.87 kPa) cavity ring-down spectra was performed testing multiple line shapes that took speed dependence and Dicke narrowing into account both separately and simultaneously. They found that the line shapes that only used Dicke narrowing were not good enough to model the line shape of the O₂ B-band lines, but a line shape that included either speed dependence or both speed dependence and Dicke narrowing produced similar quality fits, ultimately concluding that speed dependence has a larger effect than Dicke narrowing. It was noted in the study by Wójtewicz et al. (2014) that both Dicke narrowing and speed-dependent effects might simultaneously play an important role in modelling the line shape of the O₂ B-band lines. However, the speed-dependent and Dicke narrowing parameters are highly correlated at low pressures. Reducing the correlation requires either a multispectrum fit of spectra at low pressures with a high enough signal-to-noise ratio or spectra that cover a wide range of pressures (Wójtewicz et al., 2014). So, by combining the high-pressure spectra used in this study with low-pressure spectra in a multispectral fit, both the speed dependence and Dicke narrowing parameters could be retrieved. The temperature dependence of the Lorentz width coefficients of this band has never been measured before, which could have an impact on the air mass dependence of O₂. Combining high-pressure cavity ring-down absorption coefficient measurements with those for low pressures and different temperatures as done in Devi et al. (2015, 2016) for CH₄ would lead to more accurate line shape parameters for O₂.

By taking speed dependence into account for both CO₂ (in the work of Mendonca et al., 2016) and O₂ (the work presented here), we were able to significantly decrease the air mass dependence of TCCON XCO₂ and the bias between TCCON and aircraft XCO₂. XAIR calculated with the column of O₂ retrieved with the qSDV is now closer to the expected value of 1, but XAIR still has an air mass dependence which is the result of the retrieved total column of O₂ decreasing as a function of SZA at large SZA. This remaining air mass dependence could be due to neglecting effects such as Dicke narrowing and line mixing in the absorption coefficient calculations, as well as assuming a perfect instrument line shape in the retrieval algorithm. However, retrieving O₂ with the qSDV significantly decreases the air mass dependence of XCO₂. With the qSDV line shape, XCO₂ measurements taken at SZA >82^∘ no longer have to be discarded. We recommend using the full range of SZA, which would result in more XCO₂ measurement available from all TCCON sites. This is particularly important for high-latitude TCCON sites such as Eureka because measurements taken from late February to late March and from late September to mid-October are taken at SZA >82^∘. Filtering out these large SZA measurements thus limits the knowledge of the seasonal cycle of XCO₂ at high latitudes. The air mass dependence of the O₂ column not only affects XCO₂ but all trace gases measured by TCCON and in the future the air mass dependence of all XGas will be determined with these new O₂ columns.