GNSS radio occultation profiles in the neutral atmosphere from inversion of excess phase data

Long-term stability, global coverage and high resolution are characteristics that make the Global Navigation Satellite System (GNSS) radio occultation (RO) technique well-suitable to serve as anchor measurements for observing the Earth’s atmosphere. The concept of occultation soundings utilizes a receiver placed on a low Earth orbit to measure the accumulated atmospheric contribution along the limb in terms of a phase delay. A high sampling rate allows to reconstruct profiles of geophysical parameters through an inversion process of occultation signals. However, such measurements require a careful processing in order to provide accurate retrievals in the neutral atmosphere. The following development describes specific aspects in radio occultation methodology implemented in the retrieval chain from phase data to profiles of dry pressure and dry temperature. Independent retrievals from near-real time measurements are compared with occultation products provided by official processing centers to demonstrate reliability of the solution. The region within the upper troposphere and lower stratosphere (UTLS) is particularly represented by a low uncertainty being within 0.5% (K). A comparison with radiosondes shows a significant contribution of a water vapor term in the lower troposphere that comes from the dry air assumption in occultation profiles of pressure and temperature. Radiosonde measurements reproduced to refractivity profiles show very high agreement with occultation soundings, which is generally below 5%. A superior accuracy of RO refractivity is observed in the upper troposphere, where retrievals are consistent with radiosondes


INtroductIoN
The concept of radio occultation (RO) was first applied in extraterrestrial satellite missions to study the composition of planetary atmospheres from measurements collected on-board Mariner, Pioneer and Voyager spacecraft (Kliore et al. 1965;Fjeldbo et al. 1971;Fjeldbo 1973;Lindal et al. 1983). The potential of RO technique in the remote sensing of the Earth's atmosphere has been investigated in Global Positioning System/Meteorology (GPS/MET) experiment (Rocken et al. 1997). Following its success, subsequent occultation missions have been launched, which includes but is not limited to, Challenging Mini-satellite Payload (CHAMP) (Wickert et al. 2001), Gravity Recovery and Climate Experiment (GRACE) (Beyerle et al. 2005), and Meteorological Operation (MetOp) satellites (Von Engeln et al. 2009). A significant contribution to occultation soundings has been introduced with the launch of six low-Earth orbiting (LEO) satellites in U.S.-Taiwanese project: Formosa Satellite Mission 3/Constellation Observing System for Meteorology, Ionosphere and Climate (FORMOSAT-3/ COSMIC; F-3/C for brevity) (Anthes et al. 2008).
GNSS limb soundings rely on GNSS signals propagating between a pair of occulting satellites once they set or rise behind the Earth's horizon (Kursinski et al. 1997). A receiver placed on a low Earth orbiter (LEO) records signals transmitted from a satellite of Global Navigation Satellite No reference! E a r l y R e l e a s e Terr. Atmos. Ocean. Sci., Vol. 30, No. 2, 1-19, April 2019 System (GNSS) to accumulate atmosphere-induced effects. The angle at which radio waves bend is derived from excess Doppler shifts computed relative to straight-line propagation paths, given stable atomic clocks and LEO positions estimated in precise orbit determination (POD) (Hwang et al. 2009). Inversion of bending angles allows to recover refractivity profiles of the neutral atmosphere in the assumption of local spherical symmetry (Melbourne et al. 1994). However, propagating rays encounter significant disturbances in the troposphere due to increasing contribution of water vapor (Sokolovskiy 2001(Sokolovskiy , 2003. As such, the reconstruction to profiles of pressure and temperature cannot be achieved without external estimate of humidity (Poli et al. 2002) and thus, dry-air assumption applies to direct retrievals (Syndergaard 1999).
The well-established radiosonde observations, routinely launched from fixed locations for several decades, are considered a true reference for atmosphere profiling (Rao et al. 2009;Seidel et al. 2010;Liu et al. 2014;Wilgan et al. 2017). Although in-situ soundings may lack in certain aspects of a spatio-temporal resolution, due to limited coverage mostly to lands and typically two soundings daily, they serve as key observables for operational weather analyses. The radio occultation technique can be viewed as a relatively new method, but its potential Xie et al. 2012;Ho et al. 2015) and reliability (Kuo et al. 2005;Wang et al. 2013) have been already widely emphasized. Highinclination orbits allow for global coverage provided from uniformly distributed occultation events. Due to the absence of instrument-dependent systematic effects, RO soundings are capable to capture radiosonde biases (Sun et al. 2010). The high accuracy of RO profiles, especially within upper troposphere and lower stratosphere (UTLS) (Foelsche et al. 2009;He et al. 2009;Steiner et al. 2011), can be applied to correct satellite radiance measurements from nadir-viewing space-based platforms (Liu et al. 2016a) to reduce the error in the initial conditions of numerical weather prediction (NWP) models (Liu et al. 2008(Liu et al. , 2016b. The COSMIC Program Office, housed at the University Corporation for Atmospheric Research (UCAR), Boulder, USA with the COSMIC Data Analysis and Archive Center (CDAAC), serves as a leading provider of RO products at multiple data levels for available occultation missions (Kuo et al. 2004;Schreiner et al. 2011). Over time, a number of processing centers was established to focus on particular applications or mission-specific retrievals in order to ensure long-terms stability of RO measurements based on multicenter results (Steiner et al. 2013). Currently, National Space Organization (NSPO) in Taiwan establishes an independent center for processing FORMOSAT-7/COSMIC-2 data, Taiwan/TriG Radio Occultation Process System (TROPS), as a further development of the F-3/C retrieval algorithm (Chiu et al. 2008), with contributions from Global Positioning System Science and Application Research Center (GPSARC) at National Central University (NCU) and the Taiwan Analysis Center for COSMIC (TACC) at Taiwan's Central Weather Bureau (CWB).
In the following development, an implementation of the retrieval process for the inversion of radio occultation phase data is demonstrated for real-time profiling the neutral atmosphere. A joint effort of Wroclaw University of Environmental and Life Sciences (WUELS) and NCU allowed to establish a platform (ROWNP for brevity) for assessment, processing and comparison of radio occultation data. The pilot system has been set to collect regional occultation soundings using real-time F-3/C measurements and perform the inversion process to geophysical parameters. The procedure consists of pre-processor for handling the excess phase data, inversion to bending angle profiles with geometrical optics (GO) and wave optics (WO), integration to the refractive index and retrieval of atmospheric profiles of dry pressure and dry temperature. The retrieval methodology has been verified by comparing output profiles to official RO products and in-situ measurements acquired from operational radiosonde network. In the following sections, we show an architecture of the processing system, discuss a numerical implementation in the retrieval chain and present by-products at consecutive steps. Inter-comparison studies for derived atmospheric profiles provide uncertainty estimates of RO parameters and concluding remarks follow to summarize the established development.

SoftwAre ArchItecture
The RO processing software was developed to routinely process real-time F-3/C data based on initial products provided by TACC/CDAAC. The application of real-time products, which are typically distributed few hours after measurement time, makes the retrieval to be more likely affected by observation errors contrary to postprocessed products, at the expense of early data delivery. The retrieval system runs on a UNIX-type operating machine connected to Network-attached storage (NAS). The core program is written in Perl language, which is initialized by crontab and launched by shell script (sh). MATLAB serves as a computing environment for data processing and analysis.
The retrieval chain, presented in Fig. 1 and discussed in the following sections, is feed by input data from level 1b products that proceed through several processing stages. Coordinates transforms, given satellite positions and respective occultation geometry, precede the inversion process in order to define a reference frame with respect to Earth figure expressed by WGS-84 ellipsoid. Geoid undulation for the estimated tangent point (TP) is computed from EGM-96 model and applied for the altitude conversion. The excess phase data is pre-processed for the quality control of observations and demodulation of navigation message. The signal is terminated based on signal-to-noise ratios (SNR) on E a r l y R e l e a s e L1 and L2 channels. The bending angles are then computed from filtered phases and the optimized profile is integrated to the refractive index, which is further converted to geophysical parameters of dry atmosphere.

obServAtIoNAl domAIN for dAtA retrIevAl
The occultation tables available in occTab real-time product of level 1b allow to define spatio-temporal frame for data processing. Predictions of events in LEO-specific dumps of data provide locations of perigee points together with occultation time of atmospheric soundings. Routinely downloaded daily occultation files are scheduled to extract regional RO soundings limited to the area of Central Europe, with a domain centered at Poland and spanning neighboring countries. The spatially filtered data summarized in Table 1 are acquired and processed at consecutive retrieval stages.
Since radiosondes are considered key-observations for atmospheric profiling and essential data source for validations, balloon soundings from stations being a part of Global Observing System are acquired once the routine processing is initialized. Locations of three stations in the vicinity of Poland with WMO identifiers 12120 (54.75°N, 17.53°E), 12374 (52.40°N, 20.97°E), 12425 (51.13°N, 16.98°E) are presented in Fig. 2 (from north to south) together with randomly distributed occultation events observed between January and May 2015. A shell script executes Uniform Resource Locator (URL) syntax using the cURL software library to download radiosonde observations to local storage from the National Oceanic and Atmospheric Administration (NOAA) and Earth System Research Laboratory (ESRL) radiosonde database. The developed processing scenario for the atmospheric profiling with GPS-RO technique aims to complement in-situ soundings and provide estimates of retrieval uncertainty for occultation profiles. In order to satisfy this aim, selected meteorological variables measured with the ascend of radiosondes are utilized for the reconstruction to complementary profiles from the RO technique to perform a comparison once collocated soundings are identified.

INverSIoN of exceSS phASeS to GeophySIcAl pArAmeterS
The excess path lengths for a series of measurements within a single occultation event are assessed from level 1b operational product atmPhs provided by the TACC/ CDAAC. Utilized carrier phases at L1 and L2 frequencies are computed in single-difference strategy to remove satellite clock errors (Schreiner et al. 2010). The L1 frequency is modulated by coarse-acquisition (C/A) code (L1C for brevity) and the second signal is available for both C/A and precise (P) codes (L2C and L2P respectively). Improvements in the receiver software introduced with open-loop (OL) tracking (Ao et al. 2009) prevent from an early termination E a r l y R e l e a s e of the signal and allow for deep penetration of the lower troposphere. The recursive feedback in traditional closed-loop (CL) tracking applies for L1C and the receiver uses phaselocked loop (PLL) to adjust the frequency of the incoming signal based on previous measurements. The un-encrypted L2C signal, similarly to L1C, can use both PLL and OL modes, whereas L2P is tracked in semi-codeless PLL only (Sokolovskiy et al. 2014). Once certain conditions are satisfied, the receiver tracking the signal in a setting occultation switches from CL to OL and a model-aided operating mode is activated. The inversion errors in the lower troposphere related to the tracking technique (Sokolovskiy 2001;Beyerle et al. 2003Beyerle et al. , 2006 could be therefore eliminated. The particular challenge in tracking of rising occultation stems from low SNR values since the very beginning of the occultation event, which are driven by highly variable troposphere and, potentially, associated multipath conditions. Since the loss of lock issue arising from PLL no longer applies in OL mode, the number of observations substantially increased due to the ability of acquiring rising occultations. The retrieval software processes both, setting and rising occultations. The observation data are reversed for rising occultations to generate all output profiles in a setting order. Doppler shifts yield the excess path in terms of the phase data for which the amplitude data of single-to-ratio indicates the quality of the observations.

occultation Geometry
The idealized GPS-to-LEO geometry illustrated in Fig. 3 needs to be solved for a kinematic system of moving satellites at non-coplanar orbital planes and non-circular orbits. The geometrical distance S between a pair of occulting satellites is derived from LEO and GPS positions given at measurement epochs in seconds since the beginning of the occultation where the coordinates are given in geocentric inertial frame. The occultation (central) angle θ = θ L + θ G is computed from the dot product of two satellite vectors and respective radii R G and R L : product data level description gpsBit level 1b GPS navigation data modulation bits for processing open loop data in the lower troposphere atmPhs level 1b atmospheric excess phases and auxiliary data after single difference clock correction atmPrf level 2 reference observations of the bending angle and refractivity as processed by the TACC/CDAAC Table 1. Data products used in the processing and inter-comparisons. No reference! E a r l y R e l e a s e The impact parameter for a straight line distance a 0 , commonly referred also as a height of straight line (HSL), is calculated as follows: Since the Earth's oblateness introduces deviations from the perfectly spherical shape utilized in the assumption of symmetrical atmosphere, the ray tangent to the local sphere needs to be found from a set of perigee points prior to phase data inversion to atmospheric profiles. The tangent point is determined by minimizing Euclidean distances for perigee points and the radius of WGS84 ellipsoid. The tangent ray is illustrated in Fig. 4. The elliptical shape of the Earth is considered by finding a new curvature center xyz c according to Syndergaard (1998): where xyz t represents the tangent point coordinates and n e is Earth's normal unit vector as presented in Fig 4. The local radius of curvature R c is calculated from azimuth-dependent Euler formula where A is the azimuth of the occultation plane at the tangent point defined as an angle to GPS satellite with respect to North, M and N are radii of curvature along the meridian and prime vertices, respectively. Geodetic latitude and longitude of the tangent point follows from coordinates in Earth fixed reference frame (ECEF). Transformation from Earth-centered inertial (ECI) at J2000 epoch, which applies for satellite positions, is fulfilled by the rotation with respect to z-axis according to sidereal hour angle (for details see Benzon et al. 2003).

Navigation data modulation
In PLL mode, the projection of the phase ahead is reliable for tracking single-tone signals under sufficient signal-to-ratio (SNR). The phase correction based on the assumption of 50 Hz sampling allows for the internal removal of navigation data modulation (NDM) at the receiver level. When the phase rate between data samples exceeds the bandwidth (i.e., larger than π/2), the robust tracking becomes unstable and leads to loss of lock. Once the OL mode becomes activated, the tracking can be maintained due to possible prediction of atmospheric conditions, hence associated Doppler frequency shifts, based on a phase model. However, the process requires demodulation of 50 Hz navigation data through the external processing of radio occultation signals, together with resolving cycle ambiguities. Removal of NDM based on the navigation message is applied for L1C signal.
A binary sequence is transmitted in the navigation message and provided in real-time gpsBit product of level 1b. For the purpose of decoding, the samples of navigation bits have to be synchronized with L1C phase data. This is due to the fact that the phase is assigned to a receiving time, whereas the sequence of bits corresponds to a transmitting time. The half-cycle jumps observed in the excess phase rate in Fig. 5a can be related to changes in a binary stream to remove the effect of NDM. The respective phase replica is obtained by controlling the rate of change between data samples. A binary series is retrieved from the samples by the assignment of zeros (no jump) and ones (jump), respectively to the excess phase rate. Corresponding data sequences are Fig. 3. Standard occultation geometry with locations of GPS satellite (GPS), low-Earth orbit satellite (LEO), and Tangent Point (TP). For each occurrence of subscript G and L, we refer to the GPS and LEO satellites, respectively. R G and R L are satellite radii, a G and a L are impact parameters with a 0 serving as impact parameter for a straight line distance, α is the bending angle, θ G and θ L are the occultation angles, and e G and e L are the excess angles, which together with φ G and φ L yield the incidence angles I G and I L .
E a r l y R e l e a s e presented in Fig. 6 as a function of time since the beginning of the occultation. A common starting bit is obtained through cross-correlation of the signals and half-wavelength corrections apply for the excess phase data as indicated by navigation message bits. If the condition of half-cycle slip is satisfied, all following data samples are inverted.
In case the navigation message is missing or incomplete, the removal of NDM can be performed without external information about 50 Hz bits. This approach applies for demodulation of L2C signal and, if the navigation message is missing or incomplete, also for L1C. The phase corrections follow from phase replica itself. As such, the excess phase rate corresponding to half-cycle jumps needs to ex-ceed 0.03 m s -1 . The full cycle phase ambiguities of L1C and L2C, introduced by either incorrect modulation of the navigation bits or unresolved after NDM removal, are detected and resolved by exceeding the filter to the full signal wavelength, independently on the NDM removal method.
Demodulated L1C and L2C excess phases are adjusted with sliding polynomial regression based on Savitzky-Golay filter of third degree and first order in 0.5 s windows. The removal of NDM together with the application of the filter allows for a smooth derivative of the L1C phase with respect to time, as presented in Fig. 5. The NDM cannot be applied for L2P signal, which is dropped once PLL-OL is detected (approx. 54 occultation second), since a further smoothing NDM is removed from L1C raw phase data after PLL-OL transition. Strong fluctuations arise in L2P phase as no NDM is available and the signal is truncated once OL is activated for L1C.
Check! E a r l y R e l e a s e of the phase cannot provide reliable retrieval. Depending on the occultation, the tracking mode switches from PLL to OL between -10 to -20 km below the height of straight line connecting a transmitter and a receiver. Figure 7 shows an example of occultation, where both L1 and L2 signals are modulated by C/A code. In such a case, tracking in OL mode extends the applicability of the L2C signal to the troposphere. Its phase can be demodulated from binary series derived the phase replica and application of the smoothing yields consistent retrievals in terms of the cut-off altitudes.

Geometrical optics retrieval
The bending angle as a function of the impact parameter is derived from the occultation geometry given the NDMcorrected excess Doppler and orbital positions with velocities for setting/rising satellites. The applied geometrical optics (GO) approximation imposes single-ray propagation between LEO and emitting GPS satellite. In the presence of circular LEO and GPS orbits, the departure of the impact parameter (a) from the ray tangent to the straight line propagation path (a 0 ) could be derived from a simple relation of the optical path length computed from geometrical distance and the excess path (d}) and the angular rate of the occultation angle (di), so that a = a 0 + d d } i. However, realistic non-circular orbits and the effect of atmospheric refraction require an iterative approach to find the impact parameter based on the first guess value. The equation follows from a derivative of the optical path ): where i o is a derivative of the occultation angle. The angles at GPS and LEO can be therefore derived from the Snell law, since R L n L sin(I L ) = R G n G sin(I G ) = a, to yield the bending angle from geometrical relations: where I corresponds to the incidence angles at GPS and LEO satellites. Figure 8a shows the geometrical optics bending angle profile retrieved after NDM is removed and a smoothing filter is applied. The unresolved cycle slips induce significant residuals once the receiver transitions to OL mode at around 10 km above the Earth in the impact height space (i.e., impact parameter reduced by local curvature radius R C = 6384.4 km). Nevertheless, the zoomed section of the profile within the troposphere in Fig. 8b results in ambiguous determination of the bending angle, which is no longer a monotone function of the impact parameter. The geometrical optics retrieval thus enables the excess phase to be related to more than one pair of bending angles and impact parameters (Sokolovskiy 2001;Gorbunov 2002a). The effect is further resolved by applying wave optics methods and derivation of the bending angle from a raw complex signal.

wave optics retrieval
The single-valued function of the bending angle with respect to the impact parameter is determined by means of wave optics (WO), which allows for disentangling of incoming rays in multipath conditions. A number of available methods based on Fourier integral operators (Gorbunov 2002a;Gorbunov and Lauritsen 2004;Sokolovskiy et al. 2007) adopts radio-holographic solutions to enable the multipath signal propagation to be corrected, causing the impact parameter to monotonically decrease or increase for setting Check! E a r l y R e l e a s e and rising occultation, respectively. In the following software retrieval, we apply two radio-holographic methods: phase matching (PM) ) and full spectrum inversion (FSI) (Jensen et al. 2003), both of which are based on the stationary phase method (Born and Wolf 2000). The unambiguous relation of the frequency is straightforward in single ray regions. The single-valued dependence in multipath conditions is found from Doppler angular rate from a raw complex signal and the bending angle profile is derived with sub-Fresnel vertical resolution.
Within the FSI method, a transformation of the wave field with respect to a geometrical space u(t) to the impact parameter space U(a) is realized by means of the Fourier integral of the entire signal ( ) From the transformed phase ( ) z~ it can be found that its derivative leads to ( ) t-, which is a singlevalued function. In the following implementation, the time is replaced by a space parameter of the occultation angle i as the independent variable. The derivative of the Doppler angular frequency ( ) is proportional to the impact parameter k ã i = o , given a wavenumber k, assuming that any instantaneous frequency can occur only once. Hence, ( ) a i follows from the known function (Jensen et al. 2003): Check!

E a r l y R e l e a s e
The phase matching can be considered as an alternative method of the stationary phase ). The method uses the assumption that the main contribution to the phase integral occurs near the stationary point and there is only one stationary point. In order to reduce the computational load, a short width of time window is selected within the oscillating integrand and the higher order components are neglected by applying a Gaussian function. A relation between the space parameter and the partial derivative of the matching phase is expressed by The phase matching function is directly applicable to noncircular occultation geometry and the respective occultation angle for given stationary point and corresponding impact parameter ( ) a s i allows to find the refraction angle : ' 1 1 (10) Figure 9 shows resolved multipath effect in the geometrical optics profile of the bending angle with the use of two radioholographic methods. Because the phase matching can be applied for realistic orbits without approximate corrections to the complex electromagnetic field (i.e., propagation to closes circles in FSI due to application of FFT), it serves as a standard retrieval method in the troposphere. The FSI bending angle is used as a control parameter for the validation of PM. The wave optics retrievals are applied only for L1C signal. Combination of geometrical optics and wave optics methods for the inversion process to the bending angle requires both retrievals to be connected to yield a single occultation profile. The upper-air geometrical optics profile is replaced and sewed with radio-holographic retrieval at the impact height of 20 km.

calibration of bending Angles in the upper Atmosphere
In the RO retrievals for the troposphere studies, the ionosphere is considered a source of noise. The first order term of the ionospheric refraction is excluded by the linear combination of bending angles α on the L1 and L2 frequencies f according to Vorob'ev and Krasil'nikova (1994): where f L1 and f L2 signals have frequencies of 1575.42 and 1227.60 MHz, respectively. Prior to the correction, the bending angles from L1 and L2 channels, interpolated to common set of impact parameters, are smoothed in 2 second time windows to reduce the effect of residual ionospheric noise. The limited application of L2P carrier to the range within the closed-loop tracking (Fig. 5) enables the linear combination of the frequencies down to the height of approximately 10 -20 km. This makes L2P carrier applicable as an ionospheric correction provider for higher altitudes. Although L2C allows to track the signal significantly deeper, which is comparable to L1C, the retrieval software is currently set to provide the ionospheric correction down to a consistent joint altitude. The L1-L2 bending angle profile is merged to L1 profile at 20 km. Occultations with larger L2 cut-off height are discarded. Below the joint altitude, the L1-L2 profile is extrapolated down to the troposphere to provide continuous correction. The retrieval noise at higher altitude levels becomes comparable to ionospheric residuals effects. In order to avoid propagation of these errors to geophysical parameters through further processing of bending angle profiles, a statistically optimized estimate is determined by weighting observations above 40 km in terms of the impact height. The optimized bending angle at the ith impact parameter is obtained according to relation (Hocke 1997;Healy 2001): where the error of the background bending angle profile is estimated as ( ) .
and the fraction is a scaling coefficient (C) for observed o a minus background b a bending angles. The background bending angle profile is calculated by inverting a refractivity profile derived from the 1976 Standard Atmosphere model. The error of the observed bending angle, approximated as , is smoothed in 3 km windows and dynamically weighted based on SNR distribution at L1 carrier amplified by a factor of 3 to account for random errors due to ionospheric correction. L1 SNR estimated at impact parameters is normalized and smoothed in a scale length corresponding to diameter of the first Fresnel zone. The approach allows to weight measurements towards the background proportionally to ionospheric activity and accounts for small-scale irregularities. Figure 10 illustrates how scaling coefficient and observation error change with SNR slope and periodical fluctuations with respect to unweighted measurements. Since large variations in unsmoothed SNR are observed at 60 km (Fig. 10b) due to potential ionospheric effects, the weighting will affect whole occultation profile (in the statistically optimal regions) as the normalized and smoothed SNR does not start at unity. Between 50 and 60 km, where small scale E a r l y R e l e a s e E a r l y R e l e a s e variations in normalized SNR occur, the bending angle is shown to agree more with the background (Fig. 10e) with respect to unweighted approach. The scaling coefficients gradually decrease downwards 40 km and a smooth transition follows to a fixed value of 1, so that the observed profile is no longer weighted with the background. The comparison of retrieved profile converted to dry temperature with respect to official product of CDAAC shows generally a very good agreement below 40 km of mean sea level (geometric) height. The differences are mostly observed within altitudes where the statistical optimization applies.

refractivity and dry Air retrieval
The ionosphere-free and optimized retrievals are converted to profiles of geophysical parameters in the neutral atmosphere. Figure 11 shows that a statistically significant systematic error is expected after inversion of an unoptimized bending angle profile to refractivity. Below 50 km of the impact height, the error can exceed unaccounted ionospheric correction in single-carrier L1 retrieval. The systematic effect is mostly driven by oscillating characteristics in unoptimized bending angle after ionospheric correction that cannot be removed by linear combination of frequencies.
The bending angle profile at a set of impact parameters is inverted to the refractive index as a function of radius by means of the inverse Abel transform (Fjeldbo et al. 1971): where the refractive index equals n = 1 + N·10 -6 and x = rn(r) is a refractional radius. The inverted bending angle profile is extrapolated above the highest measured impact parameter so that the upper limit of the integral becomes 120 km. The definite integral is then evaluated numerically with the use of trapezoidal rule. The geometrical height with respect to the reference ellipsoid can be easily derived by subtracting the local radius of curvature from ray tangent radii of n(r) profile. Further conversions determine the atmospheric profiles of dry pressure (P d ) and dry temperature (T d ). Assuming the ideal gas law under hydrostatic equilibrium, dry pressure is directly proportional to the air density and can be computed by numerical integration of the refractivity profile: where M d is a molar mass of dry air (28.9644 g mol -1 ), k 1 is a refractivity coefficient, g is gravity acceleration, and the universal gas constant R equals 8.31432 J mol -1 K -1 . Typi-cally, the pressure value at the top of the integral is added to avoid bias. The dry temperature can be calculated by ignoring water vapor term in the refractivity expansion. The relation of dry pressure P d and refractivity N yields: The geometric altitude above mean sea level z is computed from the relation: where the undulation u is a difference between reference ellipsoid and the geoid determined with EGM96 model (Lemoine et al. 1998) for the position of the occultation point.

retrIevAl quAlIty of rAdIo occultAtIoN profIleS
The dataset with collected occultation profiles within January and May 2015 is used to assess uncertainty of the retrieval by comparison of atmospheric profiles to official CDAAC products. 347 profiles were localized in the regional domain from predictions based on real-time occultation tables. The distribution of occultation events considered in the comparison is presented in Fig. 2. The radiosonde sites are utilized for collocations with neighboring occultation events to provide independent estimates for the accuracy. Profiles of the atmosphere are compared at common altitudes expressed at mean sea level heights. The resolution of the radio occultation technique allows for high vertical sampling that in average result in 5000 perigee points per profile. Radiosoundings typically provide 50 measurements combined from mandatory pressure levels and significant levels with respect to temperature during the ascend to an altitude of 35 km. In order to preserve the high vertical resolution, altitude levels for comparisons are established based on occultation data.

comparison with cdAAc
The final occultation products from the retrieval chain are the profiles of atmospheric dry pressure and dry temperature. The water vapor contribution is therefore neglected in both CDAAC and ROWNP profiles, which predominantly applies to the troposphere. However, differences in processing scenarios at a number of stages from level 1b phase and amplitude data to level 2 products of atmospheric profiles are expected to introduce uncertainties in compared profiles. The differences in the dry temperature computed as ROWNP minus CDAAC are presented in Fig. 12. The E a r l y R e l e a s e standard deviation with a magnitude of up to 5 K indicates high variability of the differences below 5 km of mean sea level height. A larger systematic effect is observed in the lower troposphere that increases with decreasing the altitude. At the near surface levels, the bias reaches around 2 K. The section in the upper troposphere (UT, 8 to 12 km) shows another example of bias and standard deviation increments, although significantly smaller than those observed in the lower troposphere. The uncertainties reduce and become stable in the lower stratosphere (LS, 12 to 20 km). Above 20 km the standard deviation starts to grow exponentially. This suggests an influence of the increasing ionospheric re-sidual effect and the importance of the optimal weighting of occultation profiles with meteorological background. The differences induced by handling of bending angle profiles above the impact height of 40 km have been already demonstrated to play a significant role in the retrieved geophysical quantities (Fig. 10). The assessed uncertainties in the dry temperature are highly comparable to commonly observed differences in the upper atmosphere (Hajj et al. 2002;Gorbunov 2002b;Li et al. 2013). The fact that the true reference, such as in-situ observations, is difficult to assess at respective geometric altitudes makes the evaluation of results problematic. E a r l y R e l e a s e Similar trend as in the dry temperature in the upper atmosphere can be observed in the fractional differences of dry pressure. However, Fig. 13 indicates a positive bias above 40 km with a magnitude of around 2%. The corresponding standard deviation exceeds 5%. Part of the effect is due to exponentially decreasing pressure values with height, which can lead to noticeable fractional errors even for small differences. The fractional differences stabilize in the UTLS region. The retrieval becomes almost bias-free at the surface with the standard deviation below 1%.

comparison with radiosondes
In order to support the analysis of uncertainty for radio occultation retrievals, an independent quality estimate can be provided from operational in-situ soundings at three radiosonde stations. All sites use a single radiosonde type manufactured by Vaisala (RS92). Standard readings from meteorological sensors consist of total air pressure, air temperature and humidity-related variable of dew point depression. The geophysical variable corresponding to occultation retrievals in the assumption of dry atmosphere can be derived from radiosoundings for the dry pressure by subtracting the water vapor term from the observed pressure so that P d = Pe. The water vapor partial pressure is calculated based on dew point depression using the following equation: where the dew point temperature T d is calculated as the difference between the air temperature T and dew point depression T dd , both given in degree Celsius. The temperatures derived from occultations are compared to "wet" tempera-tures from radiosondes. Therefore, the differences can be expressed as a measure of water vapor. In order to reproduce profiles of refractivity, the two-term expansion with empirical coefficients according to Smith and Weintraub (1953) is applied: where k 1 = 77.6 K hPa -1 and k 3 = 3.73 × 10 5 K 2 hPa -1 . The interpolation from mandatory and significant levels to finer vertical resolution is achieved by cubic spline for the dry pressure and refractivity, and linearly for the temperature.
Prior to the comparison of occultation profiles with radiosondes, the precision of in-situ measurements has been estimated by analyzing the variability of parameters over time with respect to average values. The averages from radiosoundings at each station are calculated at predefined vertical grid to yield the differences for single measurements (i) as ( ) for the temperature, whereas for the fractional differences of dry pressure and refractivity, the observed minus average difference is divided by the observation. Figure 14 summarizes assessed anomalies represented by over 700 observations in the troposphere, which decrease to 500 at the upper levels. The de-trended parameters show altitude-specific variabilities, which depend on meteorological variables. The dry pressure displays increased fluctuations in the UT reaching 3%. The temperature shows contrasting behavior with the standard deviation increasing from 3% at 10 km toward the bounds of the UT altitudes up to 6% at 6 and 12 km. The standard deviation in the refractivity combines the effects of pressure and temperature. Its wavy structure manifests major contributions at the top of UT and in the lower troposphere. The latter one is magnified Fig. 13. As in Fig. 12 but for dry pressure in terms of fractional differences. E a r l y R e l e a s e by water vapor term introduced to refractivity by Eq. (18) but neglected in the analyzed dry pressure. The variability of the refractivity is at the order of 3% at 2 km and nearly 5% at 12 km.
The comparison of radio occultation profiles with radiosondes requires carefully defined criteria in terms of spatio-temporal framework (Sun et al. 2010). The randomly distributed occultation events need to match release times of balloon soundings evenly carried aloft from fixed station locations. This introduces uncertainty due to the variability of meteorological variables over time and space. The effect also stems from ~200 km horizontal drift of radiosonde balloon as it ascends during ~100 min flight (Seidel et al. 2011). Balloons are typically released ~40 min in advance to nominal 00 and 12 UTC. The collocation issue is minimized by 200 km search radius and ±2 hour time difference with respect to official launch times. The applied collocation mismatches affect the representativeness of the data sample relative to number of observations used to analyze the variability of radiosoundings in Fig. 14. In total, 231 sonde observations from three sites have been collected for comparisons with occultation profiles. The data sample corresponds to 37 profiles that met the collocation criteria. The accuracy estimated as RAOB minus RO differences is presented in Fig. 15 separately for CDAAC and ROWNP retrievals.
The dry pressure computed for radiosonde by removing the water vapor partial pressure leads to the overestimation of RO pressure in the lowest troposphere. The negative bias with a magnitude of nearly -5% applies to both CDAAC and ROWNP profiles and reduces at the bottom of the UT. However, slightly larger bias of -1% at 8 km remains in the ROWNP retrieval relative to CDAAC with -0.3%. The most significant systematic effect is observed at the uppermost levels with -2.7 and -1.5% for CDAAC and ROWNP respectively, which is at the expense of the increased standard deviation for CDAAC by 0.7% in the UTLS region. It should be noted that no extrapolation of profiles is applied above the burst height of radiosonde balloons, which contribute to smaller data samples above 20 km and thus, less representative quantities.
A remarkable warm bias is observed in the temperature differences below the UT. The "real" temperatures from radiosondes are larger since RO retrievals do not account for the water vapor. As the humidity decreases with increasing the altitude, the systematic error vanishes at the UT and a negative bias of around -1% develop in CDAAC. The magnitude of bias is larger by a factor of 2 in ROWNP. The introduced spike corresponds to the fluctuation of the standard deviation in CDAAC temperature of 8%, thus leading to a comparable root mean square error.
The refractivity is shown to be free of the significant systematic effect observed in the temperature differences. In fact, the RO refractivity serves as a fundamental variable for the retrieval of the (dry) pressure and temperature through the inversion process. On the contrary, mapping of radiosonde meteorological measurements to the refractivity space can be considered in terms of the forward modeling. Hence, RO uncertainties in the excess phase through the bending angle propagate to refractivity and induce differences in the dry profiles. The radiosonde minus RO differences can be therefore explained by refractivity deficits E a r l y R e l e a s e observed below 5 km in the troposphere since distributions of standard deviation are magnified in the temperature. The proportionality is also visible at the top of the UT. The larger refractivity bias, higher in the UTLS region, in ROWNP relative to CDAAC (-1.5% vs -0.7%) is related to the dry pressure bias. This is again compensated by lower standard deviation of the fractional differences by 0.7%.

SummAry ANd coNcluSIoNS
The retrieval methodology from GNSS radio occultation excess phase to atmospheric profiles of geophysical parameters implemented in the inversion software has been described and validated. The processing utilizes real-time products for FORMOSAT-3/COSMIC mission provided by CDAAC/TACC. Its capability to derive reliable estimates of occultation profiles has been demonstrated on regional domain setting in Central Europe. The inter-comparison with official occultation profiles processed by CDAAC allowed to relate differences in the retrieval chain to the uncertainty of atmospheric parameters. The comparison with in-situ measurements based on collocations of RO events with radiosonde profiles provides the accuracy of occultation soundings.
The inversion of occultation measurements is a complex process composed of multiple steps with sophisticated approaches for handling the input data. A tuning of the retrievals above the UTLS region by statistical optimization of the bending angles is a particular issue in the processing of occultation data. Since the residual ionospheric noise exceeds E a r l y R e l e a s e the atmosphere-induced variability in the upper atmosphere, observations are weighted toward the background in order to provide the most probable estimates. The emphasis in the presented implementation was put on SNR data as an additional weighting quantity, which does not dependent on the background. As expected, inconsistencies arise with respect to official CDAAC products, yet resulting in comparable uncertainty estimates typically observed in RO retrievals above the UTLS altitudes. The differences obtained by comparison of dry profiles combine other effects accumulated throughout the inversion process and cannot be certainly attributed to any particular deficits in the retrieval chain. More statistically representative datasets will be favorable for future validation studies and improvements in the software. The developed system will also serve as a processing platform for occultation data from the up-coming FORMOSAT-7/COSMIC-2 mission once level 1b products become available. Radiosondes, as a long-term and well-established source of reference for the vertical profiling of the atmosphere, were used to support the uncertainties from the inter-comparison of occultation profiles. Based on nearest collocations, a high reliability of RO refractivity could be demonstrated, which expressed in terms of the standard deviation is generally at the order of 5%. However, in the UT region the differences are as low as 1%. The accuracy of retrieved profiles generally result in a lower standard deviation relative to radiosonde-minus-CDAAC differences. However, slightly larger negative refractivity bias was observed, which further propagates to errors in the dry pressure, and to a lesser extent, to the dry temperature. Comparison of such parameters with radiosoundings is problematic in the troposphere due to neglected water vapor term in RO retrievals. The underestimated RO temperature with respect to the true causes a significant positive bias. Although the water vapor pressure was subtracted from the total pressure for radiosonde measurements, the removed humidity contribution introduced a negative bias with a magnitude of around -5% at the surface.