A New Algorithm for the Retrieval of Atmospheric Proﬁles from GNSS Radio Occultation Data in Moist Air and Comparison to 1DVar Retrievals

: The Global Navigation Satellite System (GNSS) Radio Occultation (RO) is a key technique for obtaining thermodynamic proﬁles of temperature, humidity, pressure, and density in the Earth’s troposphere. However, due to refraction e ﬀ ects of both the dry air and water vapor at low altitudes, retrieval of accurate proﬁles is challenging. Here we introduce a new moist air retrieval algorithm aiming to improve the quality of RO-retrieved proﬁles in moist air and including uncertainty estimation in a clear sequence of steps. The algorithm ﬁrst uses RO dry temperature and pressure and background temperature / humidity and their uncertainties to retrieve humidity / temperature and their uncertainties. These temperature and humidity proﬁles are then combined with their corresponding background proﬁles by optimal estimation employing inverse-variance weighting. Finally, based on the optimally estimated temperature and humidity proﬁles, pressure and density proﬁles are computed using hydrostatic and equation-of-state formulas. The input observation and background uncertainties are dynamically estimated, accounting for spatial and temporal variations. We show results from applying the algorithm on test datasets, deriving insights from both individual proﬁles and statistical ensembles, and from comparison to independent 1D-Variational (1DVar) algorithm-derived moist air retrieval results from Radio Occultation Meteorology Satellite Application Facility Copenhagen (ROM-SAF) and University Corporation for Atmospheric Research (UCAR) Boulder RO processing centers. We ﬁnd that the new scheme is comparable in its retrieval performance and features advantages in the integrated uncertainty estimation that includes both estimated random and systematic uncertainties and background bias correction. The new algorithm can therefore be used to obtain high-quality tropospheric climate data records including uncertainty estimation.


Introduction
The Global Navigation Satellite System (GNSS) Radio Occultation (RO) technique is an atmospheric sounding technique providing global-coverage, high-accuracy, and high-precision vertical profiles of Earth's atmosphere [1][2][3][4][5]. The technique uses receivers on Low Earth Orbit (LEO) satellites to receive GNSS signals after they propagated through the atmosphere in limb sounding geometry. Vertical profiling is achieved due to the satellites' orbital motions. As the signals propagate, they are bent due to the atmospheric refractivity gradients.
The accumulated bending angle can be calculated from precise orbit data and the excess phase measurements acquired by tracking the GNSS signals on a LEO satellite. The bending angle profile can in turn be converted to a refractivity profile using an Abel integration. In dry air conditions, atmospheric temperature, pressure, and density profiles can then be retrieved using a refractivity equation, hydrostatic integral, and ideal gas law [2].
In moist air conditions, however, which apply in the troposphere below about 9 km to 16 km, the refractivity is also significantly affected by moisture. In this case, there are four unknown variables, i.e., temperature, pressure, density, and humidity, but only three equations as stated above are available as constraint. This results in a temperature-humidity ambiguity problem [2,[6][7][8] that fundamentally could only be solved by way of occultation technique extension by using higher frequency signals as proposed for microwave occultation [9][10][11]. Most air retrieval algorithms, which are the focus of this study, instead solve this problem for RO by way of data processing extension including background information on tropospheric temperature and/or humidity.
In early moist air retrieval algorithm designs, scientists used a direct method to retrieve tropospheric humidity or temperature profiles, by using background profiles of either temperature or humidity [2,6]. However, this method may induce sub-optimal uncertainty from background data assumed "exactly true". As a more general alternative, the one-dimensional variational (1DVar) method [12,13], also termed optimal estimation method [14], was suggested for moist-air retrieval [15] and further investigated by several studies [7,[16][17][18][19].
The 1DVar method works by finding a maximum likelihood optimal estimate of a vertical atmospheric state profile x, given a set of observations y b and a priori knowledge on a background atmospheric state profile x b as well as the error covariance matrices of both the observation and background information. The 1DVar can be written as a minimization of the following equation [20]: where H[x] denotes a forward operator mapping the state x to the observation space y o . The matrices B and O are background and observation error covariance matrices, respectively, representing the standard uncertainties and correlations of the background data and the observation (plus forward-modeled) data. Minimizing the cost function J(x) by variation of the state x yields the retrieved state x r that minimizes the total deviation against background and observational data. The usual selection of y o in moist-air retrieval by 1DVar is the observed refractivity profile from which temperature, humidity and surface pressure are retrieved as state x r [17][18][19][20]. Currently, the RO data processing centers Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) Data Analysis and Archive Center (CDAAC), University Corporation for Atmospheric Research (UCAR) Boulder, Radio Occultation Meteorology Satellite Application Facility (ROM-SAF), Danish Meteorological Institute (DMI) Copenhagen, and National Oceanic and Atmospheric Administration (NOAA) Center for Satellite Applications and Research (STAR) Maryland, use 1DVar algorithm implementations for their (operational) moist air retrievals [20][21][22][23][24][25]. Both ROM-SAF and CDAAC moist air profiles are used for our evaluation of the new algorithm in this study.
ROM-SAF data used in this study are the latest reprocessed climate data records CDR v1.0, which are available at http://www.romsaf.org. The CDR v1.0 processing is based on ROPP 8.1 [26], with few adaptations. In its products, ROM-SAF Level 2B data provide moist-air profiles. These profiles are and evaluated in a retrieval performance validation with corresponding profile ensembles from 1DVar moist air retrievals provided by ROM-SAF and CDAAC.
The paper is structured as follows: Section 2 describes the new algorithm in terms of the algorithm basis, detailed algorithm steps for profile retrieval and uncertainty estimation, and explains the retrieval scheme and process. Section 3 presents the algorithm evaluation results in terms of its performance for individual RO events as well as of its statistical performance in different latitude bands. Finally, a summary and conclusions are given in Section 4. Appendices A-C provide complementary information on aspects of numerical implementation, background bias correction, and vertical error correlations.

Algorithm Basis
In dry air condition, the refractivity N can be expressed as N = c 1 Rρ d = c 1 (p d /T d ), where R = 287.06 J kg −1 K −1 is the dry air gas constant [2,35], c 1 = 77.60 K hPa −1 is the Smith-Weintraub refractivity formula first constant ("dry term"), and ρ d , p d , and T d are RO-retrieved dry density, dry pressure, and dry temperature, respectively. The profiles p d (z) and T d (z) are derived in RO processing by the so-called dry air retrieval step, using the hydrostatic integral and the equation of state (e.g., [2,41]), and are available as input to the moist air retrieval.
The refractivity formula embodying the dry air equation of state p d /ρ d = RT d , allows formulating the ratio of p d and T d in terms of generic refractivity at any altitude level of z: wherein T and p represent physical (moist) temperature and pressure, respectively, and V w is the water vapor volume mixing ratio. The latter relates to pressure p, water vapor partial pressure e, and specific humidity q as: where a w = 0.622 is the moist air gas constant ratio, b w = 1 − a w = 0.378 is the moist air gas constant ratio complement [35,45]. The generic refractivity on the right hand side (R.H.S.) of Equation (2) denotes any existing type of refractivity relation from Smith-Weintraub type to Thayer type, with any given coefficients [46][47][48][49][50], and can be expressed in the form: where f 0 is unity or close to unity and f 1 is close to (c 2 /T)/c 1 , where c 2 = 3.73 × 10 5 K 2 hPa −1 represents the Smith-Weintraub refractivity formula second constant ("wet term"). The exact values of f 0 and f 1 depend on which refractivity formula is used and whether ideal gas behavior is adopted. As the current OPSv5.6 and rOPS baseline, the standard Smith-Weintraub refractivity formula is used, corresponding to f 0 = 1 and f 1 = (c 2 /T)/c 1 = c T /T [K] and c T = c 2 /c 1 = 4806.7 K. These values are used later on.
Healy [48] conveys that this standard relation continues to be a very good representation and its use keeps parametric consistency with other processing chains using it as well. Kirchengast et al. [35] explain that the perturbations to f 0 and f 1 will be very small (order 10 −3 or smaller) for any more advanced refractivity formulation so that they could be readily added as "epsilon terms" within the step 1a / step 1b iteration algorithm (cf. Section 2.2) if desired. Aparicio and Laroche [50] caution that any use of an advanced refractivity formulation beyond the Smith-Weintraub form should also consistently use a correspondingly advanced equation-of-state formulation accounting for non-ideal gas behavior; an aspect that can as well be accounted for by adding "epsilon terms" in the current algorithm.  (2) as basis to explicitly express moist air profiles of T and V w , we get the following two mutually equivalent forms: Based on hydrostatic integration, the dry pressure and the moist pressure at any given altitude level z can be expressed as: , and where c w = 1/a w − 1 = 0.608 is the moist air humidity coefficient for virtual temperature [35]. By expressing the moist pressure vertical increment dlnp in terms of the dry pressure increment dlnp d , and also using Equation (3) to convert q to V w , we get: where Since the differential increments dlnp and dlnp d will be log-linearly discretized over adjacent levels, we can write , where i represents the corresponding level indices. Similarly, we can write dlnp d = ln[p d (z i )/p d (z i−1 )] for the dry pressure increment. Based on these expressions, we can then derive the expression of p at any altitude level z i as: where β(z i−1/2 ) = represents the exponent of fractional dry pressure change between levels z i and z i−1 that leads to matching this change to the fractional moist pressure change. Since T d is always smaller than T if moisture is non-zero, β is (slightly) smaller than one, expressing that p is changing less than p d , consistent with the fact that p d is always larger than p for non-zero moisture [51]. The specific formulation of β, with temperature expressed as mid-layer linear average (arithmetic mean) between the two levels, and water vapor mixing ratio as mid-layer log-linear average (geometric mean), is found helpful for high numerical accuracy at any given level spacing. Based on these general expressions of Equations (5), (6), and (11), we can either solve for T and p if q is prescribed, or for V w (and hence q via Equation (3)) and p if T is prescribed. We can do this by a simple iteration at any arbitrary altitude level z where a suitably adjacent level has been solved for p before (starting at a "tropospheric top" level with negligible moisture where p d essentially equals p). If q (and hence V w ) is prescribed, then, for any altitude level, we iterate the pairs of Equations (5) and (11) until T has converged to within a small tolerance dT tol , and p will be consistent with the converged T.
Similarly, if T is prescribed, we iterate the pair of Equations (6) and (11) until V w has converged to within a small tolerance (dV w /V w ) tol , and p will then be consistent with the converged V w . This formulation of the "direct method" of moist air retrieval is highly robust and versatile and applicable to arbitrary non-equidistant vertical grids of any level number (from minimum two levels) and vertical range from a chosen "tropospheric top" level to bottom of profile.
Given the resulting temperature and humidity profiles as well as estimates of their uncertainties and of the background profile uncertainties, we may then proceed to an optimally estimated profile for each, temperature and humidity, by combining these profiles with their corresponding background profiles in an inverse-variance-weighted manner. The following section provides more details.

Algorithm Description
The scheme and sequence of the new moist air retrieval algorithm is shown in Figure 1. The method consists of three steps. The first step includes two (formally parallel) sub-steps, which are independent from each other. Given the resulting temperature and humidity profiles as well as estimates of their uncertainties and of the background profile uncertainties, we may then proceed to an optimally estimated profile for each, temperature and humidity, by combining these profiles with their corresponding background profiles in an inverse-variance-weighted manner. The following section provides more details.

Algorithm Description
The scheme and sequence of the new moist air retrieval algorithm is shown in Figure 1. The method consists of three steps. The first step includes two (formally parallel) sub-steps, which are independent from each other. Step 1a is to retrieve temperature and pressure as well as their associated uncertainties with specific humidity and its uncertainty prescribed.
Step 1b is to retrieve specific humidity and pressure as well as their associated uncertainties with temperature and its uncertainty prescribed.
Step 2 is to combine the retrieved temperature profile from step 1a and humidity profile from step 1b with their corresponding background profiles based on inverse variance weighting in order to obtain optimally estimated temperature and humidity profiles. This core step serves to eliminate the effects from sub-optimal estimation using fixed prescribed profiles by optimally weighting retrieved and background profiles each for temperature and humidity so as to arrive at a best estimate consistent with the available input uncertainty knowledge. Step 1a is to retrieve temperature and pressure as well as their associated uncertainties with specific humidity and its uncertainty prescribed.
Step 1b is to retrieve specific humidity and pressure as well as their associated uncertainties with temperature and its uncertainty prescribed.
Step 2 is to combine the retrieved temperature profile from step 1a and humidity profile from step 1b with their corresponding background profiles based on inverse variance weighting in order to obtain optimally estimated temperature and humidity profiles. This core step serves to eliminate the effects from sub-optimal estimation using fixed prescribed profiles by optimally weighting retrieved and background profiles each for temperature and humidity so as to arrive at a best estimate consistent with the available input uncertainty knowledge.
Step 3 then calculates optimal pressure, density, water vapor volume mixing ratio, and water vapor partial pressure profiles based on the optimally estimated temperature and specific humidity profiles from step 2. It uses standard thermodynamic relations for the purpose, such as hydrostatic integration according to Equation (8), and the most air equation of state.
The algorithm, as used in the OPSv5.6 implementation and in this study, focuses on the altitude range from a "tropospheric top" level of 16 km downward and hence covers the entire range where specific humidity may be non-negligibly small [51]. ECMWF operational 24h forecast fields are used to provide prescribed background temperature and specific humidity profiles collocated to the analyzed RO events. Below we introduce these three steps in detail, with some specific but relevant aspects of numerical implementation of steps 1a and 1b described in Appendix A.
Step 1a: Retrieval of temperature and its uncertainty with specific humidity prescribed The input profiles of this step include: the prescribed background specific humidity q b (z) and its uncertainty u qb (z); the observed dry temperature T d (z) and its uncertainty u Td (z); and the observed dry pressure p d (z) and its uncertainty u pd (z). The output profiles include: temperature T q (z) and its uncertainty u Tq (z); pressure p q (z) and its uncertainty u pq (z), where the subscript q denotes variables retrieved with specific humidity prescribed. Using the prescribed q b (z), the corresponding water vapor volume mixing ratio V wb (z) can be calculated using Equation (3). Then, based on Equations (5) and (11), T q (z) and p q (z) at altitude level z i can be expressed as: where β q (z i−1/2 ) = . For each altitude z i from the initial altitude of our moist retrieval z iniMoist = 16 km to the bottom level of the profile, iteration over Equations (12) and (13) yields the profiles of T q and p q . The variance of T q (z) denoted as u 2 Tq (z) can be obtained from propagating the variance profiles u 2 Td (z) and u 2 qb (z) based on a linearized version of Equation (12), linearized with some reasonable assumptions (cf. Appendix A): where c q2T = 7727.9 K is the moist air humidity coefficient in temperature error estimation and the square-root of u 2 Tq (z) is the desired uncertainty profile. Similarly, based on a linearized version of Equation (13), uncertainty of p q (z) denoted as u pq (z) can be calculated as: where β q (z) = Step 1b: Retrieval of specific humidity and its uncertainty with background temperature prescribed The input profiles of this step include: prescribed temperature T b (z) and its uncertainty u Tb (z); observed dry temperature T d (z) and its uncertainty u Td (z); and observed dry pressure p d (z) and its uncertainty u pd (z). The output profiles include the specific humidity q T (z) and its uncertainty u qT (z) and pressure p T (z) and its uncertainty u pT (z). According to Equations (6) and (11), the corresponding water vapor volume mixing ratio profile V wT (z) and p T (z) can be expressed as: , where c T = c 2 /c 1 = 4806.7 K has been described above. V wT (z) and p T (z) can be solved by iterating level by level top-downward from z i = z iniMoist to the bottom level. After obtaining V wT (z), the corresponding specific humidity q T (z) can be calculated using an inverse version of Equation (3) and its variance u 2 qT can be propagated from the variance profiles of u 2 Tb (z) and u 2 Td (z) based on a linearized version of Equation (16) (cf. Appendix A): The square root of u 2 qT (z) is the desired uncertainty profile u qT (z). Similarly, based on a linearized version of Equation (17), the uncertainty of pressure u pT (z) can be calculated as: where Step 2: Optimal estimation of temperature and specific humidity and their uncertainties Based on the retrieved temperature profile T q (z) and its uncertainty u Tq (z) obtained in step 1a and also on the prescribed background temperature profile T b (z) and its uncertainty u Tb (z), the optimally estimated temperature profile T e (z) can be calculated by combining T q (z) and T b (z) based on inverse variance weighting at all altitude levels: Furthermore, its variance profile u 2 Te (z) can be estimated using the uncertainty propagation law as: where the square-root of u 2 Te (z) is the desired uncertainty profile u Te (z). Similarly, using the retrieved specific humidity profile q T (z) and its uncertainty u qT (z) obtained in step 1b, and also on the prescribed specific humidity profile q b (z) and its uncertainty u qb (z), the optimally estimated specific humidity profile q e (z) and its variance u 2 qe (z) can be estimated as: where again the square-root of u 2 qe is the desired uncertainty profile u qe (z).
Step 3: Optimally estimated pressure, density, water vapor volume mixing ratio, water vapor partial pressure, and their associated uncertainties Based on the optimally estimated temperature and specific humidity profiles, the optimally estimated water vapor volume mixing ratio V we (z), pressure p e (z), density ρ e (z), and water vapor partial pressure e e (z) can be calculated quite straightforwardly since the relevant retrieval operators are known. The corresponding uncertainty profiles u Vwe (z), u pe (z), u ρe (z), and u ee (z), can also be calculated using variance-based uncertainty propagation, given the state retrieval operators.
Using the optimally estimated specific humidity profile q e (z), the derived water vapor volume mixing ratio V we (z) can be calculated according to Equation (3): and its associated uncertainty can be calculated, based on linearization of Equation (24), according to the uncertainty propagation law: The optimally estimated pressure profile p e (z) can be calculated using the temperature profile T e (z) and volume mixing ratio profile V we (z) based on Equation (11): . This calculation, started at the z iniMoist level as the previous pressure retrievals, is effectively based on the hydrostatic equation (Equation (7) or (8)) (in the convenient variant available in the context of this algorithm) and provides a pressure profile hydrostatically consistent with the estimated temperature and humidity profiles. We call this a hydrostatic-equation-based closure scheme for the retrieval of pressure to emphasize that it is improved over the refractivity-equation-based closure scheme used in the OPSv5.6 approach. In the OPSv5.6 approach being part of the OPSv5.6 processing system, the pressure profile is T d (z)(1+c 2 /(c 1 T(z))·V we (z)) , which is based on the Smith-Weintraub equation and implies that pressure is consistent with refractivity, temperature, and humidity. Due to errors in the refractivity profile, this pressure profile is somewhat "noisy" against the hydrostatic pressure profile that is fully consistent with the retrieved temperature and humidity.
The uncertainty of the estimated pressure profile can be calculated using a linearized version of Equation (26) in the form: where β e (z) = T e (z)(1+2b w V we (z)) . Using Equation (3), the water vapor partial pressure profile e e (z) can be computed as: and its variance profile can be calculated as: Remote Sens. 2019, 11, 2729 10 of 34 The square root of this variance profile is the uncertainty profile u ee (z). Finally, the density profile ρ e (z) can be derived by using the equation of state in moist air: and, based on a linearized version of Equation (30), its variance profile can be calculated as: Again the square root of the variance profile is the desired uncertainty profile u ρe (z).

Modeling of Observation and Background Uncertainties
The uncertainties of observed and background / prescribed variables are key for determining their weights in the optimally estimated profiles and are therefore critical for providing accurate moist profiles and associated uncertainty estimates. In the new algorithm, we dynamically estimate the background and observation uncertainties. As evident from above, these uncertainties include the background temperature uncertainty profile u Tb (z), the background specific humidity uncertainty profile u qb (z), the observed dry temperature uncertainty profile u Td (z), and the observed dry pressure uncertainty profile u pd (z). We sequentially describe below how we estimated these uncertainties for this study.

Observation Uncertainty Modeling
The observation uncertainty of both observed dry temperature u Td and observed dry pressure u pd are modeled following the empirically derived error model developed by Scherllin-Pirscher et al. [51]. Currently, both the OPSv5.6 and the dynamic approach use this model to estimate the observation uncertainty. In the future rOPS system, the propagated individual-profile based observation uncertainties (and error correlation matrices) will be used [43].
The model structure is the same for both temperature and humidity, only with different parameter settings. We set the parameters based on our tests for moist air retrieval, close to original ones of Scherllin-Pirscher et al. [51]. The vertical structure of the model, needed here only up to the bottom of the stratosphere, is: where z Ttop is the top altitude of the troposphere domain, z Sbot is the bottom altitude of the stratosphere domain, s 0 is the standard error (uncertainty) in the upper troposphere/lower stratosphere domain, q 0 is the best-fit magnitude parameter for the tropospheric model, and p is the associated exponent parameter. The complete parameter settings are summarized in Table 1.

Background Uncertainty Modeling
The calculation follows the approach of Li et al. [52,53] and starts from the preparation of a daily updated global three-dimensional (latitude, longitude, and altitude) background uncertainty fields. The horizontal grid resolution of the uncertainty fields is 10 • latitude × 20 • longitude (center of base cell is 5 • N, 10 • E). The vertical resolution was updated to 100 m level spacing from 0.1 km to 20 km altitude. This construction yields the daily uncertainty fields at a global 18 × 18 × 200 grid. All the mean basic profiles that are required for the calculation of uncertainties are saved for each 10 • × 20 • grid cell center location and on the defined 200 vertical grid levels.
In each 10 • × 20 • grid cell and on the 200-level standard vertical grid, several types of basic variables are pre-calculated and saved for both background temperature and specific humidity. These basic variables include (the same notations of some basic variables are used for both temperature and humidity due to the same calculation method): (1) the mean analysis profile of temperature T a and specific humidity q a ; (2) the standard deviations of the ensemble of analysis values against its mean s a ; (3) bias of the mean analysis profile b a ; (4) the mean forecast profile of temperature T f and specific humidity q f ; (5) the standard deviations of the forecast-minus-analysis difference profiles s f-a ; (6) the number of values in the analysis and the forecast ensemble N a,f . These basic variables except the bias of the mean analysis profile are calculated based on statistical calculation using a large ensemble of forecast and analysis profiles in each grid cell. The details of how to extract the ensemble of profiles on the grid and how to calculate these profiles were described by Li et al. [52,53].
The bias profile b a for temperature is estimated by systematic error modeling according to Li et al. [52,53]. It is applied with no vertical variations but with latitudinal variations. The temperature biases are smallest within the ±40 • latitude band, where the values are equal to the basic mean magnitude of s 0 (0.5 K). Such values increase with the increase of latitude. Poleward of 60 • , s 0 are 20% higher than their basic mean magnitude in the summer hemisphere but twice their mean magnitude in the winter hemisphere [51,52]. The bias profile b a for specific humidity is currently adopted as a relative uncertainty value of 5% of the mean analysis humidity, an educated-guess value.
Using these variables, the background uncertainties u b (representing both for u Tb and u qb ) are estimated as: where u a and s f-a here represent the collocated values obtained from a bi-linear interpolation of their values from the four grid points surrounding the tangent-point location of the given RO event.
Preparing for this, u a at each grid point (denoted for clarity as u a_grid ) is estimated as a combination of the systematic biases and the statistical errors [52,53]: Since background temperature uncertainty u Tb between 10 km and 16 km needs to be penalized to gradually increase in uncertainty at these high tropospheric altitudes, in order to ensure that the observations always safely take increasing weight towards the stratosphere, we modified u Tb from 10 km to 16 km and used an intentional uncertainty increase of the form: where z Ttop is 10 km and H Tb is the "uncertainty scale height" set to 5 km. Background specific humidity uncertainty u qb is input in form of relative humidity values into the scheme. That is, we first use the collocated background specific humidity uncertainty divided by the collocated mean forecast humidity profile, u qb /q b , and then use this relative value to multiply it with the collocated actual background profile in order to obtain the specific humidity uncertainty in absolute values for the algorithm. Figure 2 shows the comparison between OPSv5.6 and the new (dynamic) observation and background uncertainties. It can be seen that u Td increases with decreasing altitude from 0.7 K at 16 km to more than 4 K at the surface. Dry pressure uncertainty stays around 0.2% from 16 km to 10 km and then gradually increases with decreasing altitude to about 1.5% at the surface. As noted above, the observation uncertainties are still used as global static profiles, i.e., used globally in the same way, while in the future rOPS they will be as well used dynamically such as the dynamic background uncertainties discussed next.
Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 35 and then gradually increases with decreasing altitude to about 1.5 % at the surface. As noted above, the observation uncertainties are still used as global static profiles, i.e., used globally in the same way, while in the future rOPS they will be as well used dynamically such as the dynamic background uncertainties discussed next. The third row of Figure 2 shows that the OPSv5.6 uTb is a static global profile, being 2 K at 16 km, 0.6 K at 10 km, and about 1.2 K at the surface. In comparison, the dynamic uTb exhibits latitudinal and altitudinal variations. uTb in polar regions is larger than that in non-polar regions. It is largest in the southern hemisphere polar regions, with values varying from 1.2 K close to the surface to more than Mean profiles and uncertainties of the four input variables, i.e., observed dry temperature (first row), observed dry pressure (second row), background temperature (third row), and background specific humidity (fourth row) on 15 July 2008. For the input observed profiles (first row and second row), which the OPSv5.6 and the dynamic approach share the same uncertainties, the left, middle, and right panels show for the observed mean profile, uncertainty profile as a function of altitude, and the uncertainty profile as a function of altitude and latitude, respectively. For the background profiles (third and fourth row), the left, middle, and right panels show for the background mean profile as a function of altitude, global mean static uncertainty profile of OPSv5.6 approach as a function of altitude, and dynamic background uncertainty as a function of altitude and latitude, respectively. The third row of Figure 2 shows that the OPSv5.6 u Tb is a static global profile, being 2 K at 16 km, 0.6 K at 10 km, and about 1.2 K at the surface. In comparison, the dynamic u Tb exhibits latitudinal and altitudinal variations. u Tb in polar regions is larger than that in non-polar regions. It is largest in the southern hemisphere polar regions, with values varying from 1.2 K close to the surface to more than 4 K above 12 km. In non-polar regions, the values gradually decrease from high latitudinal regions to low latitudinal regions. Furthermore, our sensitivity test results (not shown) indicate that the dynamic u Tb exhibits clear seasonal variations, with largest uncertainty in the polar winter hemisphere.
The fourth row of Figure 2 shows that relative values of OPSv5.6 uncertainty of background specific humidity u qb . It is 10% at the surface, increases to about 40% at 7 km, and then gradually decreases to about 15% at 16 km. The dynamic u qb exhibits clear latitudinal variations, with largest values (>40% from 3 km up to 10 km) in tropical regions, decreasing towards the poles.
We also correct the potential biases that exist in background profiles. This is done by subtracting a background bias profile estimated by using gridded mean forecast profiles minus analysis profiles. It is found that bias-corrected background profiles are useful for getting optimal profiles. For conciseness of the main text of this paper, Appendix B provides more details. Furthermore, we also inspected the vertical correlation structure of background and observation inputs, obtained by constructing error covariance matrices of forecast/observed minus analysis profiles. It is found that the correlations of both background and observation profiles are reasonably small. We hence disregarded to account for these error correlations in the algorithm as introduced in this study; here Appendix C provides more details.

Inspection of Intermediate Variables
In order to provide insight on the characteristics of sub-step results, we illustrate here the input, intermediate, and output variables of the new dynamic approach, using one simulated MetOp (simMetOp) event and one real COSMIC event as representative examples. The MetOp event is simulated under moderate ionospheric conditions. The observational errors represent MetOp/GRAS-type receiving system errors (precise orbit determination (POD) errors, receiver thermal noise, local multipath, clock instabilities), following the proven setting by Steiner and Kirchengast [31], also recently used by Schweitzer et al. [10] and Li et al. [52,53]. The results are shown in Figures 3 and 4.
First focusing on the top row (temperature), it can be seen that u Tb of the simMetOp event, which is located at higher latitudinal regions, is larger than that of the COSMIC event, with values that decrease from 3 K at 16 km to 1 K at 10 km and further to 0.8 K below. For both events, u Tq is smaller than u Tb above 10 km, while below 10 km, u Tq increases quickly and becomes larger than 6 K at bottom altitude levels. The optimally estimated T e is bounded between T b and T q and properly takes more weight from the profile that has ascribed less uncertainty. The differences between T e and the corresponding reference profile are generally smaller than the differences of T b and T q , indicating the effectiveness of the optimal estimation. Comparing dynamic u Tb and OPSv5.6 uncertainty u TdOPSv56 , we can see that u Tb is of similar magnitude as u TdOPSv56 , with values larger at high latitudes and smaller at low latitudes.
Next focusing on the middle row (specific humidity), and first comparing u qb and u qT , we see that u qb is larger than u qT below 8 km for the simMetOp event and below 7 km for COSMIC event. Above 7 km to 8 km, u qT increases quickly to large values. In the optimal estimation, q e takes more weight from the profile with smaller uncertainties in the optimal estimation, and its difference against the reference profile is smaller than the one between q b and q T , again indicating the effectiveness of the optimization. The OPSv5.6 specific humidity uncertainty u qbOPSv56 is a static profile globally, starting from near 20% at bottom altitude levels, increasing to about 40% at 7 km and then gradually decreasing to below 20% at 16 km.
Finally focusing on the bottom row of Figures 3 and 4 (pressure), the resulting profiles and also the uncertainties are seen basically very close to (or even identical to) each other, due to the dominating factor being the input dry pressure profile and its uncertainty (cf. Equations (A2), (A6), (A8), (A13) in Appendix A and Equations (26) and (27)  In the top row, the left panel shows temperature profiles for background temperature (blue), temperature calculated with specific humidity prescribed (green), and the optimally estimated temperature (red). The middle panel shows the estimated uncertainty profiles for the three profiles shown left and for the observed dry temperature (black) as well as the uncertainty of the background temperature from the OPSv5.6 approach (dashed blue). The right panel shows the differences between the temperature profiles shown left and the reference profiles, where the references profiles are the ECMWF co-located analysis profiles. In the three panels of the middle row, the same type of variables is shown as in the upper row, but for specific humidity q ; thus the intermediate variable here is specific humidity with temperature prescribed (subscript "T") and there is no dedicated input uncertainty profile in the middle panel (such as uTd in the upper row). Similarly, the bottom row shows the corresponding variables for pressure, whereby here the intermediate pressures from both humidity prescribed (subscript "q") and temperature prescribed (subscript "T") are shown together with the optimally estimated pressure (subscript "e"), and the middle panel also illustrates the input uncertainty profile of the dry pressure pd. In the three panels of the middle row, the same type of variables is shown as in the upper row, but for specific humidity q; thus the intermediate variable here is specific humidity with temperature prescribed (subscript "T") and there is no dedicated input uncertainty profile in the middle panel (such as u Td in the upper row). Similarly, the bottom row shows the corresponding variables for pressure, whereby here the intermediate pressures from both humidity prescribed (subscript "q") and temperature prescribed (subscript "T") are shown together with the optimally estimated pressure (subscript "e"), and the middle panel also illustrates the input uncertainty profile of the dry pressure p d .

Results-Algorithm Performance Evaluation by Comparison to 1DVar Retrievals
The new algorithm in terms of both the OPSv5.6 and dynamic approach is evaluated using simulated and real observed RO data. In addition to these two approaches, moist-air profiles retrieved at ROM-SAF and CDAAC with 1DVar approaches are used for our comparison. The basis of the comparison is to calculate difference profiles between RO retrieved profiles and their corresponding reference profiles based on which we compare the approaches, first by inspecting individual RO events and second by intercomparing statistical profile ensemble results, including systematic differences and standard deviations in different latitudinal bands.
We focus on the comparison of the retrieved profiles of temperature, specific humidity, and pressure. Other moist-air profiles such as density or water vapor pressure profiles, which can be readily derived using these three variables (see Figure 1 and related description), are found to show similar comparative characteristics and are hence for conciseness not additionally discussed here.
Following the successful basic performance evaluation approach of Li et al. [52,53], the data used for the evaluation include simMetOp and real CHAMP and COSMIC data on 15 July 2008, plus, due to the number of CHAMP RO events from one day being not sufficiently high, also CHAMP data on 14 and 16 July 2008. Co-located ECMWF operational analysis profiles are used as reference profiles for all simulated and real RO events.
The End-to-End GNSS Occultation Performance Simulation and Processing System (EGOPS) version v5.6 [54,55] was used for the forward simulation and retrieval of simMetOp data as well as for the retrieval of the CHAMP and COSMIC profiles by the OPSv5.6 and dynamic approaches. The

Results-Algorithm Performance Evaluation by Comparison to 1DVar Retrievals
The new algorithm in terms of both the OPSv5.6 and dynamic approach is evaluated using simulated and real observed RO data. In addition to these two approaches, moist-air profiles retrieved at ROM-SAF and CDAAC with 1DVar approaches are used for our comparison. The basis of the comparison is to calculate difference profiles between RO retrieved profiles and their corresponding reference profiles based on which we compare the approaches, first by inspecting individual RO events and second by intercomparing statistical profile ensemble results, including systematic differences and standard deviations in different latitudinal bands.
We focus on the comparison of the retrieved profiles of temperature, specific humidity, and pressure. Other moist-air profiles such as density or water vapor pressure profiles, which can be readily derived using these three variables (see Figure 1 and related description), are found to show similar comparative characteristics and are hence for conciseness not additionally discussed here.
Following the successful basic performance evaluation approach of Li et al. [52,53], the data used for the evaluation include simMetOp and real CHAMP and COSMIC data on 15 July 2008, plus, due to the number of CHAMP RO events from one day being not sufficiently high, also CHAMP data on 14 and 16 July 2008. Co-located ECMWF operational analysis profiles are used as reference profiles for all simulated and real RO events.
The End-to-End GNSS Occultation Performance Simulation and Processing System (EGOPS) version v5.6 [54,55] was used for the forward simulation and retrieval of simMetOp data as well as for the retrieval of the CHAMP and COSMIC profiles by the OPSv5.6 and dynamic approaches. The EGOPS software was developed by WEGC for the simulation of occultation observations and also retrieval, based on the simulated or real observed RO observations. Its retrieval subsystem for RO atmospheric profiles is also independently denoted as OPS. The version 5.6 is its most state-of-the-art version. In simulation of occultation observations, users can carefully select observational noise modeling options, antenna types, orbit accuracy settings, and many other choices. Details on the EGOPS/OPSv5.6 simulation and retrieval capabilities can be found in Fritzer et al. [54] and Schwärz et al. [55]. Figure 5 shows the differences between RO retrieved moist profiles and their corresponding reference profiles for three exemplary RO events from simMetOp, CHAMP, and COSMIC for the four approaches evaluated. The three events are intentionally selected to represent at the same time a diversity of latitudes and hence atmospheric conditions, from northern hemisphere middle latitudes (simMetOp) via southern hemisphere polar (CHAMP) to tropical region (COSMIC). EGOPS software was developed by WEGC for the simulation of occultation observations and also retrieval, based on the simulated or real observed RO observations. Its retrieval subsystem for RO atmospheric profiles is also independently denoted as OPS. The version 5.6 is its most state-of-theart version. In simulation of occultation observations, users can carefully select observational noise modeling options, antenna types, orbit accuracy settings, and many other choices. Details on the EGOPS/OPSv5.6 simulation and retrieval capabilities can be found in Fritzer et al. [54] and Schwärz et al. [55]. Figure 5 shows the differences between RO retrieved moist profiles and their corresponding reference profiles for three exemplary RO events from simMetOp, CHAMP, and COSMIC for the four approaches evaluated. The three events are intentionally selected to represent at the same time a diversity of latitudes and hence atmospheric conditions, from northern hemisphere middle latitudes (simMetOp) via southern hemisphere polar (CHAMP) to tropical region (COSMIC). From the simMetOp event we can see that the temperature and pressure differences from the dynamic approach are smaller than those of the OPSv5.6 approach, while specific humidity From the simMetOp event we can see that the temperature and pressure differences from the dynamic approach are smaller than those of the OPSv5.6 approach, while specific humidity differences from the dynamic approach are larger than those from OPSv5.6. For the CHAMP event, temperature differences of the four approaches are generally consistent. Specific humidity differences of the four approaches are generally similar, with the OPSv5.6 and dynamic approaches showing somewhat larger differences from about 4 km to 7 km and CDAAC and ROM-SAF larger ones from about 7 km to 10 km.

Insights from Individual Event Profiles
Pressure differences from OPSv5.6 are largest among the four approaches and exhibit some smaller-scale altitude variations. This is expected according to the algorithm choice (cf. Section 2.2) that optimal pressure in the OPSv5.6 approach is calculated consistent with the Smith-Weintraub formula, rather than with the hydrostatic equation, and is hence affected by the errors/fluctuations of dry temperature. Pressure differences from ROM-SAF are smallest for the CHAMP event and largest for the COSMIC event amongst the four approaches, indicating event-to-event variation in how estimated temperature and humidity play together in yielding pressure profiles.
For the COSMIC event, temperature differences for the four approaches are again rather similar. Specific humidity differences from all four approaches are generally consistent as well, with slightly larger values from the dynamic approach between 5 km to 10 km. While these inspections provide some insights to typical individual-event behavior, a more reliable comparison based on statistical results is needed as discussed below.

Statistical Ensemble Results
In order to investigate the statistical performance of the four approaches in different latitudinal regions, the error statistics in terms of the systematic differences and standard deviations of the retrieved profiles against the reference profiles are calculated in six representative latitudinal regions, comprising global total (90 • S to 90 • N) and five latitude bands (see Figure 6 and its caption).
Remote Sens. 2019, 11, x FOR PEER REVIEW 17 of 35 differences from the dynamic approach are larger than those from OPSv5.6. For the CHAMP event, temperature differences of the four approaches are generally consistent. Specific humidity differences of the four approaches are generally similar, with the OPSv5.6 and dynamic approaches showing somewhat larger differences from about 4 km to 7 km and CDAAC and ROM-SAF larger ones from about 7 km to 10 km. Pressure differences from OPSv5.6 are largest among the four approaches and exhibit some smaller-scale altitude variations. This is expected according to the algorithm choice (cf. Section 2.2) that optimal pressure in the OPSv5.6 approach is calculated consistent with the Smith-Weintraub formula, rather than with the hydrostatic equation, and is hence affected by the errors/fluctuations of dry temperature. Pressure differences from ROM-SAF are smallest for the CHAMP event and largest for the COSMIC event amongst the four approaches, indicating event-to-event variation in how estimated temperature and humidity play together in yielding pressure profiles.
For the COSMIC event, temperature differences for the four approaches are again rather similar. Specific humidity differences from all four approaches are generally consistent as well, with slightly larger values from the dynamic approach between 5 km to 10 km. While these inspections provide some insights to typical individual-event behavior, a more reliable comparison based on statistical results is needed as discussed below.

Statistical Ensemble Results
In order to investigate the statistical performance of the four approaches in different latitudinal regions, the error statistics in terms of the systematic differences and standard deviations of the retrieved profiles against the reference profiles are calculated in six representative latitudinal regions, comprising global total (90°S to 90°N) and five latitude bands (see Figure 6 and its caption).  In order to avoid outliers and ensure an identical RO event ensemble for all four approaches, RO profiles that showed bad quality in any of the processing systems, based on the quality control settings and flags in the respective data files supplied, were rejected from the joint event ensemble. In particular, the quality of retrieved profiles from the OPSv5.6 and dynamic approaches was determined by the OPSv5.6 system specifications [56] and the quality of the profiles retrieved at CDAAC and ROM-SAF by the system they used at their centers [20,26,57]. Figure 6 shows the resulting number of profiles (i.e., number of RO events) available for the joint statistical evaluation for the four approaches in the six latitudinal bands. Given our joint-events selection noted above, the number of profiles above about 5 km to 8 km altitude is the same for the four approaches; only in the lower and middle troposphere below about 8 km there is a number-of-profiles reduction depending on the specific processing systems and their specific criteria to cut the tropospheric penetration of individual moist-air profiles depending on retrieval quality.
In general, the number of profiles available deep into the troposphere from the OPSv5.6 and dynamic approach is somewhat larger than that from CDAAC and ROM-SAF. Furthermore, the number of ensemble members in the five latitude bands varies from about 50 (CHAMP in SHP) to about 500 profiles (COSMIC in SHSM), with all bands enabling reasonable statistics for this initial comparative performance evaluation study among the four approaches. Figures 7-9 illustrate the statistical results for the RO retrieved profiles of simMetOp, CHAMP, and COSMIC, respectively. Figure 7 shows for the simMetOp profiles ensemble that the systematic differences for temperature and specific humidity from the dynamic approach are smaller than those from the OPSv5.6 approach. The best relative improvements are found in tropical regions, where the temperature systematic differences of the dynamic approach are 0.2 K smaller than those of the OPSv5.6 approach, and specific humidity differences are about 15% smaller.
Standard deviations of temperature and specific humidity from the dynamic approach are smaller than or similar to those of OPSv5.6. The propagated uncertainties of temperature, u Te are larger than the statistically estimated standard deviations, which is especially related to the fact that u Te is calculated using u Td (cf. Equations (14) and (21)), which is empirically estimated for real rather than simulated data based on the model by Scherllin-Pirscher et al. [51]. That is, for simulated data, u Td is overestimated since the quality of dry temperature of our simulated data is better than real observed data. The propagated specific humidity uncertainties are of similar magnitude compared to the statistically estimated uncertainties.
The statistical differences of pressure for the dynamic approach are clearly smaller than those from OPSv5.6, which is due to the different closure-scheme of the pressure computation as discussed above (Section 2.2). The propagated pressure uncertainties are larger than the statistically estimated uncertainties, similar to temperature, which is similarly related to the overestimation of the uncertainty of dry pressure, targeted to real data, for these simulated MetOp data. Figure 8 shows for the CHAMP profiles ensemble that the temperature error statistics in terms of systematic differences and standard deviations from the dynamic, OPSv5.6, and ROM-SAF approaches are rather similar in all latitudinal bands, with systematic differences reaching around ±0.2 K and standard deviations being smaller than 1 K down to the boundary layer. Temperature statistics of CDAAC are as well rather similar to the other three approaches, with standard deviation only slightly larger below about 8 km. This slightly larger standard deviation of CDAAC is probably due to a somewhat stronger weighting of observations vs. background at low to middle troposphere levels, where observations are nosier compared to background. Furthermore, it needs to be kept in mind that the reference profiles are for mean RO event locations, while actual tangent points drift during occultation (e.g., [19]), which also contributes to enlarged deviations at lowest tropospheric levels. The propagated temperature uncertainties u Te are basically consistent with the statistically estimated uncertainty, which again indicates the reasonableness of this simplified uncertainty propagation. Remote Sens. 2019, 11, x FOR PEER REVIEW 19 of 35    The specific humidity error statistics from the four approaches are basically consistent, on a global-mean scale, with systematic differences reaching around ±10 % and standard deviations varying from 10 % to 50 %. In more detail, the specific humidity statistics reveal clear altitudinal variations, with largest values in the upper troposphere (about 5 km to 10 km), and also latitudinal variations. For example, in the NHP region, error statistics of ROM-SAF and CDAAC are found to be comparatively larger above about 10 km, which is possibly related to the treatment of very small specific humidity values, where WEGC (dynamic, OPSv5.6) places a stronger constraint towards the The specific humidity error statistics from the four approaches are basically consistent, on a global-mean scale, with systematic differences reaching around ±10% and standard deviations varying from 10% to 50%. In more detail, the specific humidity statistics reveal clear altitudinal variations, with largest values in the upper troposphere (about 5 km to 10 km), and also latitudinal variations. For example, in the NHP region, error statistics of ROM-SAF and CDAAC are found to be comparatively larger above about 10 km, which is possibly related to the treatment of very small specific humidity values, where WEGC (dynamic, OPSv5.6) places a stronger constraint towards the background. Therefore, WEGC profiles have less deviation from the (ECMWF-based) background humidity and also from the ECMWF analysis reference here.
The propagated uncertainty of the specific humidity, u qe , is generally smaller below about 10 km than the statistically estimated uncertainties (except in the SHP region, where absolute moisture content is low). This is due to the reason that observed, background, and analysis specific humidity jointly represent more (noisy) variations over the troposphere than the simplified error estimates used here do capture. In other words, part of the "representativeness error" is not captured.
Pressure error statistics below 10 km from all four approaches are basically rather consistent, with systematic differences varying within about ±0.1% and standard deviations reaching about ±0.2%. Pressure differences of the OPSv5.6 approach are somewhat larger and exhibit comparatively more fluctuations below 5 km (below 2 km for systematic differences), which is due to the different refractivity-based (rather than hydrostatic-based) closure scheme explained in the algorithm description above (Section 2.2).
The propagated uncertainties of pressure are similar to the statistically estimated uncertainties above about 5 km. Below 5 km, the propagated uncertainties are found somewhat larger than the statistically estimated uncertainty. This is due to the reason that the simplified input uncertainty estimate used here for dry pressure, which is the basis for calculating the propagated uncertainty, is rather large at low altitudes. In WEGC's future rOPS context, this input uncertainty will be more realistic as part of a full chain of uncertainty propagation [44]. Figure 9 shows for the COSMIC profiles that the comparative performances of the four approaches are similar to what is visible and was just discussed for the CHAMP profiles shown in Figure 8. We note that in the tropical lower troposphere below about 3-4 km, where moisture content is largest, the OPSv5.6 and dynamic approaches exhibit a negative systematic difference of up to around −10% in specific humidity. Similar to other differences between the approaches visible in the upper troposphere, this is likely related to different uncertainty weighting choices and points to room for further refinement in future.

Simple Observation-to-Background Weighting Ratio Profiles and Comparative Results
In order to know how much observation information was used in the retrieved moist profiles, we calculate observation-to-background uncertainty weighting ratio profiles, r obw , of temperature and specific humidity for the dynamic, OPSv5.6 and ROM-SAF approaches (CDAAC-provided moist profile files do not contain uncertainty information, hence these data are not included here). Since both our approaches and the 1DVar approach used by ROM-SAF are not (fully) linear, it is not an easy task to calculate the real r obw in a comparable manner. Hence we implemented and inspected an approximation as follows.
Considering that the calculation of retrieval-to-background uncertainty ratio (r rbu ) is straightforward and consistently possible for all four datasets, we used this ratio to calculate an approximate r obw . The retrieval-to-background uncertainty ratio is defined and calculated as r rbu = 100 u ret u b , where u ret is the uncertainty profile of the retrieved (optimally estimated) profile and u b is the corresponding background uncertainty profile. Based on the r rbu profile, we can then estimate r obw as: which implicitly assumes an inverse-variance weighted combination of observations and background in the optimal estimation, being a reasonable approximation. Figure 10, illustrating the r obw results for temperature and specific humidity, shows that the r obw for temperature from dynamic, OPSv5.6, and ROM-SAF approaches are generally consistent, with dynamic r obw comparatively largest. This indicates that, by its observational and background uncertainty choices, the dynamic approach uses more observation information in the temperature retrieval. The temperature results of the three approaches reveal clear latitudinal variations, with less observation information used in the tropics (TRO), where humidity is large, and more in polar regions (especially SHP), where humidity is small and RO likewise accurate.   The r obw profiles of specific humidity from all three approaches are generally consistent, with values from ROM-SAF approach largest (except for NHSM and NHP region), in line with how the temperature weighting is tentatively going the other way. That is, higher relative weight on temperature will generally lead to lower relative weight on humidity, and vice versa. In NHSM and NHP region, values from OPSv5.6 are largest probably due to larger background errors (cf. Figure 2) and subsequent more observation information used in the final optimal estimation.

Discussion
The results of Section 3 and Figure 5, Figure 8, and Figure 9 provide clear evidence that our new "simplified 1DVar" approach with the "direct-retrieval method" results obtained as intermediate step, both in form of the OPSv5.6 and the dynamic approaches, provide (at least) the same level of quality as the ROM-SAF and CDAAC "full 1DVar" approaches. Especially when comparing the results from the OPSv5.6 and ROM-SAF approaches, which use the most similar background uncertainties, the error statistics of these two approaches are generally very close.
In general, the background and observation uncertainties are key for determining the weights of background and observations. Hence mainly the weights of background and observations in the optimal estimation determine the statistical errors of moist profiles, which depend on RO data processing centers' evaluation of the quality of background and observation data. Comparing the OPSv5.6 and the dynamic approach, we find that the dynamic approach has reduced the systematic differences of temperature and the statistical errors of humidity as well as improved the thermodynamic consistency of the pressure results, which confirms the effectiveness of improvements of the dynamic approach on top of the OPSv5.6 approach.
The inter-comparison results here demonstrate the performance and general accuracy of the new algorithm in its basic form, without accounting for error correlations and further uncertainty propagation advancements. Schwarz [44] advanced the dynamic approach by propagating estimated random uncertainties using covariance propagation, controlled by Monte-Carlo ensemble methods. The covariance propagation, accounting for error correlations, also enabled to implement a full covariance-weighted optimal estimation. These specific most recent advancements are published elsewhere, together with a further step of performance evaluation on its added value.
Overall it is already clear from the results of this study that the "simplified 1DVar", with its special features of step-by-step transparency of state retrieval as well as systematic and random uncertainty propagation, is a viable new algorithm achieving the quality of "full 1DVar".

Conclusions
In this study, a new sequentially linearized "simplified 1DVar" algorithm was introduced that combines the so-called direct method, with temperature or humidity prescribed, with optimal estimation, for providing accurate temperature, humidity, and pressure profiles from RO in the troposphere. It was also evaluated using the "full 1DVar" algorithm implementations from the ROM-SAF and CDAAC processing centers.
While approximating the matrix inversion and iteration approach used in 1DVar algorithms in simplified form, we find the new algorithm to nevertheless effectively allow retrieving accurate optimally estimated profiles, along with systematic and random uncertainty propagation and effective observation-to-background weighting ratio tracking. The direct-method retrieval results, temperature profiles with background humidity profiles prescribed as well as humidity profiles with background temperature profiles prescribed, are available as intermediate results and can hence be considered a useful by-product.
The uncertainties of background and observational variables are dynamically estimated in the new algorithm, using statistical calculations and empirical modeling. The estimated uncertainties account for latitudinal and seasonal variations. Residual biases in background profiles (ECMWF short-range forecast profiles) are corrected for by using co-located ECMWF forecast-minus-analysis difference bias profiles and are found useful in reducing biases in resulting profiles.
The comparison of the new algorithm against the moist-air profiles provided by the current OPSv5.6 processing system and profiles from ROM-SAF and CDAAC showed that it provides robust and high quality temperature, humidity, and pressure profiles in the troposphere, comparable in performance with "full 1DVar", plus uncertainty estimates in good quality.
The new algorithm was implemented in the OPSv5.6 system with static uncertainty profiles as an initial scope, while the further advanced dynamic approach presented in this paper, is using dynamic uncertainties and the further improvements described. In future, the algorithm in a further advanced form, based on the work by   [44], will be used as part of the WEGC's new rOPS processing system. This rOPS-implemented moist-air algorithm that is built on the algorithm introduced in this study, is also used in the first large-scale reprocessing towards a tropospheric climate data record 2001-2019 by the rOPS and its integrated uncertainty propagation.  Figure A1. Differences between RO retrieved profiles and ECMWF co-located analysis profiles obtained from using the bias-corrected background profiles (red, "Bias Corr") and from using the original background profiles (black, "No Bias Corr") profiles for three exemplary RO events from simMetOp (upper), CHAMP (middle), and COSMIC (bottom) from 15 July 2008.

Bias correction effects illustrated for test-day ensemble of temperature and humidity profiles
The effects of the background bias correction scheme are investigated as well statistically, again by comparing the retrieval results between those obtained using the bias-correction retrieval results and those obtained without using the bias-correction. The uncertainties used for the retrieval examples illustrated here are the dynamic uncertainties. In Figure A2, from left to right panels, statistical results are shown for the simMetOp, CHAMP, and COSMIC missions.
In particular in tropical regions a good quality of humidity retrievals can be rather challenging so that bias correction is expected to be most helpful in such conditions. Indeed, from the results in Figure A2 we can see that the bias-correction scheme can obviously reduce the biases in retrieved moist profiles, especially for the humidity profiles in tropical regions, where the amount of moisture is significant and the humidity profiles are more readily biased. Figure A1. Differences between RO retrieved profiles and ECMWF co-located analysis profiles obtained from using the bias-corrected background profiles (red, "Bias Corr") and from using the original background profiles (black, "No Bias Corr") profiles for three exemplary RO events from simMetOp (upper), CHAMP (middle), and COSMIC (bottom) from 15 July 2008.

Bias correction effects illustrated for test-day ensemble of temperature and humidity profiles
The effects of the background bias correction scheme are investigated as well statistically, again by comparing the retrieval results between those obtained using the bias-correction retrieval results and those obtained without using the bias-correction. The uncertainties used for the retrieval examples illustrated here are the dynamic uncertainties. In Figure A2, from left to right panels, statistical results are shown for the simMetOp, CHAMP, and COSMIC missions.
In particular in tropical regions a good quality of humidity retrievals can be rather challenging so that bias correction is expected to be most helpful in such conditions. Indeed, from the results in Figure A2 we can see that the bias-correction scheme can obviously reduce the biases in retrieved moist profiles, especially for the humidity profiles in tropical regions, where the amount of moisture is significant and the humidity profiles are more readily biased.

Appendix C. Vertical correlations of observations and background errors
The vertical correlations of the input parameters were also investigated in order to understand the level of approximation if they are disregarded in the current OPSv5.6 and dynamic approach implementations of the new algorithm. The error correlation matrix was calculated by constructing a global-mean error covariance matrix using all the ensemble of difference profiles between the forecast/observed and analysis profiles. By dividing the values in the error covariance matrix by the corresponding squared-uncertainty values (variances) in the matrix diagonal, the error correlation matrix can be obtained [e.g., 52]. Figure A3 shows the correlation matrix (left), exemplary correlation function (middle), and correlation length (right) for the four input parameters, i.e., observed dry temperature, observed dry pressure, background temperature, and background specific humidity. The data shown are mainly

Appendix C Vertical Correlations of Observations and Background Errors
The vertical correlations of the input parameters were also investigated in order to understand the level of approximation if they are disregarded in the current OPSv5.6 and dynamic approach implementations of the new algorithm. The error correlation matrix was calculated by constructing a global-mean error covariance matrix using all the ensemble of difference profiles between the forecast/observed and analysis profiles. By dividing the values in the error covariance matrix by the corresponding squared-uncertainty values (variances) in the matrix diagonal, the error correlation matrix can be obtained (e.g., [52]). Figure A3 shows the correlation matrix (left), exemplary correlation function (middle), and correlation length (right) for the four input parameters, i.e., observed dry temperature, observed dry pressure, background temperature, and background specific humidity. The data shown are mainly from 15 July 2008. However, in order to show the variations of correlations with day of month, we also show the correlation functions and correlation lengths from the 5th and 25th of July. Figure A3 shows that both correlation functions and correlation lengths show little variations with day of month. Correlation functions of all the four parameters are close to Gaussian shape at the main peak. From the main peak outwards, the correlation functions of observed dry temperature, background temperature, and observed dry pressure have some negative side peaks, while the functions of background specific humidity are all positive. Except the correlation lengths of observed pressure being slightly larger, with values varying around 2 km, the correlation lengths of the other three parameters are limited to about 1 km to 1.5 km.
These results indicate that except the observed pressure, the correlations of the other three parameters are not significant. Therefore, in the new algorithm as presented here, it was considered reasonable to disregard the correlations of the input variables within the scope of this study. Further advancements that include the full covariance formulation and propagation in the algorithm are described in a separate follow on paper, based on initial descriptions in Kirchengast et al. [40] and Schwarz [44].
with day of month. Correlation functions of all the four parameters are close to Gaussian shape at the main peak. From the main peak outwards, the correlation functions of observed dry temperature, background temperature, and observed dry pressure have some negative side peaks, while the functions of background specific humidity are all positive. Except the correlation lengths of observed pressure being slightly larger, with values varying around 2 km, the correlation lengths of the other three parameters are limited to about 1 km to 1.5 km. Figure A3. Correlation matrices (left), exemplary correlation functions at three exemplary altitude levels of 11 km, 7 km, and 3 km (middle), and estimated correlation length for correlation functions (right) for the observed dry temperature uncertainty (first row), observed dry pressure uncertainty (second row), background temperature uncertainty (third row), and background specific humidity uncertainty (fourth row). The correlation matrices are shown for 15th July 2008 only, and the correlation functions and correlation lengths are shown for 5th, 15th, and 25th of July 2008. Figure A3. Correlation matrices (left), exemplary correlation functions at three exemplary altitude levels of 11 km, 7 km, and 3 km (middle), and estimated correlation length for correlation functions (right) for the observed dry temperature uncertainty (first row), observed dry pressure uncertainty (second row), background temperature uncertainty (third row), and background specific humidity uncertainty (fourth row). The correlation matrices are shown for 15th July 2008 only, and the correlation functions and correlation lengths are shown for 5th, 15th, and 25th of July 2008.