Mesospheric nitric oxide model from SCIAMACHY data
Stefan Bender, Miriam Sinnhuber, Patrick J. Espy, John P. Burrows
MMesospheric nitric oxide model from SCIAMACHY data
Stefan Bender , Miriam Sinnhuber , Patrick J. Espy , and John P. Burrows Department of Physics, Norwegian University of Science and Technology, Trondheim, Norway Institute of Meteorology and Climate Research, Karlsruhe Institute of Technology, Karlsruhe, Germany Institute of Environmental Physics, University of Bremen, Bremen, Germany
Correspondence:
Stefan Bender ([email protected])
Abstract.
We present an empirical model for nitric oxide( NO ) in the mesosphere ( ≈ –90 km) derived from SCIA-MACHY (SCanning Imaging Absorption spectroMeter forAtmospheric CHartoghraphY) limb scan data. This workcomplements and extends the NOEM (Nitric Oxide Empir-ical Model; Marsh et al., 2004) and SANOMA (SMR Ac-quired Nitric Oxide Model Atmosphere; Kiviranta et al.,2018) empirical models in the lower thermosphere. The re-gression ansatz builds on the heritage of studies by Hen-drickx et al. (2017) and the superposed epoch analysisby Sinnhuber et al. (2016) which estimate NO productionfrom 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. Weuse a non-linear regression model, incorporating a finite andseasonally varying lifetime for the geomagnetically induced NO . We estimate the parameters by finding the maximumposterior probability and calculate the parameter uncertain-ties using Markov chain Monte Carlo sampling. In additionto providing an estimate of the NO content in the meso-sphere, the regression coefficients indicate regions where cer-tain processes dominate. It has been recognized in the past decades that the meso-sphere and stratosphere are coupled in various ways (Bald-win and Dunkerton, 2001). Consequently, climate modelshave been evolving to extend to increasingly higher levelsin the atmosphere to improve the accuracy of medium- andlong-term predictions. Nowadays it is not unusual that thesemodels include the mesosphere (40–90 km) or the lower ther- mosphere (90–120 km) (Matthes et al., 2017). It is there-fore important to understand the processes in the mesosphereand lower thermosphere and to find the important driversof chemistry and dynamics in that region. The atmosphereabove the stratosphere ( (cid:38) km) is coupled to solar and ge-omagnetic activity, also known as space weather (Sinnhuberet al., 2012). Electrons and protons from the solar wind andthe radiation belts with sufficient kinetic energy enter the at-mosphere in that region. Since as charged particles they movealong the magnetic field, this precipitation occurs primarilyat high geomagnetic latitudes.Previously the role of NO in the mesosphere has beenidentified as an important free radical, and in this sense adriver of the chemistry (Kockarts, 1980; Barth, 1992, 1995;Roble, 1995; Bailey et al., 2002; Barth et al., 2009; Barth,2010), particularly during winter when it is long-lived be-cause of reduced photodissociation. NO generated in the re-gion between 90 and 120 km at auroral latitudes is stronglyinfluenced by both solar and geomagnetic activity (Marshet al., 2004; Sinnhuber et al., 2011, 2016; Hendrickx et al.,2015, 2017). At high latitudes, NO is transported down tothe upper stratosphere during winter, usually down to 50 kmand occasionally down to 30 km (Siskind et al., 2000; Ran-dall et al., 2007; Funke et al., 2005a, 2014b). At those al-titudes and also in the mesosphere, NO participates in the“odd oxygen catalytic cycle which depletes ozone” (Crutzen,1970). Additional dynamical processes also result in thestrong downward transport of mesospheric air into the up-per stratosphere, such as the strong downwelling that oftenoccurs in the recovery phase of a sudden stratospheric warm-ing (SSW) (Pérot et al., 2014; Orsolini et al., 2017). Thisdownwelling is typically associated with the formation of anelevated stratopause.Different instruments have been measuring NO in themesosphere and lower thermosphere, but at different alti- a r X i v : . [ phy s i c s . a o - ph ] F e b S. Bender et al.: SCIAMACHY mesosphere NO model tudes and at different local times. Measurements from so-lar occultation instruments such as Scisat-1/ACE-FTS orAIM/SOFIE are limited in latitude and local time (sunriseand sunset). Global observations from sun-synchronouslyorbiting satellites are available from Envisat/MIPAS be-low 70 km daily and 42–172 km every 10 days (Funkeet al., 2001, 2005b; Bermejo-Pantaleón et al., 2011); fromOdin/SMR between 45 and 115 km (Pérot et al., 2014; Kivi-ranta et al., 2018); or from Envisat/SCIAMACHY (SCanningImaging Absorption spectroMeter for Atmospheric CHar-toghraphY) 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 and18:00 (Odin) and 10:00 and 22:00 (Envisat). While MI-PAS has both daytime and night-time measurements, SCIA-MACHY provides daytime (10:00) data because of the mea-surement principle (fluorescent UV scattering; see Benderet al., 2013, 2017b). Unfortunately, Envisat stopped commu-nicating in April 2012, and therefore the data available fromMIPAS and SCIAMACHY are limited to nearly 10 yearsfrom August 2002 to April 2012. The other aforementionedinstruments are still operational and provide ongoing data aslong as satellite operations continue.Chemistry–climate models struggle to simulate the NO amounts and distributions in the mesosphere and lower ther-mosphere (see, for example, Funke et al., 2017; Randallet al., 2015; Orsolini et al., 2017; Hendrickx et al., 2018).To remedy the situation, some models constrain the NO con-tent at their top layer through observation-based parametriza-tions. For example, the next generation of climate simula-tions (CMIP6; see Matthes et al., 2017) and other recentmodel simulations (Sinnhuber et al., 2018) parametrize par-ticle effects as derived partly from Envisat/MIPAS NO mea-surements (Funke et al., 2016). NO in the mesosphere and lower thermosphere NO in the mesosphere and lower thermosphere is producedby N dissociation, N + hν → N( D) + N( S) ( λ < nm ) , (R1)followed by the reaction of the excited nitrogen atom N( D) with molecular oxygen (Solomon et al., 1982; Barth, 1992,1995): N( D) + O → NO + O . (R2)The dissociation energy of N into ground state atoms N( S) is about 9.8 eV ( λ ≈ nm) (Hendrie, 1954; Frost et al.,1956; Heays et al., 2017). This energy together with the ex-citation energy to N( D) is denoted by hν in Reaction (R1)and can be provided by a number of sources, most notablyby auroral or photoelectrons as well as by soft solar X-rays.The NO content is reduced by photodissociation, NO + hν → N + O ( λ < nm ) , (R3) by photoionization NO + hν → NO + + e − ( λ < nm ) , (R4)and by reacting with atomic nitrogen: NO + N → N + O . (R5) N O has been retrieved in the mesosphere and thermo-sphere 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 thesource of this N O to being most likely the reaction between NO and N atoms produced by particle precipitation: N + NO → N O + O . (R6)We note that photo-excitation and photolysis at 185 nm (vac-uum UV) of NO or NO mixtures in nitrogen, N , or he-lium mixtures at 1 atm leads to N O formation (Maric andBurrows, 1992). Both mechanisms explaining the produc-tion of N O involve excited states of NO . Hence these path-ways contribute to the loss of NO and potentially an addi-tional daytime source of N O in the upper atmosphere. N O acts as an intermediate reservoir at high altitudes ( (cid:39) km;see Sheese et al., 2016), reacting with O ( D) in two well-known channels to N and O as well as to 2 NO . How-ever, the largest N O abundances are located below 60 kmand originate primarily from the transport of tropospheric N O into the stratosphere through the Brewer–Dobson cir-culation (Funke et al., 2008a, b; Sheese et al., 2016) but canreach up to 70 km in geomagnetic storm conditions (Funkeet al., 2008a; Sheese et al., 2016). Both source and sink re-actions indicate that NO behaves differently in sunlit condi-tions than in dark conditions. NO is produced by particle pre-cipitation at auroral latitudes, but in dark conditions (withoutphotolysis) it is only depleted by reacting with atomic ni-trogen (Reaction (R5)). This asymmetry between productionand depletion in dark conditions results in different lifetimesof NO .Early work to parametrize NO in the lower thermosphere(100–150 km) used SNOE measurements from March 1998to September 2000 (Marsh et al., 2004). With these 2.5years of data and using empirical orthogonal functions, theso-called NOEM (Nitric Oxide Empirical Model) estimates NO in the lower thermosphere as a function of the so-lar f . cm radio flux, the solar declination angle, and theplanetary Kp index. NOEM is still used as prior input for NO retrieval, for example from MIPAS (Bermejo-Pantaleónet al., 2011; Funke et al., 2012) and SCIAMACHY (Ben-der et al., 2017b) spectra. However, 2.5 years is relativelyshort compared to the 11-year solar cycle, and the years1998 to 2000 encompass a period of elevated solar activ-ity. To address this, a longer time series from AIM/SOFIEwas used to determine the important drivers of NO inthe lower thermosphere (90–140 km) by Hendrickx et al. . Bender et al.: SCIAMACHY mesosphere NO model 3 (2017). Other recent work uses 10 years of NO data fromOdin/SMR from 85 to 115 km (Kiviranta et al., 2018). Funkeet al. (2016) derived a semi-empirical model of NO y inthe stratosphere and mesosphere from MIPAS data. Herewe use Envisat/SCIAMACHY NO data from the nominallimb mode (Bender et al., 2017b, a). Apart from providinga similarly long time series of NO data, the nominal En-visat/SCIAMACHY NO data cover the mesosphere from 60to 90 km (Bender et al., 2017b), bridging the gap between thestratosphere and lower thermosphere models.The paper is organized as follows: we present the data usedin this work in Sect. 2. The two model variants, linear andnon-linear, are described in Sect. 3. Details about the param-eter and uncertainty estimation are explained in Sect. 4, andwe present the results in Sect. 5. Finally we conclude ourfindings in Sect. 6. NO We use the SCIAMACHY nitric oxide data set ver-sion 6.2.1 (Bender et al., 2017a) retrieved from the nominallimb scan mode ( ≈ –93 km). For a detailed instrument de-scription, see Burrows et al. (1995) and Bovensmann et al.(1999), and for details of the retrieval algorithm, see Benderet al. (2013, 2017b).The data were retrieved for the whole Envisat period (Au-gust 2002–April 2012). This satellite was orbiting in a sun-synchronous orbit at around 800 km altitude, with Equatorcrossing times of 10:00 and 22:00 local time. The NO num-ber densities from the SCIAMACHY nominal mode wereretrieved from the NO gamma band emissions. Since thoseemissions are fluorescent emissions excited by solar UV,SCIAMACHY NO data are only available for the 10:00 day-side (downleg) part of the orbit. Furthermore, the retrievalwas carried out for altitudes from 60 to 160 km, but aboveapproximately 90 km, the data reflect the scaled a priori den-sities from NOEM (Bender et al., 2017b). We therefore re-strict the modelling to the mesosphere below 90 km.We averaged the individual orbital data longitudinally ona daily basis according to their geomagnetic latitude within10 ◦ bins. The geomagnetic latitude was determined accord-ing to the eccentric dipole approximation of the 12th gen-eration of the International Geomagnetic Reference Field(IGRF12) (Thébault et al., 2015). In the vertical direction theoriginal retrieval grid altitudes (2 km bins) were used. Notethat mesospheric NO concentrations are related to geomag-netically as well as geographically based processes, but dis-entangling them is beyond the scope of the paper. Follow-upstudies can build on the method presented here and study, forexample, longitudinally resolved timeseries.The measurement sensitivity is taken into account viathe averaging kernel diagonal elements, and days where its binned average was below 0.002 were excluded from thetimeseries. Considering this criterion, each bin (geomagneticlatitude and altitude) contains about 3400 data points. We use two proxies to model the NO number densi-ties, one accounting for the solar irradiance variations andone accounting for the geomagnetic activity. Various prox-ies 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 cy-cle. The NOEM (Nitric Oxide Empirical Model; Marsh et al.,2004) uses the natural logarithm of the solar 10.7 cm ra-dio flux f . . More recent work on AIM/SOFIE NO (Hen-drickx et al., 2017) uses the solar Lyman- α index becausesome of the main production and loss processes are driven byUV photons. Besides accounting for the long-term variationof NO with solar activity, the Lyman- α index also includesshort-term UV variations and the associated NO production,for example caused by solar flares. Barth et al. (1988) haveshown that the Lyman- α index directly relates to the ob-served NO at low latitudes (30 ◦ S–30 ◦ N). Thus we use it inthis work as a proxy for NO .In the same manner as for the irradiance variations, the“right” geomagnetic index to model particle-induced vari-ations of NO is a matter of opinion. Kp is the oldest andmost commonly used geomagnetic index; it was, for exam-ple, used in earlier work by Marsh et al. (2004) for mod-elling NO in the mesosphere and lower thermosphere. Kpis derived from magnetometer stations distributed at differ-ent latitudes and mostly in the Northern Hemisphere (NH).However, Hendrickx et al. (2015) found that the auroral elec-trojet index (AE) (Davis and Sugiura, 1966) correlated bet-ter with SOFIE-derived NO concentrations (Hendrickx et al.,2015, 2017); (see also Sinnhuber et al., 2016). The AE indexis derived from stations distributed almost evenly within theauroral latitude band. This distribution enables the AE in-dex to be more closely related to the energy input into theatmosphere at these latitudes. Therefore, we use the auroralelectrojet index (AE) as a proxy for geomagnetically induced NO . To account for the 10:00 satellite sampling, we averagethe hourly AE index from noon the day before to noon on themeasurement day.It should be noted that tests using Kp (or its linear equiv-alent Ap) instead of AE and using f . instead of Lyman- α suggested that the particular choice of index did not lead tosignificantly different results. Our choice of AE rather thanKp and Lyman- α over f . is physically based and moti-vated as described above. S. Bender et al.: SCIAMACHY mesosphere NO model3 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 thetime (measurement day) t : x NO ( φ, z, t ) . In the following weoften drop the subscript NO and combine the time directioninto a vector x with the i th entry denoting the density at time t i , such that x i ( φ, z ) = x ( φ, z, t i ) . In the (multi-)linear case, we relate the nitric oxide numberdensities x NO ( φ, z, t ) to the two proxies, the solar Lyman- α index (Ly α ( t )) and the geomagnetic AE index (AE( t )). Har-monic terms with ω = 2 π a − = 2 π (365 . d ) − account forannual and semi-annual variations. The linear model, includ-ing a constant offset for the background density, describesthe NO density according to Eq. (1): x NO ( φ, z, t ) = a ( φ, z ) + b ( φ, z ) · Ly α ( t ) + c ( φ, z ) · AE ( t )+ (cid:88) n =1 [ d n ( φ, z ) cos( nωt )+ e n ( φ, z ) sin( nωt )] . (1)The linear model can be written in matrix form for the n measurement times t , . . . , t n as Eq. (2), with the parametervector β given by β lin = ( a, b, c, d , e , d , e ) (cid:62) ∈ R and themodel matrix X ∈ R n × . x NO ( φ, z ) = Ly α ( t ) AE ( t ) cos( ωt ) sin( ωt ) ... Ly α ( t n ) AE ( t n ) cos( ωt n ) sin( ωt n )cos(2 ωt ) sin(2 ωt ) ... cos(2 ωt n ) sin(2 ωt n ) · abcd e d e = X · β (2)We determine the coefficients via least squares, minimizingthe squared differences of the modelled number densities tothe measured ones. In contrast to the linear model above, we modify the AE in-dex by a finite lifetime τ which varies according to season,we denote this modified version by (cid:102) AE. We then omit the har-monic parts in the model, and the non-linear model is given by Eq. (3): x NO ( φ, z, t ) = a ( φ, z )+ b ( φ, z ) · Ly α ( t )+ c ( φ, z ) · (cid:102) AE ( t ) . (3)Although this approach shifts all seasonal variations to theAE index and thus attributes them to particle-induced effects,we found that the residual traces of particle-unrelated sea-sonal effects were minor compared to the overall improve-ment of the fit. Additional harmonic terms only increase thenumber of free parameters without substantially improvingthe fit further.The lifetime-corrected (cid:102) AE is given by the sum of the pre-vious 60 days’ AE values, each multiplied by an exponentialdecay factor: (cid:102) AE ( t ) = d (cid:88) t i =0 AE ( t − t i ) · exp (cid:26) − t i τ (cid:27) . (4)The total lifetime τ is given by a constant part τ plus thenon-negative fraction of a seasonally varying part τ t : τ = τ + (cid:40) τ t , τ t ≥ , τ t < , (5) τ t = d cos( ωt ) + e sin( ωt ) , (6)where τ t accounts for the different lifetime during winter andsummer. The parameter vector for this model is given by β nonlin = ( a, b, c, τ , d, e ) (cid:62) ∈ R , and we describe how we de-termine these coefficients and their uncertainties in the nextsection. The parameters are usually estimated by maximizing thelikelihood, or, in the case of additional prior constraints, bymaximizing the posterior probability. In the linear case andin the case of independently identically distributed Gaussianmeasurement uncertainties, the maximum likelihood solu-tions are given by the usual linear least squares solutions.Estimating the parameters in the non-linear case is more in-volved. Various methods exist, for example conjugate gra-dient, random (Monte Carlo) sampling or exhaustive searchmethods. The assessment and selection of the method to es-timate the parameters in the non-linear case are given below.
Because of the complicated structure of the model func-tion in Eq. (3), in particular the lifetime parts in Eqs. (5)and (6), the usual gradient methods converge slowly, if atall. Therefore, we fit the parameters and assess their un-certainty ranges using Markov chain Monte Carlo (MCMC)sampling (Foreman-Mackey et al., 2013). This method sam-ples probability distributions, and we apply it to sample theparameter space putting emphasis on parameter values with a . Bender et al.: SCIAMACHY mesosphere NO model 5
Table 1.
Parameter search space for the non-linear model and un-certainty estimation.Parameter Lower bound Upper bound PriorformOffset ( a ) − cm − cm − flatLyman- α amplitude ( b ) − cm − cm − flatAE amplitude ( c ) cm − cm − flat τ d d exp τ cosine amplitude ( d ) − d d exp τ sine amplitude ( e ) − d d exp high posterior probability. The posterior distribution is givenin the Bayesian sense as the product of the likelihood and theprior distribution: p ( x mod | y ) ∝ p ( x mod | y , β ) p ( β ) . (7)We denote the vector of the measured densities by y andthe 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 ( x mod | y , β ) is in our case given by aGaussian distribution of the residuals, the difference of themodel to the data, given in Eq. (8): p ( x mod | y , β ) = N ( y , S y )= C exp (cid:26) −
12 ( y − x mod ( β )) (cid:62) S − y ( y − x mod ( β )) (cid:27) . (8)Note that the normalization constant C in Eq. (8) does not in-fluence the value of the maximal likelihood. The covariancematrix S y contains the squared standard errors of the dailyzonal means on the diagonal, S y = diag ( σ y ) .The prior distribution p ( β ) restricts the parameters to liewithin certain ranges, and the bounds we used for the sam-pling are listed in Table 1. Within those bounds we assumeuniform (flat) prior distributions for the offset, the geomag-netic and solar amplitudes, and in the linear case also for theannual and semi-annual harmonics. We penalize large life-times using an exponential distribution p ( τ ) ∝ exp {− τ /σ τ } for each lifetime parameter, i.e. for τ , d , and e in Eqs. (5)and (6). The scale width σ τ of this exponential distribution isfixed to 1 day. This choice of prior distributions for the life-time parameters prevents sampling of the edges of the param-eter space at places with small geomagnetic coefficients. Inthose regions the lifetime may be ambiguous and less mean-ingful. 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 de-noted by σ y , S y = diag ( σ y ) . However, the standard error ofthe mean might underestimate the true uncertainties. In ad-dition, possible correlations may occur which are not ac-counted for using a diagonal S y .Both problems can be addressed by adding a covariancekernel K to S y . Various forms of covariance kernels canbe used (Rasmussen and Williams, 2006), depending on theunderlying process leading to the measurement or residualuncertainties. Since we have no prior knowledge about thetrue correlations, we use a commonly chosen kernel of theMatérn-3/2 type (Matérn, 1960; MacKay, 2003; Rasmussenand Williams, 2006). This kernel only depends on the (time)distance between the measurements t ij = | t i − t j | and hastwo parameters, the “strength” σ and correlation length ρ : K ij = σ (cid:32) √ t ij ρ (cid:33) exp (cid:40) − √ t ij ρ (cid:41) . (9)Both parameters are estimated together with the model pa-rameter vector β . We found that using the kernel (9) in acovariance matrix S y with the entries S yij = K ij + δ ij σ y i , (10)worked best and led to stable and reliable parameter sam-pling. Note that an additional “white noise” term σ couldbe added to the covariance matrix to account for still under-estimated data uncertainties. However, this additional whitenoise term did not improve the convergence, nor did it influ-ence the fitted parameters significantly.The approximately × covariance matrix ofthe Gaussian process model for the residuals was evaluatedusing the Foreman-Mackey et al. (2017) approximationand the provided Python code (Foreman-Mackey et al.,2017). For one-dimensional data sets, this approach iscomputationally faster than the full Cholesky decomposi-tion, which is usually used to invert the covariance matrix S y . With this approximation, we achieved sensible MonteCarlo sampling times to facilitate evaluating all × latitude × altitude bins on a small cluster in about 1 day. Weused the emcee package (Foreman-Mackey et al., 2013)for the Monte Carlo sampling, set up to use 112 walk-ers and 800 samples for the initial fit of the parameters,followed by another 800 so-called burn-in samples and1400 production samples. The full code can be found at https://github.com/st-bender/sciapy (Ben-der, 2018a). We demonstrate the parameter estimates using example timeseries x NO at 70 km at 65 ◦ S, 5 ◦ N, and 65 ◦ N. NO shows dif-ferent behaviour in these regions, showing the most variation S. Bender et al.: SCIAMACHY mesosphere NO model with respect to the solar cycle and geomagnetic activity athigh latitudes. In contrast, at low latitudes the geomagneticinfluence should be reduced (Barth et al., 1988; Hendrickxet al., 2017; Kiviranta et al., 2018). We briefly only show theresults for the linear model and point out some of its short-comings. Thereafter we show the results from the non-linearmodel and continue to use that for further analysis of the co-efficients.
The fitted densities of the linear model Eq. (1) compared tothe data are shown in the upper panels of Fig 1 for the threeexample latitude bins (65 ◦ S, 5 ◦ N, 65 ◦ N) at 70 km. The linearmodel works well at high southern and low latitudes. At highnorthern latitudes and to a lesser extent at high southern lat-itudes, the linear model captures the summer NO variationswell. However, the model underestimates the high values inthe polar winter at active times (2004–2007) and overesti-mates the low winter values at quiet times (2009–2011).For the sample timeseries (65 ◦ S, 5 ◦ N, 65 ◦ N at 70 km), thefits using the non-linear model Eq. (3) are shown in the up-per panels of Fig 2. The non-linear model better capturesboth the summer NO variations as well as the high valuesin the winter, especially at high northern latitudes. However,at times of high solar activity (2003–2006) and in particu-lar at times of a strongly disturbed mesosphere (2004, 2006,2012), the residuals are still significant. At high southernand low latitudes, the improvement over the linear model isless evident. At low latitudes, the NO content is apparentlymostly related to the eleven-year solar cycle and the parti-cle influence is suppressed. Since this cycle is covered bythe Lyman- α index, both models perform similarly, but thenon-linear version has one less parameter. In both regionsthe residuals show traces of seasonal variations that are notrelated to particle effects. The linear model appears to cap-ture these variations better than the non-linear model. How-ever, by objective measures including the number of modelparameters , the non-linear version fits the data better in allbins (not shown here). At high southern latitudes, the SCIA-MACHY data are less densely sampled compared to highnorthern latitudes (see Bender et al., 2017b). In addition to Past and recent research in model selection provides a num-ber of choices on how to compare models objectively. The resultsare so-called information criteria which aim to provide a consis-tent way of how to compare models, most notably the “Akaikeinformation criterion” (AIC; Akaike, 1974), the “Bayesian infor-mation criterion” or “Schwarz criterion” (BIC or SIC; Schwarz,1978), the “deviance information criterion” (DIC; Spiegelhalteret al., 2002; Ando, 2011), or the “widely applicable informationcriterion” (WAIC; Watanabe, 2010; Vehtari et al., 2016). Alterna-tively, the “standardized mean squared error” (SMSE) or the “meanstandardized log loss” (MSLL) (Rasmussen and Williams, 2006,chap. 2) give an impression of the quality of regression models withrespect to each other. the sampling differences, geomagnetic latitudes encompassa wider geographic range in the Southern Hemisphere (SH)than in the Northern Hemisphere (NH), and the AE index isderived 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 thatimproves the fit in the NH is thus less effective in the SH. Using the non-linear model, we show the latitude–altitudedistributions of the medians of the sampled Lyman- α and ge-omagnetic index coefficients in Fig. 3. The white regions in-dicate values outside of the 95% confidence region or whosesampled distribution has a skewness larger than 0.33. TheMCMC method samples the parameter probability distribu-tions. Since we require the geomagnetic index and constantlifetime parameters to be larger than zero (see Table 1), thesesampled distributions are sometimes skewed towards zeroeven though the 95% credible region is still larger than zero.Excluding heavily skewed distributions avoids those casesbecause the “true” parameter is apparently zero.The Lyman- α parameter distribution shows that its largestinfluence is at middle and low latitudes between 65 and80 km. Another increase of the Lyman- α coefficient is in-dicated at higher altitudes above 90 km. The penetration ofLyman- α radiation decreases with decreasing altitude as aresult of scattering and absorption by air molecules. On theother hand, the concentration of air decreases with altitude.At this stage we do not have an unambiguous explanationof this behaviour, but it may be related to reaction pathwaysas laid out by Pendleton et al. (1983), which would relatethe NO concentrations to the CO and H O (or OH , re-spectively) profiles. The Lyman- α coefficients are all neg-ative below 65 km. We also observe negative values at highnorthern latitude at all altitudes and at high southern latitudesabove 85 km. These negative coefficients indicate that NO photodissociation or conversion to other species outweighsits production via UV radiation in those places. The north–south asymmetry may be related to sampling and the differ-ence in illumination with respect to geomagnetic latitudes;see Sect. 5.1.The geomagnetic influence is largest at high latitudes be-tween 50 and 75 ◦ above about 65 km. The AE coefficientspeak at around 72 km and indicate a further increase above90 km. This pattern of the geomagnetic influence matches theone found in Sinnhuber et al. (2016). Unfortunately both in-creased influences above 90 km in Lyman- α and AE cannotbe studied at higher latitudes due to a large a priori contribu-tion to the data.The latitude–altitude distributions of the lifetime parame-ters are shown in Fig. 4. All values shown are within the 95%confidence region. As for the coefficients above, we also ex-clude regions where the skewness was larger than 0.33. . Bender et al.: SCIAMACHY mesosphere NO model 7 N O d e n s i t y [ c m ] (a) datalinear model2004 2006 2008 2010 2012time [year]1000100200300 r e s i d u a l s [ c m ] (d) NO, -65°N, 70 km N O d e n s i t y [ c m ] (b) datalinear model2004 2006 2008 2010 2012time [year]10050050100150 r e s i d u a l s [ c m ] (e) NO, 5°N, 70 km N O d e n s i t y [ c m ] (c) datalinear model2004 2006 2008 2010 2012time [year]2000200400600800 r e s i d u a l s [ c m ] (f) NO, 65°N, 70 km
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) showthe data (black dots with σ error bars) and the model values (blue line). Panels (d) – (f) show the residuals as black dots with σ error bars. N O d e n s i t y [ c m ] (a) datanon-linear model2004 2006 2008 2010 2012time [year]1000100200300 r e s i d u a l s [ c m ] (d) NO, -65°N, 70 km N O d e n s i t y [ c m ] (b) datanon-linear model2004 2006 2008 2010 2012time [year]10050050100150 r e s i d u a l s [ c m ] (e) NO, 5°N, 70 km N O d e n s i t y [ c m ] (c) datanon-linear model2004 2006 2008 2010 2012time [year]2000200400600800 r e s i d u a l s [ c m ] (f) NO, 65°N, 70 km
Figure 2.
Same as Fig. 1 for the non-linear model.
The constant part of the lifetime, τ , is below 2 days inmost bins, except for exceptionally large values ( > days)at low latitudes (0–20 ◦ N) between 68 and 74 km. Althoughwe constrained the lifetime with an exponential prior distri-bution, these large values apparently resulted in a better fit tothe data. One explanation could be that because of the smallgeomagnetic influence (the AE coefficient is small in this re-gion), the lifetime is more or less irrelevant. The amplitudeof the annual variation ( | τ t | = (cid:112) τ cos + τ sin = √ d + e ; seeEq. (6)) is largest at high latitudes in the Northern Hemi-sphere and at middle latitudes in the Southern Hemisphere.This difference could be linked to the geomagnetic lati-tudes which include a wider range of geographic latitudes inthe Southern Hemisphere compared to the Northern Hemi-sphere. Therefore, the annual variation is less apparent in theSouthern Hemisphere. The amplitude also increases with de-creasing altitude below 75 km at middle and high latitudesand with increasing altitude above that. The increasing an-nual variation at low altitudes can be the result of transportprocesses that are not explicitly treated in our approach. Notethat 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 life-time with transport effects as discussed in Sinnhuber et al.(2016). For three selected latitude bins in the Northern Hemisphere(5, 35, and 65 ◦ N) we present profiles of the fitted parametersin Fig. 5. The solid line indicates the median, and the er-ror bars indicate the 95% confidence region. As indicated inFig. 3, the solar radiation influence is largest between 65 and80 km. Its influence is also up to a factor of 2 larger at lowand middle latitudes compared to high latitudes, where thecoefficient only differs significantly from zero below 65 andabove 82 km. Similarly, the geomagnetic impact decreaseswith decreasing latitude by 1 order of magnitude from highto middle latitudes and at least a further factor of 5 to lowerlatitudes. The largest impact is around 70–72 km and possi-bly above 90 km at high latitudes and is approximately con-stant between 66 and 76 km at middle and low latitudes. Notethat the scale in Fig. 5b is logarithmic. The lifetime variationshows that at high latitudes, geomagnetically affected NO persists longer during winter (the phase is close to zero for S. Bender et al.: SCIAMACHY mesosphere NO model − − −
25 0 25 50 75geomagnetic latitude [°N]60657075808590 a l t i t u d e [ k m ] (a) Lyman- α coefficients − − −
25 0 25 50 75geomagnetic latitude [°N]60657075808590 a l t i t u d e [ k m ] (b) AE coefficients − − − N O / L y α [ − c m − ( p h c m − s − ) − ] N O / A E [ c m − n T − ] Figure 3.
Latitude–altitude distributions of the fitted solar index parameter (Lyman- α , a ) and the geomagnetic index parameter (AE, b ) fromthe non-linear model. − − −
25 0 25 50 75geomagnetic latitude [°N]60657075808590 a l t i t u d e [ k m ] (a) AE constant lifetimes − − −
25 0 25 50 75geomagnetic latitude [°N]60657075808590 a l t i t u d e [ k m ] (b) AE lifetime annual amplitudes τ [ d a y s ] | τ t | [ d a y s ] Figure 4.
Latitude–altitude distributions of the fitted base lifetime τ (a) and the amplitude of the annual variation | τ t | (b) from the non-linearmodel. all altitudes at 65 ◦ N, not shown here). It persists up to 10days longer between 85 and 70 km and increasingly longerbelow, reaching 28 days at 60 km.For the same latitude bins in the Southern Hemisphere(5 ◦ S, 35 ◦ S, and 65 ◦ S) we present profiles of the fitted pa-rameters in Fig. 6. Similar to the coefficients in the North-ern Hemisphere (see Fig. 5), the solar radiation influence islargest between 65 km and 80 km and also up to a factor oftwo larger at low and middle latitudes compared to high lat-itudes. However, the Lyman- α coefficients at 65 ◦ S are sig-nificant below 82 km. Also the geomagnetic AE coefficientsshow a similar pattern in the Southern Hemisphere comparedto the Northern Hemisphere, decreasing by orders of magni-tude from high to low latitudes. Note that the AE coefficientsat high latitudes are slightly lower than in the Northern Hemi-sphere, whereas the coefficients at middle and low latitudesare slightly larger. This slight asymmetry was also found inthe study by Sinnhuber et al. (2016) and may be related toAE being derived solely from stations in the Northern Hemi-sphere (Mandea and Korte, 2011). With respect to latitude,the annual variation of the lifetime seems to be reversed com- pared to the Northern Hemisphere, with almost no variationat high latitudes and longer persisting NO at low latitudes. Afaster descent in the southern polar vortex may be responsi-ble for the short lifetime at high southern latitudes. Anotherreason may be the mixture of air from inside and outsideof the polar vertex when averaging along geomagnetic lat-itudes since the 65 ◦ S geomagnetic latitude band includes ge-ographic locations from about 45 ◦ S to 85 ◦ S. A third possibil-ity may be the exclusion of the Southern Atlantic Anomalyfrom the retrieval (Bender et al., 2013, 2017b) where pre-sumably the particle-induced impact on NO is largest. The distribution of the parameters confirms our understand-ing of the processes producing NO in the mesosphere toa large extent. The Lyman- α coefficients are related to ra-diative processes such as production by UV or soft X-rays,either directly or via intermediary of photoelectrons. Thephotons are not influenced by Earth’s magnetic field, andthe influence of these processes is largest at low latitudesand decreases towards higher latitudes. We observe nega- . Bender et al.: SCIAMACHY mesosphere NO model 9 − −
10 0 10 20NO / Ly α [10 − cm − (phcm − s − ) − ]60657075808590 a l t i t u d e [ k m ] (a) Ly α − NO / AE [10 cm − nT − ] (b) AE τ t | [d] (c) | τ t | Figure 5.
Coefficient profiles of the solar index parameter (Lyman- α , left, a ), the geomagnetic index parameter (AE, middle, b ), and theamplitude of the annual variation of the NO lifetime (right, c ) at 5 ◦ N (green), 35 ◦ N (orange), and 65 ◦ N (blue). The solid line indicates themedian, and the error bars indicate the 95% confidence region. − −
10 0 10 20NO / Ly α [10 − cm − (phcm − s − ) − ]60657075808590 a l t i t u d e [ k m ] (a) Ly α − NO / AE [10 cm − nT − ] (b) AE τ t | [d] (c) | τ t | Figure 6.
Coefficient profiles of the solar index parameter (Lyman- α , left, a ), the geomagnetic index parameter (AE, middle, b ), and theamplitude of the annual variation of the NO lifetime (right, c ) at 5 ◦ S (green), 35 ◦ S (orange), and 65 ◦ S (blue). The solid line indicates themedian, and the error bars indicate the 95% confidence region. tive Lyman- α coefficients below 65 km at all latitudes and athigh northern latitudes above 80 km. These negative Lyman- α coefficients indicate that at high solar activity, photodis-sociation by λ < nm photons, photoionization by λ < nm photons, or collisional loss and conversion to otherspecies outweigh the production from higher energy photons( < nm). At high southern latitudes these negative Lyman- α coefficients are not as pronounced as at high northern lati-tudes. As mentioned in Sect. 5.2, this north–south asymmetrymay be related to sampling and the difference in illuminationwith respect to geomagnetic latitudes, see also Sect. 5.1.The AE coefficients are largest at auroral latitudes as ex-pected for the particle nature of the associated NO pro-duction. The AE coefficient can be considered an effectiveproduction rate modulated by all short-term ( (cid:28) day) pro-cesses. To roughly estimate this production rate, we divided the coefficient of the (daily) AE by 86400 s which followsthe approach in Sinnhuber et al. (2016). We find a maxi-mum production rate of about 1 cm − nT − s − around 70–72 km. This production rate also agrees with the one esti-mated by Sinnhuber et al. (2016) by a superposed epochanalysis of summertime NO . Comparing the NO produc-tion 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 ion-ization derived from AE on these days. The AE values of105, 355, and 435 nT translate to 105, 355, and 435 NO molecules cm − s − , about 4 times larger than would be es-timated from the ionization rates in Verronen et al. (2013)but agreeing with Sinnhuber et al. (2016). The factor of 4may be related to the slightly different locations, around60 ◦ N (Verronen et al., 2013) compared to around 65 ◦ N here and in Sinnhuber et al. (2016), in which the ionization ratesmay be higher.The associated constant part τ of the lifetime ranges fromaround 1 to around 4 days, except for large τ at low latitudesaround 70 km. As already discussed in Sect. 5.2, these largelifetimes may be a side effect of the small geomagnetic coef-ficients and more or less arbitrary. The magnitude is similarto 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 north-ern latitudes with a nearly constant amplitude of 10 days be-tween 70 and 85 km. An empirical lifetime of 10 days in win-ter was used by Sinnhuber et al. (2016) to extend the NO pre-dicted by the summer analysis to the larger values in winter.Here we could confirm that 10 days is a good approximationof the NO lifetime in winter, but it varies with altitude. Thealtitude distribution agrees with the increasing photochem-ical lifetime at large solar zenith angles (Sinnhuber et al.,2016, Fig. 7b). The larger values in our study are similarlyrelated to transport and mixing effects which alter the ob-served lifetime. The small variation of the lifetime at highsouthern latitudes could be a sampling issue because SCIA-MACHY only observes small variations there in winter (seeFigs. 1 and 2). Note that the results (in particular the largeannual variation) in the northernmost latitude bin should betaken with caution because this bin is sparsely sampled bySCIAMACHY, and the large winter NO concentrations areactually absent from the data. We propose an empirical model to estimate the NO densityin the mesosphere (60–90 km) derived from measurementsfrom SCIAMACHY nominal-mode limb scans. Our modelcalculates NO number densities for geomagnetic latitudesusing the solar Lyman- α index and the geomagnetic AE in-dex. Two approaches were tested, a linear approach contain-ing annual and semi-annual harmonics and a non-linear ver-sion using a finite and variable lifetime for the geomagneti-cally induced variations. From our proposed models, the lin-ear variant only describes part of the NO variations. It candescribe the summer variations but underestimates the largenumber densities during winter. The non-linear version de-rived from the SCIAMACHY NO data describes both varia-tions using an annually varying finite lifetime for the particle-induced NO . However, in cases of dynamic disturbances ofthe 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 varia-tions are treated as statistical noise.Sinnhuber et al. (2016) use a superposed epoch analysislimited to the polar summer to model the NO data. Herewe extend that analysis to use all available SCIAMACHYnominal-mode NO data for all seasons. However, during summer the present results show comparable NO produc-tion per AE unit and similar lifetimes to the Sinnhuber et al.(2016) study.The parameter distributions indicate in which regions thedifferent processes are significant. We find that these dis-tributions match our current understanding of the processesproducing 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 ofLyman- α (or solar UV radiation in general) is largest at lowand middle latitudes, which is explained by the direct pro-duction of NO via solar UV or soft X-ray radiation (Barthet al., 1988, 2003). The geomagnetic influence is largest athigh latitudes and is best explained by the production fromcharged particles that enter the atmosphere in the polar re-gions along the magnetic field.A potential improvement would be to use actual measure-ments of precipitating particles instead of the AE index. Us-ing measured fluxes could help to confirm our current un-derstanding of how those fluxes relate to ionization (Turunenet al., 2009; Verronen et al., 2013) and subsequent NO pro-duction (Sinnhuber et al., 2016). Furthermore, including dy-namical transport, as for example in Funke et al. (2016),could improve our knowledge of the combined direct and in-direct NO production in the mesosphere. Author contributions.
SB developed the model, prepared and per-formed the data analysis and set up the manuscript, MS providedinput on the model and the idea of a variable NO lifetime. JPB andPJE contributed to the discussion and use of language. All authorscontributed to the interpretation and discussion of the method andthe results as well as to the writing of the paper. Code and data availability.
The SCIAMACHY NO data set usedin this study is available at https://zenodo.org/record/1009078 (Bender et al., 2017a). The python code to pre-pare the data (daily zonal averaging) and to perform the analy-sis is available at https://zenodo.org/record/1401370 (Bender, 2018a) or at https://github.com/st-bender/sciapy . The daily zonal mean NO data and the sampled pa-rameter distributions are available at https://zenodo.org/record/1342701 (Bender, 2018b). The solar Lyman- α indexdata were downloaded from http://lasp.colorado.edu/lisird/data/composite_lyman_alpha/ , the AE indexdata were downloaded from http://wdc.kugi.kyoto-u.ac.jp/aedir/ , and the daily mean values used in this study areavailable within the aforementioned data set. Acknowledgements.
S.B. and M.S. thank the Helmholtz-society forfunding part of this project under the grant number VH-NG-624.S.B. and P.J.E. acknowledge support from the Birkeland Centerfor Space Sciences (BCSS), supported by the Research Council ofNorway under the grant number 223252/F50. The SCIAMACHYproject was a national contribution to the ESA Envisat, funded by . Bender et al.: SCIAMACHY mesosphere NO model 11
References
Funke, B., López-Puertas, M., von Clarmann, T., Stiller, G. P.,Fischer, H., Glatthor, N., Grabowski, U., Höpfner, M., Kell-mann, S., Kiefer, M., Linden, A., Mengistu Tsidu, G., Milz, M.,Steck, T., and Wang, D. Y.: Retrieval of stratospheric NOx from5.3 and 6.2 µm nonlocal thermodynamic equilibrium emissionsmeasured by Michelson Interferometer for Passive AtmosphericSounding (MIPAS) on Envisat, J. Geophys. Res. Atmos., 110,D09 302, https://doi.org/10.1029/2004JD005225, 2005b.Funke, B., García-Comas, M., López-Puertas, M., Glatthor, N.,Stiller, G. P., von Clarmann, T., Semeniuk, K., and McConnell,J. C.: Enhancement of N O during the October–November2003 solar proton events, Atmos. Chem. Phys., 8, 3805–3815,https://doi.org/10.5194/acp-8-3805-2008, 2008a.Funke, B., López-Puertas, M., Garcia-Comas, M., Stiller, G. P.,von Clarmann, T., and Glatthor, N.: Mesospheric N y produced by energetic particle precipitationin 2002-2012, J. Geophys. Res. Atmos., 119, 13,565–13,582,https://doi.org/10.1002/2014jd022423, 2014a.Funke, B., López-Puertas, M., Stiller, G. P., and von Clarmann, T.:Mesospheric and stratospheric NOy produced by energetic par-ticle precipitation during 2002–2012, J. Geophys. Res. Atmos.,119, 4429–4446, https://doi.org/10.1002/2013JD021404, 2014b.Funke, B., López-Puertas, M., Stiller, G. P., Versick, S.,and von Clarmann, T.: A semi-empirical model formesospheric and stratospheric NO y IAGA Special Sopron Book Series ,https://doi.org/10.1007/978-90-481-9858-0, 2011.Maric, D. and Burrows, J.: Formation of N O in the photol-ysis/photoexcitation of NO, NO and air, J. Photochem.Photobiol., A, 66, 291–312, https://doi.org/10.1016/1010-6030(92)80002-d, 1992.Marsh, D. R., Solomon, S. C., and Reynolds, A. E.: Em-pirical model of nitric oxide in the lower thermo-sphere, J. Geophys. Res. Space Phys., 109, A07 301,https://doi.org/10.1029/2003JA010199, 2004.Matérn, B.: Spatial Variation: Stachastic Models and Their Ap-plication to Some Problems in Forest Surveys and Other Sam-pling Investigations, vol. 36 of Lecture Notes in Statistics . Bender et al.: SCIAMACHY mesosphere NO model 13 O production by high energy auroralelectron precipitation, J. Geophys. Res. Atmos., 113, D16 302,https://doi.org/10.1029/2007jd009690, 2008.Sheese, P. E., Walker, K. A., Boone, C. D., Bernath, P. F., andFunke, B.: Nitrous oxide in the atmosphere: First measurementsof a lower thermospheric source, Geophys. Res. Lett., 43, 2866–2872, https://doi.org/10.1002/2015gl067353, 2016.Sinnhuber, M., Kazeminejad, S., and Wissing, J. M.: Interannualvariation of NOx from the lower thermosphere to the upperstratosphere in the years 1991–2005, J. Geophys. Res. SpacePhys., 116, A02 312, https://doi.org/10.1029/2010JA015825,2011.Sinnhuber, M., Nieder, H., and Wieters, N.: EnergeticParticle Precipitation and the Chemistry of the Meso-sphere/Lower Thermosphere, Surv. Geophys., 33, 1281–1334,https://doi.org/10.1007/s10712-012-9201-3, 2012.Sinnhuber, M., Friederich, F., Bender, S., and Burrows, J. P.: Theresponse of mesospheric NO to geomagnetic forcing in 2002–2012 as seen by SCIAMACHY, J. Geophys. Res. Space Phys.,121, 3603–3620, https://doi.org/10.1002/2015JA022284, http://onlinelibrary.wiley.com/doi/10.1002/2015JA022284/abstract,2015JA022284, 2016. Sinnhuber, M., Berger, U., Funke, B., Nieder, H., Reddmann, T.,Stiller, G., Versick, S., von Clarmann, T., and Wissing, J. M.:NO yy