Retrieval of sodium number density profiles in the mesosphere and lower thermosphere from SCIAMACHY limb emission measurements

An algorithm has been developed for the retrieval of sodium atom (Na) number density on a latitude and altitude grid from SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY) limb measurements of the Na resonance fluorescence. The results are obtained between 50 and 150 km altitude and the resulting global seasonal variations of Na are analyzed. The retrieval approach is adapted from that used for the retrieval of magnesium atom (Mg) and magnesium ion (Mg) number density profiles recently reported by Langowski et al. (2014). Monthly mean values of Na are presented as a function of altitude and latitude. This data set was retrieved from the 4 years of spectroscopic limb data of the SCIAMACHY mesosphere and lower thermosphere (MLT) measurement mode (mid-2008 to early 2012). The Na layer has a nearly constant peak altitude of 90–93 km for all latitudes and seasons, and has a full width at half maximum of 5–15 km. Small but significant seasonal variations in Na are identified for latitudes less than 40, where the maximum Na number densities are 3000–4000 atomscm. At middle to high latitudes a clear seasonal variation with a winter maximum of up to 6000 atomscm is observed. The high latitudes, which are only measured in the summer hemisphere, have lower number densities, with peak densities being approximately 1000 Na atomscm. The full width at half maximum of the peak varies strongly at high latitudes and is 5 km near the polar summer mesopause, while it exceeds 10 km at lower latitudes. In summer the Na atom concentration at high latitudes and at altitudes below 88 km is significantly smaller than that at middle latitudes. The results are compared with other observations and models and there is overall a good agreement with these.


Introduction
The metal sodium, Na, was first isolated in the laboratory by Sir Humphry Davy in 1807 (Davy, 1808).This was achieved by the electrolysis of very dry molten sodium hydroxide, NaOH, with Na being collected at the cathode.More than 100 years later the Earth's atmospheric Na layer was discovered in 1929 by American astronomer Vesto Slipher (Slipher, 1929).The large scattering cross section and atmospheric column of Na in the upper atmosphere results in a relatively strong emission in the visible spectral range.Sydney Chapman, who had previously worked on explaining upper atmospheric ozone, proposed a reaction-cycle theory to explain the nightglow phenomenon and the Na emissions (see e.g., Chapman, 1938Chapman, , 1939)).
There are two possible groups of Na sources for the upper atmosphere: terrestrial sources, such as volcanic erup-tions and salt particles from the oceans, and extraterrestrial sources such as meteoroids and comet dusts.In the upper atmosphere, meteoroids are the most likely sources of Na.
Meteoroids enter the Earth's atmosphere at supersonic speed and are decelerated and frictionally heated by collisions with air molecules.These processes lead to the meteoric ablation of metals and nonmetals into the upper atmosphere.The ablated metals are transported and react with the ambient neutral atmosphere.As a result metal layers are formed that have peak densities at around 85-95 km altitude (see e.g., Plane, 2003;Plane et al., 2015 for a review).Although the metal concentrations of several thousand atoms per cubic centimeter are low, these metals are strong emitters of radiation because they have large resonance fluorescence cross sections.Therefore, the metals are readily observed by remote sensing methods.Due to their strong radiation signal and their relatively long atmospheric lifetime, metal species are used as trace species to investigate wave propagation effects and winds in the mesosphere and lower thermosphere (MLT).Furthermore, the total amount of extraterrestrial material input can be estimated from measurements of these metal layers.Additionally, metals play an important role in upper atmospheric chemistry.Their chemical transformation impacts on ozone formation and loss both in the gas phase and through metal compounds acting as nucleation nuclei for the formation of aerosols and clouds in the middle atmosphere (see e.g., Rapp and Thomas, 2006;Voigt et al., 2005;Curtius et al., 2005).A detailed understanding of the origin and the reactions of metals in the upper atmosphere is required to understand the formation and loss of ozone and particles in the upper atmosphere.Also, metal ions are the principal component of ionospheric sporadic E layers and metal ions are found throughout the ionosphere.
Na has a large number density compared to other metals in the MLT, and the lower atmosphere is nearly transparent at the wavelength of the strongest Na transitions at 589 nm.This simplifies the observation from the ground.As a result, the mesospheric Na layer is the best understood metal layer in the MLT.
As the metal layers reside at an altitude where the atmosphere is too thin for aircraft to fly, but too dense for satellites to orbit for longer time periods, in situ measurements are only possible with rockets, which can only be launched at a limited number of locations on Earth and are expensive.Remote sensing methods are used and ground observations, e.g., by lidar, yield good vertical and time-resolved results, however, only at selected locations.
In the last decade, global satellite observations of Na with long time series have been available.The first global spacebased observations were reported by Fussen et al. (2004), using the GOMOS (Global Ozone Monitoring by Occultation of Stars) instrument on the satellite Environmental Satellite (Envisat).Envisat also carries the SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHar-tographY) instrument, the measurements of which are used in this study.Other observations have been reported from OSIRIS (Optical, Spectroscopic and Infrared Remote Imaging System) on the Odin satellite (Gumbel et al., 2007;Fan et al., 2007;Hedin and Gumbel, 2011), SCIAMACHY (Casadio et al., 2007) and GOMOS (Fussen et al., 2010).Results of the Na number density retrieval from the SCIA-MACHY limb MLT measurements between 2008 and 2012 are presented in this study.These results are compared to other measurements and models.
This manuscript is structured as follows: in Sect. 2 the SCIAMACHY measurements and the Na density retrieval algorithm are explained.The results for Na number densities will be presented in Sect.3, retrieved from both Na D lines at 589 nm.In Sect. 4 differences between results from both D lines are discussed and the results are compared to other measurements and model results.The seasonal and annual changes are investigated.Finally, the findings of this study are summarized in Sect. 5.

SCIAMACHY
The limb observations of scattered solar electromagnetic radiation, observed by SCIAMACHY on board Envisat are used for this study (see Burrows et al., 1995;Bovensmann et al., 1999  In the occultation mode, SCIAMACHY observes either the Sun or the Moon through the atmosphere at sunrise and moonrise, respectively.In nadir mode, the instrument points downward towards the Earth's surface and scans the upwelling radiation at the top of the atmosphere.The nadir mode delivers a good latitudinal and longitudinal coverage.
In limb mode, the instrument points tangentially to the Earth's surface at different tangent altitudes, resulting in an adequate vertical resolution in the scanned range of altitudes, but with a poorer latitudinal and longitudinal resolution than the nadir mode.
The highest tangent altitude of the nominal limb mode, which was performed daily during the whole SCIAMACHY lifetime, is about 91 km, which is just the altitude of the Na layer peak.However, from the middle of 2008, the limb MLT mode, which scans the altitude range between 50 and 150 km in 30 steps of 3.3 km with a good vertical resolution (also around 3.3 km) at the altitude of the metal atom and ion layers, was performed for 1 day of measurements every 2 weeks.As the ion layers (e.g., Mg + , see Langowski et al., 2015) are located at slightly higher altitudes than the neutral layers, it was decided to exploit the MLT mode of SCIA-MACHY first, prior to later investigation of the standard profiling mode.However, the retrieval results from the nominal limb and limb MLT measurements should not be too different and the data set presented here will be extended later.For each MLT limb scan, there is an additional measurement of the dark signal at 350 km tangent altitude.For Na this dark signal, which is subtracted from the signal at the other tangent altitudes, is much weaker than the signal at MLT tangent altitudes.At 590 nm SCIAMACHY has a spectral resolution of 0.44 nm.This is sufficient to resolve the two Na D lines, D 1 at 589.756 nm and D 2 at 589.158 nm, which have a spectral separation of ≈ 0.6 nm.
The SCIAMACHY data set employed in this study is Level 1 data version 7.03 and 7.04 and was calibrated with ESA's calibration tool scial1c with options 1, 2, 4, 5, 7 and M-factors, which include option 3 (see Sherbakov and Lichtenberg (2008) for more details.).
The Level 1 data have been averaged using the same approach as Langowski et al. (2015) for Mg and Mg + .An average for same latitude and local time of the up to 15 orbits of SCIAMACHY data is formed before the retrieval.The multiannual monthly means of the results for the 2008-2012 data set are formed after the retrieval.Note that there is a larger latitudinal coverage for the Southern Hemisphere than the Northern Hemisphere, because the northern dayside highlatitude measurements suffer from solar straylight contaminations.This is because the Sun is partly in the field of view of the instrument (see Langowski et al., 2015;Langowski, 2015 for more details).There is also a larger coverage of high latitudes compared to the Mg/Mg + retrieval by Langowski et al. (2015), because the better signal to noise ratio of Na produced fewer edge effects in the retrieval for the outermost measurements on the altitude-latitude grid.Note that due to the high signal to noise ratio of the Na emission signal, it might be possible to also retrieve reasonable profiles from single measurements.This needs to be investigated in the future.Note that due to the nonlinear forward model, averaging data before and after the retrieval leads to different results.For this retrieval a higher variability due to statistical errors or true natural variability will increase the density, when the density is averaged after the retrieval step.This issue is discussed in Langowski et al. (2014) their Sect.4.1.,where retrievals with a high statistical error before the retrieval step, shown in their Fig. 24, show significant larger densities than the retrieval with data averaged before the retrieval step.We, however, assume this effect to be small enough to not significantly influence the retrieval result for Na.It is assumed that longitudinal variations are much smaller than latitudinal variations.

Retrieval algorithm and adaption to Na
The retrieval algorithm presented in Langowski et al. (2014), which was used for magnesium atom and ion retrievals from the SCIAMACHY limb MLT measurements (Langowski et al., 2015), is used and adjusted to the specific parameters of Na atoms.To reduce the need to refer to the original paper too often, the most important steps of the retrieval algorithm from Langowski et al. (2014) are repeated in this section.

Forward model for the single scattering algorithm
A forward model for the emission signal and absorption path of each limb measurement of an orbit is set up and inverted for the number densities of the emitting species on a 2D latitude-altitude grid.
The mathematical representation of the forward model is with I being the wavelength integrated emitted radiance, emissivity γ , density n and an absorption part -along the line of sight (LOS) and along the line from Sun (LFS) to the point of scattering into the LOS (s a stands for both absorption paths)f .The emitted radiance is the integrated product of the density n and the emissivity γ along the emission path s e , which is furthermore attenuated by selfabsorption f (see Eq. ( 6) for f ).Equation ( 1) is discretized on the 2D latitude-altitude grid and inverted for the number density n.The changes with respect to the Mg/Mg + retrieval described by Langowski et al. (2014) lie in the emissivity γ and self-absorption f calculation, while the rest of the retrieval algorithm remain unchanged, beside marginal changes (e.g., there is no correction for inelastic scattering needed for the Na retrieval, because the ratio of emission radiance to Rayleigh scattered radiance is high enough that this effect is negligibly small).The emissivity γ is calculated as follows: with π F (λ) being the solar irradiance (note that it is convention to use π F (λ) (see e.g., Chamberlain et al., 1958) and π belongs to the symbol and is not meant as a factor 3.14) and  , which agrees for both, is used in the approximate formula by McNutt and Mack (1963).The right figure shows the SAO 2010 spectrum compared to the fully resolved approximation of the D 1 and D 2 lines taken from McNutt and Mack (1963).Note that the approximate formula by McNutt and Mack (1963) is only valid close to the center of the Fraunhofer lines.The fully resolved Na Fraunhofer lines are much deeper than the lines in the SAO 2010 spectrum and the SCIAMACHY spectrum.
The process causing the emission is resonance fluorescence (A detailed derivation of Eq. (3) from basic principles as well as a discussion of the phase function P is given in the appendix of Langowski (2015).).A solar photon is absorbed by a Na atom, which is excited from the lower state i to the higher state j and is immediately re-emitted, which returns the atom to the lower state i.The two relevant transitions are from the lowermost excited states 3 2 P 1 2 for D 1 and 3 2 P 3 2 for D 2 to the ground state 3 2 S 1 2 .The relative Einstein coefficient, the probability of the resonant transition compared to all other possible transitions from the upper state to lower states, is 1 for both lines, because only the two lowest excited states are involved as upper states and the transition between the P states is highly improbable.The Na specific parameters, i.e., oscillator strength f ij and transition wavelength λ ij , are taken from the NIST (National Institute of Standards and Technology) atomic spectra database (Kramida et al., 2012).The scattering angle θ dependent phase function P is a linear combination of the phase function for Rayleigh scattering and an isotropic part: P is normalized to 4π , which is already considered in Eq. (1).The factors E 1 and E 2 depend on the change in angular momentum j and are taken from Chandrasekhar (1960) (see Table 1).The factors E 1 and E 2 are different for both D lines.The D 1 line has a purely isotropic phase function (E 1 = 0 and E 2 = 1), while the D 2 line has a mixture of both components E 1 = 0.5 and E 2 = 0.5.The wavelengthintegrated absorption cross section has to be distributed over the correct shape function of the emission line and the re-Table 1. E 1 and E 2 depend on the angular momentum j and the change of angular momentum j (from Chandrasekhar, 1960).
sulting absorption cross section profile is multiplied by the wavelength-dependent solar irradiance.This product is then integrated over all wavelengths, yielding the true combination of the second and third factor of the emissivity in Eq. ( 2).
A high-resolution spectrum of the solar irradiance in the vicinity of the Na D 1 and D 2 lines needs to be used.Figure 1 shows solar spectra for this wavelength region from SCIAMACHY, Chance and Kurucz (2010) and McNutt and Mack (1963), from which only the latter fully resolves the Na D lines.Note that the approximate formula for the spectrum measured by McNutt and Mack (1963) is only valid in the vicinity of the center of the lines.Following McNutt and Mack (1963), the solar irradiance in the vicinity of the Na D 1 and D 2 lines can be calculated for x (x is defined further below): (5) The ratio I 0 I baseline of the intensity at the line center and the baseline intensity at the edge of the Fraunhofer lines is stated in McNutt and Mack (1963).To scale this to the SCIA-MACHY spectrum, a solar irradiance of 5.44 × 10 14 photon s cm 2 nm is used for the edge of the line; for the D 2 line, the following values are used: I 0 = 0.0444 × 5.44 × 10 14 photon s cm 2 nm A = 2.16 The formula is given for wavenumbers k and the parameters x and x e are normalized to the wavenumber of the line center.
The width parameter of the line in terms of wavenumbers is denoted by σ (see McNutt and Mack, 1963).Wavenumber shifts between the solar spectrum and the absorption cross section in the mesosphere are considered.Here, a positive shift value leads to a red shift as the line center is moved toward shorter wavenumbers and, therefore, longer wavelengths.Since the width of the solar Fraunhofer line is not much larger than the width of the lines in the mesospheric absorption cross sections, these shifts have a nonnegligible influence on the emissivity.In McNutt and Mack (1963) a constant red shift for the solar lines of grav shifts k line center = 2.7 × 10 −6 is measured, which is a combination of the gravitational red shift and other smaller shifts, e.g., pressure shifts.Additionally, Doppler shifts from the rotation of Earth and the change of the Earth-Sun distance along the elliptical orbit of Earth are considered, which are similar in magnitude to the constant shift and have a combined maximum amplitude of ± 3.2 ×10 −6 .For the D 1 line the following parameters are used: I 0 = 0.0495 × 5.44 × 10 14 photon s cm 2 nm A = 2.14 x e = 12.8 × 10 −6 .
Because there are hyperfine splittings for both D lines, the line center is the weighted average of the individual line wavelengths and strengths.
Na has only one stable isotope, i.e., 23 11 Na, and therefore has no isotope effect.The stable isotope has a nuclear angular momentum of I = 3 2 , which leads to a hyperfine splitting of the energy levels.The splitting for the lower 3 2 S 1 2 state is stronger than the splitting of the upper states 3 2 P 1 2 and 3 2 P 3 2 .This can be explained phenomenologically by the smaller distance of the valence electron and the nucleus in the S state, which leads to a larger overlap of the nucleus and electron wave functions and therefore a stronger perturbation of the electron state.Due to the stronger splitting of the S state compared to the P states, the D 1 and D 2 lines each split into two groups of lines in close spectral vicinity (this is called an "s-resolved blend" in McNutt and Mack, 1963).Both are used in high resolution to calculate the emissivity γ (Eq.2) and attenuation factor f (Eq.6).
The Doppler width of the Na lines in the mesosphere is approximately 1.25 pm.The two groups of adjacent lines have a separation of about 2 pm and thus can be separated.However, the lines inside a group are too narrow to be resolved in the mesosphere.The existence of several degenerate lines, however, is important for the correct weighting of the two sresolved lines.This is, e.g., well explained in Chamberlain et al. (1958); McNutt and Mack (1963); Fricke and von Zahn (1985).
The solar spectrum as well as the mesospheric absorption cross section for the Na D 2 line are shown in Fig. 2. The Na density is large enough that a non-negligible part of the incoming solar irradiation is either absorbed along the path from the Sun to the point of resonance fluorescence, or along the line of sight, after the emission.This reduction of emissivity is considered in the self-absorption factor f : with the integrated true slant column density g = nds.Note that the integral in Eq. ( 1) contains n as a linear factor, but also has a nonlinear dependence on n because of the selfabsorption factor f .The equation is linearized, employing an iterative approach using the retrieved density of the previous step in the nonlinear part to retrieve the linear density n.As noted by Langowski et al. (2014), setting f = 1, which corresponds to no self-absorption, is a good starting step for the iteration.

Discretization of the forward model and calculation of the Jacobian
For every single limb measurement, Eq. ( 1) is discretized on a latitude-altitude grid with m elements to yield the formula: x j s LFS gc i,j .
c 1 is a constant including all parameters that are path independent in the integral in Eq. ( 1).For every grid element i along the LOS, there is the path element in this grid element s LOS i for the emission path.Furthermore, for every grid element i along the LOS, there are two matrices with path elements.One matrix contains the path elements from the satellite to the grid element with matrix elements s LOS gc i,j .The other matrix contains the path elements from the Sun to the grid element with matrix elements s LFS gc i,j .Both matrices are needed for the absorption calculations of the attenuation factor f .P i is the phase function in grid element i, and x i is the density in grid element i that should be retrieved.Figure 3 is a sketch of the considered paths.
The argument of f is summarized to To invert the forward model Kx = y, the equation system J T Jx = J T x needs to be solved, with J being the Jacobian of Kx.If K was independent of x, J simply was equal to K. The elements of the Jacobian J are calculated with the following formula: If f were independent of x, f = ∂f (g i (x)) ∂g i (x) would be zero and the second addend in Eq. ( 9) would also be zero.Then J would be J = K.In practice, omitting the second addend in Eq. ( 9) reproduces synthetic model density profiles that are forward modeled by Eq. ( 1) and inverted with the retrieval algorithm very well and this simplification J = K(x) is used.However, K still explicitly depends on x, which has to be later considered in the retrieval step.

Calculation of path matrices
Figure 4 is a sketch illustrating the crucial steps for the calculation of the path matrices.Note that this is a minimalistic description and more details are given in Langowski (2015).The path segments of the altitude grid are calculated with right-angled triangle geometry.The right-angled triangles considered for the LOS are formed by the path between the center of Earth and tangent point b, the right angle between b and the segment of the line of sight a and the path between the center of Earth and a grid element's upper altitude limit c.As b and c and the right angle are known, a can be calculated for each grid element.The path length within each altitude grid element is calculated by the difference of a for neighboring grid elements.For the LFS, the tangent point of the LFS needs to be found first and the position of the grid element with respect to the tangent point needs to be evaluated.For the calculation of the absorption part, it is crucial to separate the LOS at the tangent point and do the respective calculation of the factor f for each side of the tangent point, before adding both parts of the Jacobian together.In spherical coordinates, positions with the same latitude (north and south) can be calculated by a double cone equation.These double cone equations for the latitude boundaries of a grid element are set equal with the straight line equation r = T P + λ êLOS , where T P is the position of the tangent point, êLOS a vector in LOS direction and λ the distance between tangent point and the considered point on the line of sight.λ is determined for each latitude element.Once this is done, the latitude grid boundaries λ and the altitude grid boundaries a are sorted and the path lengths are determined by differences of neighboring boundary parameters (see Fig. 4, where red crosses are latitude changes and blue circles are altitude changes).Figure 5 shows the calculated path lengths for one limb measurement.The latitudinal separation of consecutive limb measurements is around 7 • at the equator.It is smaller near the poles at the cost of a larger variation in local time, which is considered by splitting the orbit in an ascending satellite movement and a descending satellite movement part.For daily averages, the latitudinal separation is roughly halved, as SCIAMACHY has an alternating scanning pattern of limb and nadir scans for consecutive orbits.The largest paths are found in the close vicinity of the tangent point.The emission signal strength for each grid element is roughly the product of density and path length.For the around 4 • latitudinal separation, the overlap region of consecutive limb measurements at the same altitude is less than 5 km above the tangent altitude.However, note that neighboring measurements in the averaged data come from different longitudes; thus this is not a real overlap of the same volume of air.

Inversion of the forward model
The equation system to be solved J T Jx = J T x is nonlinear as J explicitly depends on x, because of the self-absorption contributions in f and f .The equation system is linearized, by using initial values x for the calculation of J and retrieving new x values that are closer to the real values of x.This is iteratively done until convergence of x is achieved.Test retrievals showed that after roughly five iteration steps, convergence is achieved; i.e., the largest step by step changes in the retrieved densities are far less than 1 %.In practice, 20 iteration steps are used.
As J T Jx = J T x is typically an ill-posed mathematical problem, solutions oscillate and are not physically correct if no constraints are applied.Therefore, three different smoothing constraints are applied: latitudinal smoothing, which penalizes solutions with differences in densities of grid elements with neighboring latitudes; altitudinal smoothing, which penalizes solutions with differences in densities of grid elements with neighboring altitudes; as well as Tikhonov regularization (Tikhonov, 1943) with a zero a priori x a which in general favors solutions with the smallest oscillations and overall smallest distance from zero.The final equation that has to be solved is Eq. ( 10): The a priori covariance matrix S a is in fact a scalar (λ apriori ) multiplied by an identity matrix.S H and S φ are the matrices for altitudinal and latitudinal constraints (large sparse matrices with only two diagonals of non-zero values) and λ h and λ φ are the scalar weighting factors for both constraints.
x is the vector of number densities.On the right-hand side, there is the covariance matrix for the slant column densities (SCDs) (S y ), which is assumed to be diagonal: the SCDs y and the a priori x a .As there should not be any bias on the form of the profile, x a = 0 is used.
There is some arbitrariness in the choice of the constraint strength λ h , λ φ and λ apriori .We choose a ratio of these three of 10 : 2 : 1.The constraint strength should be chosen so strong that the solution is affected but not dominated by the smoothing.As will be shown in Sect.2.3.1, there is a range of 2 orders of magnitude in the choice of the constraint parameters which only results in moderate changes of the final result, which shows that the arbitrariness in the decision of which parameter is finally used does not significantly influence the result.

Extension for multiple scattering
In the previous section, the optimization and adaptation of the single scattering retrieval algorithm developed for Mg and Mg + for the retrieval of Na were described.However, the single scattering approximation for the background signal is only valid for wavelengths below about 300 nm.In the visible region, radiation may be reflected from the Earth's surface or scattered back from the lower atmosphere into the mesosphere.As a result, a part of the incoming solar irradiation may pass the grid cells more than just once and can, therefore, produce more emission.This is considered by a factor multiplied to the solar irradiance, which will be called the albedo factor in the following.
Different correction methods for Na retrievals with OSIRIS are reported by Gumbel et al. (2007) and Hedin and Gumbel (2011).In Hedin and Gumbel (2011) the background signal for the limb scan at 40 km tangent altitude is compared to a single scattering radiative transfer model considering Rayleigh scattering only.Because the lowest tangent altitude of SCIAMACHY limb MLT measurements is at 53 km, this approach cannot directly be used, because the background signal is too small at this altitude compared to straylight contamination.Because SCIAMACHY can resolve both Na D lines, which are differently sensitive to selfabsorption, another approach to determine the amount of radiation passing the Earth's atmosphere, being reflected at the Earth's surface and then being scattered into the limb field of view of the instrument, is presented in Sect.2.3.1.Unfortunately, this approach did not always yield reasonable results; therefore another approach is presented in Sect.2.3.2, which is similar to the one in Hedin and Gumbel (2011), and which was finally used.

Total to single scattering ratio estimation from direct comparison of the D 1 and the D 2 line
The approach discussed in this section is based on the expectation to retrieve the same densities from both D lines.Therefore, the albedo factor is assumed to be the factor where both density profiles are nearly the same.In practice, however, this method fails too often; hence another method is used, which is described in Sect.2.3.2.The albedo factor is determined as the factor for which both D lines yield the same Na number densities.For typical Na slant column profiles shown in Fig. 6, the identifi-  6.For albedo factors that are too large, the densities are smaller and the D 1 line shows higher densities.For a albedo factor that is too small, the D 2 line shows higher densities.For the optimum albedo factor (here 1.2), both Na lines yield the same densities.cation of the optimal albedo factor is illustrated in Fig. 7.It should be noted that the phase function (Eq.4) is changed to P new = P old +(albedo factor−1) albedo factor .This takes into account the increased part of multiply scattered radiation, which we assume is unpolarized and therefore effectively increases E 2 in Eq. ( 4).
This approach is only reasonable if the retrieval result is nearly independent of the applied constraint parameters (e.g., vertical smoothing, necessary to reduce oscillation of the result).This is not the case for the Mg/Mg + retrieval by Langowski et al. (2014).For Na the statistical errors are sufficiently small; thus the Na retrieval is much less sensitive to smoothing constraints than that of Mg.However, the dependence on the variation of the constraint parameter is unfortunately not entirely negligible, especially when the densities are large, as is the case in Fig. 8. Figure 8 shows retrievals for different constraint parameters.For moderate constraint parameter (1 × 10 −7 to 5 × 10 −10 ) the retrieved peak density is nearly independent of the choice of the constraint parameter for approximately 2 orders of magnitude.A factor of 5-10 in the constraint parameter has a similar effect to a change of the albedo factor of 0.1 in Fig. 7.The three highest constraint parameters (5 × 10 −7 -5 × 10 −5 ) show smoothing that is too strong, while the lowest constraint parameters (5×10 −9 -5×10 −10 ) show oscillations at high altitudes.Note that a stronger smoothing leads to the need of a lower albedo factor to match the D 1 and D 2 density, so a systematic error in one property is rather reduced than increased if the other one is tuned, which results in some robustness in the method.
Unfortunately, this method failed quite often for the following reasons: the densities were too small and insensitive to the albedo factor, the D 2 slant column densities were already larger than the ones for D 1 , or the differences between  6, albedo factor 1.2; see legend for constraint parameter value).
the D 1 and D 2 slant column densities initially were too large; therefore the D 2 line could yield larger densities than the D 1 even for albedo factors smaller than 1.However, although this algorithm failed the matching of densities retrieved from both Na lines, it is a good indicator for how good the calibration of the data and the radiative transfer model used work.

Total to single scattering ratio estimation from comparison of simulated single scattering limb radiance and measured limb radiance
This approach calculates the albedo factor based on the ratio of the measured radiance, compared to the simulated single scattered radiance in the vicinity of the line.This radiance should only come from Rayleigh scattering of the major atmospheric constituents (N 2 , O 2 , Ar etc.), whose concentrations in the atmosphere are well known; therefore the Rayleigh scattering is easy to simulate.However, like other similar instruments, SCIAMACHY has a straylight contamination issue at mesospheric altitudes; indeed, an altitude region must be used for this calculation, where all possible error sources are small.This altitude lies typically below 40 km and is not scanned by the SCIAMACHY limb MLT states.Nevertheless, retrieving the albedo factor from collocated nominal measurements and matching the nominal profile to the MLT measurement profile yields reasonable results.The Rayleigh scattered background radiance in the vicinity of the Na D lines can be used to calculate the total to single scattering ratio for a limb measurement.A simple approach to obtain the albedo factor is to use the ratio of the measured limb radiance and the limb radiance calculated with a radiative transfer model.The radiative transfer model SCIATRAN is able to calculate the single and total Rayleigh scattered electromagnetic radiation for known measurement geometries and atmospheric parameters.Figure 9 shows the ratio of the measured limb radiance and the simulated single scattering radiance for different tangent altitudes, as well as the ratio for the simulated total scattered simulated radiance and the single scattered radiance for different ground albedos.
As was reported by Oikarinen et al. (1999) the modeled total to single scattering ratio only shows a weak dependency on the tangent altitude.The measured limb radiance, however, has a completely different behavior and shows a nearly exponential increase in the total to single scattering ratio above 50 km altitude.The true Rayleigh scattered limb radiance is roughly proportional to the density at the tangent altitude and exponentially decreasing limb radiances with increasing tangent altitudes are expected.We assume that there is a small straylight component from lower tangent altitudes that reaches the instrument for high tangent altitudes and that this component is only weakly dependent on altitude.Above 50 km this additional straylight component is on the order of magnitude as the actual limb radiance at this tangent altitude and becomes much bigger than the actual Rayleigh scattered component.This nearly constant offset to the radiance along with the exponential decrease of the Rayleigh scattered radiance with altitude explains the nearly exponential rise of the measured to simulated single scattering ratio.In the troposphere and lower stratosphere, clouds significantly influence the radiance, and the simple approach using an albedo factor fails there.Therefore, we assume that in a region above 20 and below 45 km, there is a region where the limb measurement to simulated single scattering ratio is very close to its simulated value with the right ground albedo.
Any unwanted straylight also affects the dark signal measurement at 350 km tangent altitude, which is usually subtracted from the limb radiances at the 30 other tangent altitudes.Instead of simply subtracting the dark signal measure-   9 and the red line corresponds to the black line in Fig. 9.The green and the blue line are identical to the lines in Fig. 9.The red line is not due to the substraction of the additive component.ment, another approach is used here: we assume that a part of the total incoming radiation I inc (h, λ) is proportional to the simulated single scattering radiance I ss (h, λ), which we call the multiplicative part aI ss with the multiplicative component a, while a part of the straylight component and the actual dark radiance is an additional component b of the light: I inc (h, λ) = a(h)I ss (h, λ) + b(h) (with tangent altitude h and wavelength λ).To determine the multiplicative and additive components, the wavelength region between 650 and 660 nm is used.This spectral region is not affected strongly by atmospheric absorption in the mesosphere and upper stratosphere, and includes the H α Fraunhofer line at 656 nm, which is  a clear solar signature; hence fitting the multiplicative and additive component is not an ill posed problem (where a smaller a could be compensated by a higher b etc.).The estimated value for the total to single scattering ratio is the minimum of the multiplicative component a above 20 km altitude.Figure 10 shows the fit and Fig. 11 the results of the multiplicative and additive component for one example profile.
This estimation, however, can only be done for the nominal SCIAMACHY limb measurements measuring from ground to 90 km altitude.The MLT measurements start at around 53 km altitude; therefore the minimum between 20 and 45 km can not directly be found.However, the latitudinally and longitudinally co-located nominal measurements from the days of the same time period show very similar profile shapes, which can be fitted to MLT data to retrieve the albedo factor.Figure 12 shows the final fit of the albedo factor for an example measurement.
First the multiplicative components for all nominal and MLT measurements are found.The median for the days in the same time period (± 200 orbits were used here) of nominal limb measurements is formed for all altitudes.The albedo factor A for the median nominal measurements is found.Between 50 and 70 km the logarithms of the nominal and the MLT measurements are fitted as factor B (ln MLT = B ln nominal).Fitting the logarithm puts effectively more weight on the matching of the lower albedo factor values at lower altitudes, which considers that the perturbing effect is smallest there.The albedo factor for the MLT measurements is then given by the product AB.The resulting albedo factors show similarities to the simulated total to single scattering ratios, which are high at scattering angles at around 90 • and close to 1 for low (around 0 • ) and high (around 180 • ) scattering angles.Eichmann et al. (2015) found out that more than 90 % of SCIAMACHY limb measurements are influenced by tropo-spheric clouds.A typical spread of the albedo due to this can be seen in the spread of the nominal albedo ratios in Fig. 12.Assuming Fig. 12 represents the typical spread for this estimation using the factor A only results in an error of about 20 %, which is rather large for a number, whose value is known to be larger 1 but unlikely larger than 2.5.The additional fit of factor B should reduce this error to less than 10 %.

Monthly averaged Na densities
The monthly mean Na number densities are shown in Fig. 13 as a function of latitude and altitude.The average of both retrieval results ( 1 2 (D 1 + D 2 )) is used.The altitude of the density maximum is about 92 km and varies only by a few kilometers during the year.Na shows a seasonal cycle in number density with a winter maximum, with peak densities in the winter middle latitudes of up to 6000 cm −3 .In the summer the maximum density decreases to only ≈ 1000 cm −3 at high latitudes.The annual mean is shown in Fig. 14 and shows an average peak density of roughly 2000-4000 cm −3 apart from high latitudes, which are only measured in the summer period and, therefore, only show the small summer densities.The results for both hemispheres are nearly symmetric.
The seasonal variation of the vertical profile for low, middle and high latitudes is shown in Fig. 15, and vertical profile shapes for selected latitudes in July are shown in Fig. 16.At low latitudes a semiannual variation with maxima in March and September is observed.This variation is well correlated with the semiannual variation in temperature (see e.g., von Savigny and Lednyts'kyy, 2013), which shows a maximum during this time.The semiannual oscillation for Na was also found in model studies by Marsh et al. (2013).The vertical profiles at high latitudinal summer show a reduced width.
Figure 17 shows the seasonal variation of the Na number density for the region between 80 and 105 km altitude, as well as the seasonal variation of the normalized Na profiles at 71 • N, which is a high latitude and covered by SCIAMACHY observations for several months during boreal summer.
The peak altitude at 71 • N is about 93 km for most months.Not only is the density strongly reduced during summer, the profile also becomes thinner: from a full width at half maximum (FWHM) of 11 km in spring and autumn to only 7 km in summer.Qualitatively a similar reduction is also observed for other width-defining parameters than 50 % of the maximum value (e.g., 75, 25, 5 % of maximum).The reduction of the profile width occurs on both the lower and the upper edge of the profile.
Figure 18 shows the vertical column densities (VCDs) for different months and latitudes.The VCDs are formed by integration of the vertical profiles shown in Fig. 13.The VCDs also show the seasonal cycle with a summer minimum of slightly below 1 × 10 9 cm −2 at high latitudes and up to 1 × 10 10 cm −2 at the highest latitudes covered in the winter hemisphere.
4 Error discussion and validation

Estimation of errors in the vertical profile
Four to eight individual day measurements have been used to form the multiannual monthly averages shown in the previous section.The different number of days for each month is explained by the fact that 4 full years were not covered (hence there are fewer days available in summer) and furthermore for some individual days, the raw data are missing or the retrieval did not converge.The errors of the measured limb radiances are linearly propagated into the error of the SCDs to estimate the error of daily average SCDs.As the inversion step includes a nonlinear operation, a further linear error propagation from the SCDs to the number densities is not carried out.Instead we use the same Monte Carlo approach as Langowski et al. (2014).A random Gaussian error in the range of the typical error of a daily average SCD profile is applied to a daily average SCD profile, and the number densities are retrieved.This is repeated a large number of times (here 1000 repetitions were used).The mean and the standard deviation of the large number of Monte Carlo realizations are determined.The standard deviation determined in this way quantifies the error of the profile.The result of this method is shown in Fig. 19.
The error is smaller than the measured number density in the region between 80 and 100 km, with the D 2 line having a slightly smaller error than the D 1 line.In the maximum number density region, the relative error is roughly 10 %.In Fig. 19 both Na lines agree very well; however, this is not always the case, which will be discussed in the next section.

Differences between D 1 and D 2 retrieval results
For a comparison of the individual results based on the D 1 or the D 2 line only, Fig. 20 shows the absolute and relative differences for the VCDs of both lines.
The overall agreement of the results for the D 1 and the D 2 line is good, showing relative differences of only ±10 % for most months and latitudes.However, for the highest latitudes in southern hemispheric winter, the differences are larger, with absolute differences of up to 3 × 10 9 cm −2 which correspond to relative differences of up to 40 %.The large discrepancies occur in a region where the Na density, and therefore self-absorption, is high.
There may be different reasons for the discrepancies, e.g., remaining issues with the calibration of the data.Furthermore, small inaccuracies in the assumptions for the radiative transfer model contribute to the differences.For instance, the self-absorption approximation only considers loss along the line of sight and no contribution of multiple scattering into the line of sight.Furthermore, a constant width of the Doppler-broadened mesospheric absorption cross section is used; however, the width can slightly change with temperature.In addition to the resonance fluorescence dayglow, Na also shows a chemiluminescent nightglow (see e.g., Fussen et al., 2010;Plane et al., 2012).This nightglow shows a lower ratio of the D 2 /D 1 emission signal than 2, which is not caused by self-absorption.Thus, non-negligibly small chemiluminescence, compared to the resonance fluorescence, could explain larger retrieved densities for the D 1 line with the current method.

Comparison to independent data sets
In the following, the SCIAMACHY data set is compared to ground-and space-based measurements.Figure 21 shows a comparison to ground-based ALOMAR (Arctic Lidar Observatory for Middle Atmosphere Research) lidar measurements at Andøya, Norway, at 69 • N recorded between 2008 and 2013.The majority of these have been published by Dunker et al. (2015) and show a mean peak altitude of around 92 km, which is in good agreement with the peak altitudes found in the SCIAMACHY data.The SCIAMACHY error bars are the standard deviation of measurements in each month.The error bars of the respective daily mean Na column density measured by the lidar denote the standard deviation, which is a measure of the geophysical variation on that day.
The SCIAMACHY and the ALOMAR results are similar and on the same order of magnitude (a factor 3 for the largest differences) and overall in good agreement.In the summer months June and July, the ALOMAR and SCIAMACHY measurements of Na column density sometimes agree, but on several days, ALOMAR measures much smaller or much larger values.The ALOMAR Na lidar measures larger column density each time a sporadic Na layer appears in the lidar's observation volume.This usually happens during the night (see e.g., Heinrich et al., 2008), but not during SCIAMACHY's time of observation (11:00).The transience of sporadic Na layers results in larger geophysical variation in the red symbols of those particular nights in Fig. 21.We attribute the cases when ALOMAR observed much smaller column density than SCIAMACHY to (a) adsorption or absorption of Na atoms on noctilucent cloud (NLC) particles and smaller ice particles and (b) the temperature in the coldest phase of a gravity wave leading to a chemical balance with less atomic Na and more Na compounds.SCIAMACHY does not detect these strong variations, because the measurement volume is much larger than that of the lidar, and because SCIAMACHY averages in longitude.The ALOMAR VCDs are larger in April and September than in June and July, and thus are in general agreement with the SCIAMACHY seasonal variation.Lidar results for the winter month at Andoya are available from Tilgner and von Zahn (1988), showing VCDs of around 5×10 9 cm −2 between December and February.In their Fig. 7, Dunker et al. (2013) show VCDs between 3×10 9 and 7×10 9 cm −2 in December.Again the seasonal cycle for the two different techniques is clear.
The left panel of Fig. 22 shows a comparison to groundbased lidar measurements by She et al. (2000) at around 40 • N. The overall agreement of the seasonal variation and density is good.
The peak altitude obtained by She et al. (2000) is also 90-92 km.The FWHM of the peak is 9-10 km (not shown here), which is slightly smaller than the 11-14 km for SCIA-MACHY.However, SCIAMACHY scans a much larger volume of space, which might in part explain this difference.Also, the SCIAMACHY retrievals have a vertical resolution of about 4 km.The largest differences between SCIA-MACHY and the values presented by She et al. (2000) are about 20 % from July to November.
The data set of She et al. (2000) was also used for comparison with other global satellite measurements by Fussen et al. (2010) (their Fig. 7) and the WACCM (Whole Atmosphere Community Climate Model) model by Marsh et al. (2013) (their Fig. 10).This comparison is shown in the right panel of Fig. 22.There is a good qualitative and quantitative global agreement in seasonal variation with measurement results reported by Fussen et al. (2010) and Hedin and Gumbel (2011), as well as a model with results by Marsh et al. (2013).The reduction of the profile width in the summer high northern latitudes is tentatively explained by Marsh et al. (2013) by an increased ionization rate at the upper edge and the Na reaction into NaHCO 3 .The latter is favored by low temperatures in the polar summer mesopause region at the bottom edge.A further reduction process at the bottom edge, which is, e.g., shown by Gardner et al. (2005), Fig. 9, is the uptake of atomic Na by NLCs or smaller ice particles.This process also depends on low temperatures; therefore both reduction processes at the lower edge are correlated and hardly distinguishable.

Conclusions
The extension of a retrieval approach previously used for the retrieval of mesospheric magnesium, to the retrieval of Na number density profiles in the MLT region is presented.Monthly mean Na number densities on a latitudinal and vertical grid retrieved from the SCIAMACHY limb MLT measurements from 2008 to 2012 are presented.The Na peak has a FWHM of 5-15 km which depends on latitude and season and is smallest in summer at high latitudes.The peak density varies from 1000 to 6000 cm −3 and shows a clear seasonal cycle with a summer minimum most pronounced at high latitudes.The retrieved SCIAMACHY data set is in good agreement with other measurements and models.This data set provides a unique set of data for testing our understanding of the role of meteoroids and their release of Na coupled with the other metal and metal ions, e.g., Mg and Mg + , on upper atmospheric chemistry.
In the future, the SCIAMACHY Na data product will be extended to the period from 2002 to 2012 by also applying the retrieval algorithm to the nominal SCIAMACHY limb data, using the results of the MLT retrieval as a priori information.
of Greifswald.SCIAMACHY data were kindly provided by the European Space Agency (ESA).We would also like to thank ESA, as a part of the work was funded through the ESA MesosphEO project.Ulf-Peter Hoppe and Tim Dunker are grateful to the Research Council of Norway for funding the Na lidar measurements at ALOMAR through grants 208020/F50 and 216870/F50.For the maintenance and operation of the ALOMAR Na lidar we also wish to thank the National Science Foundation (NSF) for financial support through grant NSF AGS-1136269.
The article processing charges for this open-access publication were covered by the University of Bremen.

Figure 1 .
Figure 1.Data of the solar extraterrestrial spectrum in the region 586-594 nm measured by SCIAMACHY on 3 June 2010 and from the Smithsonian Astrophysical Observatory (SAO; Chance and Kurucz, 2010).The Fraunhofer lines are readily observed and the left figure shows both the SCIAMACHY and the SAO 2010 spectrum, from which the baseline value of 5.44 × 10 14 photon s cm 2 nm

Figure 2 .
Figure 2. Absorption cross section for the D 2 line at temperature T = 220 K (blue) and solar irradiance spectrum for different shifts.Both are used in high resolution to calculate the emissivity γ (Eq.2) and attenuation factor f (Eq.6).

Figure 3 .
Figure 3. Scheme of the altitude and latitude grid (2D projection) and the line of sight and line from Sun.Each grid element has a different line from Sun and absorption path along the line of sight.

Figure 4 .
Figure 4. Two-dimensional intersection of Earth's atmosphere, with the center of Earth, M, and the altitudes as radii of the circles.Path lengths along the line of sight for the vertical grid can be found with right-angled triangle algebra.Changes of the latitude are marked with red crosses as additional "a" sides.The path length in each grid cell is the difference of the a sides of neighboring grid cells.Note that, depending on the binning of the latitudes, it is possible that all grid cells only have one latitude (taken from Langowski et al. (2014)).

Figure 5 .
Figure 5. Path lengths in different grid cells for a typical line of sight (LOS) of a limb measurement.The biggest part of the path lies in the tangent point altitude region, but higher altitudes are also passed (taken from Langowski et al. (2014)).

Figure 6 .
Figure 6.Vertical SCD profile for both Na D lines (6 October 2009, 35 • N) for three different albedo factors (1.1-1.3).The SCDs of the D 2 line are smaller because of stronger self-absorption.

Figure 7 .
Figure7.Vertical Na number density profile retrieved from the SCD profiles shown in Fig.6.For albedo factors that are too large, the densities are smaller and the D 1 line shows higher densities.For a albedo factor that is too small, the D 2 line shows higher densities.For the optimum albedo factor (here 1.2), both Na lines yield the same densities.

Figure 8 .
Figure8.Vertical Na density profiles for the D 1 and D 2 lines and for different constraint parameters (same conditions as in Fig.6, albedo factor 1.2; see legend for constraint parameter value).

Figure 9 .
Figure9.Measured radiance to single scattering ratio (averaged between 649 and 661 nm, see Fig.10) for different tangent altitudes of a series of nominal limb measurements as well as total to single scattering ratio for different ground albedos simulated with SCIA-TRAN.

Figure 10 .
Figure 10.Fit of the multiplicative and additive radiation component around the H α line at 656 nm.

Figure 11 .
Figure11.Fit results for the multiplicative and additive component as well as the simulated values for the same ground albedo as in Fig.9.The colors are different.The blue line corresponds to the red line in Fig.9and the red line corresponds to the black line in Fig.9.The green and the blue line are identical to the lines in Fig.9.The red line is not due to the substraction of the additive component.

Figure 12 .
Figure 12.Measured radiance to simulated single scattering ratio calculation for an MLT and co-located nominal limb measurements.The albedo factor is the product of AB (AB = 2.20).

Figure 13 .Figure 14 .
Figure 13.Latitude-altitude distribution of the monthly mean Na densities.The average of the results for both Na lines is used.Note that the highest covered latitude is at 82 • N/S.

Figure 15 .
Figure 15.Seasonal variation of the vertical Na density profile for low, middle and high latitudes.

Figure 16 .
Figure 16.Normalized vertical Na profile at selected latitudes in July.The densities are normalized to the peak value.

Figure 17 .Figure 18 .
Figure 17.Left: monthly mean Na number densities at 71 • N. Right: profiles of the left figure normalized by the peak densities of each month.

Figure 19 .
Figure 19.Mean and standard deviation (error bars) of the Na D 1 and D 2 line retrieval for the equatorial SCIAMACHY measurements on 20 March 2010.

Figure 21 .
Figure 21.Comparison of the SCIAMACHY VCD profile at 69 • N with lidar observations at Andøya, Norway.

Figure 22 .
Figure 22.Left: comparison of the SCIAMACHY VCD profile at 40 • N and the results from She et al. (2000).Right: comparison of the same profiles with profiles from Fussen et al. (2010) (their Fig. 7) and Marsh et al. (2013) (their Fig. 10).