Retrieval algorithm for densities of mesospheric and lower thermospheric metal atom and ion species from satellite-borne limb emission signals

Meteoroids bombard Earth’s atmosphere during its orbit around the Sun, depositing a highly varying and significant amount of matter into the thermosphere and mesosphere. The strength of the material source needs to be characterized and its impact on atmospheric chemistry assessed. In this study an algorithm for the retrieval of metal atom and ion number densities for a two-dimensional (latitude, altitude) grid is described and explained. Dayglow emission spectra of the mesosphere and lower thermosphere are used, which are obtained by passive satellite remote sensing with the SCIAMACHY (SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY) instrument on board Envisat. The limb scans cover the tangent altitude range from 50 to 150 km. Metal atoms and ions are strong emitters in this region and form sharply peaked layers with a FWHM (full width at half maximum) of several 10 km in the mesosphere and lower thermosphere measuring peak altitudes between 90 to 110 km. The emission signal is first separated from the background signal, arising from Rayleigh and Raman scattering of solar radiation by air molecules. A forward radiative transfer model calculating the slant column density (SCD) from a given vertical distribution was developed. This nonlinear model is inverted in an iterative procedure to yield the vertical profiles for the emitting species. Several constraints are applied to the solution for numerical stability reasons and to get physically reasonable solutions. The algorithm is applied to SCIAMACHY limb-emission observations for the retrieval of Mg and Mg + using emission signatures at 285.2 and 279.6/280.4 nm, respectively. Results are presented for these three lines as well as error estimations and sensitivity tests on different constraint strength and different separation approaches for the background signal.


Introduction
Meteoroids entering Earth's atmosphere have a large metal content.The meteoroids ablate on entry as a result of frictional heating in the mesopause region at around 80 to 100 km altitude.The altitude where the ablation occurs depends on the meteoroid's mass, velocity and entry angle and the boiling point of the specific metal species (see, e.g., McNeil et al., 1998;Vondrak et al., 2008 for more details).The ablated atomic metals may also become ionized by charge exchange with ionized oxygen, ionized nitrogen and NO + .This ionization of the neutral metals leads to the formation of ion layers at slightly higher altitudes than the neutral layers.A significant amount of the ablated metals is also ionized during the ablation process by hyperthermal conditions (see, e.g., Fig. 14 in Vondrak et al., 2008).As ozone is the main reactive species in this region, the metals react quickly with ozone and further reactions lead to stable reservoir species like metal hydroxides and carbonates.In a dynamic equilibrium comprising production from meteoroid Published by Copernicus Publications on behalf of the European Geosciences Union.

M. Langowski et al.: Retrieval algorithm for MLT metal densities
ablation and loss to stable molecules, the metals form sharply peaked layers at around 85 to 95 km and the ion layers peak 5 to 15 km above the neutral layer peaks.Comprehensive laboratory studies and model simulations to obtain the loss rates of the metals have been conducted, e.g., by Plane andHelmer (1995), McNeil et al. (1998), Fritzenwallner and Kopp (1998), Plane (2003) and Plane and Whalley (2012), and detailed descriptions of the chemical processes can be found in these studies.
The densities of metal species and their ions -corresponding only to a few thousand particles per cm 3 even in the peak region -are several orders of magnitude less than the densities of the main atmospheric constituents.However, the oscillator strengths of the transitions and subsequently the absorption cross sections for metals and their ions are large, compared to the main constituents.This is not the case for most of the formed reservoir species, which have not been detected by UV-visible remote sensing (beside FeO, Evans et al., 2010).As the Mg line at 285.2 nm and the Mg + lines at 279.6 and 280.4 nm have wavelengths below 300 nm, observations from the ground are not possible because of ozone absorption, and until the mid-1990s have been carried out only sporadically through rocket and satellite measurements (e.g., Kopp, 1997;Gardner et al., 1995;Gerard and Monfils, 1978;Joiner and Aikin, 1996;Mende et al., 1985).First global measurements of Mg and Mg + with daily global coverage have been conducted with GOME (Global Ozone Monitoring Experiment aboard ERS-2) in nadir mode as total column densities and were reported in Aikin et al. (2004) and Correira et al. (2008Correira et al. ( , 2010)).The first long-term measurements with global coverage in limb mode have been available since 2002 from SCIAMACHY on board Envisat.Investigations of this data set are reported in Scharringhausen et al. (2008a, b).However, the nominal limb mode of SCIA-MACHY only extends up to 92 km and therefore only gives a low vertical resolution in the peak region.Between 2008 and the end of the Envisat mission in April 2012, a special limb mode was introduced that provides scans up to 150 km tangent altitude to better study the mesosphere and lower thermosphere emissions.These measurements are used in the current study.The algorithm presented in this paper is based on the algorithm initially developed by Scharringhausen et al. (2008a), but reflects additional key issues found and solved, which in turn have had significant influence on the results.Note that a validation of the results is not possible, because there are no collocated measurements available.
In Sect. 2 the SCIAMACHY instrument and the data used are introduced.The main focus of this study is the description of the 2-D algorithm in Sect.3, and sample results are briefly discussed in Sect. 4 to demonstrate that the algorithm works.
2 The instrument SCIAMACHY/Envisat SCIAMACHY, the SCanning Imaging Absorption spectroMeter for Atmospheric CHartographY (Burrows et al., 1995;Bovensmann et al., 1999, and references therein) is a grating spectrometer on board the European environmental research satellite Envisat.Envisat was launched on 28 February 2002 into a Sun-synchronous orbit at an altitude of about 800 km and a descending node at 10:00 LT (local time).In late October 2010, an orbit maneuver took place -lowering the satellite altitude by about 17 km -in order to reduce the risk of collisions with other spacecraft.After performing very well for a decade and without any significant gaps in SCIAMACHY data coverage, contact to Envisat was lost, unexpectedly and without warning, on 8 April 2012 and the European Space Agency (ESA) officially announced the end of the mission on 9 May 2012.
The SCIAMACHY instrument comprises 8 spectral channels ranging in wavelength from 214 to 2386 nm.The FWHM (full width at half maximum) of the slit function is different for every channel.It is roughly 0.22 nm in channel 1 (214 to 334 nm) and the sampling rate is about 2, which corresponds to a sampling every 0.11 nm in channel 1.
SCIAMACHY has three different observation modes: limb, nadir and solar/lunar occultation.In the occultation mode the instrument points directly towards the Sun or the Moon.In nadir mode the instrument's line-of-sight (LOS) is directed downwards to Earth's atmosphere and surface, while in limb mode the instrument points tangentially to the Earth's surface at different tangent altitudes.The nadir mode has the advantage of providing total column information with high latitudinal and longitudinal resolution, but vertical profile information can typically not be retrieved.Limb observations provide vertical profile information with good vertical resolution.However, only between the lowest and the highest tangent altitude of a series of consecutive limb measurements -hereafter called a limb state -the vertical resolution is good.Altitude regions below the lowest tangent altitude cannot be accessed, and although altitudes above the highest tangent point contribute to the observed limb-scatter or limbemission signal, profile information for altitudes above the highest tangent point can generally only be retrieved with a much poorer resolution than for the region between the highest and the lowest tangent point.At the beginning of the SCIAMACHY mission in 2002, limb scans were performed from about 0 km up to about 105 km tangent altitude.However, the maximum tangent altitude was reduced to about 92 km in Fall 2003, in order to improve limb-nadir-matching.This means that over the mission lifetime from 2002 to 2012, limb measurements ranging in tangent altitude from 0 to about 90 km in 3.3 km steps are available.One additional limb measurement is performed at 250 km tangent altitude.At this tangent altitude no atmospheric signal is expected to be present (see, e.g., Gerard and Monfils, 1978) and the limb radiance spectra at this tangent altitude are subtracted from the other limb spectra of the same state as the dark signal correction.For the first three or four limb states after orbital sunrise -at high northern latitudes -an additional stray-light component in the spectra at 250 km tangent altitude occurs that has a similar spectral shape as the solar irradiance spectrum.The reason for this is most likely that the Sun is still within the total clear field of view (TCFOV) of the instrument.
Since metal species and their ions and other emitting species have their density maxima at altitudes around 90 km and higher, a new mesosphere and lower thermosphere limb state was introduced in September 2008 that provides limb measurements covering the 53.5 to 150 km tangent altitude range in steps of 3.3 km.The dark signal measurement is performed at a tangent altitude of about 350 km for these states.These states are performed every two weeks for 15 consecutive orbits -which corresponds to about a full day of observations, covering nearly all longitudes.
The flux of photons is small in channel 1, where the Mg and Mg + lines used in this study occur.Therefore, the signalto-noise ratio is quite low.Because of this, the Mg and Mg + retrievals presented here are based on daily-and zonallyaveraged as well as latitudinally-binned data.Geolocations of a reference orbit were used for the binning and the averaged data is furthermore smoothed with neighboring measurements within a ±10 • latitude range and a 2 h local time range.Typically this results in an average of roughly 20 single spectra for each latitude bin.
The SCIAMACHY data set employed is level 1 data version 7.03 and 7.04 and calibrated with ESA's calibration tool scial1c with all options (0, 1, 2, 4, 5, 6, 7 and M factors, which include option 3).In the data-averaging process, measurements with very low signal (night measurements) and measurements with anomalous spectral peaks that occur especially in the Southern Atlantic Anomaly (SAA) region, which can span 1/3 of all longitudes at its largest latitudinal extent around 30 • S, were excluded.On some days inexplicably high densities at high altitudes at 20 • S were retrieved and the retrieved densities of these days where not used for further data processing.

Outline
The retrieval algorithm introduced in this study is an improved version of the algorithm developed by Scharringhausen et al. (2008b).New features that are now implemented are a Ring effect (Grainger and Ring, 1962) correction, a treatment of self-absorption, the use of high resolution solar spectra for calculating emissivities, and the separation of night and dayside measurements at the same latitudes, which increases the number of limb measurements that can be used in the retrieval at polar latitudes.
The algorithm uses all limb measurements of an orbit to simultaneously retrieve densities of emitting species on a 2-D latitudinal and vertical grid.The algorithm can be separated into two steps.In the first step the pure spectral emission signal is separated from the background signal and slant column densities (SCDs) of the relevant species are retrieved.In the second step the inverse problem is solved to retrieve vertical density profiles from the SCD profiles determined in the first step.
Assuming single scattering of light from the Sun into the LOS of the instrument (emission path s e ), for a single measurement (one limb spectrum at one altitude) the forward model for the radiance can be calculated as follows: with emissivity γ , density n and an absorption part -along the LOS and the line from Sun (LFS) (s a stands for both absorption paths)f .In the first step of the retrieval the lefthand side of Eq. ( 1) is determined, while in the second step the right-hand side of Eq. ( 1) is inverted to obtain the density n.More detailed information on γ and f is provided in Sects.3.2.3 and 3.3.4,respectively.

Determination of the background signal in the UV
In order to obtain the pure emission signal, the background signal has to be subtracted from the observed limb radiance spectra.For the UV spectral range at wavelengths below 300 nm and tangent altitudes above 70 km, the background signal is mainly formed by single Rayleigh scattering of the incoming solar radiation as most of the multiply scattered light coming from lower altitudes is absorbed by the stratospheric ozone layer.The cross section for single Rayleigh scattering for a refractive index of n = 1 is given by the following equation: with the polarizability volume α and the wavelength λ.When the observed limb spectrum is divided by the solar irradiance spectrum, the resulting ratio spectrum is -ignoring terrestrial emission features -a smooth function, which can be well approximated by a straight line within narrow spectral windows, e.g., a 2 nm window around λ = 285 nm (see Fig. 1).Ozone is a strong absorber of UV-radiation (see, e.g., Gorshelev et al., 2013;Serdyuchenko et al., 2013, for absorption cross section data).For the Mg/Mg + line wavelength region, the main influence of ozone on the limb spectrum does not originate from ozone in the mesosphere, but from stratospheric ozone.Regarding the absorption of metal emissions by ozone, the optical depth for a limb measurement is only 0 Wavelength in nm Limb radiance/solar irradiance in sr −1 Fig. 1.Limb radiance divided by solar irradiance (average 2008(average -2012, Equator, tangent altitude 89.3 km), Equator, tangent altitude 89.3 km).Division of both spectra leads to a smooth background signal in contrast to emission lines, which are clearly identified.The background signal is separated from the emission lines through a linear fit.Note that relative to the strong Mg + at 279.6 and 280.4 nm and the strong Mg line at 285.2 nm, very weak NO lines can be observed at the longer wavelength edge of the Mg + lines at 281 nm.However, this influence is negligibly small.τ = 0.01 at 70 km tangent altitude and 280 nm wavelength, i.e., absorption by ozone is negligible for tangent altitudes above 70 km.For altitudes below 75 km the uncertainty introduced by the Ring effect correction described in Sect.3.2.2 is much more important than the absorption by ozone.
Stratospheric ozone, however, has a significant effect on the radiative transfer relevant for the retrieval problem, as it essentially absorbs incoming radiation entirely below 300 nm.Multiple scattering and surface reflection are thus negligible, and only single-scattering has to be considered.Figure 2 shows that the relative difference between multiple scattering and single scattering simulated with the SCI-ATRAN radiative transfer model (Rozanov et al., 2014) is less than 1 %.This leads to significant simplifications of the required radiative transfer model.
The Rayleigh scattering cross section as well as the ozone absorption cross section vary smoothly over the typical spectral width of the metal emission lines of only a few pm, which is -at least for most space borne instruments -typically much narrower than the width of the spectrometers's slit function (for SCIAMACHY about 0.22 nm at 280 nm).
To obtain the slant column emission (SCE) signal, the following steps, also illustrated in Fig. 3, are applied.First the Ring effect correction (see Sect. 3.2.2) is applied.The limb radiance spectra are divided by the solar irradiance spectrum (also measured with SCIAMACHY) so that the emission lines in the resulting spectrum can clearly be identified as sharp peaks in the smooth spectrum.The baseline of this of simulated spectra for multiple scattering for the wavelength region of the Mg/Mg + lines.The simulations where performed with SCIA-TRAN for an equatorial scenario.This and other scenarios with, e.g., an ozone layer of 150 Dobson units, show only minor differences of less than 1 % between single and multiple scattering.Multiple scattering thus can be neglected.
spectrum is subtracted and then the spectrum is multiplied with the solar irradiance spectrum again to get the emission spectrum without background.The SCE is determined as the factor that leads to the best fit with the emission spectrum when multiplied to the normalized shape function, which is the slit function of the spectrometer.Advantageous conditions for a good separation of background and emission signal are a good signal-to-noise ratio, strong emission signals compared to the background radiation and no overlap with nearby other emission signals.Other emission signals are, e.g., the NO γ (1, 6) band close to 280 nm.However, the emission from this band is negligibly small relative to the Mg + emission.

Ring effect correction
One of the biggest challenges in retrieving upper atmospheric metal emissions from satellite remote sensing of their dayglow is that most of these metals also have absorption and/or emission signatures in the solar irradiance spectrum.The Fraunhofer lines, some of which are very deep, lead to different complications.The background signal formed by Rayleigh and Raman scattering rises with higher densities, which occur at lower altitudes.Inelastic rotational Raman scattering leads to a filling-in of the Fraunhofer lines, leading from the pure Rayleigh scattered spectrum I 0 to the Rayleigh and Raman scattered spectrum I 1 .This filling in of Fraunhofer lines is known as the Ring effect (Grainger and Ring, 1962).As the biggest part of the scattered light is resonantly Rayleigh scattered, and the spectral redistribution from Raman scattering is therefore small compared to a pure Rayleigh scattered spectrum, the changes to the limb signal due to the Ring effect can be linearized.Therefore, the filling-in is simulated a second time, leading from the measured radiance spectrum I 1 to a double filled in spectrum I 2 .The difference between I 2 and I 1 is subtracted from the original limb spectrum I 1 to correct the Ring effect and obtain I 0 : Under the assumption I R (λ) = I R2 (λ) = I R1 (λ), it follows that To test whether this linearization approach works, the Ring effect operator (see below) is applied to I 0 , which should result in I 1 again.Figure 4 shows I 0 , I 1 , I 2 and I 1 recalculated from I 0 for a sample spectrum.Both the measured I 1 and I 1 recalculated from I 0 in Fig. 4 are nearly identical.The Ring effect operator used corresponds to the convolution of the spectrum with a boxcar function with 3.3 nm Ring effect correction for the Mg 285.2 nm line for a limb measurement with a tangent altitude of 53.5 km.The Ring effect is simulated, using the measured spectrum I 1 (λ), to obtain I 2 (λ).The filling-in effect is assumed to be small and differences between the first and the second application are therefore also very small.Under this assumption, the differences are nearly the same: I 0 (λ)−I 1 (λ) ≈ I 1 (λ) − I 2 (λ).And by adding I 1 (λ) the corrected spectrum I 0 (λ) is obtained as width (see Fig. 5), corresponding to 31 wavelength grid points (15 points left and 15 points right to the center wavelength), which has an additional peak at the center of the convolution function.This additional peak accounts for the Rayleigh scattering and accounts for approximately 96 % of the area of the convolution function, while the remaining 4 % originate from adjacent wavelengths.
Instead of a boxcar for the Raman part of the spectrum, a more correct distribution, e.g., from Penney et al. (1974), could be used.However, when smoothed to the resolution of SCIAMACHY, this distribution is very close to a boxcar function.Thus, the boxcar function is probably more than adequate for this purpose.Alternatively, a more accurate description of the radiative transfer including all effects described in Vountas et al. (1998) could be used, at the cost that the radiative transfer model would become more complicated.The true percentage of Raman scattered light may vary between 3 to 6 %.The effect of different Raman scattering percentages will also briefly be discussed in Sect.4.3.A small error is introduced to the emission signal, as not only the background signal but also the emission signal is smoothed by the Ring effect correction.For example, if there was a pure emission signal and the Ring correction was applied with 4 % Raman scattering, the I 2 signal at the peak wavelength would decrease by 4 % and the corrected signal I 0 would increase by 4 % at the peak wavelength.However, even in the peak region for Mg the ratio of emission to background radiation is close to 1, and the error introduced by the Ring correction there is negligibly small compared to other error sources, while at lower altitudes the background signal dominates.
For Mg this simple model leads to significantly decreased SCDs for tangent altitudes below 85 km, which is shown in Fig. 6.Result of the Ring effect correction for a tangent altitude of 53.5 km for the ratio I /F .Shown are the uncorrected spectrum and the corrected spectra with different percentages of Raman scattering in the applied method.The emission signals become smaller with this correction and the Mg line at 285.2 nm as well as the Si line at 288 nm, which cannot be observed in the region between 90 to 100 km and originates probably purely from Ring effect, nearly vanish for the Ring correction with 4 % of inelastic scattering in the background signal.Fig. 6, while measurements above 85 km shown in Fig. 7 are only weakly affected.The Mg + lines are generally less effected by the Ring effect because due to the double line structure the filling-in mainly comes only from one edge of the lines and not from both.Furthermore, the Mg + layer lies several km above the Mg layer, so the background signal is generally lower in the peak region than it is for Mg.The Ring effect explains the apparent rise of the Mg + emission signal below 70 km.However, it is better to not apply the Ring correction at all for the Mg + retrievals, which only use the region above 70 km, because the additional errors in the peak region, which is where emission dominates, introduced by the correction outweighs the benefit of the very small corrections at the lower peak edge (see Fig. 31).
Figure 8 shows the resulting Mg SCE profiles for approaches with different fractions of Raman scattering to the total scattering.It it is concluded that Mg profiles at altitudes below 70 km cannot reliably be retrieved with the Ring effect correction schemes used here, and even between 70 and 80 km the Ring effect correction may lead to non-negligible systematic retrieval errors.The effect of different chosen parameters for the Ring correction will be shown in Sect.4.3.
Spectra simulated with the latest version of SCIATRAN (see, e.g., Rozanov et al., 2005) from March 2013, which uses a more sophisticated forward model for the background radiation for limb geometry, were used for a comparison.The Ring correction, as described in this section, was applied to SCIATRAN simulations, including inelastic scattering.Equatorial latitudes and tangent altitude from 53 to 90 km were used.The resulting spectra were compared to the SCI-ATRAN simulation without inelastic scattering, which could be well reproduced.Above 70 km altitude the differences are small enough to justify the use of the method described in this section.
In the algorithm the Ring effect correction is performed before the background subtraction described in Sect.3.2.1.

γ factors
There are several constant factors on the right-hand side of Eq. ( 1), which is repeated here: γ being the emissivity, In the case of resonance fluorescence the metal atom absorbs a photon from the incoming irradiation field πF λ 0 and becomes excited from the initial state i to state j .As the Einstein coefficients for spontaneous emissions A j i are large for strong transitions, re-emission occurs immediately.Depending on the quantum number change from i to j , the scattering is isotropic or anisotropic, and the phase function P (θ ) is a superposition of an isotropic part and an anisotropic part.It is very similar to the Rayleigh scattering phase function.The phase function -which is normalized to 4π -is given by 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 wavelength-integrated cross section depends only on the transition wavelength λ ij and the oscillator strength f ij of the transition and can be directly calculated from database values.Here, data from the NIST atomic spectra database (Kramida et al., 2012) is used.Note that the integrated cross section has to be spectrally distributed over the shape of the line for mesospheric conditions, i.e., a Doppler-broadened Gaussian line shape for each individ-Table 1. E 1 and E 2 depend on the change of angular momentum j (from Chandrasekhar, 1960).
The cross section is multiplied by the solar irradiance πF λ 0 and the product is integrated to obtain the emissivity.If the metal density is sufficiently large, the emitted radiation can be partly absorbed along the LOS by the metal itself, and there is also absorption of radiation before the emission process along the LFS.This can also be interpreted as a reduced emissivity γf with the attenuation factor f .The shifts in the cross section profile, which are 250 to 500 smaller than the resolution of SCIAMACHY and therefore cannot be resolved, are important for the calculation of f in the emission and self-absorption process, which is discussed in Sect.3.3.4(see Fig. 18).
The metals found in Earth's upper atmosphere are also abundant in the solar atmosphere.Mg even has the same isotopic ratios in the solar atmosphere as in the terrestrial atmosphere (Boyer et al., 1971).The occurrence of metal species in the solar atmosphere leads to the formation of Fraunhofer lines in the solar spectrum.As the majority of Fraunhofer lines are narrow, having line widths in the pm range, a high resolution solar spectrum is required for a proper treatment of radiative transfer effects relevant for this study.SCIAMACHY's spectral resolution of about 0.22 nm in the 280-285 nm spectral range is too poor to resolve the individual atomic lines.Ground-based instruments can have the required spectral resolution but are incapable of observing the Mg/Mg + lines because the stratosphere is optically thick below about 300 nm.Therefore, the balloon-borne measurements of Anderson and Hall (1989) and Hall and Anderson (1991) carried out in 1978 and displayed in Fig. 9 are employed.The spectral resolution of these solar irradiance measurements is 0.01 nm.
Fortunately, this spectral resolution is already sufficient to resolve the Mg and the Mg + lines.However, the Mg + emission lines show an inner minimum that is not fully resolved.The solar Fraunhofer lines are also much broader than expected from pure Doppler broadening.This broadening is mainly pressure broadening.For later considerations in Sect.3.3.4a spectrally constant solar irradiance is assumed in the considered wavelength range in the mesosphere, since the lines are much narrower there.
Another feature that is not yet considered in the solar reference spectrum is the solar cycle.Since the beginning of the space age, the MgII index -i.e., the ratio between the irradiance of the chromospheric emission lines near the center of  (1991) (data from 1978) in comparison to the SCIAMACHY solar irradiance spectrum during high solar activity ( 2002) and low solar activity (2010).There are especially strong differences between the high resolution spectrum and the SCIAMACHY spectrum in the Fraunhofer lines at 279.6, 280.4,285.2 and 288.2 nm.When smoothed and scaled, the high resolution spectrum is very similar to the SCIAMACHY spectrum (see Fig. 10).
the MgII Fraunhofer line and the wings of the MgII Fraunhofer line -is one of the most commonly used solar proxies (Skupin et al., 2004;Snow et al., 2005).There is a variability of 20 % in the 11 yr cycle and similar variability of up to 10 % in the 27 day cycle.The solar cycle variation is considered in the following way.First the high resolution spectrum has to be scaled to the SCIAMACHY spectrum with a constant factor to consider different degradation effects of the instrument.This works best when a smoothed version of the high resolution spectrum is used (see, e.g., Fig. 10) to find the best fitting factor first.Now the actual effect of the variability is investigated.Therefore, the spectrum is scaled only in the center of the emission lines by a constant factor, which simulates high or low solar activity.Note that for the inversion step from the smoothed lines to the highly resolved lines, it is assumed that the shape of the emission lines, which is not resolved in the smoothed spectrum, is the same as in the high resolution spectrum, which may not be true.The following instrumental effects of SCIAMACHY are applied to the scaled and peak-scaled high resolution spectrum.The high resolution spectrum is smoothed by convolving it with the SCIA-MACHY channel 1 slit function, i.e., a hyperbolic function with a FWHM of 0.22 nm and a shape given by Eq. ( 9): Finally, the spectrum has to be sampled as the SCIAMACHY spectrum.This is done by interpolating the smoothed spectrum to the SCIAMACHY wavelength grid.This scaling, smoothing and sampling leads to a similarly shaped spectrum as that of the SCIAMACHY spectrum and is shown in Fig. 11 for Mg and in Fig. 12 for Mg + .Smoothing and interpolation only slightly affect the integrated value over each single line.Using different scaling factors for the emission line core in order to model the solar variability results in different integrated values.This change can be described by a linear function for each individual line.To consider the variations, this linear function is inverted for the integrated values of the daily measured SCIAMACHY spectrum to obtain the scaling factor that has to be applied to the core of the emission lines.By doing this, the high resolution spectrum is used to obtain the correct line shape of the lines and the daily SCIA-MACHY spectrum is used to correct for daily variabilities.To obtain the SCDs the SCEs have to be multiplied by 4π and divided by γ .

Discretization of the forward model
Before the forward model can be inverted it first has to be formulated.This starts with the discretization of Eq. (1).All constant factors in Eq. ( 1) are combined to a factor c 1 .The detected resonance fluorescence signal corresponds to the signal integrated along the LOS.Before the solar radiation is absorbed and re-emitted it follows a line from Sun (LFS) to the center of the grid cell.Along the LFS a certain fraction of the incoming radiation is already absorbed, and after the absorption and re-emission process into the LOS the emitted radiation is also partly absorbed.Absorption by the emitting species introduces nonlinearity into the forward model.As discussed in Sect.3.2.1,absorption by ozone along the LOS  2002) solar activity in comparison with the scaled high resolution reference spectrum.When the high resolution spectrum is smoothed and sampled as the SCIAMACHY spectrum, there are only small differences left, mainly resulting from different phases of the solar cycle for the different spectra.Note that a scaling factor of 0.9 instead of 0.83 is used for Mg in the retrieval because this fits better for the line center.2002) solar activity in comparison with the scaled high resolution reference spectrum.When the high resolution spectrum is smoothed and sampled as the SCIAMACHY spectrum, there are only small differences left between both, mainly resulting from different phases of the solar cycle for the different spectra.An adjustment for the solar phase is done by scaling the core of the emission lines.The adjustment has to be done for each line individually.Possible variations of the inner structure of the emission lines are not considered.
is negligible, and only self-absorption is considered in the forward model.
Let the forward model to be inverted be Kx = y, where y represents the vector of the SCDs.If there are n individual measurements, y is a n × 1 matrix.x also has to be a vector and is a m × 1 matrix.Then K has to be a n × m matrix.The solution x for the number densities on the latitude-altitude grid has to be a vector, so technically the a × b grid matrix (with a and b being the number of latitudes and altitudes, respectively) has to be stored into a m × 1 vector.
Three different pathways have to be considered (see Fig. 13).A pathway for the emission along the whole LOS.For each grid cell along the LOS, the path from the Sun to this grid cell, as well as the path from the grid cell to the satellite, is needed for the absorption along both paths.For every grid cell i there is the pathway along the LOS s LOS i through this grid cell and the matrix (the vector) of the pathway along the LOS, starting from the grid cell to the satellite s LOS gc i,j , and the grid cell-specific pathways along the LFS s LFS gc i,j .To clarify this, the difference of s LOS i and s LOS gc i,j is that the first one includes the path length of all grid cells along the LOS from the satellite to the opposite side of the satellite, while for the second one only the grid cells j between the grid cell i and the satellite are considered as all path lengths in grid cells j along the LOS behind the grid cell i are 0. The two special cases, j = i and grid cell i that includes the tangent point are also considered.
The ray tracing for s LOS i has only to be done once, while the ray tracing for s LFS gc i,j and s LOS gc i,j has to be done once for every grid point.The discretized formula is as follows: x j s LFSgc i,j .(10) The argument of f is summarized to The forward model for the measured radiance can be expressed through a matrix equation Kx = y, with the state vector x containing the metal species densities to be retrieved for all latitudes and altitudes of the model atmosphere.The measurement vector y contains the measured radiances for all limb scans of one orbit, at the different tangent altitudes and latitudes.K is the matrix representing the forward model.For a linear problem K is independent of x and therefore K is identical to the Jacobian J (J = K) of Kx.However, only the emission part is linear; the absorption part is nonlinear.To find the least squares of the residuum between model and measurements, the equation system J T Jx = J T y has to be solved and therefore J is required.For a single wavelength the absorption part is exponential.However, several different wavelengths are involved so the exact calculation is a mixture of different exponentials.Therefore, the exact function f should be used, which depends on the optical depth along the LOS and LFS.The function f and its first derivative f with respect to optical depth can be derived numerically (shown in Sect.3.3.4).Both are quite smooth and monotonic.To get J the derivative of Eq. (10) has to be formed.

Formation of the Jacobian J
Considerations start with how one element of the Jacobian J is calculated, as J is built element-wise from this.To calculate one element the derivative of Eq. ( 10) has to be formed.This is shown in the Appendix and leads to Eq. ( 12): Note that if f would be independent of x, then f = ∂f (g i (x)) would vanish and with it the second term.This finally would lead to K = J.However, as f depends on x, a slightly different matrix has to be formed.However, for this special mathematical problem, omitting the second addend in Eq. ( 12) reproduces synthetic model density profiles that are forward modeled with the forward model and then again inverted with the retrieval algorithm very well.This is easier and faster than solving the full equation.However, K still depends on x.This and further consequences of the dependency of f on x are discussed in Sect.3.3.5.Now J can be built element-wise with Eq. ( 12).The biggest part of the path lies in the tangent point altitude region, but higher altitudes are also passed.Note that for the contribution of a grid element's emission to the total emission along the LOS, the path length in the grid interval as well as the density of the emitters in the grid interval and the absorption of the emitted light along the LOS have to be considered.For measurements with tangent altitudes far below the metal layer peak, where the density is far smaller than in the peak region, the emission signal mainly comes from the region above the tangent altitude.

Calculation of path matrices
To solve Eq. ( 12) the path matrices s LOS i , s LOS gc i,j and s LFS gc i,j introduced in Sect.3.3.1 are needed.These matrices are built once during the initialization of the program and are later used together with the initial density profile (introduced in Sect.3.3.5) of the last iteration step to calculate g i (x) in Eq. ( 12).
The paths are calculated on a 2-D grid of altitudes and latitudes.To overcome ambiguity at the poles, additional 82 • of latitudes are added at each pole to separate dayside and nightside latitudes.Latitudes higher than 82 • are ignored in both hemispheres as this is the highest latitude covered by SCIAMACHY measurements.
An example plot for path lengths in different grid cells for a typical LOS is shown in Fig. 14.All altitude intervals above the tangent altitude are passed twice by the LOS (from the satellite's point of view, downwards and upwards).To overcome ambiguities of different non-connected parts of the LOS that are, e.g., several 100 km apart, but are within the same altitude and latitude interval, the LOS is separated at the tangent point into 2 parts.This separation is necessary, because when there is a high density region between both parts of the LOS within the same grid cell, then the optical depth and with it the absorption specific function f and f are quite different for both parts.So the absorption relevant matrices are built up separately for each side of the tangent point of the LOS.However, in the end after building up the separate Jacobians for both parts of the LOS, both Jacobians are added.
For each side on the LOS, s LOS i has to be computed once for each measurement and s LOS gc i,j (only the elements of s LOS i that are closer to the satellite than the grid cell i) and can be quickly built from s LOS i by copying the valid nonzero elements.For the far side of the LOS matrices, s LOS i for both parts of the LOS are needed.
There are only minor differences in the calculation of s LFS gc i,j to the calculation of s LOS i , e.g., the additional finding of the tangent point of the straight line, which contains the LFS, and to find out whether this tangent point is also a part of the LFS.For SZA < 75 • , LFS contributions are negligible.For SZA ≈ 90 • , the LFS contribution is as strong as the LOS contribution.Also, measurements with a too high SZA and with a tangent point of the LFS in the lower atmosphere have to be excluded as the model does not include lower atmospheric contributions.
For a vertical grid, right-angled triangle algebra is applied, using the maximum altitude of the grid cell as hypotenuse c and the tangent height as adjacent side b to derive the path lengths in each vertical grid level, as illustrated in Fig. 15.Changes of latitude within a vertical grid level can be considered by combining the straight line equation for the LOS with the cone equation of the latitude cones.More details for this will be shown in Langowski (2014).The initial incoming light is partly absorbed along the line from Sun (LFS) before the absorption and re-emission process into the lineof-sight (LOS).After this it is further absorbed along the LOS.Only a small part of the light on the LFS is absorbed and emitted into the LOS.However, this just leads to a constant scaling for the remaining light before and after the resonance fluorescence process and the shape of the spectrum stays the same.Therefore, the parts of g from the LOS and the LFS may not be treated separately.The initially as constant approximated irradiance profile is absorbed where the cross section is high.This leads to a spectral hole burning.

Calculation of self-absorption contribution f and f
As pointed out in Sect.3.2.3, the emission depends on the product of the incoming solar irradiation and the absorption cross section.The same applies for the absorption.However, while for the emission the spectrally integrated values of both quantities are sufficient, proper consideration of the spectral variations is crucial for the treatment of the absorption part.
For monochromatic radiation, the exponential Beer-Lambert absorption law I = I 0 e −σ nds can be applied with a path independent absorption cross section σ in the exponent.The density n is integrated along the absorption path s.However, because of the relatively low spectral resolution of SCIA-MACHY, the spectral radiance measured is integrated over a certain spectral range, which is essentially determined by the width of the instrument function.Therefore, the Beer-Lambert Law is simulated for monochromatic radiation and integrated for different g = nds (Eq.11).The emission is large for wavelengths with large absorption cross sections.However, a larger absorption cross section also means more absorption.Therefore, the metallic layer becomes optically thick for lower g in the center of the line compared to the wings of the line.This leads to a "hole-burning" effect in  16).To obtain the emissivity of the spectral line, this spectrum has to be integrated over all wavelengths.
For strong absorption the light comes mainly from the edges of the line instead of the line center.This reduces the effective cross section for the total integrated profile.
the product of cross section σ , irradiance π F and absorption attenuation e −σ nds .This is illustrated in Figs.16 and 17.
For small nds one obtains an approximate Beer-Lambert Law with a constant effective cross section.Hunten (1956) showed that, for a Gaussian-shaped cross section profile, the effective cross section is 1 √ 2 times the cross section at the maximum of the Gaussian.For larger nds, the effective cross section becomes smaller as the remaining light comes more and more from the edges of the absorption cross section profile.f is determined numerically using Eq. ( 13): f is obtained by numerical differentiation of f .For Mg and Mg + , isotopic shifts reduce the cross section in the center of the line (Fig. 18) and therefore lead to reduced selfabsorption.Results of the numerical calculation of f are shown in Fig. 19.For a ray starting from the satellite along the LOS through a strongly absorbing layer, the f factors along the LOS decrease rapidly.For the total emission from the LOS, the single emissions of each line segment have to be integrated.Because f decreases more rapidly for higher densities, this leads to a compression, and for a monochromatic line even to a saturation effect for the conversion of true SCD g = nds to the measured SCD ( f nds).This is illustrated in Fig. 20.For typically obtained values of SCDs in the order With isotope shifts Without isotope shifts Fig. 18.Absorption cross section for the Mg 285.2 nm line with and without considering isotopic shifts.The integrated cross section is the same for both cases.However, through the isotopic shifts for the heavier isotopes 25 Mg and 26 Mg, a part of the cross section is shifted from the center of the line to the edge of the line, which effectively reduces the absorption cross section for the integrated line.
of 0.5×10 11 cm −2 for Mg, this can lead to an issue with measurement noise (see also Fig. 24).

Linearization and iteration
The system of linear equations J T Jx = J T y has to be solved.This equation system is not linear as J = J(x) because of the self-absorption contributions in f and f .This problem is overcome by an iterative approach using J = J(x i−1 ) to calculate x i with starting values x 0 for the first iteration step.It is assumed that the metal layer is not opaque because then the SCD profile would be viewing angle independent and therefore almost constant below a certain tangent altitude, which is not the case (see, e.g., Fig. 8).Under this assumption, starting with x 0 = 0 is a good choice because when choosing x 0 = 0 the Jacobian is the same, as if self-absorption was not considered at all, and a profile shape is obtained that is already close to the final solution after the first step.This first step solution can be interpreted as the actual starting condition, which is close enough to the final result so that the iteration converges.
In order to investigate how many iterations typically are required to achieve convergence, up to 50 iterations were run for typical profiles.Convergence is already achieved quickly.After 20 iterations, which are used to obtain the results, the largest difference for an individual grid cell from step 19 to step 20 is less than 1 %.The changes in most of the other grid cells are at least a factor 10 to 100 smaller.

Constraints
Minimizing only the sum of the squares of the differences between the forward model Kx and the SCD profiles y typically leads to strongly oscillating solutions x.Therefore, constraints are introduced to stabilize the retrieval and to obtain more realistic results by reducing the oscillations.An a priori constraint (this is also called Tikhonov regularization Tikhonov, 1943) and smoothing constraints for neighboring latitude and altitude grid cells are used.The constraints have to be chosen as strong as required to stabilize the solution, but as weak as possible to not influence the solution too strongly.The final equation that has to be solved is Eq. ( 14): with the Jacobian J of the forward model Kx.The a priori covariance matrix S a is in fact a scalar (λ apriori ) multiplied with an identity matrix.S H and S φ are the matrices for altitudinal and latitudinal constraints (large sparse matrices with only 2 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 SCDs (S y ), which is assumed to be diagonal: the SCDs y and the a priori solution x a .As there should not be any bias on the form of the profile, x a = 0 is used.The effect of choosing different vertical smoothing constraints is demonstrated in Sect.4.2.

Results
In this section sample results are briefly presented to demonstrate that the algorithm works well.As most of the discussed parameters and effects only have an influence on the vertical profile, only results for vertical profiles in a limited equatorial latitude range are shown.

Equatorial vertical retrieval results and error estimations
Figures 21-23 show the results of vertical profile retrievals of Mg at 285.2 nm and Mg + at 279.6/280.4nm for the equatorial region (10 • S to 10 • N).The red curve in these figures shows the retrieval result using the mean equatorial SCD profile from 2008-2012 as a test profile.The Mg profile peaks at 90 km and has a FWHM of ≈ 15 km.The Mg + profile peaks at 96 km and has a FWHM of ≈ 12 km.Both Mg + lines show similar peak values.However, for small peak values the more weakly absorbing 280.4 nm line shows slightly smaller peak values, while for regions with high peak values (3000-4000 cm −3 ), the stronger absorbing 279.6 nm line shows the higher values.However, the largest differences are smaller than 25 % in the peak region.Attenuation f f with isotope shift exponential approximation for small g f without isotope shift f with isotope shift for T=170 K f with isotope shift for T=250 K Fig. 19.Calculated attenuation f as a function of g for the Mg 285.2 nm line.Taking isotopic shifts into account reduces the absorption effect and f is higher for this case.The black lines show variations in the absorption behavior in the range of typical temperatures in the 90 km to 100 km altitude range, resulting in a different Doppler broadening of the absorption cross section profiles.
To calculate the Mg/Mg + density errors the mean error of the SCE of the line is taken, corresponding to roughly 1 × 10 8 ph cm −2 s −1 sr −1 for single measurements.It is half as large as the highest SCE for Mg and 1/4 (280.4 nm) and 1/8 (280.4 nm) as large as the SCE for Mg + in the peak region.To obtain the SCDs the SCEs have to be multiplied by 4π and divided by γ .The same applies to the errors.For roughly 20 measurements per day (see Sect. 2), the single measurement error is divided by √ 20 ≈ 4.5 to derive the daily error.Note that the daily error for Mg is roughly as large as the single measurement error for Mg + at 279.6 nm.
A Monte Carlo method is used to propagate the radiance errors to the retrieved Mg/Mg + density profiles.The Gaussian error is applied to the mean SCD profile, and then the density profile is retrieved.This is repeated often enough (1000 times) so that the mean and the standard deviation of the results converge, and the standard deviation is interpreted as the error of the retrieval.The blue line in Figs.21-24 shows the mean of these Monte Carlo runs with 1σ error bars.Note that the error of the mean is the standard deviation divided by a factor √ 1000 ≈ 32.For Mg the relative errors for daily averages in the peak region are 20-30 %.The mean of the retrievals with additional errors is up to 100 cm −3 higher than the retrieval without errors.This is explained by the nonlinearity of the forward model (see Sect. 3.3.4). Figure 24 shows the errors for the Mg retrieval using single measurements, and for this case the mean is shifted significantly between both methods.This leads to a systematic difference between averaging the SCDs before the retrieval and averaging densities after the retrieval.Attenuation using f Exp.approx.f. small g Fig. 20.Plot of attenuated slant column densities (SCDs) versus the true SCDs for the Mg 285.2 nm line.The attenuated SCDs are determined using f , as well as the exponential approximation of f for small true SCDs.Typical measured SCDs (y axis) in the peak region for Mg are 6-8×10 10 cm −2 .The exponential approximation cannot be used in the retrieval, since it is already saturated below 8 × 10 10 cm −2 .For Mg the self-absorption effect is so strong that the measured SCDs are a factor 3 to 4 smaller than the true SCDs in the Peak region.For Mg + (not in the plot) the conversion factor is close to 1.
Therefore, an averaging must be applied to the SCDs before the retrieval.
On the other hand the true natural variability has to be taken into account, which is still there, even if perfect measurements without errors are used.Like the variability resulting from errors, the true variability leads to an increase of the mean values.Furthermore, the longer the time span of averaged SCEs is, the more the boundary conditions, like scattering angles and solar zenith angles, which are needed for the retrieval, change.Therefore, a compromise between averaging enough data before the retrieval step so that shifts due to errors are excluded, and averaging enough individual results after the retrieval step to account for the natural variability, is made to get the mean density profile at the Equator.As there is one day of measurements available roughly every 14 days, daily averages of the spectra are formed, and densities are retrieved from these daily averages.These densities are further averaged, e.g., to monthly averaged results, to reduce the errors.
The Mg density errors are only small in the peak region, while below 80 km and above 100 km the relative error is bigger than 100 %.However, for the 4 yr average (average of all spectra to obtain one SCD profile before retrieving densities), the errors are negligibly small, so that the small peaks at 113 km and between 130-140 km may also be real and, e.g., could originate from sporadic layers.Furthermore, Mg + also shows small peaks at the same altitudes.To estimate the error, the retrieval is repeated 1000 times with a typical random Gaussian error for daily averaged data applied to the initial slant column densities.The blue line shows the mean result of the 1000 runs, while the error bars show the standard deviation from the mean values of the 1000 runs.For daily averaged data (which are used here) the differences between both methods are small.For Mg + the relative errors in the peak region are less than 20 %.The 1σ error is close to 100 % below 85 km, but although the error is also large above 105 km, the density is still larger than 50 cm −3 within the error limits and not zero at the top altitude of 150 km.
Although the same constraints are applied for both Mg + lines, the results for the 280.4 nm line oscillate more strongly above the peak.And while the peak at 113 km occurs for both lines, the higher peaks are at different altitudes.However,   Mg (285.2 nm).The same methods as in Fig. 21 are applied.However, here the single measurements error is used.As for high relative errors the error propagation cannot be approximated to be linear and nonlinearities in the forward model lead to a shift of the mean to higher values.Therefore, a certain averaging of the data is needed before applying the retrieval algorithm.Note that the error bars estimate the error of single measurements while the error for the blue line is a factor √ 1000 ≈ 32 smaller.
both Mg + profiles are in good agreement within the error range.

Influences of different constraint strengths
As mentioned in Sect.3.2.2,different constraints are used.
To be effective these constraints should influence but not dominate the retrieval.The a priori and vertical-smoothnessconstraint have a similar effect so both must be tuned   (285.2 nm) in the equatorial region from 15 • S to 15 • N.For too low constraint strength (λ h < 5 × 10 −9 ) the profile is strongly oscillating, while for too high constraint strength (λ h > 5 × 10 −7 ) the profile is too smooth.Acceptable solutions between these two extrema show peak densities from 500 to 2000 cm −3 .However, at least the peak altitude is nearly independent of λ h .The solution with λ h = 5 × 10 −8 is favored, considering, e.g., seasonal variations etc. together.To simplify the search for an optimal constraint, a fixed ratio of λ h : λ φ : λ apriori (see Eq. 14, where λ apriori influences S a ) of 10 : 2 : 1 is used, as this was empirically found to yield realistic solutions.To investigate the sensitivity of the retrieval to the constraint parameters, many different combinations with the fixed ratio were tested, ranging from solutions that are obviously too heavily smoothed, to solutions that are obviously too strongly oscillating.show results for the vertical equatorial profile for Mg/Mg + for one example day.From this, as well as from latitude-altitude density plots, and considering other boundary conditions, which cannot be well quantified (e.g., similarity of results for same months, estimation of the difference between retrieving densities from the 4 yr averaged data and averaging results from single days, seasonal variations, smoothness of the result compared with the resolution of at least 3.3 km steps, etc.), the range of acceptable parameters for λ h for both species is estimated to be between 5 × 10 −7 and 5×10 −9 .For further studies 5×10 −8 will be used for Mg and 1 × 10 −7 will be used for Mg + the optimal solution.Within the range of acceptable parameters the uncertainty of the peak value can be estimated to be ≈ 30 % for Mg + and ≈ 50 % for Mg.

Altitude in km
One of the reasons for the fuzzy transition between the too smooth and too oscillating solutions can once again be found in the errors.Figure 28 shows the sensitivity of the retrieval to the constraint strength using the 4 yr averaged data, while Fig. 29 shows the results if an error is applied to the same data the same way as discussed in Sect.  in the equatorial region from 15 • S to 15 • N. The same approach to find the optimal constraint strength for Mg (see Fig. 25) is applied for Mg + .However, the variation of the peak density with constraint strength, especially for low strength, is much weaker, so errors from a wrong choice of the constraint strength are smaller.For Mg + the solution with λ h = 1 × 10 −7 will be used for further investigations.figures show results for Mg).For simplicity a 1-D pure vertical model has been used for both figures, which is very similar to the 2-D model.

Altitude in km
The error-free solution changes from the strongly oscillating case to the optimal solution, within a narrow constraint parameter range.The vertical constraint strength can be changed by more than 10 orders of magnitude without significantly changing the retrieval result.The solution is smoothed only if the smoothing constraint becomes even stronger.Therefore, it is easy to find the optimal constraint strength with error-free measurements.If an error is applied

Altitude in km
Fig. 28.Retrieval results with a 1-D model for the error free 4 yr averaged equatorial profile for Mg (comparable with the red line solution in the plots in Sect.4.1).Different vertical constraints are chosen and λ h = 10×λ apriori .With increasing constraint strength there is a sharp transition from a totally oscillating solution to a smooth solution.This result is stable within 12 orders of magnitude in the constraints strength until the result is smoothed too strongly.Note that the numerical value of the constraints is slightly different here than for the 2-D model.
-here errors for daily averaged data are used again -the situation becomes more complicated.The over-smoothed solutions are also retrieved for the same strong constraint as in the case without errors.However, the solutions stop oscillating only for much larger constraint strength compared to the error-free case.And furthermore the result depends more sensitively on the constraint strength than in the error-free case.However, note that this averaged result still looks much smoother than the single day profiles, as not only the constraints but also the averaging process in the results lead to smoother profiles.

Influence of different Raman scattering percentage on the profile retrievals
As discussed in Sect.3.2.2, the Ring effect has significant influence on the Mg profile retrievals for altitudes below 80 km.
Figures 30 and 31 show the retrieved vertical profiles for Mg and Mg + with a different contribution of Raman scattering to the total scattering (Raman and Rayleigh).
The Raman scattering percentage is assumed to be ≈ 4 %.The results for Mg still show seasonal and latitudinal variations for the densities at low altitudes, which are not expected to be there (not shown here).However, the densities at low altitudes are zero (within the error limit) for Raman scattering contribution of 3-6 %.The profiles with no or only very weak Ring effect correction still show significant densities at 70 km and only the constraints force the profile near the lower peak edge not to diverge.The densities

Altitude in km
Fig. 29.Retrieval results with a 1-D model when an error is applied for Mg with the same method as for the blue line in Fig. 21.The constraint strength is varied.Oscillating results start at a much higher constraint strength than for the error-free case shown in Fig. 28.Furthermore, the result varies much stronger with changing constraint strength, so that it is once again difficult to find the constraint strength for the optimal solution.For the case of very strong smoothing, the same profiles as in the error-free case are retrieved for the same constraint strength.
above 95 km are nearly independent of the Ring correction.The peak value is higher for lower Raman-scattering contribution.However, the variations between the most extreme expected contributions of 3-6 % are less than 100 cm −3 , so changing the Raman contribution by 1 % changes the peak value by only 2 %.The Mg peak altitude is slightly lower for lower Raman-scattering contributions.However, within the range of expected Raman-scattering contributions the peak altitude does not change by more than 1 km.Below the peak the differences are bigger, but down to a certain altitude the descent is parallel for all of the considerable Raman scattering contributions.The differences are largest below 80 km, as different Ring corrections lead to different residual densities at the lowermost altitude.
For Mg + the Ring contribution above 70 km is small, as the double line structure leads only to a filling-in from one side of the absorption line.Furthermore, the Mg + lines are not as deep as the Mg line (see Fig. 9) and the Mg + lines have emission lines in the core.The lower and upper Peak edges are similar.The only differences for retrievals with different Ring effect corrections is the peak value, which is higher for more Raman scattering contribution, and the densities at the lowermost altitudes that are slightly higher for lower Raman scattering contribution.
The Mg + peak is at a higher altitude than the Mg peak, so that the signal in the peak region is mainly dominated by pure emission.As discussed in Sect.3.2.2, this leads to slightly larger values in the peak region.And since the corrections above 70 km are negligibly small, it is better not to use the Ring effect correction at all for Mg + .

Conclusions and outlook
An algorithm to retrieve 2-D (latitude-altitude) number density fields for Mg and Mg + from satellite limb measurements of the scattered electromagnetic radiation in the wavelength region 275-290 nm in the mesosphere and lower thermosphere between 70 and 150 km altitude was presented.The approach and all relevant issues were described and explained in detail.The derivation of errors and sensitivity studies with respect to critical parameters were shown.The performance of this algorithm under realistic conditions was demonstrated and different issues that can be investigated to further improve the algorithm were pointed out.
This algorithm is used to retrieve Mg and Mg + from the SCIAMACHY data set of limb mesosphere-thermosphere scans from 2008 to 2012.As there is more data available from SCIAMACHY for nominal scans up to 90 km, the 2008 to 2012 mesospheric and lower thermospheric (MLT) data set will be used as a priori information for the 2008 to 2012 nominal data set to fill the temporal gaps between days with MLT scans, or it can be used to reduce the noise in the signals below 90 km through averaging.Furthermore, this a priori can be used to expand the data set to the full time span of SCIA-MACHY measurements from 2002 to 2012 to create a data set which will be used to study the geophysical behavior of the metals and their ions over the past decade.Additionally, for different percentages of Raman scattering in the Ring effect correction.For altitudes above 80 km the Ring effect correction is negligible and it is very small above 70 km.In the peak region the emission dominates the background signal.Because of this the peak value rises with higher Raman scattering percentage (as discussed in Sect.3.2.2).
the algorithm can be applied to data from similar satellite instruments as well as to other UV emission lines.The algorithm can also be applied to emission lines in the visible region, if further assumptions are made.This has already been tested for, e.g., the Na D-lines, and has led to promising results.
It was demonstrated that this algorithm works well with the available data set and there will be a comprehensive analysis of the available data set in a separate publication.

Formation of the Jacobi matrix
The Jacobian J is built element-wise.To calculate one element the derivative of Eq. ( 10) has to be formed.
Here the sum rule is applied to pull the derivation into the sum: The constants are pulled before the derivation: The product rule is applied: (A7) As f is derived numerically and f is a smoothing function, ∂f (g i (x)) ∂g i (x) can also be derived numerically and for simplicity f = ∂f (g i (x)) ∂g i (x) is introduced.(A8) For the derivative of g the sum rule is applied to pull the derivation into the sums and get δ distributions.Then the sums vanish: P i s LOS i x i f (g i (x)) ( s LOSgc i,k + s LFSgc i,k ).(A10)

Fig. 3 .
Fig. 3. Slant column emission (SCE) determination.First the Ring effect correction is applied to the limb radiance I (top left panel).The limb radiance is divided by the solar irradiance F and the wavelength region left and right to the line are used for a linear fit of the background radiation (top right panel).The linear fit is subtracted from the ratio I /F (bottom left panel) and the spectrum is again multiplied with the solar irradiance F (bottom right panel).Finally the area of the emission line is fitted with the slit function of the instrument.
Fig. 4.Ring effect correction for the Mg 285.2 nm line for a limb measurement with a tangent altitude of 53.5 km.The Ring effect is simulated, using the measured spectrum I 1 (λ), to obtain I 2 (λ).The filling-in effect is assumed to be small and differences between the first and the second application are therefore also very small.Under this assumption, the differences are nearly the same: I 0 (λ)−I 1 (λ) ≈ I 1 (λ) − I 2 (λ).And by adding I 1 (λ) the corrected spectrum I 0 (λ) is obtained as I 0 (λ) = 2I 1 (λ) − I 2 (λ).

Fig. 5 .
Fig. 5. Ring effect smoothing function.96 % of the background light is Rayleigh scattered.The Raman scattered part is approximated with a boxcar for the remaining 4 %.

Fig. 7 .
Fig. 7. Result of the Ring effect correction for a tangent altitude of 90 km for the ratio I /F .The background signal is low at this altitude and the Ring effect correction is only very small.

Fig. 9 .
Fig.9.High resolution solar spectrum fromHall and Anderson  (1991) (data from 1978)  in comparison to the SCIAMACHY solar irradiance spectrum during high solar activity (2002) and low solar activity (2010).There are especially strong differences between the high resolution spectrum and the SCIAMACHY spectrum in the Fraunhofer lines at 279.6, 280.4,285.2 and 288.2 nm.When smoothed and scaled, the high resolution spectrum is very similar to the SCIAMACHY spectrum (see Fig.10).

Fig. 10 .
Fig.10.The high resolution spectrum from Fig.9is smoothed and scaled for the SCIAMACHY spectrum by a factor of 0.83, so that the peak edges of the spectra are in agreement.
Fig. 11.SCIAMACHY solar spectrum in the Mg 285.2 nm line region for low (2010) and high (2002) solar activity in comparison with the scaled high resolution reference spectrum.When the high resolution spectrum is smoothed and sampled as the SCIAMACHY spectrum, there are only small differences left, mainly resulting from different phases of the solar cycle for the different spectra.Note that a scaling factor of 0.9 instead of 0.83 is used for Mg in the retrieval because this fits better for the line center.
Fig. 12. SCIAMACHY solar spectrum in the Mg + double line region for low (2010) and high (2002) solar activity in comparison with the scaled high resolution reference spectrum.When the high resolution spectrum is smoothed and sampled as the SCIAMACHY spectrum, there are only small differences left between both, mainly resulting from different phases of the solar cycle for the different spectra.An adjustment for the solar phase is done by scaling the core of the emission lines.The adjustment has to be done for each line individually.Possible variations of the inner structure of the emission lines are not considered.

Fig. 13 .
Fig. 13.Scheme of the altitude (radius) and latitude (sectors) grid and the line-of-sight (LOS) and line from Sun (LFS).Each grid element may have a different LFS.Note that this figure shows only a specific 2-D projection of the atmosphere along a constant meridian so that latitude sectors are equally sized.Each grid cell on the LOS has its own matrices for the LFS from the Sun to the grid cell and the LOS from the grid cell to the satellite, which are needed for absorption calculations.

Fig. 14 .
Fig. 14.Path lengths in different grid cells for a typical line-ofsight (LOS) of a limb measurement.The biggest part of the path lies in the tangent point altitude region, but higher altitudes are also passed.Note that for the contribution of a grid element's emission to the total emission along the LOS, the path length in the grid interval as well as the density of the emitters in the grid interval and the absorption of the emitted light along the LOS have to be considered.For measurements with tangent altitudes far below the metal layer peak, where the density is far smaller than in the peak region, the emission signal mainly comes from the region above the tangent altitude.
Fig. 15.2-D 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 in a vertical grid cell are added as additional sides a (red crosses).The path length in each grid cell is the difference of the sides a of neighboring grid cells.Note that, depending on the binning of the latitudes, it is possible that all grid cells only have one latitude.

Fig. 16 .
Fig.16.Absorption cross section for the Mg 285.2 nm line (black line and left-hand ordinate), which has to be multiplied with the attenuated solar irradiance (other lines and right-hand ordinate) to obtain the wavelength-specific emissivities shown in Fig.17.The initial incoming light is partly absorbed along the line from Sun (LFS) before the absorption and re-emission process into the lineof-sight (LOS).After this it is further absorbed along the LOS.Only a small part of the light on the LFS is absorbed and emitted into the LOS.However, this just leads to a constant scaling for the remaining light before and after the resonance fluorescence process and the shape of the spectrum stays the same.Therefore, the parts of g from the LOS and the LFS may not be treated separately.The initially as constant approximated irradiance profile is absorbed where the cross section is high.This leads to a spectral hole burning.
Fig.16.Absorption cross section for the Mg 285.2 nm line (black line and left-hand ordinate), which has to be multiplied with the attenuated solar irradiance (other lines and right-hand ordinate) to obtain the wavelength-specific emissivities shown in Fig.17.The initial incoming light is partly absorbed along the line from Sun (LFS) before the absorption and re-emission process into the lineof-sight (LOS).After this it is further absorbed along the LOS.Only a small part of the light on the LFS is absorbed and emitted into the LOS.However, this just leads to a constant scaling for the remaining light before and after the resonance fluorescence process and the shape of the spectrum stays the same.Therefore, the parts of g from the LOS and the LFS may not be treated separately.The initially as constant approximated irradiance profile is absorbed where the cross section is high.This leads to a spectral hole burning.

Fig. 17 .
Fig.17.Calculation of the wavelength-specific emissivity for the Mg 285.2 nm line as a function of wavelength difference from the center of line.It is the product of the solar irradiance, attenuated along the line from Sun and the line-of-sight, and the absorption cross section (both shown in Fig.16).To obtain the emissivity of the spectral line, this spectrum has to be integrated over all wavelengths.For strong absorption the light comes mainly from the edges of the line instead of the line center.This reduces the effective cross section for the total integrated profile.

Fig. 21 .
Fig.21.Plot of retrieved Mg (285.2 nm) for the 4 yr averaged equatorial data set (red).To estimate the error, the retrieval is repeated 1000 times with a typical random Gaussian error for daily averaged data applied to the initial slant column densities.The blue line shows the mean result of the 1000 runs, while the error bars show the standard deviation from the mean values of the 1000 runs.For daily averaged data (which are used here) the differences between both methods are small.

Fig. 22 .
Fig. 22. Mg + (279.6 nm) retrieval result.The same methods as in Fig. 21 are applied.The errors for Mg + are smaller than for Mg and furthermore nonlinear self-absorption effects are smaller for Mg + .Therefore, both methods lead to the same mean results.

Fig. 23 .
Fig. 23.Mg + (280.4 nm) retrieval result.The same methods as in Fig. 21 are applied.The errors for Mg + are smaller than for Mg and furthermore nonlinear self-absorption effects are smaller for Mg + .Therefore, both methods lead to the same mean results.

Fig. 24 .
Fig. 24.Retrieval ofMg (285.2 nm).The same methods as in Fig.21are applied.However, here the single measurements error is used.As for high relative errors the error propagation cannot be approximated to be linear and nonlinearities in the forward model lead to a shift of the mean to higher values.Therefore, a certain averaging of the data is needed before applying the retrieval algorithm.Note that the error bars estimate the error of single measurements while the error for the blue line is a factor

Fig. 25 .
Fig. 25.Single day (7 September 2009) vertical profile for Mg(285.2nm)  in the equatorial region from 15 • S to 15 • N.For too low constraint strength (λ h < 5 × 10 −9 ) the profile is strongly oscillating, while for too high constraint strength (λ h > 5 × 10 −7 ) the profile is too smooth.Acceptable solutions between these two extrema show peak densities from 500 to 2000 cm −3 .However, at least the peak altitude is nearly independent of λ h .The solution with λ h = 5 × 10 −8 is favored, considering, e.g., seasonal variations etc.

Fig. 26 .
Fig. 26.Single day (7 September 2009) result for Mg + (279.6 nm)in the equatorial region from 15 • S to 15 • N. The same approach to find the optimal constraint strength for Mg (see Fig.25) is applied for Mg + .However, the variation of the peak density with constraint strength, especially for low strength, is much weaker, so errors from a wrong choice of the constraint strength are smaller.For Mg + the solution with λ h = 1 × 10 −7 will be used for further investigations.

Fig. 30 .
Fig.30.Plot of Mg average profile of all single day results of the 2-D retrieval in the equatorial region (15 • S to 15 • N) for different percentages of Raman scattering in the Ring effect correction.The different corrections mainly affect the peak density and the lower peak edge.The peak altitude is nearly unaffected.Differences for variations of 1 % in the Raman scattering percentage vary the lower peak edge altitude by ≈ 1 km for altitudes above 80 km.For altitudes below 80 km, the differences drastically rise.

Fig. 31 .
Fig.31.Plot of Mg + (279.6 nm) average profile of all single day results of the 2-D retrieval in the equatorial region (15 • S to 15 • N) for different percentages of Raman scattering in the Ring effect correction.For altitudes above 80 km the Ring effect correction is negligible and it is very small above 70 km.In the peak region the emission dominates the background signal.Because of this the peak value rises with higher Raman scattering percentage (as discussed in Sect.3.2.2).
∂g i (x) ∂x k = s LOS gc i,k + s LFS gc i,k .
Vertical slant column emission (SCE) profile with different percentages of inelastic scattering in the Ring effect correction.The difference between the corrected and uncorrected profile increases with decreasing tangent altitude below 90 km.Note that the SCE of Mg at 53.5 km is still in the order of 50 % of the SCE at 90 km for typical recalculations of the forward model from retrieval results.
In the first part d dx k x i is δ ik and the first sum over i disappears:d dx k I = c 1 P k s LOS k f (g k (x))Now the first part is in its final state and the second part is evaluated.For the derivation of f the chain rule is applied to getd dx k I = c 1 P k s LOS k f (g k (x))