An optimal estimation-based retrieval of upper atmospheric oxygen airglow and temperature from SCIAMACHY limb observations

An optimal estimation-based algorithm is developed to retrieve the number density of excited oxygen (O₂) molecules that generate airglow emissions near 0.76 µm ($b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$ or A band) and 1.27 µm (a¹Δ_g or ¹Δ band) in the upper atmosphere. Both oxygen bands are important for the remote sensing of greenhouse gases. The algorithm is applied to the limb spectra observed by the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) instrument in both the nominal (tangent heights below ∼ 90 km) and mesosphere–lower thermosphere (MLT) modes (tangent heights spanning 50–150 km). The number densities of emitting O₂ in the a¹Δ_g band are retrieved in an altitude range of 25–100 km near-daily in 2010, providing a climatology of O₂ a¹Δ_g-band airglow emission. This climatology will help disentangle the airglow from backscattered light in nadir remote sensing of the a¹Δ_g band. The global monthly distributions of the vertical column density of emitting O₂ in a¹Δ_g state show mainly latitudinal dependence without other discernible geographical patterns. Temperature profiles are retrieved simultaneously from the spectral shapes of the a¹Δ_g-band airglow emission in the nominal limb mode (valid altitude range of 40–100 km) and from both a¹Δ_g- and $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$-band airglow emissions in the MLT mode (valid range of 60–105 km). The temperature retrievals from both airglow bands are consistent internally and in agreement with independent observations from the Atmospheric Chemistry Experiment Fourier transform spectrometer (ACE-FTS) and the Michelson Interferometer for Passive Atmospheric Sounding (MIPAS), with the absolute mean bias near or below 5 K and root mean squared error (RMSE) near or below 10 K. The retrieved emitting O₂ number density and temperature provide a unique dataset for the remote sensing of greenhouse gases and constraining the chemical and physical processes in the upper atmosphere.

1 Introduction

This study of upper atmospheric airglow from oxygen (O₂) is driven by the need to measure O₂ simultaneously with methane (CH₄) and carbon dioxide (CO₂) in satellite remote sensing. Methane and CO₂ are two of the most important anthropogenic greenhouse gases, but the spatiotemporal variations in their sources and sinks are poorly understood, leading to significant uncertainties in projections of future climate trends (Miller et al., 2007; Turner et al., 2019; Friedlingstein et al., 2020; Saunois et al., 2020). Methane, in particular, has a global warming potential of about 56–105 times higher than that of CO₂ for a 20-year time period (Howarth, 2014). Reducing methane emissions is among the most impactful actions that can be taken to reduce the global warming, which requires the identification of emission sources with a high degree of accuracy. Spaceborne observations offer a powerful tool in quantifying the spatiotemporal distributions of methane and CO₂ and inferring their sources and sinks due to the extensive spatial and temporal coverage of satellites (Eldering et al., 2017; Turner et al., 2018; Lorente et al., 2021). To separate the sources and sinks of methane and CO₂ from variations in surface pressure and specific humidity, the abundances of CO₂ and methane observed by satellites are usually represented by the column-averaged dry mole fractions ($X_{{\mathrm{CO}}_{\mathrm{2}}}$ and $X_{{\mathrm{CH}}_{\mathrm{4}}}$). Simultaneous observation of O₂ is often necessary for deriving $X_{{\mathrm{CO}}_{\mathrm{2}}}$ and $X_{{\mathrm{CH}}_{\mathrm{4}}}$ and accounting for contamination of aerosol and clouds (Butz et al., 2011).

The two widely used O₂ absorption bands in greenhouse gas remote sensing are the $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ band near 0.76 µm (O₂ A band hereafter) and the $a^{\mathrm{1}}{\mathrm{\Delta }}_{\mathrm{g}}\leftarrow X^{\mathrm{3}}{\mathrm{\Sigma }}_{\mathrm{g}}^{-}$ band near 1.27 µm (O₂ ¹Δ band hereafter). The O₂ A band is commonly used in existing and planned spaceborne CO₂ and methane sensors (Bovensmann et al., 1999; Crisp et al., 2017; Veefkind et al., 2012; Moore et al., 2018; Buchwitz et al., 2013). The O₂ absorption feature in the A band is much stronger than the shortwave infrared CO₂ and methane bands and separated by significant spectral gaps, which challenges the spectral fitting and radiative transfer modeling. The O₂ A band also overlaps with strong terrestrial solar-induced fluorescence that provides valuable information on plant photosynthesis but perturbs the O₂ absorption features and, if not properly accounted for, leads to systematically biased greenhouse gas retrievals (Frankenberg et al., 2012).

The O₂ ¹Δ band plays an instrumental role in ground-based greenhouse gas remote sensing (Wunch et al., 2011; Fu et al., 2014; Frey et al., 2019) and is adopted by the MicroCarb (Bertaux et al., 2020) and MethaneSAT (Staebell et al., 2021) satellite missions. Both O₂ A and ¹Δ bands feature upper-atmospheric airglow emissions (Zarboo et al., 2018) that overlap with the backscattered solar light containing O₂ absorption signals. The A-band airglow, which peaks in the mesopause region, has previously been neglected in nadir remote sensing due to its low intensity (Sioris, 2003), but there is a lack of quantitative assessment of how the A-band airglow impacts greenhouse gas retrieval. The ¹Δ-band airglow emitted from the upper stratosphere and mesosphere is much stronger due to photochemically generated O₂ in a¹Δ_g state from ozone photolysis and was the main reason for choosing the A band over the ¹Δ band by spaceborne sensors (Kuang et al., 2002). However, recent advances in modeling the O₂ ¹Δ band airglow spectra enable the disentangling of the airglow from backscattered light (Sun et al., 2018a; Bertaux et al., 2020) and open up the opportunity for spaceborne nadir remote sensing of the O₂ ¹Δ band.

In order to account for the airglow contributions in the nadir-observed O₂ spectra, it is crucial to have an accurate understanding of the spatial, temporal, and spectral distribution of airglow emissions. Nevertheless, observation-based studies of O₂ A- and ¹Δ-band airglow are rare and generally lack the information needed for nadir spaceborne remote sensing of greenhouse gases. Zarboo et al. (2018) retrieved volume emission rates (VERs) for the O₂ A- and ¹Δ-band airglow over the vertical range of 50–150 km using the mesosphere and lower thermosphere (MLT) limb observation mode of the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (SCIAMACHY) instrument. These MLT observations could not capture a significant portion of the O₂ ¹Δ-band airglow, which peaks near the stratopause, and the linear inversion applied by Zarboo et al. (2018) is subject to systematic biases in the retrieved airglow VERs below 90 km for the O₂ A band and below 60 km for the O₂ ¹Δ band due to O₂ self-absorption (Sun et al., 2018a). Li et al. (2020) conducted high vertical and along-track resolution retrieval of the O₂ ¹Δ-band airglow for the purpose of retrieving mesospheric ozone using the Optical Spectrograph and InfraRed Imaging System (OSIRIS) limb radiance, which could not reveal the necessary horizontal distributions and spectral shape of the airglow. Sun et al. (2018a) and Bertaux et al. (2020) fit SCIAMACHY limb spectra in the O₂ ¹Δ band with onion-peeling algorithms and retrieved airglow from selected orbits, although the spectral variation in airglow emission lines due to temperature was not fully incorporated.

The distribution of temperature in the atmospheric regions where airglow is emitted is important for reconstructing the emission spectra. The spectral band shape has been used to retrieve upper atmospheric temperature using CO₂ (Boone et al., 2005; Marshall et al., 2011; García-Comas et al., 2014), O₂ A-band absorption (Nowlan et al., 2007), and the O₂ A-band airglow (Sheese et al., 2010; Yang et al., 2021). Using the O₂ ¹Δ-band airglow spectra, Sun et al. (2018a) was able to retrieve mesospheric temperature above 60 km, limited by excessive uncertainties below it. The mass spectrometer incoherent scatter (MSIS) model provides temperature estimates throughout the atmosphere and is often used when observations are unavailable (Picone et al., 2002).

To this end, we develop an optimal estimation-based algorithm to retrieve the airglow emissions and upper atmospheric temperature from limb-viewed radiance spectra. Compared with the onion-peeling algorithm in Sun et al. (2018a) that fits limb spectra from high to low tangent heights progressively, this algorithm combines spectra over a range of tangent heights and simultaneously retrieves the vertical profiles of local emissions and temperature in a consistent manner. The use of Bayesian inversion enables the incorporation of a priori knowledge, balancing of measurement error and prior error, and detailed posterior error analysis, including the averaging kernel matrix and degrees of freedom for signal (DOFS; Rodgers, 2000; Brasseur and Jacob, 2017). The airglow emission spectra are simulated based on the spectral model for the O₂ ¹Δ band proposed by Sun et al. (2018a), which we demonstrate can be extended to the O₂ A band with a simple generalization. We compare the retrieved temperatures in the upper stratosphere, mesosphere, and lower thermosphere with independent measurements. The algorithm is applied to 1 year of SCIAMACHY limb observations, including the MLT mode, to construct a climatology of O₂ airglow and upper atmospheric temperature at 10:00 local solar time (LST).

2 Data

The SCIAMACHY limb spectra and MSIS model outputs are necessary for the retrieval algorithm. The retrieved temperature profiles from SCIAMACHY are compared with observations from Atmospheric Chemistry Experiment Fourier transform spectroscopy (ACE-FTS) and Michelson Interferometer for Passive Atmospheric Sounding (MIPAS) instruments.

2.1 SCIAMACHY

The SCIAMACHY instrument is an eight-channel grating spectrometer that measures radiation that is backscattered, reflected, transmitted, or emitted by the Earth's atmosphere and surface in limb and nadir geometries from 240 to 2380 nm. The instrument was launched on board the Envisat satellite, which was operational on a sun-synchronous orbit with an Equator-crossing time in the descending node of 10:00 LST from March 2002 until April 2012. In this study, only limb scattering measurements in spectral channels 4 (597–789 nm) and 6 (990–1750 nm) were used. In the nominal limb mode, SCIAMACHY observed the atmosphere from the surface up to 93 km in 2010. It also observed the MLT region, covering 50–150 km altitudes in the limb view over 2 d in each month, from July 2008 until April 2012. We conduct spectral and radiometric calibrations of SCIAMACHY spectra using version 3.2.6 of the SciaL1C command line tool software package. We then use a customized algorithm to convert the Level 1 data to the NetCDF4 file format. In contrast to previous studies (Bender et al., 2017; Sun et al., 2018a; Zarboo et al., 2018; Bertaux et al., 2020), where SCIAMACHY spectra were averaged across the track, we keep the eight across-track positions separate to enhance the horizontal space resolution. Limb observations over different altitudes at the same across-track position (usually there are 15 observing altitudes per across-track position for both nominal and MLT modes) are grouped as a single vertical sounding in order to retrieve atmospheric profiles. The latitude and longitude of the tangent point near 50 km altitude is used to represent the location of the entire vertical sounding, which gives a location ambiguity of less than ∼ 20 km. The average vertical distance between adjacent tangent heights within the same across-track position is about 6.6 km for both nominal and MLT modes. We have noticed that the eight vertical soundings across the track can always be grouped into four pairs. The soundings within each pair are close in latitude and longitude but zigzag in tangent heights. Therefore, the eight vertical soundings can be combined into four vertical soundings with a doubled vertical sampling of 3.3 km. We retrieve profiles at eight native across-track positions and only pair the retrieved profiles into four when calculating the column number density of emitting O₂.

To cover the A band, radiance spectra in the 759 to 772 nm range are obtained from channel 4. In addition to airglow emissions, the observed radiance spectra contain photons from Rayleigh scattering and multiple scattering by the atmosphere and the surface. The scattering signal can be approximated by limb views at maximum tangent heights and contains the O₂ absorption feature for the A band. No O₂ absorption is observed for the ¹Δ band due to much weaker atmospheric scattering and lower O₂ absorption. To account for the scattered light in the A band, a background signal consisting of averaged A-band spectra at high tangent heights (130–150 km) from the same sounding is subtracted from each limb spectrum (Zarboo et al., 2018). Before each subtraction, the background signal is scaled to match its out-of-band radiance with the out-of-band radiance of the limb spectrum to be corrected for. This step assures the out-of-band radiance centers at zero after correction. This correction assumes the spectral shape of scattered light is the same for all tangent heights and may leads to systematic errors at low tangent heights from the A band, where airglow emission is low, and the scattering path may differ significantly from the thermosphere. Following Zarboo et al. (2018), we consider 750–759 and 767–780 nm to be out of band. The O₂ ¹Δ-band spectra are taken from 1240 to 1300 nm in channel 6 with two bad pixels being removed. The multiple scattering is considered negligible in the O₂ ¹Δ band, and only a linear background fitted using the out-of-band radiance at 1210–1240 and 1300–1340 nm is subtracted from each limb spectrum.

2.2 MSIS model

The temperature, pressure, and the ground state O₂ number densities are sampled from the NRLMSISE-00 model, a MSIS model developed at the U.S. Naval Research Laboratory (Picone et al., 2002) using a Python package by Hirsch and Kastinen (2021). The MSIS temperature profile is used as a priori values in the optimal estimation. The pressure profile is fixed in the algorithm, whereas the relative changes in O₂ number densities from MSIS values are retrieved, which affects the O₂ self-absorption.

2.3 ACE-FTS

The ACE-FTS instrument was launched aboard the Canadian SCISAT-1 satellite in August 2003 and is still active. It operates from 750 to 4400 cm⁻¹ (2.2 to 13.3 µm) at a high spectral resolution (0.02 cm⁻¹). A sequence of atmospheric transmission spectra in the limb geometry is recorded using the Sun as an infrared source during sunrise and sunset (i.e., solar occultation). The temperature and pressure retrievals from the ACE-FTS were performed at 15–125 km by fitting a set of spectral windows containing CO₂ lines with the inverse of temperature as a free parameter (Boone et al., 2020). The temperature retrieval process was done by dividing the atmosphere into two altitude regimes with a crossover at about 50 km. The temperature retrieval was performed assuming a fixed CO₂ mixing ratio up to about 50 km, while above it an empirical function was used to describe CO₂ mixing ratio in order to force the results to exhibit smooth behavior. We use the temperature profiles from version 4.1 of ACE-FTS Level 2 data for the comparison of temperature retrieved from SCIAMACHY airglow. The ACE-FTS temperature profiles have been extensively used in the validation of profiles from other instruments such as SOFIE (Solar Occultation for Ice Experiment; Marshall et al., 2011), MIPAS (García-Comas et al., 2014), and OSIRIS (Sheese et al., 2012).

2.4 MIPAS

The MIPAS instrument was a Fourier transform spectrometer for the detection of limb emission spectra in the middle and upper atmosphere on board the Envisat satellite since March 2002. It observed a wide spectral interval throughout the mid-infrared with a high spectral resolution of 0.025 cm⁻¹. Operating in a wavelength range from 4.15 to 14.6 µm, MIPAS detected and spectrally resolved a large number of emission features of atmospheric minor constituents playing a major role in atmospheric chemistry. The retrieval of temperature is done from measurements of the CO₂ atmospheric radiance at 15 µm for each MIPAS single limb scan. Because MIPAS was on the same platform as SCIAMACHY, there were abundant colocated observations, enabling monthly and zonally resolved temperature intercomparison. The measurement modes of MIPAS used in this study include the nominal measurement mode with an altitude coverage of roughly 6–70 km, the middle atmosphere (MA) mode covering 18–102 km, the upper atmospheric (UA) mode covering 42–172 km, and the noctilucent cloud (NLC) mode covering 39–102 km. The nominal measurement mode makes up the bulk of MIPAS measurements, whereas the MA and UA modes were available at least every 10 d, and the NLC mode only happened on a few days in 2010. We use the nominal temperature profiles from version 8 of MIPAS Level 2 data retrieved by the European Space Agency (ESA; Dinelli et al., 2021). Version 8 data from the other modes are obtained through the Institute of Meteorology and Climate Research in cooperation with the Instituto de Astrofísica de Andalucía (IMK/IAA) retrieval algorithm (García-Comas et al., 2012; Kiefer et al., 2021). The typical total errors are 0.5–2 K below 70 km and 2–7 K above (for MA, UA, and NLC modes). The typical vertical resolutions in the comparison range of this study are 3–7 km.

3 Methodology

We assume the atmosphere as homogeneous layers and calculate airglow emission and O₂ self-absorption using HITRAN (high-resolution transmission molecular absorption database; Rothman, 2021). The local emission and absorption from these layers are calculated and integrated along the line-of-sight to simulate the limb-viewed radiance. A Bayesian inversion is applied to retrieve airglow emission and temperature profiles by minimizing the difference between simulated and observed limb radiance spectra with a priori regularization.

3.1 Airglow emission from a single layer

An atmospheric layer bounded by two tangent heights of SCIAMACHY limb observations is the basic spatial resolving unit of this study. We create an additional layer above the outermost tangent height by assuming a layer thickness equal to the average difference between adjacent tangent heights. The atmospheric properties are assumed to be uniform within a layer, and the layer height is represented by its middle altitude when sampling meteorological parameters from the MSIS model and comparing with other observations. The excited O₂ molecules that give rise to A and ¹Δ band airglow are denoted as O${}_{\mathrm{2}}\left(b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}\right)$ and O₂(a¹Δ_g), respectively, and generalized as O${}_{\mathrm{2}}^{\ast }$. Their number densities ($\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$) within the layer are proportional to the VERs via the band-integrated Einstein A coefficients a as follows:

$$\begin{array}{}\text{(1)}& \mathrm{VER}=\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]a.\end{array}$$

Here VER is measured by photons per cubic centimeter per second (hereafter photons cm⁻³ s⁻¹), $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ is measured by molecules per cubic centimeter (hereafter molec. cm⁻³), and Einstein A coefficient is $\mathrm{2.27}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ for the ¹Δ band (Sun et al., 2018a) and 0.08693 s⁻¹ for the A band (Long et al., 2010). It is worth pointing out that $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ is the number density that directly contributes to the airglow at the wavelengths of the O₂ A band or ¹Δ band and does not necessarily include all excited O₂ in the $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$ or a¹Δ_g state. For example, the emission channel $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$–a¹Δ_g, which corresponds to the so-called Noxon band at 1.91 µm (Noxon, 1961), is ignored. Nevertheless, this does not change the derived parameters.

The absorption and emission spectra in a layer can be resolved by either wavelength (λ) or wavenumber ($\mathit{\nu }=\mathrm{1}/\mathit{\lambda }$). We will use the subscript λ to denote a parameter resolved by wavelength and the subscript ν to denote a parameter resolved by wavenumber. We calculate the monochromatic absorption cross section of O₂ using the Python library HITRAN Application Programming Interface (HAPI; Kochanov et al., 2016) at a grid space of 0.0002 nm for the O₂ A band and 0.001 nm for the ¹Δ band, assuming the Voigt line shape profile. Note the monochromatic absorption cross section in centimeters squared per molecule (hereafter cm² molec.⁻¹) is the same at wavelength and wavenumber coordinates, i.e., σ_λ=σ_ν, so we may simplify it as σ. We have tested both HITRAN2016 (Gordon et al., 2017) and HITRAN2020 (Gordon et al., 2022) spectroscopic databases and provide results using HITRAN2020 unless otherwise noted. The Jacobian of σ with respect to temperature T, $\partial \mathit{\sigma }/\partial T$, is calculated by finite difference of 0.01 K using HAPI.

The monochromatic VER (i.e., emissivity) for airglow in wavenumber space is denoted as ε_ν with a unit of photons per cubic centimeter per second per centimeter inverse (hereafter photons cm⁻³ s⁻¹ (cm⁻¹)⁻¹). Integrating ε_ν across the band should recover the VER, as follows:

$$\begin{array}{}\text{(2)}& \mathrm{VER}=\underset{\mathit{\nu }}{\int }{\mathit{\epsilon }}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }.\end{array}$$

According to the spectroscopic airglow model described in Sun et al. (2018a) and Bertaux et al. (2020), the O₂ airglow emission spectrum ε_ν is related to its absorption spectrum σ by the following:

$$\begin{array}{}\text{(3)}& {\mathit{\epsilon }}_{\mathit{\nu }}={\displaystyle \frac{c_{\mathrm{0}}\mathit{\sigma }{\mathit{\nu }}^{\mathrm{2}}}{\exp (c_{\mathrm{2}}\mathit{\nu }/T)-\mathrm{1}}},\end{array}$$

where c₀ is a scaling constant, and c₂ is the second radiation constant in Planck's law with a value of 1.4387769 cm K. Combining Eqs. (1), (2), and (3), we can solve for c₀ as follows:

$$\begin{array}{}\text{(4)}& c_{\mathrm{0}}={\displaystyle \frac{\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]a}{{\int }_{\mathit{\nu }}\frac{\mathit{\sigma }{\mathit{\nu }}^{\mathrm{2}}}{\exp (c_{\mathrm{2}}\mathit{\nu }/T)-\mathrm{1}}\mathrm{d}\mathit{\nu }}}.\end{array}$$

Plugging in Eq. (4) into Eq. (3), we obtain the following:

$$\begin{array}{}\text{(5)}& {\mathit{\epsilon }}_{\mathit{\nu }}=\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]a{\displaystyle \frac{t_{\mathit{\nu }}}{{\int }_{\mathit{\nu }}{t}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }}},\end{array}$$

where

$$\begin{array}{}\text{(6)}& t_{\mathit{\nu }}={\displaystyle \frac{\mathit{\sigma }{\mathit{\nu }}^{\mathrm{2}}}{\exp (c_{\mathrm{2}}\mathit{\nu }/T)-\mathrm{1}}}\end{array}$$

is an intermediate function to facilitate Jacobian derivation in the following step.

The SCIAMACHY retrieval works in wavelength space, so it is more convenient to define airglow emissivity as a function of wavelength (i.e., using ε_λ instead of ε_ν). The unit of ε_λ is photons per cubic centimeter per second per nanometer (hereafter photons cm⁻³ s⁻¹ nm⁻¹). Integrating ε_λ across the band in wavelength space should equally give the VER, as follows:

$$\begin{array}{}\text{(7)}& \mathrm{VER}=\underset{\mathit{\lambda }}{\int }{\mathit{\epsilon }}_{\mathit{\lambda }}\mathrm{d}\mathit{\lambda }.\end{array}$$

ε_λ and ε_ν are, hence, related via the following:

$$\begin{array}{}\text{(8)}& {\mathit{\epsilon }}_{\mathit{\lambda }}={\mathit{\epsilon }}_{\mathit{\nu }}{\displaystyle \frac{\mathit{\nu }}{\mathit{\lambda }}},\end{array}$$

and similar relationships apply for the Jacobians of ε_λ and ε_ν. The Jacobian of airglow emissivity with respect to $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ is as follows:

$$\begin{array}{}\text{(9)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\lambda }}}{\partial \left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]}}& ={\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\nu }}}{\partial \left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]}}{\displaystyle \frac{\mathit{\nu }}{\mathit{\lambda }}}\\ & ={\displaystyle \frac{a\mathit{\nu }{t}_{\mathit{\nu }}}{\mathit{\lambda }{\int }_{\mathit{\nu }}{t}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }}}.\end{array}\end{array}$$

The Jacobian of ε_λ with respect to temperature is more complicated and can be derived by differentiating Eq. (5) as follows:

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\lambda }}}{\partial T}}& ={\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\nu }}}{\partial T}}{\displaystyle \frac{\mathit{\nu }}{\mathit{\lambda }}}\\ & =\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]a{\displaystyle \frac{\frac{\partial t_{\mathit{\nu }}}{\partial T}{\int }_{\mathit{\nu }}{t}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }-t_{\mathit{\nu }}{\int }_{\mathit{\nu }}\frac{\partial t_{\mathit{\nu }}}{\partial T}\mathrm{d}\mathit{\nu }}{{\left({\int }_{\mathit{\nu }}{t}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }\right)}^{\mathrm{2}}}}{\displaystyle \frac{\mathit{\nu }}{\mathit{\lambda }}}.\end{array}\end{array}$$

Here, we leverage the following relationship:

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial {\int }_{\mathit{\nu }}{t}_{\mathit{\nu }}\mathrm{d}\mathit{\nu }}{\partial T}}=\underset{\mathit{\nu }}{\int }{\displaystyle \frac{\partial t_{\mathit{\nu }}}{\partial T}}\mathrm{d}\mathit{\nu }.\end{array}$$

Equation (10) is completed by plugging in the following:

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial t_{\mathit{\nu }}}{\partial T}}={\displaystyle \frac{\partial \mathit{\sigma }}{\partial T}}{\displaystyle \frac{{\mathit{\nu }}^{\mathrm{2}}}{\exp (c_{\mathrm{2}}\mathit{\nu }/T)-\mathrm{1}}}+\mathit{\sigma }{\displaystyle \frac{c_{\mathrm{2}}{\mathit{\nu }}^{\mathrm{3}}\exp (c_{\mathrm{2}}\mathit{\nu }/T)}{{\left(T\exp (c_{\mathrm{2}}\mathit{\nu }/T)-\mathrm{1}\right)}^{\mathrm{2}}}},\end{array}$$

which is the derivative of Eq. (6) with respect to T.

3.2 Forward model for limb-viewed airglow spectra

Assuming that we include N limb observations in a vertical sounding, limited by an altitude range, the number of atmospheric layers will also be N. We use index i to denote each limb view in a collection of limb observations, and hence, i ranges from 1 to N. Figure 1 illustrates the line of sight at tangent layer i=2 when there are N=7 layers (and equivalently, tangent heights) under consideration. We use index j to denote layers penetrated by the line of sight of the limb view i, and j ranges from i to N. The line of sight slices through those layers twice under the assumption of homogeneous layers. To uniquely identify each segment of the line of sight, another index l is introduced to denote those segments, from near to far, as shown in Fig. 1. Namely, a segment is a slant path of the line of sight through a layer. The segment index l ranges from 1 to 2N−i. A layer j can be uniquely identified for any segment l, as follows:

$$\begin{array}{}\text{(13)}& j=\left|l-N+{\displaystyle \frac{i-\mathrm{1}}{\mathrm{2}}}\right|+{\displaystyle \frac{i+\mathrm{1}}{\mathrm{2}}}.\end{array}$$

However, a given layer j will correspond to two segments (i.e., each layer is penetrated twice), as follows:

$$\begin{array}{}\text{(14)}& l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}.\end{array}$$

Figure 1 Illustration of the airglow limb-viewing line of sight. The observer is located at the left end of the horizontal line. The layer thicknesses are greatly exaggerated relative to the Earth's radius for visualization. The dashed vertical line indicates the tangent point and the thickness of each layer.

We denote the native resolution radiance spectrum at limb view i as r_λ,i, where the subscript λ implies it is resolved in wavelength space. Its unit is photons per centimeter squared per second per nanometer per steradian (hereafter photons cm⁻² s⁻¹ nm⁻¹ sr⁻¹). To model this radiance, we need to sum emissions of all segments along the line of sight while considering absorption by downstream ground state O₂ molecules, as follows:

$$\begin{array}{}\text{(15)}& r_{\mathit{\lambda },i}=\sum _{l=\mathrm{1}}^{\mathrm{2}N-i}\left({\displaystyle \frac{L_{ij}{\mathit{\epsilon }}_{\mathit{\lambda },j}}{\mathrm{4}\mathit{\pi }}}\exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}-\sum _{l^{\prime }=\mathrm{1}}^{l-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}\right)\right).\end{array}$$

Here the outermost summation loops from the nearest (l=1) to the farthest ($l=\mathrm{2}N-i$) segment, and each segment l maps uniquely to a layer j through Eq. (13). L_ij is the path length of the tangent view i through layer j in units of centimeters. $L_{ij}{\mathit{\epsilon }}_{\mathit{\lambda },j}/\left(\mathrm{4}\mathit{\pi }\right)$ is the unattenuated radiance emitted from such a segment. If there were no absorption, then the exponential term in Eq. (15) would be absent, and the observed limb radiance is simply the sum of radiances of all segments.

The exponential term in Eq. (15) represents the O₂ self-absorption between the emitting segment l and the observer. The summation within the exponential term gives the total absorption optical depth of the downstream segments ranging from $l^{\prime }=l-\mathrm{1}$, the segment immediately downstream, to $l^{\prime }=\mathrm{1}$, the segment closest to the observer. Each downstream segment l^′ again maps uniquely to a layer j^′, whose optical depth is given by the following:

$$\begin{array}{}\text{(16)}& {\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}=[{\mathrm{O}}_{\mathrm{2}}{]}_{j^{\prime }}{\mathit{\sigma }}_{\mathit{\lambda },j^{\prime }}{L}_{i{j}^{\prime }}.\end{array}$$

Here $[{\mathrm{O}}_{\mathrm{2}}{]}_{j^{\prime }}$ is the O₂ number density (in ground state, not to be confused with the emitting O₂ or O${}_{\mathrm{2}}^{\ast }$) in layer j^′, and ${\mathit{\sigma }}_{\mathit{\lambda },j^{\prime }}$ is the monochromatic absorption cross section of O₂ at layer j^′ (the subscript λ is added back to emphasize that σ is resolved in wavelength space).

The ${\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}$ term in Eq. (15) represents the self-absorption happening within the emitting segment. It has been neglected in the onion-peeling algorithm in Sun et al. (2018a), leading to large errors in spectral fitting at lower tangent heights where self-absorption is significant. We treat it here as an effective optical depth and as if all emitters are concentrated at one end of the segment. The following can be proven:

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}=-\ln \left({\displaystyle \frac{\mathrm{1}-\exp \left(-[{\mathrm{O}}_{\mathrm{2}}{]}_j{\mathit{\sigma }}_{\mathit{\lambda },j}{L}_{ij}\right)}{[{\mathrm{O}}_{\mathrm{2}}{]}_j{\mathit{\sigma }}_{\mathit{\lambda },j}{L}_{ij}}}\right).\end{array}$$

Moreover, such an effective optical depth of the emitting segment approaches half of its optical depth as a transmitting-only segment (i.e., $\left(\right[{\mathrm{O}}_{\mathrm{2}}{]}_j{\mathit{\sigma }}_{\mathit{\lambda },j}{L}_{ij})/\mathrm{2}$) when it is optically thin. The proof is given in Appendix A.

Equation (15) serves as the forward model for limb-viewed radiance, and the Jacobians, with respect to parameters to be retrieved (i.e., the state vector elements), can be obtained by finite difference. To reduce the computing cost, we derive the analytical Jacobians when possible. The derivations of analytical Jacobians with respect to temperature, emitting O₂ number density, and the relative change in the ground state O₂ number density are provided in Appendix B. The analytical Jacobians are also validated with finite difference. In addition, we include a squeeze factor of the SCIAMACHY instrument line shape (ILS) and a wavelength shift in the state vector. Before each inversion iteration, the radiance and Jacobians at all tangent heights are concatenated, convolved with the ILS, and sampled at the shifted SCIAMACHY wavelength grid. The Jacobians with respect to all state vector elements are then assembled into the full Jacobian matrix to be used in the inversion.

3.3 Optimal estimation-based inversion

After background subtraction, the limb-viewed radiance from tangent heights within the retrieval range is concatenated as the measurement vector r. The forward model described in the previous section maps the state vector x into the measurement space and is denoted in the following as f(x):

$$\begin{array}{}\text{(18)}& \mathit{r}=\mathit{f}\left(\mathit{x}\right)+\mathit{e},\end{array}$$

where e is the measurement error representing SCIAMACHY detector noise and any other inadequacy of the retrieval system in reproducing the limb-observed radiance. We denote the error covariance matrix of e as S_r and assume no correlation between individual measurements. The diagonal elements of S_r are calculated as the sum of a scaled radiance and the square of the readout noise. The scaled radiance term approximates the detector shot noise and model–data discrepancy, whose variances are assumed to scale linearly with radiance. The scaling factors applied to the radiance are tuned to balance the model–data discrepancy and retrieval–prior discrepancy. Their values are chosen to be 5×10⁸ for the O₂ ¹Δ band and 1×10⁷ for the O₂ A band. The readout noise is approximated by the standard deviation of the out-of-band radiance for each limb-viewed spectrum.

The retrieval algorithm minimizes the cost function as follows:

$$\begin{array}{}\text{(19)}& \begin{array}{rl}{\mathit{\chi }}^{\mathrm{2}}& ={\left(\mathit{r}-\mathit{f}\left(\mathit{x}\right)\right)}^T{\mathbf{S}}_{\mathrm{r}}^{-\mathrm{1}}\left(\mathit{r}-\mathit{f}\left(\mathit{x}\right)\right)\\ & +{\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right)}^T{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\left(\mathit{x}-{\mathit{x}}_{\mathrm{a}}\right),\end{array}\end{array}$$

where x_a is the prior vector, and S_a is the prior error covariance matrix. A smaller prior error leads to heavier influence by the prior values and vice versa. For each retrieval, we first conduct a linear inversion of the spectra and use the vertical mean value of the inverted $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ profile as the prior values for the $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ profile. Note that the linear inversion ignores self-absorption and leads to inaccurate profile shape (Sun et al., 2018a), but the purpose here is only to obtain a very rough estimate of $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$. The corresponding prior error is set to be 100 times the prior, which is a constant value for all altitudes. This effectively gives no prior constraint to the $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ profile and assures its information all comes from observations through near-unity DOFS of retrieved $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$. The a priori temperature profile is sampled from the MSIS model. Retrieval of temperature becomes difficult for the O₂ ¹Δ band in the stratosphere because of strong self-absorption and interference from upper layers. To stabilize the temperature retrieval, we apply a tighter prior error of 10 K below 50 km and relax to 30 K above 50 km. A logistic function with a steepness scale of 2.5 km is applied to smooth the transition at 50 km. Above 90 km, which is only observed in the MLT mode, the temperature prior error is further relaxed to 60 K to account for the potentially large observation–prior difference in the thermosphere. In the state vector, we also include changes in the O₂ profile relative to the a priori from MSIS, which is equivalent to retrieving the natural log of [O₂]. The prior values are set to be zero, i.e., no change to [O₂] from MSIS, and the prior errors are set at 50 %. The O₂ number density in the fitting is not constrained with temperature and pressure through the ideal gas law or hydrostatic balance. This lends some freedom to the O₂ number density in deviating from the MSIS model and slightly enhances the goodness of fit in layers with large self-absorption. The DOFS of the retrieved ln ([O₂]) also serves as a good qualitative indicator for the emergence of O₂ self-absorption when going from high to low altitudes. The prior errors of each individual profile ($\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$, temperature, and ln ([O₂])) are correlated internally with a length scale of 1 scale height.

The forward model is nonlinear, so the cost function has to be minimized iteratively. We adopt the Levenberg–Marquardt modification of the Gauss–Newton method. On each iteration i, we solve for the state vector update dx_i+1, as follows:

$$\begin{array}{}\text{(20)}& \begin{array}{rl}d{\mathit{x}}_{i+\mathrm{1}}=& {\left((\mathrm{1}+\mathit{\gamma }){\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}+{\mathbf{K}}_i^T{\mathbf{S}}_{\mathrm{r}}^{-\mathrm{1}}{\mathbf{K}}_i\right)}^{-\mathrm{1}}\\ & \left({\mathbf{K}}_i^T{\mathbf{S}}_{\mathrm{r}}^{-\mathrm{1}}(\mathit{r}-\mathit{f}({\mathit{x}}_i\left)\right)-{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}({\mathit{x}}_i-{\mathit{x}}_{\mathrm{a}})\right),\end{array}\end{array}$$

where K_i is the full Jacobian matrix (i.e., $\partial \mathit{r}/\partial \mathit{x}$) at iteration i, and the state vector is initialized by the a priori, i.e., x₀=x_a. γ is the Levenberg–Marquardt parameter that helps stabilize the iteration compared with the standard Gauss–Newton method. We initialize and update the value of γ and determine convergence following the OCO-2 retrieval algorithm (Crisp et al., 2021).

After convergence is achieved, the posterior error covariance matrix can be calculated as follows:

$$\begin{array}{}\text{(21)}& \widehat{\mathbf{S}}={\left({\mathbf{K}}^T{\mathbf{S}}_{\mathrm{r}}^{-\mathrm{1}}\mathbf{K}+{\mathbf{S}}_{\mathrm{a}}^{-\mathrm{1}}\right)}^{-\mathrm{1}},\end{array}$$

where K is the Jacobian matrix in the final iteration, and the averaging kernel matrix is as follows:

$$\begin{array}{}\text{(22)}& \mathbf{A}=\widehat{\mathbf{S}}{\mathbf{K}}^T{\mathbf{S}}_{\mathrm{r}}^{-\mathrm{1}}\mathbf{K}.\end{array}$$

In this study, the DOFS of each state vector element refers to the corresponding diagonal element in the averaging kernel (Liu et al., 2010).

4 Results

The algorithm described in the previous section is applied to retrieve O₂ ¹Δ-band airglow from SCIAMACHY nominal limb observations and both ¹Δ- and A-band airglow from SCIAMACHY MLT mode limb observations. For nominal limb observations, the measurement data are limited to tangent heights above 25 km and below 100 km (usually the uppermost tangent height for nominal limb scans is at ∼ 90 km). When the solar zenith angle (SZA) is above 70^∘, the lower height bound is gradually increased to 40 km as the airglow vanishes at lower altitude. The A-band airglow retrieval is not attempted using the nominal limb observations due to limited coverage.

4.1 Retrieval demonstration using individual vertical soundings

Figure 2 demonstrates the spectral fitting of a vertical sounding of SCIAMACHY O₂ ¹Δ-band nominal limb spectra on 3 January 2010 at 28.0^∘ N, 99.5^∘ E. In total, 10 limb observations (i.e., N=10) with tangent heights ranging from 28.4 to 87.4 km are included. Figure 2a shows the concatenated radiance spectra, and Fig. 2b shows the fitting residuals. The predicted 95 % confidence intervals, approximated by twice the measurement error, are also overlaid with the residuals. The measurement errors at the ¹Δ band are dominated by the portion that scales with radiance, whereas readout noise that is constant across the band is relatively insignificant, except at the highest two tangent heights. This highlights the importance of appropriately considering the portion of measurement errors that changes with radiance. The measurement error would be grossly underestimated if only the out-of-band variance were considered. The goodness of fit is indicated by the χ² value of 1.74 at convergence. Retrieval using the HITRAN2016 line list gives a slightly higher χ² of 1.79 (result not shown), indicating a more accurate O₂ ¹Δ-band spectroscopy in HITRAN2020. Indeed, the latest edition of the database makes use of very accurate experimental work that was recently carried out in Grenoble (Konefał et al., 2020; Tran et al., 2020) and at National Institute of Standards and Technology (NIST; Fleurbaey et al., 2021). These measurements allowed improved spectroscopic models that efficiently decouple contributions of the magnetic dipole and much weaker but not negligible electric quadrupole transitions. The theoretical background is provided in Gordon et al. (2010) and Mishra et al. (2011).

Figure 2 (a) SCIAMACHY limb-viewed spectra concatenated over 10 tangent heights within one vertical sounding (black), simulated concatenated spectra using the prior values of the state vector (blue), and simulated concatenated spectra using the posterior values of the state vector (red). Each limb-viewed ¹Δ-band spectrum contains 77 data points, so the measurement vector length is 770. Altitudes of each tangent height are labeled above the corresponding spectra. The radiance unit is photons cm⁻² s⁻¹ nm⁻¹ sr⁻¹. (b) Residuals from the fit in panel (a) and the 95 % confidence intervals of residuals are shown.

The retrieved atmospheric profiles from the vertical sounding in Fig. 2 are shown in Fig. 3. The retrieved O₂(a¹Δ_g) number density is displayed in Fig. 3a, with error bars indicating twice the posterior error or approximately the 95 % confidence interval. As demonstrated in previous studies, O₂ ¹Δ-band airglow concentrates in the lower mesosphere and upper stratosphere and peaks above the stratopause. The retrieved [O₂(a¹Δ_g)] becomes increasingly uncertain down to the stratosphere because (1) the lower the layer, the smaller the number of tangent views that can detect them; (2) those limb views include the contributions by emissions from all layers above, and (3) O₂ self-absorption becomes significant. The nominal limb observations that stop at ∼ 90 km capture most of the upper tail of O₂ ¹Δ-band airglow, as supported by comparison to the results of the MLT mode retrieval shown in Figs. 4 and 5 that extend beyond 120 km. The DOFS for the retrieved [O₂(a¹Δ_g)] is effectively unity at all retrieval altitudes, as shown in Fig. 3c. The retrieved temperature profile with posterior uncertainty is shown in Fig. 3b together with MSIS temperature that is used for prior values (green) and a colocated MIPAS temperature profile (black). The MIPAS observation is located at 28.7^∘ N, 99.1^∘ E with a 91 km spatial separation and a 15 min temporal separation from SCIAMACHY. In the mesosphere, the retrieved SCIAMACHY temperature deviates from prior values and matches closely with independent MIPAS observations. This lends confidence to the spectral modeling of the O₂ ¹Δ-band airglow and the forward model for limb-viewed spectra. The DOFS for the retrieved temperature is close to unity above 50 km and quickly drops below 50 km (Fig. 3c), indicating heavier prior influence for the stratospheric temperature. This is by design as it is challenging to extract observational information for the stratosphere due to O₂ self-absorption, whereas the prior stratospheric temperature is more trustworthy. Also shown in Fig. 3c is the DOFS for the relative changes in the O₂ number density or, equivalently, ln ([O₂]). The influence of O₂ self-absorption is felt by the retrieval below 80 km and becomes significant below 60 km. Below 40 km, the DOFS drops due to the general loss of observational constraint. The DOFS for [O₂(a¹Δ_g)] does not appear to drop because its prior constraint is much lower.

Figure 3 (a) Retrieved [O₂(a¹Δ_g)] profile from the spectral fitting in Fig. 2. (b) Retrieved temperature profile simultaneously with the [O₂(a¹Δ_g)] (pink), colocated MSIS model temperature (green), and colocated MIPAS observation (black). Error bars indicate twice the posterior error in panels (a) and (b). (c) DOFS for the retrieved [O₂(a¹Δ_g)], temperature, and ln ([O₂]) profiles. Note that the profiles are plotted at the layer center altitude that is ∼ 3.3 km above each tangent height.

Figure 4 Similar to Fig. 2 but using 10 limb views in an MLT vertical sounding.

Figure 5 Similar to Fig. 3 but using an MLT vertical sounding.

The same O₂ ¹Δ-band retrieval is applied to an MLT mode vertical sounding on 4 January 2010 at 55.8^∘ N, 92.0^∘ E. The spectral fitting is shown in Fig. 4. The lowest tangent height of this particular MLT sounding is at 57 km, and the upper limit is set to below 120 km, above which the airglow is negligible. The retrieval χ² of 0.70 (0.74 if using HITRAN2016) is significantly lower than the nominal vertical sounding, indicating a better fitting quality. This is consistent with the fact that most challenges of the ¹Δ-band retrieval are in the stratosphere, which is not observed in the MLT mode.

Figure 5 shows profiles retrieved from the spectral fit in Fig. 4. The portions below 90 km are qualitatively similar to the nominal mode retrieval in Fig. 3, although the MLT mode reveals further information at the mesopause region. There is a minor peak of [O₂(a¹Δ_g)] at near 90 km, whose maximum is ∼ 20 times lower than the major peak near the stratopause (comparing Figs. 5a and 3a). The nominal vertical soundings miss the upper part of the minor peak. Figure 5b compares retrieved MLT temperature with MSIS temperature (green) and colocated ACE-FTS temperature (black). The ACE-FTS observation is located at 55.5^∘ N, 92.7^∘ E, with a spatial separation of 56 km and a temporal separation of 1:57 h from SCIAMACHY. The retrieved temperature from SCIAMACHY ¹Δ-band airglow agrees well with ACE-FTS in the upper mesosphere but becomes more uncertain above 95 km, where the ¹Δ-band airglow becomes too dim to adequately reflect temperature from its spectral shape. The behaviors of DOFS in Fig. 5c are consistent with the nominal mode retrieval where they overlap and show the loss of observational constraint on temperature above 95 km due to fading airglow.

Figure 6 Similar to Fig. 4 but showing the O₂ A-band fitting using 12 limb views from the same MLT vertical sounding. Each limb-viewed A-band spectrum contains 62 data points, so the entire measurement vector has 744 elements.

The O₂ A-band retrieval is applied to the same MLT mode vertical sounding and shown in Fig. 6. The upper limit of the tangent height is increased to 130 km as the A-band airglow extends further into the thermosphere. The airglow radiance observed at the A band is 2 orders of magnitude lower than the ¹Δ band and closer to the readout noise level, which appears to dominate the total measurement uncertainty. The residual does not show significantly larger variance at high radiance, in contrast to the ¹Δ band. Overall, the A-band airglow spectra shapes are well captured. The χ² value for this fit is 0.77, with little difference between HITRAN2016 and HITRAN2020.

Figure 7 shows retrieved profiles from the A-band spectral fitting in Fig. 6. The emitting O₂ number density ($\left[{\mathrm{O}}_{\mathrm{2}}\right(b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}\left)\right]$) can be retrieved with high precision above ∼ 90 km. The retrieval error quickly grows below the mesopause region (below 80 km) due to O₂ self-absorption and loss of observational constraint. The temperature profiles from both A-band and ¹Δ-band airglow are shown together with MSIS and ACE-FTS temperature profiles in Fig. 7b. Because the A-band airglow extends further into the thermosphere, the A-band-derived temperature shows less uncertainty than the ¹Δ-band-derived temperature above the mesopause. Below the mesopause, the ¹Δ-band-derived temperature shows better precision. Both SCIAMACHY-based temperature profiles capture the double mesopause feature observed by ACE-FTS and agree well between each other at 80–95 km, where both airglow are strong and less interfered. These indicate that the airglow spectral model can be reliably applied to both O₂ A and ¹Δ bands.

Figure 7 Similar to Fig. 5 but in the A band. The A-band-derived emission of the O₂ number density and temperature are shown in blue. In panel (b), the ¹Δ-band-derived temperature, the same one shown in Fig. 5b, is overlaid in pink.

4.2 Climatology of O₂ ¹Δ-band airglow and temperature retrieved from nominal limb mode observations

The [O₂(a¹Δ_g)] profiles are retrieved at eight vertical soundings across track. A numerical artifact appears when integrating the [O₂(a¹Δ_g)] profiles vertically, as four vertical soundings are defined at staggered altitudes compared to the other four. Therefore, we interweave the [O₂(a¹Δ_g)] profiles retrieved from the adjacent two across-track positions into one single profile and integrate the interwoven profile to obtain the column number density of O₂(a¹Δ_g). As a result, the O₂(a¹Δ_g) column number density is defined at four across-track positions instead of eight. The latitude and longitude of the paired across-track positions are averaged from the ones of the original vertical soundings. Figure 8 exemplifies the column number density of O₂(a¹Δ_g) calculated from 1 d of retrieval on 3 January 2010 using the SCIAMACHY nominal limb observations. Out of the 14 orbits in Fig. 8, one (orbit no. 41015) is highlighted to show the spatial extent of a single orbit. Some missing sounding locations are identifiable and due to failure to converge. These missing points are generally located at high latitudes and high solar zenith angles. In these transition regions between daytime and nighttime, the horizontal variation in airglow intensity is significant, which violates the homogeneous layer assumption for the retrieval algorithm. Retrieved data are often available for part of the ascending phase of the orbit at the summer hemisphere (most valid data are in the descending phase), leading to some repeated observations at the same latitude, although at different SZAs and potentially in the nighttime. To eliminate such a latitudinal ambiguity and nighttime data, we remove the ascending portion when averaging over multiple days and limit the SZA to within 100^∘. The temperature sounding locations of MIPAS on the same day are shown in Fig. 8 as small gray and black symbols. The MIPAS temperature will be used for comparison in Sect. 5.

Figure 8 Large symbols are colored by the O₂(a¹Δ_g) column number density retrieved from SCIAMACHY nominal limb orbits (nos. 41010–41023) on 3 January 2010. A single SCIAMACHY orbit (no. 41015) is highlighted in opaque colors with a black edge. Small gray symbols mark the temperature sounding locations from MIPAS that share orbits with SCIAMACHY. A single MIPAS orbit (no. 41023) is highlighted in black.

Figure 9 Monthly average of O₂(a¹Δ_g) column number density retrieved from SCIAMACHY nominal limb orbits.

Figure 10 Vertical–latitudinal distribution of [O₂(a¹Δ_g)]. Each data point in the retrieved profiles is binned to 3D grid boxes with sizes of 3${}^{\circ }\times$3${}^{\circ }\times$3.2 km and then averaged zonally.

We then apply the retrieval algorithm to the O₂ ¹Δ-band spectra from all SCIAMACHY nominal limb observations over the year 2010. The use of HITRAN2020 improves the goodness of fit in mid- to high latitude but not in the tropical region, especially over land. This is not fully understood, as HITRAN2020 clearly reduces some structural residuals through more accurate line intensities. Figure 9 illustrates the global maps of monthly averaged O₂(a¹Δ_g) column number density in 2010. The sounding locations are binned to a 3 × 3^∘ grid using the drop-in-the-box method (Sun et al., 2018b). Although the single-day coverage is sparse (Fig. 8), continuous coverage is achievable at a monthly scale. In general, the value in each grid box is averaged from 5–10 sounding profiles. The maximum abundance of emitting O₂ follows the amount of solar radiation closely and generally colocates with the subsolar latitude. The zonal variation is insignificant, and no other spatial pattern of O₂ ¹Δ-band airglow can be convincingly identified. Figure 10 displays the vertical and latitudinal distribution of O₂ ¹Δ-band airglow, represented by the zonal mean of retrieved [O₂(a¹Δ_g)] profiles for each month in 2010. The seasonal shift of airglow following the subsolar latitude is consistent with the spatial distribution of O₂(a¹Δ_g) column number density in Fig. 9. The vertical peak locations of [O₂(a¹Δ_g)] at 45–50 km does not appear to change significantly over the low and midlatitudes or at different times of the year. The spatial and temporal distributions of O₂ ¹Δ-band airglow are in line with previous reports using SCIAMACHY and OSIRIS limb observations (Wiensz, 2005; Bertaux et al., 2020; Li et al., 2020), while providing a more complete characterization.

The O₂ ¹Δ-band airglow is generated by several mechanisms (Zarboo et al., 2018; He et al., 2019; Bertaux et al., 2020), and the most import mechanism of O₂(a¹Δ_g) production is the solar UV photolysis of ozone. The primary ozone layer exists in the stratosphere, resulting from the photolysis via the Herzberg continuum. A secondary ozone layer, generated by photolysis via the Schumann–Runge continuum, has been observed in the mesopause region (Smith and Marsh, 2005; Li et al., 2020). The O₂ ¹Δ-band airglow also features two peaks. The declining ozone concentration and increasing UV with altitude give rise to the primary peak of [O₂(a¹Δ_g)] at 45–50 km at low and midlatitudes. Figure 11 gives a closer look at the dependency of [O₂(a¹Δ_g)] profiles on the SZA. The retrieved [O₂(a¹Δ_g)] profiles north of the subsolar latitude over January 2010 are averaged over 2^∘ SZA bins. The [O₂(a¹Δ_g)] profile shape shows little variation for SZA < 50^∘, whereas the primary airglow peak shifts to higher altitudes as the SZA further increases. At the maximum observed SZA of 96^∘ (left-most line in Fig. 11), the primary peak occurs at around 60 km. This is because, at large SZAs, the solar UV radiation no longer penetrates into the stratosphere and is mostly absorbed by ozone. The secondary [O₂(a¹Δ_g)] peak corresponding to the mesopause ozone layer is vaguely identifiable at 80–100 km and does not vary much over the SZA range shown in Fig. 11.

Figure 11 Binned profiles of [O₂(a¹Δ_g)] for SZAs between 32–96^∘ at 2^∘ intervals during January 2010. Only data north of the subsolar latitude are used.

Figure 12 Retrieved temperature using the ¹Δ-band airglow observed by nominal limb orbits. Only temperature points with DOFS greater than 0.5 are included in the averaging, and only grid boxes with more than five temperature points are shown.

Vertical and latitudinal distributions of temperature retrieved from O₂ ¹Δ airglow in the nominal limb mode is shown in Fig. 12 for each month in 2010. Only temperature retrievals with DOFS greater than 0.5 are included in the averaging, and only grid cells with more than five averaged measurements are shown, which limits the temperature profiles to above ∼ 40 km. The retrieved temperature profiles peak at the stratopause region (40–60 km) due to the ozone absorption of the solar UV radiation. The temperature minima in the mesopause region (80–100 km) are also largely captured. Temperature in the mesopause region is determined by a combination of radiative effects, including the absorptive heating of solar UV radiation by O₂ and ozone (Mlynczak and Solomon, 1993) and radiative cooling of CO₂ infrared emissions (Rodgers et al., 1992), chemical effects of exothermic reactions of odd hydrogen and odd oxygen (Mlynczak and Solomon, 1993), and the dynamics of stratosphere–mesosphere, including turbulence heating (Fritts and Vanzandt, 1993; Lübken et al., 1993) and the vertical transfer of heat by the upwind and gravity waves over the summer mesosphere. As shown in Fig. 12, the mesopause temperature is cold in polar summer and relatively warm in polar winter, driven by dynamical effects. In the polar summer, stronger radiative heating from ozone and O₂ leads to upward air motion in the mesosphere, which in turn causes an adiabatic expansion of rising air and cools the mesopause. The upward motion that drives the cooling process is amplified by breaking upward gravity waves in the polar summer caused by the stratospheric easterly winds. An opposite process in the winter pole causes an adiabatic compression that results in a warm mesosphere (Björn, 1984; Smith, 2004). Another interesting feature in Fig. 12 is the implausibly high temperature in the northern winter at a high latitude and large SZA (generally larger than 90^∘). A closer examination of the individual fits at those sounding locations reveals large residuals at lower tangent heights than 60 km, abnormally high χ², and slow (if ever) convergence. The same region also has a large warm bias relative to ACE-FTS and MIPAS (to be shown in Sect. 5). We believe those are retrieval artifacts potentially due to the horizontal gradient of airglow intensity at sunrise or sunset conditions.

5 Comparison of temperature profiles

Here we focus on comparing the retrieved temperature profiles with other instruments, as temperature is the most challenging part of the state vector to be retrieved. Reliable temperature retrievals will support the validity of the algorithm and product. An accurate upper atmospheric temperature climatology will also help better simulate the airglow spectral shapes.

5.1 Internal temperature comparison using MLT mode orbits

For the MLT mode, it is possible to intercompare temperature retrievals at the same observing locations using both the O₂ ¹Δ- and A-band airglow. Figure 13 compares temperature profiles retrieved by both bands in an MLT orbit (orbit no. 41252 on 19 January 2010) and colored by the altitude of retrieved layer. Figure 13a shows all data with DOFS larger than 0.1 over the full vertical range of profiles (roughly 50–120 km). The reddish points above 300 K are generally above 100 km and in the thermosphere. The other points below 300 K can be further classified into the mesopause region (80–100 km) and the mesosphere (below 80 km). Figure 13b magnifies the mesopause region with a much higher DOFS threshold of 0.8 and shows a tight correlation (Pearson correlation coefficient r=0.91) between ¹Δ- and A-band temperatures. As illustrated in the retrieved profiles in Figs. 5 and 7, 80–100 km is the sweet spot where both ¹Δ- and A-band airglow are strong and generally free from self-absorption. Although the ¹Δ- and A-band retrievals use the same MSIS temperature as the prior, the information in the mesopause temperature is mostly from observations with high DOFS values. Out of the 504 temperature pairs with DOFS larger than 0.1 for both bands, 449 pairs have DOFS larger than 0.8. Below the mesopause region, the A-band-derived temperature becomes unreliable due to interference from self-absorption and loss of observational constraint. Over the full vertical range, the differences between those two temperatures are within the 95 % confidence intervals labeled by the horizontal and vertical error bars in Fig. 13, which indicates that the posterior error generated by the optimal-estimation-based algorithm adequately characterizes the uncertainty.

Figure 13 (a) Temperatures retrieved from both ¹Δ and A bands in an MLT mode orbit on 19 January 2010. Gray error bars indicate the 95 % confidence intervals approximated by twice the posterior error standard deviation. Data are filtered by DOFS greater than 0.1 for both bands and colored by altitude. The red dashed line is 1 : 1. Panel (b) is similar to panel (a) but only shows data between 80 and 100 km. Data with DOFS greater than 0.1 but less than 0.8 are in gray (N=504). Data with DOFS greater than 0.8 are colored by altitude (N=449). The Pearson correlation coefficient (r) and root mean squared error (RMSE) are also labeled.

Figure 14 Comparisons between temperatures retrieved from ¹Δ- and A-band airglow over the mesopause region. Each panel includes all MLT orbits in each month of 2010.

Figure 14 extends the temperature comparison between the O₂ ¹Δ and A bands over the mesopause region to all MLT orbits in each month of 2010. The MLT mode observations are conducted on 2 d in each month, so there are ∼ 30 MLT orbits per month and 352 MLT orbits in 2010. The total number of data pairs, Pearson correlation coefficient, RMSE, and mean bias are labeled in the corresponding panel of each month. The mesopause temperature spans broader ranges in both polar winter and summer (January and December, as well as June and July) than in other months due to the dynamics-driven cold summer polar mesopause. Overall, the correlation coefficients of mesopause temperatures from the two bands ranges from 0.75 to 0.90, with the A-band temperature colder than the ¹Δ-band temperature by 5–8 K. The low bias in the A-band temperature is likely caused by errors propagated from the lower altitude, where retrieving temperature from A-band airglow becomes challenging due to strong self-absorption. The RMSE between the two bands is 10–12 K.

5.2 MLT and nominal limb temperature comparison with ACE-FTS

The ACE-FTS temperature retrieval covers the vertical ranges of airglow-derived temperature in both MLT and nominal modes of SCIAMACHY and can be compared with all temperature retrievals in this study. With coincidence criteria of 500 km and 2 h, we identify ∼ 170 SCIAMACHY MLT vertical soundings and ∼ 1800 nominal vertical soundings that colocate with ACE-FTS soundings. The colocated soundings are sparse and skew toward mid- to high latitudes (70–35^∘ S and 42–70^∘ N). Figure 15a and b display the mean bias and RMSE within 3.3 km vertical bins between the temperature profiles retrieved from ¹Δ (pink) and A (blue) bands relative to the colocated temperature profiles from ACE-FTS. The statistics for the prior temperature from MSIS are shown in green. The retrieved temperature data are limited to DOFS greater than 0.5. Compared with the prior MSIS temperature, both ¹Δ and A-band-derived posterior temperatures demonstrate smaller RMSE relative to ACE-FTS, except at the high end (above 100 km for the A band and above 95 km for the ¹Δ band), where the airglow weakens, and at the low end (below 72 km) for the A band, where the O₂ self-absorption strongly reduces observational constraints (see Fig. 7). The mean bias for the ¹Δ-band temperature stays generally less than ±5 K below 90 km. In the mesopause region (80–100 km), the A-band temperature outperforms the ¹Δ-band temperature with lower absolute mean bias and RMSE relative to ACE-FTS. However, the quality of A-band temperature quickly degrades below 80 km. The relative mean biases are consistent with the internal airglow temperature comparisons in Sect. 5.1. To merge these two temperature retrievals from the MLT observations, one may use the ¹Δ-band-derived temperature below 80 km, use the A-band-derived temperature above 90 km, and average the two with linearly varying weights between 80 and 90 km.

Figure 15 The mean bias (a) and RMSE (b) of ¹Δ band (pink) and A band (blue) temperature retrieved from the MLT orbits relative to the colocated ACE-FTS profiles are shown. Panels (c) and (d) show similar statistics between the ¹Δ-band temperature retrieved from nominal limb orbits and the colocated ACE-FTS profiles. Retrievals using HITRAN2016 spectroscopic data are shown as dotted lines, and retrievals using HITRAN2020 are shown as solid lines with the same color.

A similar comparison is made between the ¹Δ-band-derived temperature in the nominal limb orbits and ACE-FTS in Fig. 15c–d. The overall shapes of the mean biases of the ¹Δ-band temperature are consistent between Fig. 15a and c, with the absolute mean bias in the nominal limb orbits slightly lower. The ¹Δ-band-derived temperature agrees well with the ACE-FTS temperature (mean bias $<\pm$5 K and significant reduction of RMSE relative to the prior) from 55 to 90 km. The substantial warm bias relative to ACE-FTS seen below 55 km originates from the retrieval artifact in northern high latitude winter and is overly amplified because the SCIAMACHY vs. ACE-FTS colocations coincidentally cluster in northern high latitude winter. A more clear picture will be depicted in the latitudinally resolved comparison with MIPAS in the following subsection. Results using HITRAN2016 spectroscopic data are shown as dotted lines in Fig. 15, while results using HITRAN2020 are in solid lines. Overall, the ¹Δ-band retrievals become warmer and show more pronounced artifacts near the northern high latitude stratopause when using HITRAN2020.

5.3 Nominal limb temperature comparison with MIPAS

Both MIPAS and SCIAMACHY were aboard the Envisat platform. Consequently, the MIPAS temperature soundings are spatiotemporally close to the SCIAMACHY soundings for all available orbits (see an example of 1 d SCIAMACHY and MIPAS orbits in Fig. 8). Here we focus on the MIPAS comparison to ¹Δ-band-derived temperature for the nominal limb orbits. Figure 16 shows the latitudinally and vertically resolved mean bias between O₂ ¹Δ-band-derived temperature and MIPAS nominal mode temperature. The mean bias is calculated by subtracting the interpolated, same-orbit MIPAS temperature from each SCIAMACHY temperature profile and then regridding the resultant bias profile in the same way that is shown in Fig. 12. To avoid interpolation artifacts, we use the extended profiles in MIPAS Level 2 data, which fill the space above the highest retrieval level with a seasonally and diurnally varying climatology. As a result, the comparison should be limited to below the top MIPAS nominal tangent height at 70 km. The mean bias is in general within ± 5 K in the mesosphere and upper stratosphere throughout the year, and the airglow-derived temperature tends to be warmer than MIPAS. García-Comas et al. (2014) compared MIPAS with a range of satellite and ground-based temperature observations and found that MIPAS temperature differs from others by 2 K at 50–80 km in spring, autumn, and winter at all latitudes and in summer at low- to midlatitudes. Differences between MIPAS and the other instruments in the summer high latitudes are typically smaller than 2 K at 50–65 km and 5 K at 65–80 km. MIPAS in general shows colder temperatures in the mid-mesosphere.

Figure 16 Monthly mean biases between temperature retrieved from O₂ ¹Δ airglow in nominal limb orbits and the MIPAS nominal mode temperature. MIPAS temperature in the same orbit is first interpolated to the latitude and altitude of SCIAMACHY layers. The longitudinal difference is neglected. The comparison should be limited to below 70 km, which is the maximum height for MIPAS nominal mode.

Figure 17 Monthly RMSE between temperature retrieved from O₂ ¹Δ airglow in nominal limb orbits and the MIPAS nominal mode temperature product.

Figure 18 Similar to Fig. 16 but using MIPAS MA, UA, and NLC observation modes instead of the nominal mode.

A consistent cold bias up to 10 K exists in the Southern Hemisphere stratopause region. As the DOFS for ¹Δ-band-derived temperature drops quickly below 50 km (Fig. 3), this bias stems mainly from a cold bias of the MSIS prior relative to MIPAS. The implausibly warm winter polar stratopause seen in Fig. 12 is more prominent in Fig. 16, with warm biases of over 20 K compared with MIPAS. Although it is clearly a retrieval artifact, we choose to keep those soundings because the temperature retrieval in the upper portion is still valid, and the retrieved emitting O₂ number densities are valuable for information about the airglow distributions at large SZA. Figure 17 shows the RMSE between the ¹Δ-band temperature and MIPAS temperature. The RMSE values are largely below 10 K throughout the upper stratosphere and mesosphere with significant portions below 5 K. The exceptions are the aforementioned southern stratopause and northern polar winter stratopause biases. The fact that the RMSE values are similar to the absolute mean bias values indicates low variances in the retrieved ¹Δ-band temperature.

In addition, we compare the O₂ ¹Δ-band-derived temperature with the MIPAS IMK/IAA temperature product for the MA, UA, and NLC modes, as shown in Fig. 18. The numbers of coincidence with SCIAMACHY for these modes are ∼ 20 % of the MIPAS nominal mode, but they provide coverage above 70 km through the top of SCIAMACHY nominal mode retrieval. The SCIAMACHY–MIPAS temperature difference is consistent with the MIPAS nominal mode, as in Fig. 16. The absolute temperature difference in the 70–100 km vertical range is generally within ±5 K, except the summer mesopause at northern high latitudes.

6 Conclusions and discussion

We develop an optimal-estimation-based algorithm to retrieve O₂ ¹Δ- and A-band airglow, as well as temperature ranging from the upper stratosphere to lower thermosphere, using limb-viewed spectra by the SCIAMACHY instrument. The algorithm is applied to nominal limb orbits that are available near-daily to retrieve the ¹Δ-band airglow and MLT limb orbits that appear on 2 d per month to retrieve both ¹Δ- and A-band airglow in the year 2010. We demonstrate the monthly climatology for O₂ ¹Δ-band airglow and temperature retrieved from nominal limb observations, which will provide crucial information for the consideration of airglow for future remote sensing of the O₂ ¹Δ band. The global monthly distributions of the vertical column density of emitting O₂ in the a¹Δ_g state show mainly latitudinal dependence without other discernible geographical patterns. The O₂ ¹Δ-band-derived temperature agrees well with MIPAS Level 2 temperature with mean bias generally near or within ±5 K and RMSE below 10 K. Notable discrepancies include the stratospheric cold bias attributable to a priori influence, strong warm bias at winter polar stratopause region likely due to the horizontal gradient of the airglow intensity that violates the homogeneous layer assumption for the retrieval algorithm, and disagreements at the summer mesopause in northern high latitudes. The reliable retrieval of temperature from airglow indicates that we can confidently reproduce the spectral shape of airglow emissions.

The successful retrieval of the O₂ A-band airglow and the associated temperature demonstrates the generalizability of the airglow spectral model and the algorithm. The A-band- and ¹Δ-band-derived temperature profiles are in good agreement in the mesopause region, with the A-band-derived temperature lower by 5–8 K. Intercomparisons with ACE-FTS suggest that the A-band-derived temperature is of higher quality above ∼ 90 km due to the sharper decline of ¹Δ-band airglow with altitude. At lower altitudes, the ¹Δ-band may provide more reliable temperature retrieval. Driven by the need for nadir observations of the ¹Δ band by MethaneSAT, this study focuses on the ¹Δ-band retrieval, whereas the synergy between the ¹Δ- and A-band airglow may be further explored by future studies to improve understanding of the chemistry, radiation, and dynamics of the MLT region. It would also be interesting to include the Noxon band (Noxon, 1961) at 1.91 µm in future analyses. This band corresponds to the emission from the $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$ to the a¹Δ_g excited electronic states and is another avenue for depopulation of the photochemically produced oxygen molecules in the $b^{\mathrm{1}}{\mathrm{\Sigma }}_{\mathrm{g}}^{+}$ state. It is allowed through the electric quadrupole mechanism (Gordon et al., 2010).

The SCIAMACHY instrument made 10 years of near-daily nominal limb observations and 4 years of MLT-mode limb observations. It is a natural future step to process the entire SCIAMACHY lifetime to generate long-term records of O₂ airglow emissions and upper atmospheric temperature. Comparison with MIPAS indicates that the current a priori temperature from the MSIS model may be inadequate, so it may be desirable to use more accurate reanalyzed temperature profiles, where available, and stitch the MSIS temperature above. A recent development improved the MSIS model especially at the MLT regions (Emmert et al., 2021), and future adoption of the new MSIS profiles will improve the a priori estimates. In contrast to previous studies using the SCIAMACHY limb observations, we do not aggregate across the track, which leads to the options of keeping eight soundings across the track or interweaving each pair of adjacent soundings into four soundings across the track. The retrievals in this study are conducted at the native across-track resolution and, hence, a vertical resolution of 6.6 km. Future efforts will test the second option, which enhances the vertical resolution to 3.3 km, in line with previous SCIAMACHY limb studies (Bender et al., 2017; Zarboo et al., 2018).

Appendix A: Derivation of emitting layer optical depth

Radiance at one end of the emitting segment, with light path length L and uniform absorber/emitter concentrations, is as follows:

$$\begin{array}{}\text{(A1)}& \begin{array}{rl}{r}_{\mathit{\lambda }}& ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathit{\lambda }}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{L}{\int }}\exp \left(-\left[{\mathrm{O}}_{\mathrm{2}}\right]{\mathit{\sigma }}_{\mathit{\lambda }}(L-L^{\prime })\right)\mathrm{d}{L}^{\prime }\\ & ={\displaystyle \frac{{\mathit{\epsilon }}_{\mathit{\lambda }}L}{\mathrm{4}\mathit{\pi }}}\left({\displaystyle \frac{\mathrm{1}-\exp (-[{\mathrm{O}}_{\mathrm{2}}\left]{\mathit{\sigma }}_{\mathit{\lambda }}L\right)}{\left[{\mathrm{O}}_{\mathrm{2}}\right]{\mathit{\sigma }}_{\mathit{\lambda }}L}}\right).\end{array}\end{array}$$

The effective optical depth of the emitting segment (${\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}$), assuming all emitters concentrated at the far end, should satisfy the following:

$$\begin{array}{}\text{(A2)}& r_{\mathit{\lambda }}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathit{\lambda }}L}{\mathrm{4}\mathit{\pi }}}\exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}\right).\end{array}$$

The optical depth of this segment as a transmitting layer is simply as follows:

$$\begin{array}{}\text{(A3)}& {\mathit{\tau }}_{\mathit{\lambda }}=\left[{\mathrm{O}}_{\mathrm{2}}\right]{\mathit{\sigma }}_{\mathit{\lambda }}L.\end{array}$$

Combining Eqs. (A1), (A2), and (A3) results in the following:

$$\begin{array}{}\text{(A4)}& {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}=-\ln \left({\displaystyle \frac{\mathrm{1}-\exp (-{\mathit{\tau }}_{\mathit{\lambda }})}{{\mathit{\tau }}_{\mathit{\lambda }}}}\right).\end{array}$$

This is equivalent to Eq. (17). At the optical thin limit, applying L'Hôpital's rule, leads to the following:

$$\begin{array}{}\text{(A5)}& \begin{array}{rl}\underset{{\mathit{\tau }}_{\mathit{\lambda }}\to \mathrm{0}}{lim}{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}}{{\mathit{\tau }}_{\mathit{\lambda }}}}& =\underset{{\mathit{\tau }}_{\mathit{\lambda }}\to \mathrm{0}}{lim}{\displaystyle \frac{-\ln \left(\frac{\mathrm{1}-\exp (-{\mathit{\tau }}_{\mathit{\lambda }})}{{\mathit{\tau }}_{\mathit{\lambda }}}\right)}{{\mathit{\tau }}_{\mathit{\lambda }}}}\\ & =\underset{{\mathit{\tau }}_{\mathit{\lambda }}\to \mathrm{0}}{lim}{\displaystyle \frac{{\mathit{\tau }}_{\mathit{\lambda }}\exp (-{\mathit{\tau }}_{\mathit{\lambda }})+\exp (-{\mathit{\tau }}_{\mathit{\lambda }})-\mathrm{1}}{{\mathit{\tau }}_{\mathit{\lambda }}\left(\exp \right(-{\mathit{\tau }}_{\mathit{\lambda }})-\mathrm{1})}}\\ & =\underset{{\mathit{\tau }}_{\mathit{\lambda }}\to \mathrm{0}}{lim}{\displaystyle \frac{-\exp (-{\mathit{\tau }}_{\mathit{\lambda }})+{\mathit{\tau }}_{\mathit{\lambda }}\exp (-{\mathit{\tau }}_{\mathit{\lambda }})}{-\mathrm{2}\exp (-{\mathit{\tau }}_{\mathit{\lambda }})+{\mathit{\tau }}_{\mathit{\lambda }}\exp (-{\mathit{\tau }}_{\mathit{\lambda }})}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}.\end{array}\end{array}$$

Namely, the effective optical depth of the emitting segment due to the self-absorption is half of its optical depth as a transmitting segment at optical thin limit, which makes intuitive sense. At the optical thick limit, the effective optical depth of emitting segment due to self-absorption approaches the natural logarithm of its optical depth as a transmitting segment, as follows:

$$\begin{array}{}\text{(A6)}& \underset{{\mathit{\tau }}_{\mathit{\lambda }}\to \mathrm{\infty }}{lim}{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}=\ln {\mathit{\tau }}_{\mathit{\lambda }}.\end{array}$$

Appendix B: Jacobians of limb-viewed radiance with respect to temperature, $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$, and ln ([O₂])

The Jacobian of radiance observed at limb view i with respect to the temperature of layer j (j≥i, otherwise the Jacobian values are all zero) is as follows:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{\displaystyle \frac{\partial r_{\mathit{\lambda },i}}{\partial T_j}}& =\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\left({\displaystyle \frac{L_{ij}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\lambda },j}}{\partial T_j}}\right.\\ & \left.\times \exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}-\sum _{l^{\prime }=\mathrm{1}}^{l-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}\right)\right)\\ & +\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\left({\displaystyle \frac{L_{ij}{\mathit{\epsilon }}_{\mathit{\lambda },j}}{\mathrm{4}\mathit{\pi }}}\left(-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}}{\partial T_j}}\right)\right.\\ & \left.\times \exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}-\sum _{l^{\prime }=\mathrm{1}}^{l-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}\right)\right)\\ & +\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\sum _{l^{\prime }=l+\mathrm{1}}^{\mathrm{2}N-i}\left({\displaystyle \frac{L_{i{j}^{\prime }}{\mathit{\epsilon }}_{\mathit{\lambda },j^{\prime }}}{\mathrm{4}\mathit{\pi }}}\right.\\ & \left.\times \left(-{\displaystyle \frac{\partial {\mathit{\tau }}_{\mathit{\lambda },ij}}{\partial T_j}}\right)\exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },i{j}^{\prime }}-\sum _{l^{\prime \prime }=\mathrm{1}}^{l^{\prime }-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime \prime }}\right)\right).\end{array}\end{array}$$

Here the first term on the right-hand side reflects the sensitivity of airglow emission to the temperature of layer j (with two corresponding segments as in Eq. 14), the second term reflects the temperature sensitivity of self-absorption of layer j as an emitting layer, and the third term reflects the temperature sensitivity of layer j as a transmitting layer, where temperature alters its absorption from upstream emissions.

The temperature sensitivities of the transmitting optical depth ($\partial {\mathit{\tau }}_{\mathit{\lambda }}/\partial T$) and its effective optical depth as an emitting layer ($\partial {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}/\partial T$) for both layers can be related to the temperature sensitivity of the absorption cross section ($\partial {\mathit{\sigma }}_{\mathit{\lambda }}/\partial T$, obtained by finite difference using HAPI) by applying the chain rule to Eq. 16 and 17. We assume that the ground state O₂ number density [O₂] does not depend on the temperature and retrieval [O₂] profile independently to account for errors in the a priori O₂ profile. The temperature sensitivity of emissivity ($\partial {\mathit{\epsilon }}_{\mathit{\lambda }}/\partial T$) is given by Eq. 10. The derived temperature Jacobian in Eq. B1 is confirmed with numerical finite difference derivative of Eq. 15 at the instrument resolution and sampling grid. The comparison is shown in Fig. B1. The analytical Jacobians are consistent with finite difference within 10⁻⁶ for all tangent heights.

The Jacobians of emitting O₂ number density are derived similarly, but only involving the emitting segment, as follows:

$$\begin{array}{}\text{(B2)}& \begin{array}{rl}{\displaystyle \frac{\partial r_{\mathit{\lambda },i}}{\partial [{\mathrm{O}}_{\mathrm{2}}^{\ast }{]}_j}}& =\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\left({\displaystyle \frac{L_{ij}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{\partial {\mathit{\epsilon }}_{\mathit{\lambda },j}}{\partial [{\mathrm{O}}_{\mathrm{2}}^{\ast }{]}_j}}\right.\\ & \left.\times \exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}-\sum _{l^{\prime }=\mathrm{1}}^{l-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}\right)\right),\end{array}\end{array}$$

where $\partial {\mathit{\epsilon }}_{\mathit{\lambda },j}/\partial [{\mathrm{O}}_{\mathrm{2}}^{\ast }{]}_j$ is the Jacobian of airglow emissivity with respect to $\left[{\mathrm{O}}_{\mathrm{2}}^{\ast }\right]$ for layer j and given by Eq. (9).

Besides, we also include relative changes to the ground state O₂ number density in the state vector, which is equivalent to retrieving the natural logarithm of ground state O₂ number density, ln ([O₂]). The corresponding Jacobians are given by the following:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\displaystyle \frac{\partial r_{\mathit{\lambda },i}}{\partial \ln \left(\right[{\mathrm{O}}_{\mathrm{2}}{]}_j)}}=[{\mathrm{O}}_{\mathrm{2}}{]}_j{\displaystyle \frac{\partial r_{\mathit{\lambda },i}}{\partial [{\mathrm{O}}_{\mathrm{2}}{]}_j}}\\ & =[{\mathrm{O}}_{\mathrm{2}}{]}_j\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\left({\displaystyle \frac{L_{ij}{\mathit{\epsilon }}_{\mathit{\lambda },j}}{\mathrm{4}\mathit{\pi }}}\right.\\ & \left.\times \left(-{\displaystyle \frac{\partial {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}}{\partial [{\mathrm{O}}_{\mathrm{2}}{]}_j}}\right)\exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },ij}-\sum _{l^{\prime }=\mathrm{1}}^{l-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime }}\right)\right)\\ & +[{\mathrm{O}}_{\mathrm{2}}{]}_j\sum _{l\in \mathit{\{}N+j-i,N-j+\mathrm{1}\mathit{\}}}\sum _{l^{\prime }=l+\mathrm{1}}^{\mathrm{2}N-i}\left({\displaystyle \frac{L_{i{j}^{\prime }}{\mathit{\epsilon }}_{\mathit{\lambda },j^{\prime }}}{\mathrm{4}\mathit{\pi }}}\right.\\ & \left.\times \left(-{\displaystyle \frac{\partial {\mathit{\tau }}_{\mathit{\lambda },ij}}{\partial [{\mathrm{O}}_{\mathrm{2}}{]}_j}}\right)\exp \left(-{\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda },i{j}^{\prime }}-\sum _{l^{\prime \prime }=\mathrm{1}}^{l^{\prime }-\mathrm{1}}{\mathit{\tau }}_{\mathit{\lambda },i{j}^{\prime \prime }}\right)\right).\end{array}\end{array}$$

Here the first term on the right-hand side of Eq. (B3) reflects the sensitivity of the effective optical depth of layer j to its own O₂ number density, and the second term reflects the sensitivity of its optical depth to its O₂ number density as a transmitting layer. Both $\partial {\mathit{\tau }}_{\mathit{\lambda }}/\partial \left[{\mathrm{O}}_{\mathrm{2}}\right]$) and $\partial {\stackrel{\mathrm{̃}}{\mathit{\tau }}}_{\mathit{\lambda }}/\partial \left[{\mathrm{O}}_{\mathrm{2}}\right]$ can be readily calculated by differentiating Eqs. (16) and (17).

Figure B1 Validation of analytical Jacobians of radiance with respect to temperature, based on Eq. (B1) (black lines), with finite difference (white dots). Each vertical column represents one limb view (fixed i) over all layers (different j). In total, 10 tangent views (i.e., N=10) are included using a single vertical sounding at the O₂ ¹Δ band from SCIAMACHY orbit no. 41011 on 3 January 2010. Only Jacobians of layers above the tangent point are non-zero.