Annales Geophysicae (2001) 19: 71–81 c ○ European Geophysical Society 2001

Abstract. An error analysis for mesospheric profiles retrieved from absorptive occultation data has been performed, starting with realistic error assumptions as would apply to intensity data collected by available high-precision UV photodiode sensors. Propagation of statistical errors was investigated through the complete retrieval chain from measured intensity profiles to atmospheric density, pressure, and temperature profiles. We assumed unbiased errors as the occultation method is essentially self-calibrating and straight-line propagation of occulted signals as we focus on heights of 50–100 km, where refractive bending of the sensed radiation is negligible. Throughout the analysis the errors were characterized at each retrieval step by their mean profile, their covariance matrix and their probability density function (pdf). This furnishes, compared to a variance-only estimation, a much improved insight into the error propagation mechanism. We applied the procedure to a baseline analysis of the performance of a recently proposed solar UV occultation sensor (SMAS – Sun Monitor and Atmospheric Sounder) and provide, using a reasonable exponential atmospheric model as background, results on error standard deviations and error correlation functions of density, pressure, and temperature profiles. Two different sensor photodiode assumptions are discussed, respectively, diamond diodes (DD) with 0.03% and silicon diodes (SD) with 0.1% (unattenuated intensity) measurement noise at 10 Hz sampling rate. A factor-of-2 margin was applied to these noise values in order to roughly account for unmodeled cross section uncertainties. Within the entire height domain (50–100 km) we find temperature to be retrieved to better than 0.3 K (DD) / 1 K (SD) accuracy, respectively, at 2 km height resolution. The results indicate that absorptive occultations acquired by a SMAS-type sensor could provide mesospheric profiles of fundamental variables such as temperature with unprecedented accuracy and vertical resolution. A major part of the error analysis also applies to refractive (e.g., Global Navigation Satellite System based) occultations as well as to any temperature profile retrieval based on air density or major species density measurements (e.g., from Rayleigh lidar or falling sphere techniques). Key words. Atmospheric composition and structure (pressure, density, and temperature; instruments and techniques) – Radio science (remote sensing)


Introduction
Absorptive occultation data have a great capacity to provide high-accuracy profiles of atmospheric key variables like temperature and trace gases such as ozone (see, e.g., the review of Smith and Hunten, 1990, and for a detailed modern occultation experiment overview Russell et al., 1993). The limb geometry pertaining to occultation sounding enables high vertical resolution, the self-calibrating nature of the method, which rests on normalized rather than absolute intensity data, allows high accuracy, and the frequency-selective absorptive properties of the different species of gases in the atmospheric medium render it possible to simultaneously measure different species densities with multiple-channel sensors. From the profile of a major species like molecular oxygen, pressure and temperature in turn can be deduced by exploiting the law of hydrostatic equilibrium and the equation of state. Depending on ray tangent height, which is the height of closest approach to Earth's surface of a limb sounding ray, refractive bending of occultation rays needs to be accounted for or discounted; for heights above the stratopause (∼50 km) relevant in our mesospheric profiling context refraction effects can be neglected (if desired small residual bending up to ∼70 km is readily accounted for) and the concept of straight-linepropagated radiation applies (see Smith and Hunten, 1990). Furthermore, we assume all sensor channels at wavelengths with strong absorption by some species so that also the contribution of scattering to the intensity data is negligible (or readily rendered negligible by climatological correction).
The problem of retrieving species density profiles from the measured intensity data can be divided into two distinct steps (e.g., Kyrölä et al., 1993;Fussen et al., 1997;Fussen, 1998): First, transmission data from a number of wavelength channels are used in a "spectral inversion" in order to estimate the composition of the atmospheric gas for each individual sounding ray, i.e., for each tangent height. Having subsequently performed this for all rays (data points) of an occultation event, this yields along-ray columnar content profiles for all species included. Second, a "spatial inversion" step is performed for retrieving the number density profiles of the species from the columnar-content profiles, the relation being essentially described by an Abelian integral equation.
Since these two steps are not independent of each other in reality (e.g., the "spectral inversion" depends on the temperature profile which can be actually known after the "spatial inversion" only), a few iterations of the procedure may be required in practice. Alternatively, it may be feasible to unify these two steps in a single consistent formulation of the retrieval process. Subsequent work is planned to explore this latter option.
In this study we restrict ourselves to the inversion problem of one species, where the spectral inversion step is degenerated into rearrangement of Beer-Lambert's transmission law in scalar form. This restriction is sensible for the present focus on a baseline analysis of the temperature profiling performance of the "Sun Monitor and Atmospheric Sounder" (SMAS) absorptive occultation sensor proposed by Kirchengast et al. (1998), which uses molecular oxygen (O 2 ) absorption in the mid-UV range between 180 and 210 nm for providing temperature information. Generally, the assumption of one absorbing species is reasonable in wavelength ranges, where the main absorbing species significantly dominates the total absorption compared to the contribution of all other species. In the case of the SMAS sensor, O 2 is the main absorber in most of the selected wavelength range (Schumann-Runge bands), at >200 nm ozone (O 3 ) absorption (Hartley band) increasingly dominates (see, e.g., Brasseur and Solomon, 1986). We can reasonably ignore the O 3 influence for the present error analysis, however, since the SMAS sensor with its eight mid-UV channels in total (including also channels >210 nm at 224 nm and 246 nm, respectively) allows very effective decomposition into O 2 and O 3 absorption contributions (Rehrl, 2000).
Generally either stars (stellar occultation) or the Sun (solar occultation) are used as radiation source for absorptive occultations (Smith and Hunten, 1990). Utilization of stars can in principle furnish a large number of occultation events due to the multitude of stars available, but is limited by the problem of rather weak source intensities. This excludes them as option for wavelengths <250 nm as used by the SMAS sensor, while the Sun provides sufficient intensity for the targeted high-accuracy mesospheric temperature and ozone profiling. On the other hand, the Sun as a single source leads to a rather small number of occultation events given a single sensor (e.g., Russell et al., 1993). Using the mission proposed by Kirchengast et al. (1998), which is based on a small constellation of six micro-satellites each equipped with a SMAS sensor (besides a refractive occultation sensor using Global Navigation Satellite System signals), would nevertheless allow us to acquire about 150 occultation events per day and thus lead to a database of globally distributed mesospheric temperature and ozone measurements of unprecedented quality.
The most consistent approach for dealing with statistical errors, which we adopt for this study, is to describe all random variables based on their probability density function (pdf). Gaussian errors are fully described by a pdf involving a mean and a covariance matrix, respectively, corresponding to the first and second moment of a Gaussian pdf, for which all higher moments vanish (e.g., Anderson, 1984). Errors of non-Gaussian shape are fully described by their pdf as well, but non-zero higher moments are involved. Furthermore, we perform a purely algebraic propagation of errors, i.e., no statistical estimation methods and no supporting a priori information are involved. This ensures that all pdfs obtained are rigorous and essentially unique solutions for the retrieval errors of interest ("essentially" since some approximations such as Taylor expansions are involved, all of which are robust and fully reasonable in the given context, however).
The work describes the error analysis by proceeding along the following four steps: from intensity measurements to columnar content profiles (Sect. 2), from columnar content to density profiles (Sect. 3), from density to pressure profiles (Sect. 4), and from density and pressure to temperature profiles (Sect. 5), respectively. The main conclusions of the study are detailed in Sect. 6.

From intensity measurements to columnar content profiles
An occultation event is acquired by a satellite-borne solar occultation sensor each time the Sun sets or rises from the perspective of the sensor behind the atmosphere at Earth's limb. At each sampling time during such an event, typically 10 times per second, the sensor's optical detectors, in case of the SMAS sensor either UV diamond (Di) diodes or UV silicon (Si) diodes, simultaneously observe intensities in all wavelength channels λ k . By geometry, each sampling time corresponds to a specific tangent height z i (in short "height" hereafter). For the mid-UV SMAS channels of interest here, the sensor's accuracy is limited by thermal noise (intrinsic detector noise) so that the measured intensity values, denoted I ki (channel λ k and height z i ), are expected to be gathered around the correct intensity I c,ki = I ki following a Gaussian distribution of errors, where · means averaging over a large ensemble. Each intensity I ki can thus be modeled as a Gaussian random variable, which is (as every random variable) fully described by the relevant probability density function (pdf). Suppressing for convenience the channel and height indices and the normalization factor, the pdf reads From the performance characteristics of the optical detectors the sensor's noise level is usually well known. We express it in form of a reference precision, termed γ hereafter. In case of the SMAS sensor, the expected precision γ of unattenuated solar intensity measurements I 0 at 180 nm, at a sampling rate of 10 Hz, is 0.03% (γ = 3 · 10 −4 ) for Di diodes (N. Fox, National Physics Lab, Teddington/UK, private communication, 1999) and 0.1% (γ = 1 · 10 −3 ) for Si diodes (G. Schmidtke, Fraunhofer-Institute for Physical Measurement Techniques, Freiburg/Germany, private communication, 1998, 2000), respectively. Since the accuracy is thermal noise limited, the precision of an arbitrary intensity measurement I can be expressed as i.e., as the ratio of a constant root-mean-square noise value γ I 0 to the correct value I c . Inserting this into Eq. (1) yields The absorptive occultation method is said to be self-calibrating, since the absolute intensity values do not matter (γ needs to be sufficiently small, though) but rather the intensities I ki normalized by the intensity I k0 estimated "above the atmosphere", i.e., at heights where virtually no attenuation occurs. Thus we are only interested in the transmission Tr ki , the ratio between the attenuated and the "vacuum" signal, which reads for each channel λ k and height z i Tr ki = I ki I k0 . ( We assume that I k0 can be measured very precisely and thus set I k0 = I c0 . This is justified in this baseline error analysis for two reasons: Firstly, I k0 can be repeatedly measured many times "above the atmosphere" and averaged thereafter so that its statistical error is suppressed according to the inverse square root of the number of averaging ensemble members. Secondly, an eventual instrumental bias in I k0 is likely to be accompanied by virtually the same bias in the attenuated intensities I ki given that a typical mesospheric occultation event lasts less than one minute; such bias is thus largely cancelled in Tr ki . Given I k0 non-random and known, Eq. (3) provides a relation for linearly transforming Eq. (2). Since a Gaussian pdf does not change (except for normalization) by linear transformations (e.g., Anderson, 1984), the pdf for Tr ki can be expressed, by inserting Eq. (3) into Eq. (2), similar to P (I ) as (4) For the present baseline error analysis for the SMAS sensor we selected five channels located in the O 2 Schumann-Runge bands at wavelengths from 185 nm to 205 nm (185, 191, 195, 198, and 205 nm), representative of the five SMAS "O 2 /temperature" channels within 180 and 210 nm (the remaining three "O 3 " channels are centered at 210, 224, and 246 nm; see Rehrl, 2000). To be conservative, the precision γ at 180 nm introduced above was baselined for all five channels, though γ would somewhat improve in practice with increasing wavelength due to increasing solar flux. Figure 1 illustrates the transmission profiles for the five channels, which were computed based on an exponential atmosphere with a constant scale height of 7 km and assuming at the surface a total pressure of 1013.25 mbar, a temperature of 288 K, and an O 2 volume mixing ratio of 20.948%. The O 2 absorption cross section values were taken from Nicolet et al. (1989).
The profiles represent error-free "correct" profiles (Tr = Tr c ), used as starting point for the error analysis, for which we assume two different scenarios with Gaussian sensor errors according to Eq. (4); one with γ = 6 · 10 −4 (Di diodes) and one with γ = 2 · 10 −3 (Si diodes), respectively. These values involve an increase of the detector noise introduced already by a factor of 2, which was applied in order to add a reasonable error margin in compensation for a simplified absorption cross section treatment (discussed later). The profiles have been computed on a 200 m height grid, approximately mimicking the 10 Hz sampling rate of the SMAS sensor, which will experience a typical vertical scan rate of ∼2 km/sec. As Fig. 1 shows, a sensible choice of five channels is enough to probe the full temperature profile from 50 km (stratopause region) up to about 100 km (mesopause region). (The actual SMAS design includes eight channels for accurate decomposition of O 2 and O 3 contributions and for providing, in addition, a mesospheric ozone profile.) Based on Eq. (4), the next step is to relate the transmission to the optical thickness, τ ki , utilizing Beer-Lambert's transmission law, Since the transmission Tr ki is valued between 0 and 1 by definition (Eq. 3), the optical thickness τ ki , a logarithmic projection, is valued between ∞ and 0. The transformation from Tr ki to τ ki changes the pdf substantially. Combining Eq. (4) (multiplied in the exponent by (Tr c /Tr c ) −2 ) and Eq.
(5) yields the pdf for τ ki as P (τ ) in this form is clearly not of Gaussian shape (c.f. Kyrölä et al., 1993), a property, which turns out to be of importance to keep the subsequent analysis quasi-analytically tractable.
Favorably, however, since τ/τ c ≈ 1 due to the precision of the optical detectors at ∼1% or better, an expansion to first order of the inner exponential function is valid allowing us to linearize the problem. The expansion is very accurate when τ ki is not too high, i.e., for transmissions Tr ki > 0.1. This is neatly compliant with the fact that the transmission range actually foreseen to be exploited for each SMAS sensor channel is bounded to 0.1 < Tr ki < 0.9 anyway, since beyond this range adequate precision of the transmission data will be difficult to maintain in practice (e.g. Smith and Hunten, 1990). The linearization simplifies Eq. (6) to which shows that the optical thickness τ can be reasonably described by a Gaussian pdf with the correct optical thickness τ c as mean and a variance (γ exp[τ c ]) 2 . The last step to a columnar content profile is to compute it based on the available optical thicknesses. Letting σ k denote the absorption cross section (extinction cross section, in general) in channel λ k , the optical thickness and the columnar content, respectively, read where n(z) is the number density profile, ds i is the raypath for the sounding ray with tangent height z i , and d ki is the columnar content value at height z i based on the optical thickness in channel λ k , respectively. According to Eq. (8) a linear relation between d ki and τ ki is assumed. Furthermore, we assume σ k to be non-random and known for each λ k in this baseline analysis. These assumptions do not strictly apply to the 5 nm wide SMAS channels situated in the Schumann-Runge bands, which exhibit highly structured and temperature-sensitive O 2 cross sections (e.g., Nicolet et al., 1989). They reflect, however, that the non-linearities ignored by Eq. (8) are presumably of minor importance (and can be largely modeled), that in practice the effective cross sections for the SMAS channels can be determined fairly accurately (e.g., with the aid of pre-flight laboratory measurements), and that the temperature dependence can be coped with by using a reasonable guessed temperature profile and an iterative approach (as addressed in Sect. 1). In order to roughly account for potential errors due to residual cross section uncertainties given the present simplified formulation, we increased the detector noise by a factor of 2. Future work will investigate this error component with quantitative rigor.
Given that Eq. (8) with σ k non-random is a linear transformation, the pdf for d ki , P (d), is of the same Gaussian form as The dependence of the columnar content variance on σ is important, because the cross sections as a function of wavelength span a range of about 3 orders of magnitude as indicated in Fig. 1. Due to the channels with smallest σ probing the lowest heights, this leads to a dramatic increase of the variance of d with decreasing height. At the same time, however, also the absolute d values increase over about 3 orders of magnitude so that the precision of d is roughly conserved. For an individual channel, the variance increases with decreasing height following the transmission decrease (note that γ exp(σ d c )/σ = γ /(σ Tr c )). Given overlap in height of the exploitable transmission ranges of channels, those furnishing smaller variance at a given height should receive higher weight. This is best implemented by calculating the total variance of d(z i ) at any height z i as the joint optimal estimate (inverse-variance sum) of the variances of all channels fulfilling 0.1 < Tr ki < 0.9 at z i . Overall, the variance properties discussed lead to the relative error in d(z i ) being roughly constant with height.
In practice, at this point the "spectral inversion" part is finished and a columnar content vector d = d(z i ) associated approximately with a variance as derived is available.

From columnar content to density profiles
We start with the columnar content profile d from the "spectral inversion" (O 2 columnar content in the present context) and its covariance matrix S d . Profile d is computed from the transmission profiles illustrated in and described along with Fig. 1. To this end the d k profiles for the five channels at 200 m sampling are ingested into an optimal (inversevariance weighted) combination in order to obtain a single best-estimate d profile at all heights. In this baseline analysis, where we start with the correct transmission Tr c,ki from Fig. 1, the d obtained is the correct mean profile d c . The corresponding total covariance matrix, S d , assembles the total variance derived in Sect. 2 in its diagonal and is otherwise zero-valued, since each d(z i ) value is computed independently during "spectral inversion" and thus uncorrelated with all others. Formally it reads where S d,k are the diagonal covariance matrices of the five channels (with variance values for heights outside the bounds 0.1 < Tr ki < 0.9 set to a very high number mimicking infinity). S d is produced for both the γ = 6 · 10 −4 /Di diodes ("DD") scenario and the γ = 2 · 10 −3 /Si diodes ("SD") scenario.
Next, d is downsampled from 200 m to 2 km (mimicking from 10 Hz to 1 Hz), by averaging over independent groups of 10 values, in order to smooth the data to the natural resolution foreseen for the SMAS sensor, which is ∼2 km (Suncentered field-of-view vertical width ∼1/30 deg). The variance values in S d , for both DD and SD scenarios, are similarly averaged and divided by 10 in order to account for the suppression of statistical noise by the averaging. All subsequent scenario processing is done at the 2 km sampling grid and starts with (d = d c ; S d ) prepared as described.
According to Eq. (8), the columnar content, d(z i ), is an integral over the number density profile, n(z), along the path of the sounding ray (O 2 number density in our context). Adopting the reasonable assumption of spherical symmetry of the atmosphere over a few hundred kilometers about the mean tangent point location of the occultation event, the integral can be analytically cast into a pure radial form which is found to be of Abelian type (e.g., Smith and Hunten, 1990), where r i is the local radius of curvature corresponding to the tangent height z i (e.g. Steiner et al., 1999; well approximated by r i = z i + 6371 km in our context), and r is the radial coordinate in general corresponding to the height coordinate z (r = z + 6371 km in our context). Evaluation of Eq. (10) for all heights z i of interest yields the full d(z i ) profile.
To deduce the number density profile n = n(z j ) from d(z i ) we need the inverse of this integral, which analytically exists as inverse Abel transform (e.g., Smith and Hunten, 1990). For this study, we require the integral in matrix form, however, thus we perform a discretization of Eq. (10) into matrix notation, which we then can use both for inverting d to n and for analyzing the error propagation. Following a method widely called "onion peeling", Eq. (10) is discretized from top (z = 120 km) downwards into height shells (r → r j ), with shell boundaries r j and a shell depth r j −1 −r j of 2 km in order to match the sampling of (d = d c ; S d ).
Densities are for convenience assumed constant in each shell, but in order to sufficiently suppress discretization errors, the scheme is enhanced (following Smith and Hunten, 1990) by accounting for the quasi-exponential density variation in each shell by attributing the obtained density to a height level z j = r j + 1/3(r j −1 − r j ) − 6371 km (rather than to the average height in the shell).
In matrix notation, the Abel integral, Eq. (10), reads where A is the geometry matrix representing the forward Abel transform (a lower triangular matrix), d = d(r i ), n = n(r j ), and i, j = 0, . . . , N . For the special upper boundary domain beyond the top shell (z > 120 km; r > r 0 ), A ij is adopted in an analytical form (Swider, 1964) involving the complementary error function and the scale height H (set to 7 km), Equation (12) corresponds to the assumption of a spherically symmetric, isothermal (constant scale height) atmosphere beyond 120 km, which is readily replaced for real data by a more elaborated formulation. While Eq. (11) maps a density profile n to a columnar content profile d, its inversion in order to retrieve n from d as needed in this study is straightforward, where A −1 denotes the inverse of matrix A. Since the sensor channels sample well with their 0.1 < Tr ki < 0.9 measurement ranges the entire height range of interest, the mapping expressed in Eq. (13) is well conditioned, A is robustly of full rank and the inversion works stable. Thus no inclusion of indomain a priori information is required in this case such as for example in related refractive occultation techniques (see, e.g., Steiner et al., 1999), where a Bayesian optimal estimation approach (e.g., Rodgers, 2000) is very helpful. In our scenario processing, Eq. (13) yields n ∼ = n c , as we insert d = d c and since discretization and boundary errors in A are negligible in the context of this baseline analysis. (The linear operator A, or A −1 , is readily implemented in a more accurate approximation if desired; see, e.g., Syndergaard (1999), where more levels and a piecewise linear discretization instead of constant value per shell were adopted.) The covariance matrix of the retrieved density profile, S n , is retrieved from the covariance matrix of the columnar content profile, S d , by where A T is the transpose of matrix A, and the whole quadra- (see, e.g., Anderson, 1984). Formally, the validity of Eq. (14) is tied to the assumption of a Gaussian pdf around a mean state d c , which we found valid as expressed by P (d) in Eq. (9). Given that Eq. (13) is a linear transform, also the pdf for n, P (n), is Gaussian, with n c being the correct mean profile estimated via Eq. (13) and S n the covariance estimated via Eq. (14), respectively. In general, covariance matrices are clear indications of a Gaussian pdf assumption. While S −1 d , following Eq. (9), is simply a diagonal matrix with values ∝ σ 2 /γ 2 , this does not hold for S −1 n and S n , respectively. Nevertheless, the magnitude of the values of S n is scaled by S d so that the magnitude of the retrieval error in the density profile n is dominated by the errors specified in S d . Hence, a proper quantification of the covariance matrix to start with in an error propagation analysis is a step of key importance for the credibility of results further down in the retrieval chain.
Since, in addition to S d being diagonal, the elements of A are relatively simple geometric expressions (see Eq. 11) and the matrix itself is of special form (A lower triangular, and hence A T upper triangular), we can gain deeper insight into the generation of S n by closer inspection of Eq. (14). Although the quadratic form introduces covariances between all pairwise elements of n, its inversion ensures the columnar content errors to be transported from the top shell, and higher shells in general, down to lower shells, while keeping the errors of lower shells irrelevant to the density profile above. This is in line with the Abel transform pair fundamentally being limited to the half-space above any height z i of interest. Utilizing this property of S −1 n and that S d is diagonal together with the need for symmetry between rows and columns, Eq. (14) can be reformulated to which instructively highlights the structure of S −1 n and thus also provides a more explicit insight into the shape of the density pdf, P (n) (Eq. 15). Unfortunately no similarly simple expression exists for S n directly, but the determinant of S −1 n can be written as which gives a hint of the magnitude of density retrieval errors as the elements of S n are inversely proportional to det(S −1 n ). Fig. (2) illustrates, for both the DD and SD scenario, the estimates for the absolute and relative standard deviations ( S n,ii and 100 × S n,ii /n c,i ), respectively, and for the correlation functions (rows S n,ij / S n,ii S n,jj ) for the density profile retrieval from columnar content. Also the mean profile n ∼ = n c itself is shown for reference. We note that while we use O 2 number density in this section, this can be readily converted to the total number density of air by computing n air = n/V O 2 , where V O 2 = 0.20948 is the O 2 volume mixing ratio.
The standard deviations (Fig. 2, left panel) show the magnitude of statistical retrieval errors. The relative errors are found to be roughly constant with height and smaller than 0.15% / 0.5% for the DD / SD scenario. The "wavelike" shape of the standard deviation profiles (best visible in the relative error) reflects the fact that different errors ∝ γ /(σ k Tr c,ki ) are associated with the different sensor channels, each contributing over a limited height range. Comparing Figs. 1 and 2 reveals that "valleys" (minimum error) occur where transmission ranges of adjacent channels well overlap and thus the total variance in d(z i ) is small, while "spikes" (maximum error) occur at heights with least overlap, i.e., within the small height intervals where only a single channel fulfils 0.1 < Tr ki < 0.9.
In addition to the magnitude of errors, the covariances amongst the density profile elements (Fig. 2, right panel) are relevant, not least for subsequent retrieval of pressure and temperature profiles. The covariances of density values of adjacent shells amount to about 30% of the variances, so they have to be accounted for when interpreting the results. The covariances decrease with increasing distance from the reference shells, the (anti-)correlation being below 10% from the second closest shell outwards.
The covariances, in terms of their order of magnitude and their structure, reveal several properties of the retrieval of density profiles from column densities. First, the entropy of the retrieval error is appreciably lower, and the information content of the measurement is thus appreciably higher, than it might be believed from considering the variances only (e.g., Rodgers, 2000). In other words, the pair-wise correlation between shells contains additional information not included in the variance (standard deviation) and should thus be accounted for in the exploitation of the data, for example, when fused with other data by some "optimal combination". Next, all covariances are negative (see Fig. 2, right panel), which is a direct implication of the onion-peeling scheme: The densities are calculated, starting from the top shell, through determining the density of the actual shell by subtracting the known densities of the higher shells; thus a positive error in a higher shell leads to a negative error in a lower shell or vice versa. Another related feature of interest is that a columnar content error in a particular shell predominantly affects the number density error in this shell, but the downward propagation of the error is damped rather quickly as can be verified with numerical examples using synthetic "impulse response" cases. This damping mechanism is an implication of the described onion-peeling scheme, which gives rise to anti-correlation between adjacent shells.
A particularly important property of the covariances is their beneficial significance with respect to the pressure profile retrieval to follow. As this step requires summation of erroneous quantities when evaluating the hydrostatic integral, a negative correlation means a damping and a partial compensation of the errors within this sum. This property is of high relevance for the quality of retrieved pressure profiles, ignoring the importance of the covariances would lead to significant underestimation of their accuracy and, in turn, of the temperature accuracy.

From density to pressure profiles
The retrieval of a total air pressure profile p = p(z i ) from a number density profile n = n(z i ) of a major species (O 2 in our context) requires height integration over the total air mass density profile ρ(z i ), which is related to the species num-ber density profile n(z i ) from Sect. 3 via ρ(z i ) = m eff n(z i ), where m eff is the effective molecular mass for converting species number density to total mass density (m eff =M/(N A V O 2 ) = 2.2960 · 10 −25 kg in our context, whereM is the mean molar mass of air, N A is the Avogadro constant, and V O 2 is the volume mixing ratio of O 2 ). p(z i ) is then obtained by integrating the hydrostatic law over all heights above the pressure level of interest, where g(z) denotes the local acceleration of gravity (g(z) = g mesosph = 9.6 ms −2 in our processing). The densities ρ(z i ) given for discrete shells only, we approximate Eq. (18) by the sum where z j is the depth of the shell at level z j and p 120km is the pressure at the boundary of the topmost shell (in our context z j = 2 km throughout and p 120km is the pressure at 120 km of the exponential atmospheric model used). Using j = 1 in the lower bound of the sum implies that the density (covariance) values related to z > 120 km (i, j = 0) are dropped from this point on. The "−1/3" in the upper bound of the sum expresses that the lowest shell is included only down to height z i (i.e., as g(z j )ρ(z j )(2 z j /3)), a notation also used in Eq. (20). It is reasonable to assume that p 120km is non-random and the errors due to its systematic uncertainty are negligible more than ∼ 3 scale heights lower in the height range of interest <100 km. In the scenario processing, Eq. (19) yields p ∼ = p c , as we insert n ∼ = n c into ρ(z i ) and since discretization and boundary errors are negligible in this baseline analysis. The summation of density values with their errors represented by S n leads via Eq. (19) to top-downward accumulation of errors in the pressure profile p. Invoking the "Gaussian pdf linear transformation" theorem (Anderson, 1984), the latter errors can be represented by a pressure error covariance matrix S p = B · S n · B T , where B is a lower triangular matrix of full rank such as A defined via p = B · n + p 120km , i.e., it is the linear operator in Eq. (19), 1 st r.h.s. term, which maps n to p. S p can be explicitly expressed in the form where B ik = (I i ) · (m eff g(z k ) z k ) for k < i (I i is a vector with all values unity), B ik = m eff g(z k )(2 z k /3) for k = i, B ik = 0 for k > i, and B j l = B ik since both i/k and j/ l are indices of rows / columns of B. Equation (20) indicates that in S p all elements of S n down to the height of interest are summed up. The inclusion of all co-variance values besides the variance values in the sum is essential, since the anticorrelation of the density error at a given level with the other levels (discussed in Sect. 3) suppresses the standard deviation in p compared to a variance-only sum by a factor of ∼4.
Given that Eq. (19) is a strictly linear transform, the shape of P (n) (Eq. 15) remains unchanged and the pdf obtained for p, P (p), is again Gaussian, with the mean pressure p c estimated via Eq. (19), and the covariance S p estimated via Eq. (20), respectively. Fig. 3 illustrates, for both the DD and SD scenario, the error analysis results (standard deviations, mean, and correlations) for pressure retrievals in exactly the same format as Fig. 2 did for density retrievals.
The standard deviations (Fig. 3, left panel) for the pressure profiles are found to be smaller than 0.04% / 0.12% for the DD / SD scenario. Comparing Figs. 3 and 2 shows that height dependence and shape of the standard deviations in pressure are quite similar to those in density, due to the linearity of the retrieval step from density to pressure. The relative error is about a factor of 4 smaller in pressure, however, owing to the effective suppression of density errors via Eq. (20) as discussed already. The correlation functions (Fig. 3, right panel) have significantly broadened in S p compared to S n because of the hydrostatic integration involved in Eq. (20). This indicates that significant information is contained in the co-variance values, i.e., that the information content of the pressure data is substantially higher than suggested by considering only the variances (left panel).
Accounting for S p in its complete form is, similar as for S n , of high relevance for the proper subsequent estimation of the temperature error statistics, which involves both pressure and density via the equation of state.

From density and pressure to temperature profiles
The retrieval of a temperature profile T = T (z i ) from a number density profile n = n(z i ) obtained in Sect. 3 and a pressure profile p = p(z i ) obtained in Sect. 4 is straightforward as we can locally apply, in each shell, the equation of state (ideal gas law), where K = k B /V species , with k B the Boltzmann constant and V species the volume mixing ratio of the major species utilized (V species = V O 2 = 0.20948 in our context; V species would be unity if total air density n air were used). Equation (22) yields T ∼ = T c if we insert p ∼ = p c estimated via Eq. (19) and n ∼ = n c estimated via Eq. (13). With the exponential model used (see Sect. 2), T c is retrieved via Eq. (22) constant with height at a value consistent with the assumed scale height of 7 km, i.e., close to 235 K, which is a typical mesospheric temperature. For illustrating the relative temperature errors in Fig. 4 we will thus use T c =T = 235 K. Regarding the error statistics, Eq. (22) indicates that we need to consider the ratio of two multivariate Gaussian random variables with pdfs P (p) (Eq. 21) and P (n) (Eq. 15), respectively, in order to derive the pdf for T , P (T ). Clearly P (T ) will be non-Gaussian, as the ratio represents a grossly non-linear transform. According to the "Gaussian pdf invariance" theorem (Anderson, 1984), if the multivariate pdfs P (p) and P (n) are Gaussian, also the local univariate pdfs P (p(z i )) and P (n(z i )) are Gaussian of the form P (x(z i )) ∼ exp[−1/2(x(z i ) − x c (z i )) 2 S −1 x,ii ], with x = p or n. Thus, in principle, at least the variances of P (T ) could be evaluated locally based on univariate Gaussian pdfs, which may deem more simple. Unfavorably, it can be proven, however, that no (finite) moments exist, both in the multivari- ate and the univariate problem, for the distribution of a ratio X/Y of two real-valued Gaussian random variables X and Y , except in the trivial case of perfect correlation between X and Y (W. Schachermayer et al., Inst. for Statistics, Probability Theory, and Actuarial Mathematics, Techn.Univ. of Vienna, private communication, 1999; see also Kendall and Stuart, 1969, who show that for independent univariate X and Y a Cauchy pdf is obtained for which no moments exist). Rigorous analytical computation of temperature (co)-variances is thus basically impossible. A study on solving the problem with different other methods will be presented in a forthcoming paper.
In this study we step back to a linearization of the problem in order to obtain an approximate temperature error covariance matrix S T . This will still be a very accurate approximate solution given that the relative errors of profiles in our context do not exceed the 1% level. We follow an approach used by Syndergaard (1999) in an error analysis for refractive occultations. A Taylor expansion of T (z i ) about T c (z i ) up to first order is performed which yields for δT = T − T c , where δn(z) = n(z) − n c (z) and B ij corresponds to B ik defined for Eq. (20). Equation (23) is a linear transform for retrieving δT , and hence T = T c + δT , from δn = δn(z i ) known from Sect. 3, where δ ij is the Kronecker Delta. Given that n, and hence δn, follows a Gaussian pdf (Eq. 15) and given the linearity and full rank of matrix operator C (a lower triangular matrix such as A and B) we can obtain S T analogously to S p (Eq. 20) as The pdf for T , P (T ), is thus well approximated by a Gaussian pdf, which can be written as with the mean temperature T c estimated via Eq. (22), and the covariance S T estimated via Eq. (25), respectively. Figure 4 illustrates, for both the DD and SD scenario, the error analysis results (standard deviations and correlations) for temperature retrievals in the same format as Fig. 2 did for density retrievals, except that the mean profile is not thought worth being shown in this case (given that T c =T = 235 K).
The standard deviations (Fig. 4, left panel) for the temperature profiles are found to be smaller than 0.3 K / 1 K for the DD / SD scenario. In shape the relative errors are quite similar to those in pressure and density, reflecting the linearity of the retrieval chain, which is besides the self-calibrating nature of the occultations a very favorable property for use of profiles in climate applications (Kirchengast et al., 1998). Given the isothermal mean temperature profile, also the absolute errors are roughly constant with height (rather than decreasing as density and pressure errors do, following their exponential mean profile). This weak height dependence will essentially hold for realistic temperature profiles as well, since the actual temperature at any mesospheric height will not deviate more than up to ∼20% from the mean mesospheric value adopted here.
The correlations of temperature errors (Fig. 4, right panel) are very close to the ones for density (see Fig. 2), due to the fact that the operator C in Eq. (25) propagates S n to S T , i.e., there is no direct involvement of S p . Inspection of Eq. (24) reveals that in C the diagonal values dominate and that the shape is smooth, governed by C ij ∝ 1/n c (z i ) (if j ≤ i, else C ij = 0). This has the favorable implication on the temperature error structure that C does not appreciably change the correlation structure in S n , so that the temperature errors are, in essence, no more correlated than the density errors, although the intermediate pressure errors exhibit much more correlation. Weak correlation, ideally no correlation, is useful in optimal data fusion (e.g., in data assimilation into numerical weather prediction models), since it simplifies the specification of reasonable covariance matrices. The correlation obtained for T errors is, though weak, not negligible so that the information in the co-variances should be accounted for when utilizing the data (see the discussion on S n in Sect. 3).
Finally, it is worth noting in this context that the often used simplified way to estimate S T ,ii by linearizing the equation of state and exploiting only the standard deviations of p and n, i.e., using at each height δT = (∂T /∂p)δp+(∂T /∂n)δn = (nδp − pδn)/(Kn 2 ), is overly conservative: It overestimates S T ,ii by a factor of ∼1.5.

Summary and conclusions
An analysis of statistical errors involving the full probability density function (pdf) of the relevant random variables was performed for mesospheric temperature profiling by absorptive occultation sensors. The analysis is capable of explaining the origin and propagation of errors within the retrieval chain. Starting with Gaussian-distributed intensity measurements by the sensor, each step of calculation down to temperature profiles has been assigned a pdf with its describing parameters (mean profile and covariance matrix). The analysis chain was applied to a baseline quantification of the temperature profiling performance of a solar UV sensor of the type proposed by Kirchengast et al. (1998) ("Sun Monitor and Atmospheric Sounder" -SMAS). Two different photodiode detector scenarios were considered for the SMAS-type sensor, respectively, a diamond diode (DD) with 0.03% and a silicon diode (SD) with 0.1% unattenuated intensity measurement noise at 10 Hz sampling rate. These noise values were increased by a factor of 2 in order to roughly account for unmodeled potential errors due to cross section uncertainties.
The following main conclusions can be drawn.
1. The absolute error of the retrieval products (density, pressure, temperature) at a given height is predominantly determined by the ratio of the precision of the detectors to the product of absorption cross section and transmission of the channel furnishing the highest value for this product at that height. For SMAS, the absolute temperature error as well as all relative errors were found to be roughly constant with height, indicating a well balanced channel selection.
2. The retrieval error covariance matrices of density, pres-sure, and temperature are dominated by the variances but, nevertheless, contain co-variances which carry significant information. These co-variances thus need to be accounted for when propagating errors further through the chain (particularly important for density covariances), when interpreting the errors in the result profiles (e.g., pressure error correlations extend over more than a scale height before decreasing to <10%), and when exploiting them in data fusion setups like assimilation into atmospheric analysis and prediction systems.
3. The correlation lengths of temperature errors are equally small as those of density errors and thus not appreciably affected by the broad pressure correlation functions. This is somewhat counter-intuitive, but from the point of view of data exploitation favorable, finding is rooted in the fact that temperature variances are directly dependent on density rather than pressure variances as the error analysis instructively shows.

4.
A major part of the analysis applies also to refractive occultations such as those using Global Navigation Satellite System (e.g., Global Positioning System) signals. In that context a bending angle profile is to be processed instead of a columnar content profile, otherwise the basic chain is the same. Specifics like the need for a priori information due to large high-altitude errors need special consideration, however, a matter which will be studied in a forthcoming paper focusing on a theoretical error analysis of refractive occultations.
5. The analysis from density onwards applies to any temperature profile retrieval based on density profiles, i.e., not only to absorptive or refractive occultations but also e.g., to Rayleigh lidar and falling sphere techniques. While some error analyses have been performed for the latter methods in the past (e.g., Hauchecorne and Chanin, 1980;Schmidlin et al., 1991), the treatment in this work is significantly more rigorous than previous ones.
6. Concerning the temperature profiling performance of the two SMAS-type sensor scenarios investigated, the temperature is found to be retrieved throughout the mesosphere to better than 0.3 K (DD) and 1 K (SD) accuracy, respectively, at 2 km height resolution. The other retrieval products (from earlier steps in the processing chain) are of corresponding quality.
In summary, though several aspects of interest have been loosely addressed only in this baseline analysis (uncertain absorption cross sections, finite channel widths, top-of-atmosphere normalization uncertainty, etc.), the results indicate that absorptive occultations acquired by a SMAS sensor could provide mesospheric profiles of fundamental variables such as temperature with unprecedented accuracy and vertical resolution. A next step currently being worked out is a study of the performance of the full SMAS sensor for simultaneous sensing of mesospheric temperature and ozone profiles. It is hoped that SMAS sensors are built and placed in orbit in the near future.