High-resolution temperature profiles (HRTP) retrieved from bi-chromatic stellar scintillation measurements by GOMOS/Envisat

Abstract. In this paper, we describe the inversion algorithm for retrievals of high vertical resolution temperature profiles using bi-chromatic stellar scintillation measurements in the occultation geometry. This retrieval algorithm has been improved with respect to nominal ESA processing and applied to the measurements by Global Ozone Monitoring by Occultation of Stars (GOMOS) operated on board Envisat in 2002–2012. The retrieval method exploits the chromatic refraction in the Earth's atmosphere. The bi-chromatic scintillations allow the determination of the refractive angle, which is proportional to the time delay between the photometer signals. The paper discusses the basic principle and detailed inversion algorithm for reconstruction of high resolution density, pressure and temperature profiles (HRTP) in the stratosphere from scintillation measurements. The HRTP profiles are retrieved with very good vertical resolution of ~200 m and high accuracy of ~1–3 K for altitudes of 15–32 km and with a global coverage. The best accuracy is achieved in in-orbital-plane occultations, and the accuracy weakly depends on star brightness. The whole GOMOS dataset has been processed with the improved HRTP inversion algorithm using the FMI's Scientific Processor; and the dataset (HRTP FSP v1) is in open access. The validation of small-scale fluctuations in the retrieved HRTP profiles is performed via comparison of vertical wavenumber spectra of temperature fluctuations in HRTP and in collocated radiosonde data. We found that the spectral features of temperature fluctuations are very similar in HRTP and collocated radiosonde temperature profiles. HRTP can be assimilated into atmospheric models, used in studies of stratospheric clouds and in analysis of internal gravity waves activity. As an example of geophysical applications, gravity wave potential energy has been estimated using the HRTP dataset. The obtained spatio-temporal distributions of gravity wave energy are in good agreement with the previous analyses using other measurements.



Introduction
This paper is dedicated to the description of a unique method for high-resolution temperature and density profiling using bichromatic satellite stellar scintillation measurements and to assessment of the retrieved temperature profiles. The bi-30 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas (http://envisat.esa.int/instruments/gomos; Bertaux et al., 2010). Before the description of the measurements and inversion algorithm, we would like to define precisely what atmospheric parameter is retrieved (or "measured").

1.1
What is a high-resolution temperature profile? 5 In the case of the HRTP (high-resolution temperature profile), the underlying atmospheric parameter that we aim to characterize is the temperature field. However, the temperature is a four-dimensional scalar field, which is defined for each time moment over a three-dimensional (3D) space. Due to very active dynamical processes in the Earth's atmosphere, this 3D field can contain significant fluctuations down to the viscous scale, which is typically smaller (and sometimes much smaller) than one meter within the stratosphere and the troposphere. 10 During occultations, the velocity of the sounding ray within the atmosphere is much larger than the velocity of any atmospheric motion (for GOMOS/Envisat, it is more than 3000 m/s), therefore the "frozen-field" approximation during measurement time can safely be considered (Tatarskii, 1971). Regarding the spatial variation, the trace of the line of sight within the atmosphere during an occultation defines a 2D surface, which differs only slightly from a plane because of refractive effects. The only temperature fluctuations able to affect the GOMOS measurements lie in this surface . 15 Furthermore, since the signal recorded by a detector is intrinsically one-dimensional, the retrieved parameter (temperature or density) is also one-dimensional.
Due to stable stratification of the stratosphere, most of the variations of meteorological parameters, such as the temperature, occur along the vertical direction. The field of temperature within the atmosphere is strongly anisotropic and the direction of its gradient is close to the vertical. Consequently, a measurement of temperature variations along the vertical direction, 20 known as a "temperature profile", describes most of the field variation within the considered region. However, the validity of the above statements is strongly dependent on the considered scale. Large-scale temperature fluctuations are strongly stratified (anisotropic), they contain the largest fraction of potential energy (or equivalently temperature variance). Most of the energetic dynamical processes, including meteorological flow and gravity waves correspond to this anisotropic part. The characteristic ratio between dominant horizontal and vertical scales is typically equal to the ratio of maximal and minimal 25 intrinsic frequencies of the gravity waves field. In the stratosphere, the ratio of buoyancy (Brunt-Väisälä) frequency N to the Coriolis parameter f, N/f, is typically larger than 100 (Fritts and Alexander, 2003).
Small-scale fluctuations, mostly turbulence and, more generally, instable and dissipative processes, are much more isotropic (and so are the temperature gradients associated with such small-scale processes). For this kind of fluctuations, the concept of a (vertical) profile is essentially meaningless. The transition between strongly anisotropic and roughly isotropic 30 fluctuations occurs within the scale range separating the domains of waves and turbulence. For the stratosphere, it covers roughly a decade between 10 and 100 meters (of vertical scale) (Gurvich and Kan, 2003;Nastrom et al., 1997).

3
These general considerations about the structure of the atmospheric temperature field indicate that a one-dimensional "vertical profile" is only meaningful (for remote sensing measurements) for vertical scales larger than 30-100 m. The detailed characteristics of the measurement process must also be considered in order to grasp the real meaning of the retrieved profile and to understand better its relationship with the atmospheric temperature field. In case of GOMOS, the concept of vertical profiles is adequate, as high-resolution temperature profiles, which will be discussed in our paper, have 5 the vertical resolution of ~200 m.

Bi-chromatic scintillation measurements by GOMOS and previous works on HRTP
For retrievals of high-resolution temperature profiles, we use bi-chromatic scintillation measurements by the GOMOS fast photometers, which record the stellar flux with the sampling frequency of 1 kHz at blue (475-525 nm) and red (650-700 nm) wavelengths synchronously, as a star sets behind the Earth limb. 10 Two fast GOMOS photometers on board Envisat recorded the intensity fluctuations induced on the star's light by the refractivity fluctuations encountered within the atmosphere, at two wavelengths. A short description of the inversion algorithm for retrievals of high-resolution temperature profiles from bi-chromatic scintillation can be found in (Dalaudier et al., 2006;Sofieva et al., 2009c); it is presented also below in our paper.
The main advantage of HRTP is its vertical resolution, which is ~200-250 m. Such resolution allows probing gravity 15 wave (GW) spectra. The validation of the small-scale structure of HRTP is therefore an important issue before using the data in GW research. The validation of the small-scale structure is a challenging task, because temperature fluctuations are rapidly varying due to gravity wave activity. Sofieva et al. (2008Sofieva et al. ( , 2009c proposed to use spectral analysis for validation of small-scale structure in temperature profiles, as this approach allows using measurements separated by several hundreds of kilometers and by several hours. The previous validation has been performed on HRTP processed by the FMI scientific 20 processor (analogous to the ESA IPF v6 algorithm) using collocated radiosonde data. It has shown that the small-scale fluctuations in HRTP have similar rms as in collocated radiosonde profiles, for vertical (in orbital plane) occultations of bright stars (Sofieva et al., 2009c). In case of oblique occultations or dim stars, the HRTP fluctuations are of larger amplitude than those of in-situ measurements. An analogous study with the method (Sofieva et al., 2009c) but applied to much larger datasets of GOMOS HRTP IPF v6 and collocated radiosonde profiles has shown that fluctuations in HRTP v6 are nearly 25 always larger than in the collocated radiosonde data (for both vertical and oblique occultations). Therefore, IPF v6 HRTP data were not recommended for gravity wave analyses. This will be illustrated in Section 4 of this paper.
In this paper, we introduce an improved version of the GOMOS HRTP algorithm and present the whole GOMOS HRTP dataset processed using the FMI's Scientific Processor, HRTP FSP v1.

1.3
The paper structure 30 The basic principle of HRTP retrievals is described in Sect. 2. Section 3 is dedicated to a detailed description of the retrieval algorithm. Examples of retrieved HRTP profiles, their characterization and validation of the small-scale fluctuations are Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 24 September 2018 c Author(s) 2018. CC BY 4.0 License. 4 shown in Sect.4. Illustrations of using HRTP for gravity wave analyses is presented in Section 5. The information about the GOMOS HRTP dataset and data access is presented in Section 6. Summary and discussion conclude the paper (Sect. 7) 2 Basic principle of HRTP retrieval As discussed in Sofieva et al., 2009b), the light intensity transmitted through the atmosphere is affected not only by absorption and scattering, but also by refraction and diffraction. The scintillations (or large intensity fluctuations) 5 observed at the satellite level are the result of the interaction of stellar light and atmospheric air density irregularities, which are generated mainly by internal gravity waves and turbulence.
The retrieval of HRTP is based on the chromatic refraction in the atmosphere (Dalaudier et al., 2006). The refraction angle α depends on wavelength, due to the optical dispersion of air. For two rays of different color having the same impact parameter p (Figure 1), the blue ray bends more than the red one ( Figure 1, R B    ) and will consequently be observed later by 10 GOMOS (only star settings are used in GOMOS observations). The scintillation spikes produced by atmospheric density fluctuations are observed by both photometers with a time delay BR tt  (Figure 1), which is proportional to the refraction The idea of such measurements of refractive angle has been first proposed by Sokolovskiy (1991, 1992).

15
Figure 1 A scheme of chromatic refraction and the principle of refraction angle measurement by GOMOS. Both the refraction angles and the effect of dispersion are strongly exaggerated. However, the width of each beam relative to the angle difference Δ is realistic. The impact parameter p is the geometric distance of the ray from the Earth's center. The vertical separation of blue and red rays at ray perigee 30 km is ~10 m ).

20
Using the accurate knowledge of the direction to the star and of the ENVISAT orbit, it is possible to convert the measured time delay into an angle difference Δ, and then into the refraction angle  at the reference wavelength (for GOMOS, the reference wavelength is 500 nm, i.e. =B). The conversion factor is equal to 94 ) ( n0(λ) -1 is the standard refractivity (n0 is the refractive index) given for dry air at standard pressure and temperature in Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 24 September 2018 c Author(s) 2018. CC BY 4.0 License. (Edlen, 1966), retractivities 0 B  and 0 R  correspond to central wavelengths of GOMOS photometers. After that, the method is similar to that used in radio occultation (e.g., Kursinski et al., 1997). Assuming local spherical symmetry of the atmosphere, the refractive index profile n(p) can be retrieved from the refraction angle profile (p) using the Abel transform (Kursinski et al., 1997;Tatarskiy, 1968):

5
The tangent (or minimal) radius r is related to the impact parameter p through the refractive index n : ( ( )) p r n r n p r  . The corresponding pressure profile is reconstructed by integrating the hydrostatic equation, as it is done in radio-occultation and lidar measurements (e.g., Hauchecorne and Chanin, 1980;10 Kursinski et al., 1997). Finally, the temperature profile () Tr is obtained from the state equation of a perfect gas.

From simplified theory to real experiment: HRTP processing algorithm
The main steps of processing the red and blue photometer signals to high-resolution temperature profiles can be outlined as follows:  Estimation of the chromatic time delay as the position of the maximum of the cross-correlation function of blue and red 15 photometer signals after smoothing the red one. Since refractive angle is proportional to chromatic time delay, the profile of the refractive angle is obtained  Determination of the refractivity profile from the refractive angle profile via the Abel integral inversion. For upper limit initialization, an atmospheric model is used.
 The density profile is obtained from the refractivity profile using the Edlen's formula. 20  From the density profile, the pressure profile is calculated using the hydrostatic equation.
 Finally, the temperature profile is determined from these data using the state equation of a perfect gas. This basic algorithm is given for a highly simplified situation. Below we present the detailed description of the most important inversion steps and discuss various effects, which have occurred in the GOMOS experiment.

From photometer signals to the profile of time delay
The new HRTP processing starts at the altitude 32 km, where time delay is larger than 1 ms and the scintillation are of large 5 amplitude. In the previous retrievals, the upper altitude, where the HRTP processing starts, depends on strength of scintillation and value of a priori time delay a  (estimated using ECMWF data). Since it is impossible to retrieve accurately time delay (and temperature profile) in the range where time delay is smaller than the sampling rate of the photometers, the upper limit is set to 32 km in the new retrievals.
The bandwidth of GOMOS photometers is nearly the same in wavelength, but the refractivity bandwidths of photometer 10 optical filters are significantly different ( Figure 3).

7
As a result, scintillation features are more smoothed in the blue signal than in the red one (Figs. 1, 2). In order to make the chromatic smoothing similar, the signal of red photometer is convolved with a Gaussian window G(t). The width of the smoothing window WG is defined as a differential width of the blue and red signals for a Dirac perturbation in the refractive angle, It increases proportionally to the refractive angle as the line of sight deepens into the atmosphere and can be approximated as: 5 where 0 B   is the variation of the standard refractivity corresponding to the spectral width B   of the blue photometer, R 0   is defined analogously for the red photometer, dhd /dt is the vertical velocity of the line of sight and L is the distance from tangent point to the satellite.
It is evident that the determination of time delay (as a function of time) by visual recognition of characteristic scintillation 10 structures is not feasible. Therefore, the time delay is estimated via calculation of the cross-correlation function (CCF), followed by a determination of the position of its maximum. For computation of the CCF, photometer signals are cut into ~50 % overlapping sections. The length of the sections should be chosen in order to contain a "sufficient" number of structures (scintillation spikes) while preserving the best available resolution for the angle profile. We found that the optimal length of segments window  corresponds to a vertical displacement of the line of sight within the atmosphere varying from 15 250 m at 32 km to 500 m at 5 km. The smoothed red signal is pre-shifted by the (smooth) a priori time delay a  rounded to the nearest millisecond (which will be hereafter referred to as a pre-shifted time delay) in order to best align with the blue signal. For computation of a  , we use ECMWF density data at the occultation location: The CCF is calculated with 1 ms resolution and the position of its maximum is searched around zero delay, in the range 20 8 where C   is the second derivative of the cross-correlation function of photometer signals at the point of its maximum, C is the error (standard deviation) of the CCF, and t is the discretization step (1 ms). The derivation of the formula (4) uses Taylor expansion of the condition for the CCF maximum is 0 ) (    C at the vicinity of the maximum 0 and the subsequent Gaussian error propagation. As follows from Eq. (4), uncertainty in determination of CCF maximum is proportional to uncertainty of CCF values and depends on the shape of CCF: the broader the CCF, the larger error. In 5 assumption of a Gaussian distribution of intensity fluctuations, the cross-correlation coefficient C has an asymptotically Gaussian distribution with the standard deviation where n is the size of samples participating in calculation of the cross-correlation coefficient. This approximation is valid for large samples. Finally, the error of time delay estimated can be written as 10  In (6), Cmax is the maximum of the cross-correlation function, or the cross-correlation coefficient; C   can be calculated from parameters of the parabolic fit. Due to assumptions made in deriving Eq. (6), it is clear that this formula gives only an approximate estimate for the error of time delay reconstruction.
The uncertainty of time delay increases at lower altitudes ( Figure 5, panel C, red line) due to broadening of cross-correlation function (Figure 4, right). At some levels, correlation between photometer signals can be low due to presence of turbulence, 5 which leads to very uncertain time delay.

3.2
Regularization of time delay estimation

Motivation: influence of isotropic turbulence
The atmosphere contains small-scale turbulence producing nearly isotropic fluctuations of density. These fluctuations also 10 produce scintillation during stellar occultation. The contribution of turbulence to the observed scintillations can be significant and sometimes even dominant (Sofieva et al., 2007a(Sofieva et al., , 2007b. The cross-correlation of bi-chromatic scintillations caused by isotropic turbulence is significant only when the chromatic separation distance of the ray trajectories does not exceed the Fresnel scale, which is ~1 m for GOMOS (for illustration and more details, see Fig. 4 and the corresponding text in (Sofieva et al., 2009b)). The chromatic de-correlation for nearly vertical (in orbital plane) occultations is always small, 15 while it can be significant in case of oblique (off orbital plane) occultations (Kan, 2004). On the other hand, the smoothing induced by the finite optical band of the photometers will selectively damp the fluctuations associated with turbulence in the vertical direction because of their smaller size. As a result, the cross-correlation between the two photometers has a minimum at some altitude depending on the obliquity of the occultation (hereafter we define the obliquity angle β as is the angle between the direction of the apparent motion of the observed star and the local vertical at the ray perigee point, 20 (Gurvich and Brekhovskikh, 2001;Sofieva et al., 2007b)). This is illustrated in Figure 4 (left) for the oblique occultation R07673/S001 with the obliquity angle =23. Strong turbulence is observed at upper altitudes, resulting in the drop of crosscorrelation at 30-45 km. A more quantitative consideration of this effect is given in (Kan, 2004) and (Gurvich et al., 2005).
In some situations, the correlation between photometers signals is too low for an accurate determination of the time delay (Figure 4 left). These situations are handled through regularization applied to the time delay profile, which is described 25 below.

Regularization algorithm
In the case of low correlation between the recordings of the fast photometers In the V6 algorithm, the data points corresponding to low correlation between photometer signals (with cross-correlation coefficient CCC<0.7) are replaced by the weighted mean of "measured" (obtained from cross-correlation) m eas  and a priori (computed external data source, e.g. ECMWF analysis data) The weights used in (7) are inversely proportional to the uncertainties of time delay 2   (defined by Eq. (6)) and the a priori 5 profile 2 a  . Hereafter, we will refer to this regularization as to the optimal filtration method. In the HRTP IPF v6 algorithm, the uncertainty of the a priori time delay is computed assuming that a priori air density has an uncertainty of 2.5% below 25 km, 5 % at 35-50 km with the linear transition between these two altitude regions. The effects of optimal filtration on time delay, its uncertainty, and the resulting temperature profile in illustrated in Figure 5 by blue lines. As observed in Figure 5, such filtration handles exceptional values, but the resulting temperature profile has enhanced amplitude of temperature 10 fluctuations at altitudes below 17 km compared to collocated sonde data. Validation of HRTP profiles, which have been processed with optimal filtration, against collocated radiosonde data (Sofieva et al., 2009c) has shown that the amplitude of temperature fluctuations in HRTP is realistic for vertical occultations of bright stars (not affected by turbulence), but often excessive in oblique occultations. In new algorithm, we use the statistical optimization (Bayesian approach). It acts as a linear operator applied to the differences between measured and a priori time delays: where Ca is the covariance matrix of a priori time delay uncertainties and Cmeas is the covariance matrix of measured time delay uncertainties. This formulation corresponds to the Bayesian estimator (maximum a posteriori method) provided 5 measurement errors and a priori uncertainties have Gaussian distribution. If both matrices Ca and Cmeas are chosen to be diagonal, the Bayesian approach coincides with the optimal filtration (7). The diagonal elements of Ca and Cmeas are "a priori" and "measured" uncertainties (variance of the corresponding errors), while off-diagonal elements characterize the correlation length scale of measurements and a priori profiles. In the new algorithm, the correlation length for measurements is equal to the window used for computation of cross-correlation function, while the correlation length of a priori profile is 10 set as twice larger. The diagonal elements of Ca and Cmeas are In the stratosphere, it is close to 1, while it decreases at lower altitudes. 20

From time delay to refractive angle
Refractive angle is proportional to time delay between photometer signals. The proportionality coefficient depends on difference in refractivity corresponding to the central wavelengths of photometers, distance from ray perigee point to the satellite, and the satellite velocity (Eq. (3)). However, this simple relation is complicated by the fact that the stellar spectrum is modified by absorption and scattering during an occultation. As a result, the effective wavelength of the photometer signal where min and max correspond to the wavelength range of the optical filter.
Then the refractive angle  at the reference wavelength, which will be used in the further processing, can be computed as: Since the refractive angle is proportional to time delay, its uncertainty can be easily obtained by multiplication of time delay 5 uncertainty by the corresponding factor.

From refractive angle to refractivity profile
By purely geometrical considerations and because the refractive angle is small, the impact parameter p for a given wavelength λ is given by: Applying the inverse Abel transform (Eq. 1), we can obtain the profile of refractive index n(p). The application of the Abel transform assumes local spherical symmetry of the atmosphere. This assumption is also used in retrievals from radiooccultation measurements. The error due to horizontal gradients of the refractive index at right angles to the direction of light propagation has been estimated in (Healy, 2001;; it is less than 1 % for altitudes. The integration of Eq. (1) can be carried out numerically using any standard quadrature method. The weak singularity of the integrand at the lower 15 limit does not cause problems for a numerical realization: the singularity can be estimated or the midpoint product integration method can be applied. The upper limit should be chosen high enough (~120 km), therefore the refractive angle profile calculated using ECMWF&MSIS data is used at altitudes above HRTP range in the processing. Application of Abel integration requires monotonous impact parameter. This requirement can be violated, because the impact parameter is computed using measured (noisy) refractive angle. In the current implementation, the impact parameter is computed using 20 the smoothed refractive angle and its monotonicity is checked.
Real geometric (tangent) altitudes can be determined as where R is the local radius of curvature of the Earth surface.

13
The error of refractivity reconstruction can be estimated using the matrix of the discretized Abel transform. Due to the fact that the Abel integral acts as a linear operator connecting refractive angle and refractivity, the covariance matrix of refractivity uncertainty  C can be estimated using the classical error propagation formula: T   C AC A (15) where A is the matrix of the discretized Abel transform (see e.g.  and C is the covariance matrix 5 of refractive angle uncertainties.

From refractivity to density, pressure and temperature
The density profile can be obtained from the refractivity profile using Edlen's formula. By using the hydrostatic equation we can calculate the pressure P at the altitude h as where R=8.3144 J/mol/K is the universal gas constant and M is the molar mass of dry air.
The covariance matrix of air density errors C is proportional to the covariance matrix of refractivity errors C (relative 15 errors are equal).
Two main terms contributing to the error in temperature are the error in local density (small scale structures in density and temperature are anti-correlated and of equal relative amplitudes) and the error in pressure at the top of the high resolution profile Ptop (error of upper limit initialization). They are added quadratically, thus giving the uncertainty of HRTP: (18) 20 In HRTP retrievals, the vertical resolution is defined mainly by the length of scintillation records that are used for calculation of cross-correlation function, which is ~ 250 m. Due to using overlapping samples, the actual vertical resolution is somewhat smaller. The regularization on time delay, as well as other inversion steps from time delay to temperature profile, can slightly degrade the vertical resolution in case of oblique occultations, thus overall vertical resolution of HRTP is expected to be close to 250 m. 25 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. uncertainties (shaded area). In Fig.6, HRTP profiles are for bright stars in vertical occultations (the best data quality), while in Fig. 7 other occultations are illustrated (oblique or of not bright stars). These temperature profiles are collocated with high-resolution radiosonde data from the SPARC data center (http://www.sparc.sunysb.edu/html/hres.html). The collocated 5 temperature profiles are shown by blue lines in the left panels of Figs. 6 and 7, and the information about the spatio-temporal difference is provided in the figure. We would like to note that the fine structure in the HRTP and radiosonde profiles are not expected to coincide, because of the evolution of the gravity wave field in the space-time window. For similarity of temperature profiles, including their small-scale fluctuation, the horizontal separation should be ideally less than 20 km and the time difference should not exceed 2-3 h, as discussed in Sofieva et al. (2008Sofieva et al. ( , 2009c. The time separation results in 10 additional spatial separation in the atmosphere caused by advection of air masses. The relatively long measurement time of temperature profiles by radiosondes during balloon flights (it takes ∼1 hour for balloon to raise to 10-30 km) has a similar effect (the quantitative estimates of these effects can be found in Sofieva et al. (2009c). In the left panels of Figures 6 and 7, the HRTP and collocated radiosonde profiles are similar, but not fully coinciding, as expected. The temperature profiles for vertical occultation of bright stars are of similar quality as in other occultations, as follows from comparison of Figures 6 and  15 7.
Despite differences in small-scale temperature fluctuations, we can expect similar spectral properties of the temperature field at locations not far from each other (e.g., less than 500 km) during some time period (a few hours), as shown in (Sofieva et al., 2008). The power spectra density of relative temperature fluctuations in HRTP and collocated radiosonde profiles are shown in the right panels of Figures 6 and 7. For the spectral analysis, the collocated profiles were 20 interpolated to an equidistant altitude grid with a 30 m resolution in the altitude range 18-30 km. Hanning filtering with a 3 km cut-off scale was used to obtain the smooth component. One can notice a good agreement of wavenumber spectra of HRTP and radiosonde temperature fluctuations, for all the considered GOMOS occultations. In contrast to the previous HRTP validation results, the spectra of temperature fluctuations are similar also in case of oblique occultation or not bright stars. 25 The HRTP wavenumber spectra in Figures 6 and 7 have visible cut-offs corresponding to scales ~150-250 m. This is the experimental conformation of the vertical resolution of HRTP ~200 m. This agrees with the theoretical estimates of the HRTP retrievals (it is defined mainly by the lengths of the scintillation records used for evaluation of time delay).    Typical estimated (in the retrieval algorithm) uncertainties of HRTP retrievals are shown in Figure 8, for occultations of different types. The HRTP temperature profiles are considered in the equatorial region (20S-20N, tropopause is ~18 km) in 2004. As seen in Figure 8, the estimated HRTP uncertainty is 1-3 K in the stratosphere, at altitudes from ~ 2 km above the tropopause to ~30 km. The best quality is achieved in vertical occultations. The uncertainty of the retrievals depends weakly 5 on star brightness, it is noticeably larger only for occultations of dim stars with visual magnitude mv>2.5.

10
In this paper, we focus on the validation of small-scale fluctuations in the HRTP, as the HRTP vertical resolution allows probing gravity waves. The validation of small-scale structure is a challenging task, because temperature fluctuations are rapidly varying due to gravity wave activity. Following Sofieva et al. (2008Sofieva et al. ( , 2009c, we use spectral analysis for validation of small-scale structure in temperature profiles, as this approach allows using measurements separated by several hundreds of kilometers and by several hours. This study applies the method by Sofieva et al. (2009c) to much larger datasets 15 of HRTP and collocated radiosonde profiles.
We have used the radiosonde data from the US high-resolution radiosonde archive available at http://www.sparc.sunysb.edu/html/hres.html (SPARC data center). The US high-resolution radiosonde archive contains data from 93 US operated stations from years 1998-2011. The stations are located across the mainland US, Alaska, Pacific islands and the Caribbean. Most of the data is in 6-second temporal resolution (the vertical resolution is ~30 m), but in recent years, 20 Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. The rms of temperature fluctuations in collocated HRTP and radiosonde data, hrtp  and sonde  , respectively, are presented as scatter plots (Figure 9). The colored markers in Figure 9 correspond to different versions of HRTP processing: IPF v6 (green), v6 algorithm implemented in FMI scientific processor (blue), and the new HRTP algorithm presented in this paper (red). We found that rms of HRTP v6 fluctuations is overall larger than that in collocated radiosonde profiles for all 15 occultation types, despite the vertical resolution is finer for radiosonde data (and thus the opposite behavior is expected).
However, the v6 algorithm implemented in FMI HRTP scientific processor shows the behavior consistent with the previous analysis of ( Sofieva et al., 2009c) 2x and y = (1/1.2)x, solid black lines: y=1.4 x and y=(1/1.4)

x.
Several examples of wavenumber spectra of relative temperature fluctuations are shown in Figures 6 and 7, which show 5 quite typical behavior: the spectra are similar for HRTP and collocated radiosonde profiles. In this section, we show illustrations of applications of HRTP for analyses of gravity waves. The spatio-temporal distributions are presented only for illustrations that new HRTP dataset provide valuable geophysical information, which is in agreement with analyses using other datasets.
The gravity wave potential energy per unit mass is defined as 5 where s T T  are relative temperature fluctuations with respect to the smooth (background) profile s T , and g is acceleration of gravity.
In our analysis, smooth profiles s T are obtained by smoothing HRTP down to 4 km resolution. The Brunt-Väisälä frequency N 2 is estimated using the smoothed HRTP profile. The GW potential energy has been evaluated for each temperature profile 10 in the altitude range 20-30 km.  Figure 10 are in very good agreement with previous estimates of this parameter from radiosonde, lidar and GPS radio-occultation measurements, both qualitatively and quantitatively (Allen and Vincent, 1995;15 de la Torre et al., 2006;Sofieva et al., 2009a;Tsuda et al., 1991Tsuda et al., , 2000. The overall distribution reproduces the known features: strong GW activity close to the edge of polar vortex, a peak near Antarctic Peninsula in local winter. Enhancements in equatorial regions are clearly observed, analogous to those found in global analyses of radio-occultation data (de la Torre et al., 2006;Tsuda et al., 2000).
Temporal evolution of GW potential energy, for different latitudes, is shown in Figure 11. This time series is very similar to 20 that shown in (de la Torre et al., 2006, Figure 1) obtained from radio-occultation data. Enhancements at polar and midlatitudes in winter are observed; they are larger in the Southern Hemisphere and follow the evolution of polar night jet (de la Torre et al., 2006;Sofieva et al., 2009a). The enhancements in the equatorial region is also observed, which seem to be annual but might be also modulated by quasi-biennial oscillations. Similar equatorial enhancements are observed by de la Torre et al. (2006). 25 We would like to note that, despite similarities of the GW energy morphology presented in our paper with the previous studies, there are expectedly some differences and peculiar features, because also contributions of small-scale gravity waves (down to 250 m vertical scales) are present in HRTP profiles. Detailed analyses of gravity wave distributions using HRTP might be the subject of future works and publications.    Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 24 September 2018 c Author(s) 2018. CC BY 4.0 License. Summary and discussion We described the improved algorithm for high-resolution temperature and air density profiling using the bi-chromatic scintillation measurements. This method allows temperature profiling in the altitude range 10-35 km with the vertical resolution of ~ 250 m and accuracy in the stratosphere of 1-2 K. The retrieval algorithm is applied to the whole 5 GOMOS/Envisat dataset, and the new GOMOS HRTP FSP v1 data set is now available. The best accuracy is achieved for vertical occultations bright stars. The uncertainty of HRTP retrievals depends weakly on stellar brightness.
The spectral analysis of HRTP and collocated radiosonde profiles has been applied for validation of small-scale fluctuations.
It has been shown that the HRTP fluctuations are realistic (in terms of their 1D vertical spectra).
The main factors limiting accuracy of HRTP retrievals are due to instrumental properties in combination with the 10 specifics of refraction in the Earth atmosphere. The upper limit of HRTP is defined mainly by the sampling frequency of the photometers: the detectable time delay should be larger than the photometer integration time, 1 ms; it is usually below 35-37 km. For faster photometers and larger wavelength separation, the upper limit can be higher. The presence of uncorrelated scintillations generated by locally isotropic turbulence reduces the useful information content in the photometer data. At lower altitudes, the influence of isotropic turbulence is low due to selective filtering by the photometer optical filters. 15 However, lower signal-to-noise ratio at lower altitudes due to influence of absorption and broadening of scintillation peaks due to chromatic smoothing degrade accuracy of HRTP retrievals at altitudes below 15-17 km. Narrower optical filters would allow slightly better retrievals an lower altitudes. The physical model for HRTP retrievals is adapted for vertical and moderately oblique occultations (for which tan(β) is smaller than anisotropy of air density irregularities). Such occultations constitute absolute majority of GOMOS measurements. The crossing of rays (strong scintillation) at low altitudes is not 20 taken into account by the model. However, its influence is reduced due to the spatio-temporal averaging by the GOMOS optical filters , thus allowing acceptable temperature retrievals from GOMOS scintillation measurements also at altitudes below 25 km.
HRTP can be assimilated into atmospheric models, used in studies of stratospheric clouds and in analysis of internal gravity waves activity. As an illustration of application of HRTP for gravity wave research, GW potential energy has been 25 evaluated using the GOMOS HRTP dataset. The obtained spatio-temporal distributions of GW potential energy are in good agreement with previous analyses using other datasets. This paper is dedicated mainly to the retrievals, validation and geophysical assessment of small-scale fluctuations in the retrieved GOMOS high-resolution temperature profiles. However, HRTP can be smoothed down to lower resolution, and used in other analyses, including analyses of temperature trends, in combination of 10-years GOMOS HRTP data with other 30 limb profile temperature measurements. This can be a subject of future research. Atmos. Meas. Tech. Discuss., https://doi.org/10.5194/amt-2018-270 Manuscript under review for journal Atmos. Meas. Tech. Discussion started: 24 September 2018 c Author(s) 2018. CC BY 4.0 License.

Data availability
The HRTP dataset is available from http://ikaweb.fmi.fi. The HRTP profiles presented in the dataset are interpolated to a common altitude grid from 10 to 32 km with 50 m spacing and stored in yearly netcdf-4 files. The README document provides the information about the parameters included in the data files. 5