Mesospheric nitric oxide model from SCIAMACHY data

We present an empirical model for nitric oxide (NO) in the mesosphere (≈60–90 km) derived from SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartoghraphY) limb scan data. This work complements and extends the NOEM (Nitric Oxide Empirical Model; Marsh et al., 2004) and SANOMA (SMR Acquired Nitric Oxide Model Atmosphere; Kiviranta et al., 2018) empirical models in the lower thermosphere. The regression ansatz builds on the heritage of studies by Hendrickx et al. (2017) and the superposed epoch analysis by Sinnhuber et al. (2016) which estimate NO production from particle precipitation.

Our model relates the daily (longitudinally) averaged NO number densities from SCIAMACHY (Bender et al., 2017b, a) as a function of geomagnetic latitude to the solar Lyman-α and the geomagnetic AE (auroral electrojet) indices. We use a non-linear regression model, incorporating a finite and seasonally varying lifetime for the geomagnetically induced NO. We estimate the parameters by finding the maximum posterior probability and calculate the parameter uncertainties using Markov chain Monte Carlo sampling. In addition to providing an estimate of the NO content in the mesosphere, the regression coefficients indicate regions where certain processes dominate.

1 Introduction

It has been recognized in the past decades that the mesosphere and stratosphere are coupled in various ways (Baldwin and Dunkerton, 2001). Consequently, climate models have been evolving to extend to increasingly higher levels in the atmosphere to improve the accuracy of medium- and long-term predictions. Nowadays it is not unusual that these models include the mesosphere (40–90 km) or the lower thermosphere (90–120 km) (Matthes et al., 2017). It is therefore important to understand the processes in the mesosphere and lower thermosphere and to find the important drivers of chemistry and dynamics in that region. The atmosphere above the stratosphere (≳40 km) is coupled to solar and geomagnetic activity, also known as space weather (Sinnhuber et al., 2012). Electrons and protons from the solar wind and the radiation belts with sufficient kinetic energy enter the atmosphere in that region. Since as charged particles they move along the magnetic field, this precipitation occurs primarily at high geomagnetic latitudes.

Previously the role of NO in the mesosphere has been identified as an important free radical and in this sense a driver of the chemistry (Kockarts, 1980; Barth, 1992, 1995, 2010; Barth et al., 2009; Roble, 1995; Bailey et al., 2002), particularly during winter when it is long-lived because of reduced photodissociation. NO generated in the region between 90 and 120 km at auroral latitudes is strongly influenced by both solar and geomagnetic activity (Marsh et al., 2004; Sinnhuber et al., 2011, 2016; Hendrickx et al., 2015, 2017). At high latitudes, NO is transported down to the upper stratosphere during winter, usually down to 50 km and occasionally down to 30 km (Siskind et al., 2000; Randall et al., 2007; Funke et al., 2005a, 2014b). At those altitudes and also in the mesosphere, NO participates in the “odd oxygen catalytic cycle which depletes ozone” (Crutzen, 1970). Additional dynamical processes also result in the strong downward transport of mesospheric air into the upper stratosphere, such as the strong downwelling that often occurs in the recovery phase of a sudden stratospheric warming (SSW) (Pérot et al., 2014; Orsolini et al., 2017). This downwelling is typically associated with the formation of an elevated stratopause.

Different instruments have been measuring NO in the mesosphere and lower thermosphere, but at different altitudes and at different local times. Measurements from solar occultation instruments such as Scisat-1/ACE-FTS or AIM/SOFIE are limited in latitude and local time (sunrise and sunset). Global observations from sun-synchronously orbiting satellites are available from Envisat/MIPAS below 70 km daily and 42–172 km every 10 days (Funke et al., 2001, 2005b; Bermejo-Pantaleón et al., 2011); from Odin/SMR between 45 and 115 km (Pérot et al., 2014; Kiviranta et al., 2018); or from Envisat/SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartoghraphY) between 60 and 90 km daily (Bender et al., 2017b) and 60–160 km every 15 days (Bender et al., 2013). Because the Odin and Envisat orbits are sun-synchronous, the measurement local times are fixed to around 06:00 and 18:00 (Odin) and 10:00 and 22:00 (Envisat). While MIPAS has both daytime and night-time measurements, SCIAMACHY provides daytime (10:00) data because of the measurement principle (fluorescent UV scattering; see Bender et al., 2013, 2017b). Unfortunately, Envisat stopped communicating in April 2012 ,and therefore the data available from MIPAS and SCIAMACHY are limited to nearly 10 years from August 2002 to April 2012. The other aforementioned instruments are still operational and provide ongoing data as long as satellite operations continue.

Chemistry–climate models struggle to simulate the NO amounts and distributions in the mesosphere and lower thermosphere (see, for example, Funke et al., 2017; Randall et al., 2015; Orsolini et al., 2017; Hendrickx et al., 2018). To remedy the situation, some models constrain the NO content at their top layer through observation-based parametrizations. For example, the next generation of climate simulations (CMIP6; see Matthes et al., 2017) and other recent model simulations (Sinnhuber et al., 2018) parametrizes particle effects as derived partly from Envisat/MIPAS NO measurements (Funke et al., 2016).

NO in the mesosphere and lower thermosphere

NO in the mesosphere and lower thermosphere is produced by N₂ dissociation,

$$\begin{array}{}\text{(R1)}& {\displaystyle }{\mathrm{N}}_{\mathrm{2}}+h\mathit{\nu }\to \mathrm{N}{(}^{\mathrm{2}}\mathrm{D})+\mathrm{N}{(}^{\mathrm{4}}\mathrm{S})(\mathit{\lambda }<\mathrm{102}\mathrm{nm}),\end{array}$$

followed by the reaction of the excited nitrogen atom N(²D) with molecular oxygen (Solomon et al., 1982; Barth, 1992, 1995):

$$\begin{array}{}\text{(R2)}& {\displaystyle }\mathrm{N}{(}^{\mathrm{2}}\mathrm{D})+{\mathrm{O}}_{\mathrm{2}}\to \mathrm{NO}+\mathrm{O}.\end{array}$$

The dissociation energy of N₂ into ground state atoms N(⁴S) is about 9.8 eV (λ≈127 nm) (Hendrie, 1954; Frost et al., 1956; Heays et al., 2017). This energy together with the excitation energy to N(²D) is denoted by hν in Reaction (R1) and can be provided by a number of sources, most notably by auroral or photoelectrons as well as by soft solar X-rays.

The NO content is reduced by photodissociation,

$$\begin{array}{}\text{(R3)}& {\displaystyle }\mathrm{NO}+h\mathit{\nu }\to \mathrm{N}+\mathrm{O}(\mathit{\lambda }<\mathrm{191}\mathrm{nm}),\end{array}$$

by photoionization,

$$\begin{array}{}\text{(R4)}& {\displaystyle }\mathrm{NO}+h\mathit{\nu }\to {\mathrm{NO}}^{+}+e^{-}(\mathit{\lambda }<\mathrm{134}\mathrm{nm}),\end{array}$$

and by reacting with atomic nitrogen,

$$\begin{array}{}\text{(R5)}& {\displaystyle }\mathrm{NO}+\mathrm{N}\to {\mathrm{N}}_{\mathrm{2}}+\mathrm{O}.\end{array}$$

N₂O has been retrieved in the mesosphere and thermosphere from MIPAS (see, e.g. Funke et al., 2008b, a) and from Scisat-1/ACE-FTS (Sheese et al., 2016). Model–measurement studies by Semeniuk et al. (2008) attributed the source of this N₂O to being most likely the reaction between NO₂ and N atoms produced by particle precipitation:

$$\begin{array}{}\text{(R6)}& {\displaystyle }\mathrm{N}+{\mathrm{NO}}_{\mathrm{2}}\to {\mathrm{N}}_{\mathrm{2}}\mathrm{O}+\mathrm{O}.\end{array}$$

We note that photo-excitation and photolysis at 185 nm (vacuum UV) of NO or NO₂ mixtures in nitrogen, N₂, or helium mixtures at 1 atm leads to N₂O formation (Maric and Burrows, 1992). Both mechanisms explaining the production of N₂O involve excited states of NO. Hence these pathways contribute to the loss of NO and potentially an additional daytime source of N₂O in the upper atmosphere. N₂O acts as an intermediate reservoir at high altitudes (⪆90 km; see Sheese et al., 2016), reacting with O(¹D) in two well-known channels to N₂ and O₂ as well as to 2NO. However, the largest N₂O abundances are located below 60 km and originate primarily from the transport of tropospheric N₂O into the stratosphere through the Brewer–Dobson circulation (Funke et al., 2008a, b; Sheese et al., 2016) but can reach up to 70 km in geomagnetic storm conditions (Funke et al., 2008a; Sheese et al., 2016). Both source and sink reactions indicate that NO behaves differently in sunlit conditions than in dark conditions. NO is produced by particle precipitation at auroral latitudes, but in dark conditions (without photolysis) it is only depleted by reacting with atomic nitrogen (Reaction R5). This asymmetry between production and depletion in dark conditions results in different lifetimes of NO.

Early work to parametrize NO in the lower thermosphere (100–150 km) used SNOE measurements from March 1998 to September 2000 (Marsh et al., 2004). With these 2.5 years of data and using empirical orthogonal functions, the so-called NOEM (Nitric Oxide Empirical Model) estimates NO in the lower thermosphere as a function of the solar f_10.7 cm radio flux, the solar declination angle, and the planetary Kp index. NOEM is still used as prior input for NO retrieval, for example, from MIPAS (Bermejo-Pantaleón et al., 2011; Funke et al., 2012) and SCIAMACHY (Bender et al., 2017b) spectra. However, 2.5 years is relatively short compared to the 11-year solar cycle, and the years 1998 to 2000 encompass a period of elevated solar activity. To address this, a longer time series from AIM/SOFIE was used to determine the important drivers of NO in the lower thermosphere (90–140 km) by Hendrickx et al. (2017). Other recent work uses 10 years of NO data from Odin/SMR from 85 to 115 km (Kiviranta et al., 2018). Funke et al. (2016) derived a semi-empirical model of NO_y in the stratosphere and mesosphere from MIPAS data. Here we use Envisat/SCIAMACHY NO data from the nominal limb mode (Bender et al., 2017b, a). Apart from providing a similarly long time series of NO data, the nominal Envisat/SCIAMACHY NO data cover the mesosphere from 60 to 90 km (Bender et al., 2017b), bridging the gap between the stratosphere and lower thermosphere models.

The paper is organized as follows: we present the data used in this work in Sect. 2. The two model variants, linear and non-linear, are described in Sect. 3. Details about the parameter and uncertainty estimation are explained in Sect. 4, and we present the results in Sect. 5. Finally we conclude our findings in Sect. 6.

2 Data

2.1 SCIAMACHY NO

We use the SCIAMACHY nitric oxide data set version 6.2.1 (Bender et al., 2017a) retrieved from the nominal limb scan mode (≈0–93 km). For a detailed instrument description, see Burrows et al. (1995) and Bovensmann et al. (1999), and for details of the retrieval algorithm, see Bender et al. (2013, 2017b).

The data were retrieved for the whole Envisat period (August 2002–April 2012). This satellite was orbiting in a sun-synchronous orbit at around 800 km altitude, with Equator crossing times of 10:00 and 22:00 local time. The NO number densities from the SCIAMACHY nominal mode were retrieved from the NO gamma band emissions. Since those emissions are fluorescent emissions excited by solar UV, SCIAMACHY NO data are only available for the 10:00 dayside (downleg) part of the orbit. Furthermore, the retrieval was carried out for altitudes from 60 to 160 km, but above approximately 90 km, the data reflect the scaled a priori densities from NOEM (Bender et al., 2017b). We therefore restrict the modelling to the mesosphere below 90 km.

We averaged the individual orbital data longitudinally on a daily basis according to their geomagnetic latitude within 10^∘ bins. The geomagnetic latitude was determined according to the eccentric dipole approximation of the 12th generation of the International Geomagnetic Reference Field (IGRF12) (Thébault et al., 2015). In the vertical direction the original retrieval grid altitudes (2 km bins) were used. Note that mesospheric NO concentrations are related to geomagnetically as well as geographically based processes, but disentangling them is beyond the scope of the paper. Follow-up studies can build on the method presented here and study, for example, longitudinally resolved time series.

The measurement sensitivity is taken into account via the averaging kernel diagonal elements, and days where its binned average was below 0.002 were excluded from the time series. Considering this criterion, each bin (geomagnetic latitude and altitude) contains about 3400 data points.

2.2 Proxies

We use two proxies to model the NO number densities, one accounting for the solar irradiance variations and one accounting for the geomagnetic activity. Various proxies have been used or proposed to account for the solar-irradiance-induced variations in mesospheric–thermospheric NO, which are in particular related to the 11-year solar cycle. The NOEM (Nitric Oxide Empirical Model; Marsh et al., 2004) uses the natural logarithm of the solar 10.7 cm radio flux f_10.7. More recent work on AIM/SOFIE NO (Hendrickx et al., 2017) uses the solar Lyman-α index because some of the main production and loss processes are driven by UV photons. Besides accounting for the long-term variation of NO with solar activity, the Lyman-α index also includes short-term UV variations and the associated NO production, for example, caused by solar flares. Barth et al. (1988) have shown that the Lyman-α index directly relates to the observed NO at low latitudes (30^∘ S–30^∘ N). Thus we use it in this work as a proxy for NO.

In the same manner as for the irradiance variations, the “right” geomagnetic index to model particle-induced variations of NO is a matter of opinion. Kp is the oldest and most commonly used geomagnetic index; it was, for example, used in earlier work by Marsh et al. (2004) for modelling NO in the mesosphere and lower thermosphere. Kp is derived from magnetometer stations distributed at different latitudes and mostly in the Northern Hemisphere (NH). However, Hendrickx et al. (2015) found that the auroral electrojet index (AE) (Davis and Sugiura, 1966) correlated better with SOFIE-derived NO concentrations (Hendrickx et al., 2015, 2017; see also Sinnhuber et al., 2016). The AE index is derived from stations distributed almost evenly within the auroral latitude band. This distribution enables the AE index to be more closely related to the energy input into the atmosphere at these latitudes. Therefore, we use the auroral electrojet index (AE) as a proxy for geomagnetically induced NO. To account for the 10:00 satellite sampling, we average the hourly AE index from noon the day before to noon on the measurement day.

It should be noted that tests using Kp (or its linear equivalent, Ap) instead of AE and using f_10.7 instead of Lyman-α suggested that the particular choice of index did not lead to significantly different results. Our choice of AE rather than Kp and Lyman-α over f_10.7 is physically based and motivated as described above.

3 Regression model

We denote the number density by x_NO as a function of the (geomagnetic; see Sect. 2.1) latitude ϕ, the altitude z, and the time (measurement day) t: $x_{\mathrm{NO}}(\mathit{\varphi },z,t)$. In the following we often drop the subscript NO and combine the time direction into a vector x, with the ith entry denoting the density at time t_i, such that $x_i(\mathit{\varphi },z)=x(\mathit{\varphi },z,t_i)$.

3.1 Linear model

In the (multi-)linear case, we relate the nitric oxide number densities $x_{\mathrm{NO}}(\mathit{\varphi },z,t)$ to the two proxies, the solar Lyman-α index (Lyα(t)) and the geomagnetic AE index (AE(t)). Harmonic terms with $\mathit{\omega }=\mathrm{2}\mathit{\pi }{\mathrm{a}}^{-\mathrm{1}}=\mathrm{2}\mathit{\pi }(\mathrm{365.25}\mathrm{d}{)}^{-\mathrm{1}}$ account for annual and semi-annual variations. The linear model, including a constant offset for the background density, describes the NO density according to Eq. (1):

$$\begin{array}{ll}{\displaystyle }{x}_{\mathrm{NO}}(\mathit{\varphi },z,t)& {\displaystyle }=a(\mathit{\varphi },z)+b(\mathit{\varphi },z)\cdot \mathrm{Ly}\mathit{\alpha }\left(t\right)+c(\mathit{\varphi },z)\cdot \mathrm{AE}\left(t\right)\\ {\displaystyle }& {\displaystyle }+\sum _{n=\mathrm{1}}^{\mathrm{2}}\left[d_n(\mathit{\varphi },z)\cos \left(n\mathit{\omega }t\right)\right.\\ \text{(1)}& {\displaystyle }& {\displaystyle }+e_n(\mathit{\varphi },z)\left.\sin \left(n\mathit{\omega }t\right)\right].\end{array}$$

The linear model can be written in matrix form for the n measurement times $t_{\mathrm{1}},\mathrm{\dots },t_n$ as Eq. (2), with the parameter vector β given by ${\mathit{\beta }}_{\mathrm{lin}}=(a,b,c,d_{\mathrm{1}},e_{\mathrm{1}},d_{\mathrm{2}},e_{\mathrm{2}}{)}^{\top }\in {\mathbb{R}}^{\mathrm{7}}$ and the model matrix $\mathbf{X}\in {\mathbb{R}}^{n\times \mathrm{7}}$.

$$\begin{array}{}\text{(2)}& \begin{array}{rl}{\displaystyle }{\mathit{x}}_{\mathrm{NO}}(\mathit{\varphi },z)& {\displaystyle }=\left(\begin{array}{ccccc}\mathrm{1}& \mathrm{Ly}\mathit{\alpha }\left(t_{\mathrm{1}}\right)& \mathrm{AE}\left(t_{\mathrm{1}}\right)& \cos \left(\mathit{\omega }{t}_{\mathrm{1}}\right)& \sin \left(\mathit{\omega }{t}_{\mathrm{1}}\right)\\ \mathrm{⋮}\\ \mathrm{1}& \mathrm{Ly}\mathit{\alpha }\left(t_n\right)& \mathrm{AE}\left(t_n\right)& \cos \left(\mathit{\omega }{t}_n\right)& \sin \left(\mathit{\omega }{t}_n\right)\end{array}\right.\\ {\displaystyle }& {\displaystyle }\left.\begin{array}{cc}\cos \left(\mathrm{2}\mathit{\omega }{t}_{\mathrm{1}}\right)& \sin \left(\mathrm{2}\mathit{\omega }{t}_{\mathrm{1}}\right)\\ & \mathrm{⋮}\\ \cos \left(\mathrm{2}\mathit{\omega }{t}_n\right)& \sin \left(\mathrm{2}\mathit{\omega }{t}_n\right)\end{array}\right)\cdot \left(\begin{array}{c}a\\ b\\ c\\ d_{\mathrm{1}}\\ e_{\mathrm{1}}\\ d_{\mathrm{2}}\\ e_{\mathrm{2}}\end{array}\right)\\ {\displaystyle }& {\displaystyle }=\mathbf{X}\cdot \mathit{\beta }\end{array}\end{array}$$

We determine the coefficients via least squares, minimizing the squared differences of the modelled number densities to the measured ones.

3.2 Non-linear model

In contrast to the linear model above, we modify the AE index by a finite lifetime τ, which varies according to season; we denote this modified version by $\stackrel{\mathrm{̃}}{\mathrm{AE}}$. We then omit the harmonic parts in the model, and the non-linear model is given by Eq. (3):

$$\begin{array}{}\text{(3)}& x_{\mathrm{NO}}(\mathit{\varphi },z,t)=a(\mathit{\varphi },z)+b(\mathit{\varphi },z)\cdot \mathrm{Ly}\mathit{\alpha }\left(t\right)+c(\mathit{\varphi },z)\cdot \stackrel{\mathrm{̃}}{\mathrm{AE}}\left(t\right).\end{array}$$

Although this approach shifts all seasonal variations to the AE index and thus attributes them to particle-induced effects, we found that the residual traces of particle-unrelated seasonal effects were minor compared to the overall improvement of the fit. Additional harmonic terms only increase the number of free parameters without substantially improving the fit further.

The lifetime-corrected $\stackrel{\mathrm{̃}}{\mathrm{AE}}$ is given by the sum of the previous 60 days' AE values, each multiplied by an exponential decay factor:

$$\begin{array}{}\text{(4)}& \stackrel{\mathrm{̃}}{\mathrm{AE}}\left(t\right)=\sum _{t_i=\mathrm{0}}^{\mathrm{60}\mathrm{d}}\mathrm{AE}(t-t_i)\cdot \exp \left\{-{\displaystyle \frac{t_i}{\mathit{\tau }}}\right\}.\end{array}$$

The total lifetime τ is given by a constant part τ₀ plus the non-negative fraction of a seasonally varying part τ_t:

$$\begin{array}{}\text{(5)}& {\displaystyle }\mathit{\tau }& {\displaystyle }={\mathit{\tau }}_{\mathrm{0}}+\left\{\begin{array}{ll}{\mathit{\tau }}_t,& {\mathit{\tau }}_t\ge \mathrm{0}\\ \mathrm{0},& {\mathit{\tau }}_t<\mathrm{0}\end{array}\right.,\text{(6)}& {\displaystyle }{\mathit{\tau }}_t& {\displaystyle }=d\cos \left(\mathit{\omega }t\right)+e\sin \left(\mathit{\omega }t\right),\end{array}$$

where τ_t accounts for the different lifetime during winter and summer. The parameter vector for this model is given by ${\mathit{\beta }}_{\mathrm{nonlin}}=(a,b,c,{\mathit{\tau }}_{\mathrm{0}},d,e{)}^{\top }\in {\mathbb{R}}^{\mathrm{6}}$, and we describe how we determine these coefficients and their uncertainties in the next section.

4 Parameter and uncertainty estimation

The parameters are usually estimated by maximizing the likelihood, or, in the case of additional prior constraints, by maximizing the posterior probability. In the linear case and in the case of independently identically distributed Gaussian measurement uncertainties, the maximum likelihood solutions are given by the usual linear least squares solutions. Estimating the parameters in the non-linear case is more involved. Various methods exist, for example, conjugate gradient, random (Monte Carlo) sampling, or exhaustive search methods. The assessment and selection of the method to estimate the parameters in the non-linear case are given below.

4.1 Maximum posterior probability

Because of the complicated structure of the model function in Eq. (3), in particular the lifetime parts in Eqs. (5) and (6), the usual gradient methods converge slowly, if at all. Therefore, we fit the parameters and assess their uncertainty ranges using Markov chain Monte Carlo (MCMC) sampling (Foreman-Mackey et al., 2013). This method samples probability distributions, and we apply it to sample the parameter space, putting emphasis on parameter values with a high posterior probability. The posterior distribution is given in the Bayesian sense as the product of the likelihood and the prior distribution:

$$\begin{array}{}\text{(7)}& p\left({\mathit{x}}_{\mathrm{mod}}\right|\mathit{y})\propto p({\mathit{x}}_{\mathrm{mod}}|\mathit{y},\mathit{\beta })p\left(\mathit{\beta }\right).\end{array}$$

We denote the vector of the measured densities by y and the modelled densities by x_mod, similar to Eqs. (1) and (3). To find the best parameters β for the model, we maximize log p(x_mod|y).

The likelihood $p\left({\mathit{x}}_{\mathrm{mod}}\right|\mathit{y},\mathit{\beta })$ is in our case given by a Gaussian distribution of the residuals, the difference of the model to the data, given in Eq. (8):

$$\begin{array}{ll}{\displaystyle }p\left({\mathit{x}}_{\mathrm{mod}}\right|\mathit{y},\mathit{\beta })& {\displaystyle }=\mathcal{N}\left(\mathit{y},{\mathbf{S}}_y\right)\\ {\displaystyle }& {\displaystyle }=C\exp \left\{-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left(\mathit{y}-{\mathit{x}}_{\mathrm{mod}}\left(\mathit{\beta }\right)\right)}^{\top }\right.\\ \text{(8)}& {\displaystyle }& {\displaystyle }\left.{\mathbf{S}}_y^{-\mathrm{1}}\left(\mathit{y}-{\mathit{x}}_{\mathrm{mod}}\left(\mathit{\beta }\right)\right)\right\}.\end{array}$$

Note that the normalization constant C in Eq. (8) does not influence the value of the maximal likelihood. The covariance matrix S_y contains the squared standard errors of the daily zonal means on the diagonal, ${\mathbf{S}}_y=\mathrm{diag}\left({\mathit{\sigma }}_y^{\mathrm{2}}\right)$.

The prior distribution p(β) restricts the parameters to lie within certain ranges, and the bounds we used for the sampling are listed in Table 1. Within those bounds we assume uniform (flat) prior distributions for the offset, the geomagnetic and solar amplitudes, and in the linear case also for the annual and semi-annual harmonics. We penalize large lifetimes using an exponential distribution $p\left(\mathit{\tau }\right)\propto \exp \mathit{\{}-\mathit{\tau }/{\mathit{\sigma }}_{\mathit{\tau }}\mathit{\}}$ for each lifetime parameter, i.e. for τ₀, d, and e in Eqs. (5) and (6). The scale width σ_τ of this exponential distribution is fixed to 1 day. This choice of prior distributions for the lifetime parameters prevents sampling of the edges of the parameter space at places with small geomagnetic coefficients. In those regions the lifetime may be ambiguous and less meaningful.

Table 1 Parameter search space for the non-linear model and uncertainty estimation.

4.2 Correlations

In the simple case, the measurement covariance matrix S_y contains the measurement uncertainties on the diagonal, in our case the (squared) standard error of the zonal means denoted by σ_y, ${\mathbf{S}}_y=\mathrm{diag}\left({\mathit{\sigma }}_y^{\mathrm{2}}\right)$. However, the standard error of the mean might underestimate the true uncertainties. In addition, possible correlations may occur which are not accounted for using a diagonal S_y.

Both problems can be addressed by adding a covariance kernel K to S_y. Various forms of covariance kernels can be used (Rasmussen and Williams, 2006), depending on the underlying process leading to the measurement or residual uncertainties. Since we have no prior knowledge about the true correlations, we use a commonly chosen kernel of the Matérn 3/2 type (Matérn, 1960; MacKay, 2003; Rasmussen and Williams, 2006). This kernel only depends on the (time) distance between the measurements $t_{ij}=|t_i-t_j|$ and has two parameters, the “strength” σ and correlation length ρ:

$$\begin{array}{}\text{(9)}& K_{ij}={\mathit{\sigma }}^{\mathrm{2}}\left(\mathrm{1}+{\displaystyle \frac{\sqrt{\mathrm{3}}{t}_{ij}}{\mathit{\rho }}}\right)\exp \left\{-{\displaystyle \frac{\sqrt{\mathrm{3}}{t}_{ij}}{\mathit{\rho }}}\right\}.\end{array}$$

Both parameters are estimated together with the model parameter vector β. We found that using the kernel (9) in a covariance matrix S_y with the entries

$$\begin{array}{}\text{(10)}& {S_y}_{ij}=K_{ij}+{\mathit{\delta }}_{ij}{{\mathit{\sigma }}_y}_i^{\mathrm{2}}\end{array}$$

worked best and led to stable and reliable parameter sampling. Note that an additional “white noise” term σ²I could be added to the covariance matrix to account for data uncertainties that are still underestimated. However, this additional white noise term did not improve the convergence, nor did it influence the fitted parameters significantly.

The approximately 3000×3000 covariance matrix of the Gaussian process model for the residuals was evaluated using the Foreman-Mackey et al. (2017a) approximation and the provided Python code (Foreman-Mackey et al., 2017b). For one-dimensional data sets, this approach is computationally faster than the full Cholesky decomposition, which is usually used to invert the covariance matrix S_y. With this approximation, we achieved sensible Monte Carlo sampling times to facilitate evaluating all 18×16 latitude × altitude bins on a small cluster in about 1 day. We used the emcee package (Foreman-Mackey et al., 2013) for the Monte Carlo sampling, set up to use 112 walkers and 800 samples for the initial fit of the parameters, followed by another 800 so-called burn-in samples and 1400 production samples. The full code can be found in Bender (2018a).

5 Results

We demonstrate the parameter estimates using example time series x_NO at 70 km at 65^∘ S, 5^∘ N, and 65^∘ N. NO shows different behaviour in these regions, showing the most variation with respect to the solar cycle and geomagnetic activity at high latitudes. In contrast, at low latitudes the geomagnetic influence should be reduced (Barth et al., 1988; Hendrickx et al., 2017; Kiviranta et al., 2018). We briefly only show the results for the linear model and point out some of its shortcomings. Thereafter we show the results from the non-linear model and continue to use that for further analysis of the coefficients.

5.1 Time series fits

The fitted densities of the linear model Eq. (1) compared to the data are shown in Fig. 1a, b, and c for the three example latitude bins (65^∘ S, 5^∘ N, 65^∘ N) at 70 km. The linear model works well at high southern and low latitudes. At high northern latitudes and to a lesser extent at high southern latitudes, the linear model captures the summer NO variations well. However, the model underestimates the high values in the polar winter at active times (2004–2007) and overestimates the low winter values at quiet times (2009–2011).

Figure 1 Time series data and linear model values and residuals at 70 km for 65^∘ S (a, d), 5^∘ N (b, e), and 65^∘ N (c, f). Panels (a)–(c) show the data (black dots with 2σ error bars) and the model values (blue line). Panels (d)–(f) show the residuals as black dots with 2σ error bars.

For the sample time series (65^∘ S, 5^∘ N, 65^∘ N at 70 km), the fits using the non-linear model Eq. (3) are shown in Fig. 2a, b, and c. The non-linear model better captures both the summer NO variations as well as the high values in the winter, especially at high northern latitudes. However, at times of high solar activity (2003–2006) and in particular at times of a strongly disturbed mesosphere (2004, 2006, 2012), the residuals are still significant. At high southern and low latitudes, the improvement over the linear model is less evident. At low latitudes, the NO content is apparently mostly related to the 11-year solar cycle, and the particle influence is suppressed. Since this cycle is covered by the Lyman-α index, both models perform similarly, but the non-linear version has one less parameter. In both regions the residuals show traces of seasonal variations that are not related to particle effects. The linear model appears to capture these variations better than the non-linear model. However, by objective measures including the number of model parameters¹, the non-linear version fits the data better in all bins (not shown here). At high southern latitudes, the SCIAMACHY data are less densely sampled compared to high northern latitudes (see Bender et al., 2017b). In addition to the sampling differences, geomagnetic latitudes encompass a wider geographic range in the Southern Hemisphere (SH) than in the Northern Hemisphere (NH), and the AE index is derived from stations in the NH. Both effects can lower the NO concentrations that SCIAMACHY observes in the SH, particularly at the winter maxima. The lifetime variation that improves the fit in the NH is thus less effective in the SH.

Figure 2 Same as Fig. 1 but for the non-linear model.

5.2 Parameter morphologies

Using the non-linear model, we show the latitude–altitude distributions of the medians of the sampled Lyman-α and geomagnetic index coefficients in Fig. 3. The white regions indicate values outside of the 95 % confidence region or whose sampled distribution has a skewness larger than 0.33. The MCMC method samples the parameter probability distributions. Since we require the geomagnetic index and constant lifetime parameters to be larger than zero (see Table 1), these sampled distributions are sometimes skewed towards zero, even though the 95 % credible region is still larger than zero. Excluding heavily skewed distributions avoids those cases because the “true” parameter is apparently zero.

Figure 3 Latitude–altitude distributions of the fitted solar index parameter (Lyman-α, a) and the geomagnetic index parameter (AE, b) from the non-linear model.

The Lyman-α parameter distribution shows that its largest influence is at middle and low latitudes between 65 and 80 km. Another increase of the Lyman-α coefficient is indicated at higher altitudes above 90 km. The penetration of Lyman-α radiation decreases with decreasing altitude as a result of scattering and absorption by air molecules. On the other hand, the concentration of air decreases with altitude. At this stage we do not have an unambiguous explanation of this behaviour, but it may be related to reaction pathways as laid out by Pendleton et al. (1983), which would relate the NO concentrations to the CO₂ and H₂O (or OH, respectively) profiles. The Lyman-α coefficients are all negative below 65 km. We also observe negative values at high northern latitude at all altitudes and at high southern latitudes above 85 km. These negative coefficients indicate that NO photodissociation or conversion to other species outweighs its production via UV radiation in those places. The north–south asymmetry may be related to sampling and the difference in illumination with respect to geomagnetic latitudes; see Sect. 5.1.

The geomagnetic influence is largest at high latitudes between 50 and 75^∘ above about 65 km. The AE coefficients peak at around 72 km and indicate a further increase above 90 km. This pattern of the geomagnetic influence matches the one found in Sinnhuber et al. (2016). Unfortunately both increased influences above 90 km in Lyman-α and AE cannot be studied at higher latitudes due to a large a priori contribution to the data.

The latitude–altitude distributions of the lifetime parameters are shown in Fig. 4. All values shown are within the 95 % confidence region. As for the coefficients above, we also exclude regions where the skewness was larger than 0.33.

Figure 4 Latitude–altitude distributions of the fitted base lifetime τ₀ (a) and the amplitude of the annual variation $\left|{\mathit{\tau }}_t\right|$ (b) from the non-linear model.

The constant part of the lifetime, τ₀, is below 2 days in most bins, except for exceptionally large values (>10 days) at low latitudes (0–20^∘ N) between 68 and 74 km. Although we constrained the lifetime with an exponential prior distribution, these large values apparently resulted in a better fit to the data. One explanation could be that because of the small geomagnetic influence (the AE coefficient is small in this region), the lifetime is more or less irrelevant. The amplitude of the annual variation ($\left|{\mathit{\tau }}_{\mathrm{t}}\right|=\sqrt{{\mathit{\tau }}_{\cos }^{\mathrm{2}}+{\mathit{\tau }}_{\sin }^{\mathrm{2}}}=\sqrt{d^{\mathrm{2}}+e^{\mathrm{2}}}$; see Eq. 6) is largest at high latitudes in the Northern Hemisphere and at middle latitudes in the Southern Hemisphere. This difference could be linked to the geomagnetic latitudes which include a wider range of geographic latitudes in the Southern Hemisphere compared to the Northern Hemisphere. Therefore, the annual variation is less apparent in the Southern Hemisphere. The amplitude also increases with decreasing altitude below 75 km at middle and high latitudes and with increasing altitude above that. The increasing annual variation at low altitudes can be the result of transport processes that are not explicitly treated in our approach. Note that the term lifetime is not a pure (photo)chemical lifetime; rather it indicates how long the AE signal persists in the NO densities. In that sense it combines the (photo)chemical lifetime with transport effects as discussed in Sinnhuber et al. (2016).

5.3 Parameter profiles

For three selected latitude bins in the Northern Hemisphere (5, 35, and 65^∘ N) we present profiles of the fitted parameters in Fig. 5. The solid line indicates the median, and the error bars indicate the 95 % confidence region. As indicated in Fig. 3, the solar radiation influence is largest between 65 and 80 km. Its influence is also up to a factor of 2 larger at low and middle latitudes compared to high latitudes, where the coefficient only differs significantly from zero below 65 and above 82 km. Similarly, the geomagnetic impact decreases with decreasing latitude by 1 order of magnitude from high to middle latitudes and at least a further factor of 5 to lower latitudes. The largest impact is around 70–72 km and possibly above 90 km at high latitudes and is approximately constant between 66 and 76 km at middle and low latitudes. Note that the scale in Fig. 5b is logarithmic. The lifetime variation shows that at high latitudes, geomagnetically affected NO persists longer during winter (the phase is close to zero for all altitudes at 65^∘ N, not shown here). It persists up to 10 days longer between 85 and 70 km and increasingly longer below, reaching 28 days at 60 km.

Figure 5 Coefficient profiles of the solar index parameter (Lyman-α, left, a), the geomagnetic index parameter (AE, middle, b), and the amplitude of the annual variation of the NO lifetime (right, c) at 5 (green), 35 (orange), and 65^∘ N (blue). The solid line indicates the median, and the error bars indicate the 95 % confidence region.

For the same latitude bins in the Southern Hemisphere (5, 35, and 65^∘ S) we present profiles of the fitted parameters in Fig. 6. Similar to the coefficients in the Northern Hemisphere (see Fig. 5), the solar radiation influence is largest between 65 and 80 km and also up to a factor of 2 larger at low and middle latitudes compared to high latitudes. However, the Lyman-α coefficients at 65^∘ S are significant below 82 km. Also the geomagnetic AE coefficients show a similar pattern in the Southern Hemisphere compared to the Northern Hemisphere, decreasing by orders of magnitude from high to low latitudes. Note that the AE coefficients at high latitudes are slightly lower than in the Northern Hemisphere, whereas the coefficients at middle and low latitudes are slightly larger. This slight asymmetry was also found in the study by Sinnhuber et al. (2016) and may be related to AE being derived solely from stations in the Northern Hemisphere (Mandea and Korte, 2011). With respect to latitude, the annual variation of the lifetime seems to be reversed compared to the Northern Hemisphere, with almost no variation at high latitudes and longer persisting NO at low latitudes. A faster descent in the southern polar vortex may be responsible for the short lifetime at high southern latitudes. Another reason may be the mixture of air from the inside and outside of the polar vertex when averaging along geomagnetic latitudes since the 65^∘ S geomagnetic latitude band includes geographic locations from about 45 to 85^∘ S. A third possibility may be the exclusion of the South Atlantic Anomaly from the retrieval (Bender et al., 2013, 2017b), for which presumably the particle-induced impact on NO is largest.

Figure 6 Coefficient profiles of the solar index parameter (Lyman-α, left, a), the geomagnetic index parameter (AE, middle, b), and the amplitude of the annual variation of the NO lifetime (right, c) at 5 (green), 35 (orange), and 65^∘ S (blue). The solid line indicates the median, and the error bars indicate the 95 % confidence region.

5.4 Discussion

The distribution of the parameters confirms our understanding of the processes producing NO in the mesosphere to a large extent. The Lyman-α coefficients are related to radiative processes such as production by UV or soft X-rays, either directly or via intermediary of photoelectrons. The photons are not influenced by Earth's magnetic field, and the influence of these processes is largest at low latitudes and decreases towards higher latitudes. We observe negative Lyman-α coefficients below 65 km at all latitudes and at high northern latitudes above 80 km. These negative Lyman-α coefficients indicate that at high solar activity, photodissociation by λ<191 nm photons, photoionization by λ<134 nm photons, or collisional loss and conversion to other species outweigh the production from higher energy photons (<40 nm). At high southern latitudes these negative Lyman-α coefficients are not as pronounced as at high northern latitudes. As mentioned in Sect. 5.2, this north–south asymmetry may be related to sampling and the difference in illumination with respect to geomagnetic latitudes; see also Sect. 5.1.

The AE coefficients are largest at auroral latitudes, as expected for the particle nature of the associated NO production. The AE coefficient can be considered an effective production rate modulated by all short-term (≪1 day) processes. To roughly estimate this production rate, we divided the coefficient of the (daily) AE by 86 400 s, which follows the approach in Sinnhuber et al. (2016). We find a maximum production rate of about 1 cm⁻³ nT⁻¹ s⁻¹ around 70–72 km. This production rate also agrees with the one estimated by Sinnhuber et al. (2016) through a superposed epoch analysis of summertime NO. Comparing the NO production to the ionization rates from Verronen et al. (2013) from 1 to 3 January 2005 (assuming approximately 1 NO molecule per ion pair), our model overestimates the ionization derived from AE on these days. The AE values of 105, 355, and 435 nT translate to 105, 355, and 435 NO molecules cm⁻³ s⁻¹, about 4 times larger than would be estimated from the ionization rates in Verronen et al. (2013), but agreeing with Sinnhuber et al. (2016). The factor of 4 may be related to the slightly different locations, around 60^∘ N (Verronen et al., 2013) compared to around 65^∘ N here and in Sinnhuber et al. (2016), in which the ionization rates may be higher.

The associated constant part τ₀ of the lifetime ranges from around 1 to around 4 days, except for large τ₀ at low latitudes around 70 km. As already discussed in Sect. 5.2, these large lifetimes may be a side effect of the small geomagnetic coefficients and more or less arbitrary. The magnitude is similar to what was found in the study by Sinnhuber et al. (2016) using only the summer data.

The annual variation of the lifetime is largest at high northern latitudes with a nearly constant amplitude of 10 days between 70 and 85 km. An empirical lifetime of 10 days in winter was used by Sinnhuber et al. (2016) to extend the NO predicted by the summer analysis to the larger values in winter. Here we could confirm that 10 days is a good approximation of the NO lifetime in winter, but it varies with altitude. The altitude distribution agrees with the increasing photochemical lifetime at large solar zenith angles (Sinnhuber et al., 2016, Fig. 7b). The larger values in our study are similarly related to transport and mixing effects which alter the observed lifetime. The small variation of the lifetime at high southern latitudes could be a sampling issue because SCIAMACHY only observes small variations there in winter (see Figs. 1 and 2). Note that the results (in particular the large annual variation) in the northernmost latitude bin should be taken with caution because this bin is sparsely sampled by SCIAMACHY, and the large winter NO concentrations are actually absent from the data.

6 Conclusions

We propose an empirical model to estimate the NO density in the mesosphere (60–90 km) derived from measurements from SCIAMACHY nominal-mode limb scans. Our model calculates NO number densities for geomagnetic latitudes using the solar Lyman-α index and the geomagnetic AE index. Two approaches were tested, a linear approach containing annual and semi-annual harmonics and a non-linear version using a finite and variable lifetime for the geomagnetically induced variations. From our proposed models, the linear variant only describes part of the NO variations. It can describe the summer variations but underestimates the large number densities during winter. The non-linear version derived from the SCIAMACHY NO data describes both variations using an annually varying finite lifetime for the particle-induced NO. However, in cases of dynamic disturbances of the mesosphere, for example, in early 2004 or in early 2006, the indirectly enhanced NO (see, for example, Randall et al., 2007) is not captured by the model. These remaining variations are treated as statistical noise.

Sinnhuber et al. (2016) use a superposed epoch analysis limited to the polar summer to model the NO data. Here we extend that analysis to use all available SCIAMACHY nominal-mode NO data for all seasons. However, during summer the present results show comparable NO production per AE unit and similar lifetimes to the Sinnhuber et al. (2016) study.

The parameter distributions indicate in which regions the different processes are significant. We find that these distributions match our current understanding of the processes producing and depleting NO in the mesosphere (Funke et al., 2014a, b, 2016; Sinnhuber et al., 2016; Hendrickx et al., 2017; Kiviranta et al., 2018). In particular, the influence of Lyman-α (or solar UV radiation in general) is largest at low and middle latitudes, which is explained by the direct production of NO via solar UV or soft X-ray radiation (Barth et al., 1988, 2003). The geomagnetic influence is largest at high latitudes and is best explained by the production from charged particles that enter the atmosphere in the polar regions along the magnetic field.

A potential improvement would be to use actual measurements of precipitating particles instead of the AE index. Using measured fluxes could help to confirm our current understanding of how those fluxes relate to ionization (Turunen et al., 2009; Verronen et al., 2013) and subsequent NO production (Sinnhuber et al., 2016). Furthermore, including dynamical transport, as for example in Funke et al. (2016), could improve our knowledge of the combined direct and indirect NO production in the mesosphere.