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

Abstract

The Global Navigation Satellite System (GNSS) Radio Occultation (RO) is a key technique for obtaining thermodynamic profiles of temperature, humidity, pressure, and density in the Earth’s troposphere. However, due to refraction effects of both the dry air and water vapor at low altitudes, retrieval of accurate profiles is challenging. Here we introduce a new moist air retrieval algorithm aiming to improve the quality of RO-retrieved profiles in moist air and including uncertainty estimation in a clear sequence of steps. The algorithm first uses RO dry temperature and pressure and background temperature/humidity and their uncertainties to retrieve humidity/temperature and their uncertainties. These temperature and humidity profiles are then combined with their corresponding background profiles by optimal estimation employing inverse-variance weighting. Finally, based on the optimally estimated temperature and humidity profiles, pressure and density profiles 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 profiles 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 find 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.

1. 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]:

$$J\left(x\right)=\frac{1}{2}{\left(x-x_{\mathrm{b}}\right)}^{\mathrm{T}}{B}^{-1}\left(x-x_{\mathrm{b}}\right)+\frac{1}{2}{\left(y_{\mathrm{o}}-H\left[x\right]\right)}^{\mathrm{T}}{O}^{-1}\left(y_{\mathrm{o}}-H\left[x\right]\right),$$

(1)

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 retrieved by a 1DVar algorithm that uses retrieved refractivity and geometric altitude [27] together with background data from ERA-Interim (ERA-I) [28] as input. For each occultation event, the background temperature, specific humidity, and surface pressure are interpolated in time and space at 60 model levels from the ERA-I forecast available at a 3 hour and 1° × 1° grid. The 1DVar configuration is defined by a few choices described in detail in the ROM-SAF 2018 report on Level 2B and 2C 1DVar products [20].

Refractivity uncertainty is parameterized as function of height. In the troposphere it is a straight line fixed at 0.2% at the tropopause and 2% at surface. Above the tropopause it is 0.2% but never below 0.02 N-units. Refractivity vertical correlations are assumed exponential with a 3 km correlation length. Background uncertainty is taken directly from the ERA-I error of first guess estimate provided by ECMWF. The surface pressure uncertainty has been inflated in order to adapt to an evident pressure difference in ERA-I forecast and analysis. The ERA-I vertical uncertainty profiles are averaged into 5° latitude bins while the vertical correlations are provided by ECMWF as a fixed correlation matrix.

CDAAC-provided moist air profiles used in this study are the reprocessed wetPrf data records of the Challenging Minisatellite Payload CHAMP and COSMIC RO missions (CHAMP data: 2016.2430 version, available online via CDAAC website; COSMIC: update of 2013.3520 version, available online via CDAAC in future). At CDAAC, background profiles are taken from ECMWF gridded low-resolution analysis data collocated to RO locations [29]. The observation uncertainties include both systematic and random components, which are latitude-dependent and are estimated from the statistics of innovation vectors for a reasonably long period such as one month, and the information was updated regularly from the statistics for a recent period to yield best-possible performance [30]. The estimation of the background uncertainties at CDAAC is similar to the estimation of observation uncertainties. The correlation matrix was estimated using a fifth-order correlation function, such as used by Steiner and Kirchengast [31], which is similar to a Gaussian function in shape but compactly supported (in a mathematical sense).

The formulation of B and O is critical for the moist air retrieval, since it determines the weights of background and observed data that lead to the formally optimal profiles according to Equation (1). The 1DVar method is successful and retrieved moist profiles have been used in several climate and weather studies and good results were obtained [32,33,34].

In this study, we introduce a “linearized 1DVar” moist air retrieval algorithm as a robust 1DVar alternative that sequentially combines the direct method with optimal estimation. The new method is designed to derive tropospheric temperature, humidity, and pressure profiles at the same quality as 1DVar, and provides a robust linear non-iterative propagation chain, including “direct method” humidity and temperature profile retrievals as interim results, and transparent and comprehensive uncertainty estimates. It includes empirical models of background and observation uncertainties, to optimally determine the weights of background and observed data. The new algorithm was initially motivated, designed, and theoretically derived by Kirchengast et al. [35]. It is introduced here in detail in its current updated form, together with the formulation of its input ingredients, including the uncertainty formulations involved.

The new scheme is implemented since 2013 already—in line with its initial design with refractivity-equation closure for pressure retrieval [35] and in a basic form with static input uncertainty profiles—in the Wegener Center for Climate and Global Change (WEGC) current Occultation Processing System version 5.6 (OPS v5.6). It has shown reliable results for entire climate records in several studies [24,25,36,37,38].

In this study we denote this initial implementation using static uncertainties the “OPSv5.6 approach”, while we denote the updated form that uses dynamic input uncertainties, an equation-of-state closure for pressure retrieval, and forecast-minus-analysis bias correction of background profiles the “dynamic approach”. The advanced inclusions of propagating full covariance matrices as well as estimated systematic uncertainty and observation-to-background weighting ratio profiles, as implemented in the new rOPS system [39,40,41,42,43,44], are beyond the scope of this study and will be introduced in a separate follow-on paper.

Both the OPSv5.6 and the dynamic approach are tested using exemplary ensembles of simulated Meteorological Operational (MetOp) satellite data as well as real-observed CHAMP and COSMIC data 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. Appendix A, Appendix B and Appendix C provide complementary information on aspects of numerical implementation, background bias correction, and vertical error correlations.

2. Methodology—The New Moist Air Retrieval Algorithm

2.1. Algorithm Basis

In dry air condition, the refractivity N can be expressed as N = c₁Rρ_d = c₁(p_d/T_d), where R = 287.06 J kg⁻¹ K⁻¹ is the dry air gas constant [2,35], c₁ = 77.60 K hPa⁻¹ 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:

$$c_1\frac{p_{\mathrm{d}}\left(z\right)}{T_{\mathrm{d}}\left(z\right)}=N\left(T\left(z\right),V_{\mathrm{w}}\left(z\right),p\left(z\right)\right),$$

(2)

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:

$$V_{\mathrm{w}}\left(z\right)=\frac{e\left(z\right)}{p\left(z\right)}=\frac{q\left(z\right)}{a_{\mathrm{w}}+b_{\mathrm{w}}q\left(z\right)},$$

(3)

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:

$$N\left(T\left(z\right),V_{\mathrm{w}}\left(z\right),p\left(z\right)\right)=N\left(z\right)=c_1\frac{p\left(z\right)}{T\left(z\right)}\left(f_0+f_1·V_{\mathrm{w}}\left(z\right)\right),$$

(4)

where f₀ is unity or close to unity and f₁ is close to (c₂/T)/c₁, where c₂ = 3.73 × 10⁵ K² hPa⁻¹ represents the Smith-Weintraub refractivity formula second constant (“wet term”). The exact values of f₀ and f₁ 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₀ = 1 and f₁ = (c₂/T)/c₁ = c_T/T [K] and c_T = c₂/c₁ = 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₀ and f₁ will be very small (order 10⁻³ 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.

Using the R.H.S. of Equation (4) equated with the left-hand-side (L.H.S.) of Equation (2) as basis to explicitly express moist air profiles of T and V_w, we get the following two mutually equivalent forms:

$$T\left(z\right)=T_{\mathrm{d}}\left(z\right)\frac{p\left(z\right)}{p_{\mathrm{d}}\left(z\right)}\left(f_0+f_1·V_{\mathrm{w}}\left(z\right)\right),$$

(5)

$$V_{\mathrm{w}}\left(z\right)=\frac{\frac{p_{\mathrm{d}}\left(z\right)}{p\left(z\right)}T\left(z\right)-f_0·T_{\mathrm{d}}\left(z\right)}{f_1·T_{\mathrm{d}}\left(z\right)}.$$

(6)

Based on hydrostatic integration, the dry pressure and the moist pressure at any given altitude level z can be expressed as:

$$\frac{\mathrm{d}\ln p_{\mathrm{d}}\left(z\right)}{\mathrm{d}z}=-\frac{g\left(z\right)}{R{T}_{\mathrm{d}}\left(z\right)},\mathrm{and}$$

(7)

$$\frac{\mathrm{d}\ln p\left(z\right)}{\mathrm{d}z}=-\frac{g\left(z\right)}{RT\left(z\right)\left(1+c_{\mathrm{w}}q\left(z\right)\right)},$$

(8)

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:

$$d\ln p\left(z\right)=\beta \left(z\right)·d\ln p_{\mathrm{d}}\left(z\right),$$

(9)

$$\mathrm{where} \beta \left(z\right)=\frac{T_{\mathrm{d}}\left(z\right)\left(1+b_{\mathrm{w}}{V}_{\mathrm{w}}\left(z\right)\right)}{T\left(z\right)\left(1+2b_{\mathrm{w}}{V}_{\mathrm{w}}\left(z\right)\right)}.$$

(10)

Since the differential increments dlnp and dlnp_d will be log-linearly discretized over adjacent levels, we can write dlnp = dlnp(z_i) − dlnp(z_i−1) = ln[p(z_i)/p(z_i−1)], 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:

$$p\left(z_i\right)=p\left(z_{i-1}\right){\left(\frac{p_{\mathrm{d}}\left(z_i\right)}{p_{\mathrm{d}}\left(z_{i-1}\right)}\right)}^{\beta \left(z_{i-1/2}\right)},$$

(11)

where $\beta \left(z_{i-1/2}\right)=\frac{T_{\mathrm{d}}\left(z_i\right)+T_{\mathrm{d}}\left(z_{i-1}\right)}{T\left(z_i\right)+T\left(z_{i-1}\right)}·\frac{1+b_{\mathrm{w}}\sqrt{V_{\mathrm{w}}\left(z_i\right)V_{\mathrm{w}}\left(z_{i-1}\right)}}{1+2b_{\mathrm{w}}\sqrt{V_{\mathrm{w}}\left(z_i\right)V_{\mathrm{w}}\left(z_{i-1}\right)}}$ 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.

2.2. 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 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:

$$T_q\left(z_i\right)=T_{\mathrm{d}}\left(z_i\right)\frac{p_q\left(z_i\right)}{p_{\mathrm{d}}\left(z_i\right)}\left(1+\frac{c_T}{T_q\left(z_i\right)}·V_{\mathrm{wb}}\left(z_i\right)\right),\mathrm{and}$$

(12)

$$p_q\left(z_i\right)=p_q\left(z_{i-1}\right){\left(\frac{p_{\mathrm{d}}\left(z_i\right)}{p_{\mathrm{d}}\left(z_{i-1}\right)}\right)}^{{\beta }_q\left(z_{i-1/2}\right)},$$

(13)

where ${\beta }_q\left(z_{i-1/2}\right)=\frac{T_{\mathrm{d}}\left(z_i\right)+T_{\mathrm{d}}\left(z_{i-1}\right)}{T_q\left(z_i\right)+T_q\left(z_{i-1}\right)}·\frac{1+b_{\mathrm{w}}\sqrt{V_{\mathrm{wb}}\left(z_i\right)V_{\mathrm{wb}}\left(z_{i-1}\right)}}{1+2b_{\mathrm{w}}\sqrt{V_{\mathrm{wb}}\left(z_i\right)V_{\mathrm{wb}}\left(z_{i-1}\right)}}$. 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_{Tq}^2\left(z\right)$ can be obtained from propagating the variance profiles $u_{T\mathrm{d}}^2\left(z\right)$ and $u_{q\mathrm{b}}^2\left(z\right)$ based on a linearized version of Equation (12), linearized with some reasonable assumptions (cf. Appendix A):

$$u_{Tq}^2\left(z\right)={\left(\frac{p_q\left(z\right)}{p_{\mathrm{d}}\left(z\right)}\right)}^2{u}_{T\mathrm{d}}^2\left(z\right)+{\left(\frac{p_q\left(z\right)}{p_{\mathrm{d}}\left(z\right)}\frac{T_{\mathrm{d}}\left(z\right)}{T_q\left(z\right)}{c}_{q2T}\right)}^2{u}_{q\mathrm{b}}^2\left(z\right),$$

(14)

where c_q2T = 7727.9 K is the moist air humidity coefficient in temperature error estimation and the square-root of $u_{Tq}^2\left(z\right)$ is the desired uncertainty profile.

Similarly, based on a linearized version of Equation (13), uncertainty of $p_q\left(z\right)$ denoted as $u_{pq}\left(z\right)$ can be calculated as:

$$u_{pq}\left(z\right)={\beta }_q\left(z\right)\left(\frac{p_q\left(z\right)}{p_{\mathrm{d}}\left(z\right)}\right)u_{p\mathrm{d}}\left(z\right),$$

(15)

where ${\beta }_q\left(z\right)=\frac{T_{\mathrm{d}}\left(z\right)\left(1+b_{\mathrm{w}}{V}_{\mathrm{wb}}\left(z\right)\right)}{T_q\left(z\right)\left(1+2b_{\mathrm{w}}{V}_{\mathrm{wb}}\left(z\right)\right)}$.

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:

$$V_{\mathrm{w}T}\left(z_i\right)=\frac{\frac{p_{\mathrm{d}}\left(z_i\right)}{p_T\left(z_i\right)}{T}_{\mathrm{b}}\left(z_i\right)-T_{\mathrm{d}}\left(z_i\right)}{c_T·\frac{T_{\mathrm{d}}\left(z_i\right)}{T_{\mathrm{b}}\left(z_i\right)}},\mathrm{and}$$

(16)

$$p_T\left(z_i\right)=p_T\left(z_{i-1}\right){\left(\frac{p_{\mathrm{d}}\left(z_i\right)}{p_{\mathrm{d}}\left(z_{i-1}\right)}\right)}^{{\beta }_T\left(z_{i-1/2}\right)},$$

(17)

where ${\beta }_T\left(z_{i-1/2}\right)=\frac{T_{\mathrm{d}}\left(z_i\right)+T_{\mathrm{d}}\left(z_{i-1}\right)}{T_{\mathrm{b}}\left(z_i\right)+T_{\mathrm{b}}\left(z_{i-1}\right)}·\frac{1+b_{\mathrm{w}}\sqrt{V_{\mathrm{w}T}\left(z_i\right)V_{\mathrm{w}T}\left(z_{i-1}\right)}}{1+2b_{\mathrm{w}}\sqrt{V_{\mathrm{w}T}\left(z_i\right)V_{\mathrm{w}T}\left(z_{i-1}\right)}}$, where c_T = c₂/c₁ = 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_{qT}^2$ can be propagated from the variance profiles of $u_{T\mathrm{b}}^2\left(z\right)$ and $u_{T\mathrm{d}}^2\left(z\right)$ based on a linearized version of Equation (16) (cf. Appendix A):

$$u_{qT}^2\left(z\right)={\left(\frac{2\frac{p_{\mathrm{d}}\left(z\right)}{p_T\left(z\right)}{T}_{\mathrm{b}}\left(z\right)-T_{\mathrm{d}}\left(z\right)}{T_{\mathrm{d}}\left(z\right)}{c}_{T2q}\right)}^2{u}_{T\mathrm{b}}^2\left(z\right)+{\left(\frac{\frac{p_{\mathrm{d}}\left(z\right)}{p_T\left(z\right)}{T}_{\mathrm{b}}^2\left(z\right)}{T_{\mathrm{d}}^2\left(z\right)}{c}_{T2q}\right)}^2{u}_{T\mathrm{d}}^2\left(z\right).$$

(18)

The square root of $u_{qT}^2\left(z\right)$ 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:

$$u_{pT}\left(z\right)={\beta }_T\left(z\right)\left(\frac{p_T\left(z\right)}{p_{\mathrm{d}}\left(z\right)}\right)u_{p\mathrm{d}}\left(z\right),$$

(19)

where ${\beta }_T\left(z\right)=\frac{T_{\mathrm{d}}\left(z\right)\left(1+b_{\mathrm{w}}{V}_{\mathrm{w}T}\left(z\right)\right)}{T_{\mathrm{b}}\left(z\right)\left(1+2b_{\mathrm{w}}{V}_{\mathrm{w}T}\left(z\right)\right)}$.

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:

$$T_{\mathrm{e}}\left(z\right)=\left(\frac{u_{T\mathrm{b}}^2\left(z\right)}{u_{Tq}^2\left(z\right)+u_{T\mathrm{b}}^2\left(z\right)}\right)T_q\left(z\right)+\left(\frac{u_{Tq}^2\left(z\right)}{u_{Tq}^2\left(z\right)+u_{T\mathrm{b}}^2\left(z\right)}\right)T_{\mathrm{b}}\left(z\right).$$

(20)

Furthermore, its variance profile $u_{T\mathrm{e}}^2\left(z\right)$ can be estimated using the uncertainty propagation law as:

$$u_{T\mathrm{e}}^2\left(z\right)={\left(\frac{1}{u_{Tq}^2\left(z\right)}+\frac{1}{u_{T\mathrm{b}}^2\left(z\right)}\right)}^{-1}=\frac{u_{Tq}^2\left(z\right)u_{T\mathrm{b}}^2\left(z\right)}{u_{Tq}^2\left(z\right)+u_{T\mathrm{b}}^2\left(z\right)},$$

(21)

where the square-root of $u_{T\mathrm{e}}^2\left(z\right)$ 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_{q\mathrm{e}}^2$(z) can be estimated as:

$$q_{\mathrm{e}}\left(z\right)=\left(\frac{u_{q\mathrm{b}}^2\left(z\right)}{u_{qT}^2\left(z\right)+u_{q\mathrm{b}}^2\left(z\right)}\right)q_T\left(z\right)+\left(\frac{u_{qT}^2\left(z\right)}{u_{qT}^2\left(z\right)+u_{q\mathrm{b}}^2\left(z\right)}\right)q_{\mathrm{b}}\left(z\right),\mathrm{and}$$

(22)

$$u_{q\mathrm{e}}^2\left(z\right)={\left(\frac{1}{u_{qT}^2\left(z\right)}+\frac{1}{u_{q\mathrm{b}}^2\left(z\right)}\right)}^{-1}=\frac{u_{qT}^2\left(z\right)u_{q\mathrm{b}}^2\left(z\right)}{u_{qT}^2\left(z\right)+u_{q\mathrm{b}}^2\left(z\right)},$$

(23)

where again the square-root of $u_{q\mathrm{e}}^2$ 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):

$$V_{\mathrm{we}}\left(z\right)=\frac{q_{\mathrm{e}}\left(z\right)}{a_{\mathrm{w}}+b_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)},$$

(24)

and its associated uncertainty can be calculated, based on linearization of Equation (24), according to the uncertainty propagation law:

$$u_{V\mathrm{we}}\left(z\right)=\frac{a_{\mathrm{w}}}{{\left(a_{\mathrm{w}}+b_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)\right)}^2}{u}_{q\mathrm{e}}\left(z\right).$$

(25)

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):

$$p_{\mathrm{e}}\left(z_i\right)=p_{\mathrm{e}}\left(z_{i-1}\right){\left(\frac{p_{\mathrm{d}}\left(z_i\right)}{p_{\mathrm{d}}\left(z_{i-1}\right)}\right)}^{{\beta }_{\mathrm{e}}\left(z_{i-1/2}\right)},$$

(26)

with ${\beta }_{\mathrm{e}}\left(z_{i-1/2}\right)=\frac{T_{\mathrm{d}}\left(z_i\right)+T_{\mathrm{d}}\left(z_{i-1}\right)}{T_{\mathrm{e}}\left(z_i\right)+T_{\mathrm{e}}\left(z_{i-1}\right)}·\frac{1+b_{\mathrm{w}}\sqrt{V_{\mathrm{we}}\left(z_i\right)V_{\mathrm{we}}\left(z_{i-1}\right)}}{1+2b_{\mathrm{w}}\sqrt{V_{\mathrm{we}}\left(z_i\right)V_{\mathrm{we}}\left(z_{i-1}\right)}}$. 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 derived as $p_{\mathrm{e}}\left(z\right)=p_{\mathrm{d}}\left(z\right)\frac{T_{\mathrm{e}}\left(z\right)}{T_{\mathrm{d}}\left(z\right)\left(1+c_2/\left(c_{1}T\left(\mathrm{z}\right)\right)·V_{\mathrm{we}}\left(z\right)\right)}$, 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:

$$u_{p\mathrm{e}}\left(z\right)={\beta }_{\mathrm{e}}\left(z\right)\frac{p_{\mathrm{e}}\left(z\right)}{p_{\mathrm{d}}\left(z\right)}{u}_{p\mathrm{d}}\left(z\right),$$

(27)

where ${\beta }_{\mathrm{e}}\left(z\right)=\frac{T_{\mathrm{d}}\left(z\right)\left(1+b_{\mathrm{w}}{V}_{\mathrm{we}}\left(z\right)\right)}{T_{\mathrm{e}}\left(z\right)\left(1+2b_{\mathrm{w}}{V}_{\mathrm{we}}\left(z\right)\right)}$.

Using Equation (3), the water vapor partial pressure profile e_e(z) can be computed as:

$$e_{\mathrm{e}}\left(z\right)=V_{\mathrm{we}}\left(z\right)p_{\mathrm{e}}\left(z\right),$$

(28)

and its variance profile can be calculated as:

$$u_{\mathrm{e}e}^2\left(z\right)=p_{\mathrm{e}}^2\left(z\right)u_{V\mathrm{we}}^2\left(z\right)+V_{we}^2{u}_{p\mathrm{e}}^2\left(z\right).$$

(29)

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:

$${\rho }_{\mathrm{e}}\left(z\right)=\frac{p_{\mathrm{e}}\left(z\right)}{R{T}_{\mathrm{e}}\left(z\right)\left(1+c_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)\right)},$$

(30)

and, based on a linearized version of Equation (30), its variance profile can be calculated as:

$$\begin{array}{cc}{u}_{\rho \mathrm{e}}^2\left(z\right)={\left(\frac{1}{R{T}_{\mathrm{e}}\left(z\right)\left(1+c_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)\right)}\right)}^2& u_{p\mathrm{e}}^2\left(z\right)+{\left(\frac{p_{\mathrm{e}}\left(z\right)}{R{T}_{\mathrm{e}}^2\left(z\right)\left(1+c_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)\right)}\right)}^2{u}_{T\mathrm{e}}^2\left(z\right)\\ & +{\left(\frac{c_{\mathrm{w}}{p}_{\mathrm{e}}}{R{T}_{\mathrm{e}}\left(\mathrm{z}\right){\left(1+c_{\mathrm{w}}{q}_{\mathrm{e}}\left(z\right)\right)}^2}\right)}^2{u}_{q\mathrm{e}}^2\left(z\right)\hfill \end{array}$$

(31)

Again the square root of the variance profile is the desired uncertainty profile u_ρe(z).

2.3. 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.

2.3.1. 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:

$$s_{\mathrm{model}}\left(z\right)=\{\begin{array}{ll}{s}_0+q_0\left[\frac{1}{z^p}-\frac{1}{z_{\mathrm{Ttop}}^p}\right]& \mathrm{for}z\le z_{\mathrm{Ttop}}\\ s_0& \mathrm{for}{z}_{\mathrm{Ttop}}\le z<z_{\mathrm{Sbot}}\end{array},$$

(32)

where z_Ttop is the top altitude of the troposphere domain, z_Sbot is the bottom altitude of the stratosphere domain, s₀ is the standard error (uncertainty) in the upper troposphere/lower stratosphere domain, q₀ 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. Parameter settings for the observational uncertainty model for dry temperature and dry pressure.

	${\mathit{z}}_{\mathbf{Ttop}}$	${\mathit{z}}_{\mathbf{Sbot}}$	${\mathit{s}}_0$	${\mathit{q}}_0$	p
$u_{T\mathrm{d}}$	10.0 km	20.0 km	0.7 K	3 K kmp	0.5
$u_{p\mathrm{d}}$	10.0 km	17.0 km	0.15%	0.7 %kmp	0.5

.

2.3.2. 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 ${\overline{T}}_{\mathrm{a}}$ and specific humidity ${\overline{q}}_{\mathrm{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 ${\overline{T}}_{\mathrm{f}}$ and specific humidity ${\overline{q}}_{\mathrm{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.5 K). Such values increase with the increase of latitude. Poleward of 60°, s₀ 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:

$$u_{\mathrm{b}}\left(z\right)={\left[{\left(u_{\mathrm{a}}\left(z\right)\right)}^2+{\left(s_{\mathrm{f}-\mathrm{a}}\left(z\right)\right)}^2\right]}^{1/2},$$

(33)

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]:

$$u_{\mathrm{a}\_\mathrm{grid}}={\left[b_{\mathrm{a}}^2+\left(s_{\mathrm{a}}^2/N_{\mathrm{a},\mathrm{f}}\right)\right]}^{1/2}.$$

(34)

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:

$$u_{T\mathrm{b}}\left(z\right)=u_{T\mathrm{b}}\left(z_{T\mathrm{top}}\right)·{\mathrm{e}}^{\frac{z-z_{T\mathrm{top}}}{H_{T\mathrm{b}}}},$$

(35)

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_{q\mathrm{b}}/{\overline{q}}_{\mathrm{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.

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.

2.4. 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 Figure 3 and Figure 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 Figure 3 and Figure 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) above). This makes transparent that the pressure chain of computations is quite simple in terms of uncertainty setup in this version, and quite robust in terms of co-estimating the pressure together with temperature and humidity in a manner consistent with the hydrostatic equation and the equation of state.

3. 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].

3.1. Insights from Individual Event Profiles

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 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.

3.2. 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.

Figure 7, Figure 8 and Figure 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.

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.

3.3. 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_{\mathrm{rbu}}=100\frac{u_{\mathrm{ret}}}{u_{\mathrm{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:

$$r_{\mathrm{obw}}=100\left(1-{\left(\frac{r_{\mathrm{rbu}}}{100}\right)}^2\right)=100\left(1-\frac{u_{\mathrm{ret}}^2}{u_{\mathrm{b}}^2}\right),$$

(36)

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_{\mathrm{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.

4. 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”.

5. 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 Schwarz (2018) [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.