The COSMIC/FORMOSAT-3 RO products are obtained from the University Corporation for Atmospheric Research (UCAR) COSMIC Data Analysis and Archive Center (CDAAC). The observed RO refractivity is obtained from the Level-2 wetPrf products, along with the retrieved temperature, pressure, and water vapor partial pressure at every 0.1 km of altitude, between the surface level and 20 km. The observed refractivity included in the wetPrf files is the same product as in the atmPrf files, provided here at the same height levels as the other thermodynamic products. These observations are collocated with the European Center for Medium range Weather Forecast (ECMWF) ERA-Interim (e.g., Dee et al., 2011), the ECMWF high-resolution operational analysis, and the National Centers for Environmental Prediction (NCEP) operational analysis, the Global Forecast System (GFS) (NOAA/NCEP, 2003). These collocated profiles are obtained also at the CDAAC in the Level-2 eraPrf, echPrf, and gfsPrf products. The RO products are interpolated into the analyses height levels when the comparisons are performed.

Data from the TRMM and GPM precipitation missions are obtained from the NASA Goddard Earth Sciences Data and Information Services Center (GES DISC). The TRMM data used here are the Level-2 orbital 2B31 products, which provide vertical structure information of precipitation and have a limited swath coverage. The GPM data used are the final run of the Integrated Multi-satellitE Retrievals for GPM (IMERG) products (Huffman et al., 2018), which provide a surface rain rate for the region comprised between ±60^∘ latitude. In order to assess the precipitation intensity and structural characteristics, data from the vertically profiling TRMM radar are used, while the GPM IMERG data are used to increase the statistics. The TRMM 2B31 products provide precipitation information for the region sensed by the TRMM Precipitation Radar (PR), such as rain rate, with a swath of approximately 250 km, a horizontal resolution of 5×5 km, and a vertical resolution of 250 m. The IMERG product provides an estimate of the global surface precipitation every 30 min with a horizontal resolution of 0.1^∘ latitude × 0.1^∘ longitude.

For this study, the precipitation information comes solely from the TRMM and GPM-retrieved products, and no precipitation information is used nor assumed from the analyses and reanalyses. Therefore, the analyses and the reanalysis products might or might not be associated with different precipitation conditions when their products are generated. However, it is not the aim of this work to assess the ability of analyses to reproduce precipitation but to evaluate and compare the RO products with their provided thermodynamic fields for a given location and time in the actual presence or absence of precipitation. Nevertheless, we consider that this could solely lead to minor effects, since the water vapor field is spatiotemporally smoother than the cloud water content field.

2.1 Collocations with the TRMM PR

The COSMIC/FORMOSAT-3 RO products between 2006 and 2015 were collocated with TRMM orbital products. A total of 16 881 COSMIC RO soundings are identified to be within the swath of the TRMM precipitation measurements (250 km) and within ±15 min. After 2013, the number of COSMIC/FORMOSAT-3 RO observations dropped significantly. However, the quality and distribution of the observations was not affected. These events were then classified depending on the presence or not of precipitation and its intensity. Henceforth, each event is linked to the number of pixels of the TRMM radar with a reflectivity (Z) larger than 30 dBZ, used as a proxy for heavy precipitation events, in the surroundings (100 km) of the occultation location within the radar swath.

For each event with evidence of precipitation in its surroundings, the approximated RO ray trajectories have been simulated using ray-tracing techniques and geo-located together with the radar-retrieved three-dimensional reflectivity observations. Therefore, it is possible to interpolate the precipitation information into the set of RO ray trajectories. An example of such interpolation is shown in Fig. 1. We can therefore estimate the amount of precipitation crossed by each of the rays, estimate its effect, and compare it with the actual RO observables such as the excess phase (or the Doppler shift), the signal-to-noise ratio (SNR), or the atmospheric vertical retrievals. We use this information to assess the impact of precipitation on the RO signal propagation and its retrievals, as described in Sect. 4.

Figure 1Three-dimensional collocation of an RO event with a TRMM PR precipitation measurement. It corresponds to the coincidence between the C004.2006.329.22.20.G19 RO event and the 2B31.20061125.51450 TRMM PR product. Here, the precipitation structure is shown in a three-dimensional grid, along with the set of RO ray trajectories (in gray). Black stars indicate the tangent point of the rays. Only a few rays are shown for illustration purposes. The orange dashed lines indicate the edges of the TRMM PR swath. The precipitation information (rain rate) interpolated into the RO plane is shown in the two-dimensional projection in the latitude–height plane.

Download

It is worth mentioning that in this study, we focus on the effect of rain and hydrometeors large enough to exhibit a significant reflectivity signature in the TRMM radar retrievals (working at Ku band), which in its turn are the ones we expect to have the largest impact on RO retrievals. The scattering effects of smaller particles, especially above the melting layer, have larger uncertainty and must be treated carefully.

A clear positive refractivity bias is observed between ∼3 and ∼8 km of altitude when precipitation is present in the occultation position, with respect to the refractivity from weather analyses and reanalyses. In Fig. 2 the bias is shown, for the comparison between the GPM IMERG collocated RO products and the three different analyses and reanalyses introduced in Sect. 2. In this case, the data are separated according to the amount of rain in the surroundings: events with no rain (no-rain profiles) and events where 〈R〉>10 (mm h⁻¹) in the $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$ surrounding area.

Figure 2Fractional difference between the RO-observed refractivity and that from (blue) ERA-Interim, (red) ECMWF high-resolution analysis, and (green) NCEP GFS operational analysis. The RO profiles are divided into no rain (solid) and heavy rain (dashed) based on the collocated GPM IMERG precipitation measurements (see Sect. 2.2). These data correspond to 2016.

Download

While the bias is clearly seen for the two analyses and one reanalysis used in the comparison, their performance within heavy precipitation is also different. When precipitation is not present close by the RO sounding, the RO refractivity and that of analyses and reanalyses agree with each other (i.e., no significant bias) as well as among themselves.

In Fig. 3 we show the regional dependence of the bias, at a height of 6 km. Through the paper we focus on the height range around 6 km because this is where the bias is maximal, as seen in Fig. 2. Here, the globe is divided in hexagons with a diameter of approximately 20^∘, and the events are separated according to their 〈R〉: 〈R〉=0 mm h⁻¹; $\mathrm{0}<\langle R\rangle <\mathrm{2.5}$ mm h⁻¹; and 〈R〉>2.5 mm h⁻¹. This separation is shown at each different column, while the rows separate the analyses or reanalyses used in the comparison. The size of the hexagons is chosen so that all of them contain a significant number of observations, and spatial patterns are clearly seen. Only those bins with a minimum of 25 observations inside them are shown, and the typical range of observations inside the bins is between 1000 and 7000 observations per bin for the no-rain scenarios, between 200 and 1600 observations per bin in the low-rain regime, and between 25 and 150 observations per bin in the heavy-rain regime. This figure shows how the positive bias is present globally under heavy precipitation, although it is larger in certain regions and it depends on the analysis in use. Common features for all three reanalyses are, for example, the positive bias under heavy precipitation that is present in the west Pacific warm pool, the eastern part of the Pacific, the Indian Ocean, the equatorial part of the Atlantic, and over South America and central Africa. These regions are associated with extreme precipitation features (Liu and Zipser, 2015), either to large extension precipitation events or to precipitation systems with a high deep convective cores.

Figure 3Regional averaged fractional difference between the RO-observed refractivity and that from (a, b, c) ERA-Interim: (d, e, f) ECMWF high-resolution analysis and (g, h, i) NCEP GFS operational analysis for a height of 6 km. The compared profiles are classified as belonging to no-rain profiles (a, d, g; 〈R〉=0 mm h⁻¹), low and moderate precipitation (b, e, h; $\mathrm{0}<\langle R\rangle <\mathrm{2.5}$ mm h⁻¹), and heavy-rain profiles (c, f, i; 〈R〉>2.5 mm h⁻¹). The grid corresponds to hexagons with a diameter of about 20^∘. Only those with a minimum of 25 observations inside them are shown.

Download

Besides the positive bias in the region above an altitude of 4 km, a negative bias is also clearly observed below 3 km, both for the rainy and no-rain events. This bias is not assessed here, since it has already been discussed previously in other studies (e.g., Ao et al., 2003; Sokolovskiy, 2003; Xie et al., 2006, 2012; Wang et al., 2017). Similarly, other potential sources of bias have been checked, for example, the angle of incidence of the occultation ray to the receiver, with respect to the transmitter position. The larger the angle, the larger the tangent point drift. This implies that the theoretical spherically symmetric atmosphere could depart from a realistic approximation and induce errors in the retrievals (Foelsche et al., 2011). Also, large incident angles correspond to low SNRs, which could be introducing positive biases (Sokolovskiy et al., 2010). Therefore, the positive bias has been checked grouping the occultation events according to its azimuth angle, in addition to rain variables. The results (not shown) reveal no significant changes to the positive N bias and confirm that the RO observation geometry is not a contributing factor to the positive bias.

Once other known observational issues are discarded as plausible sources of the bias, the influence of the scattering term in Eq. (1) is assessed. In order to further investigate its importance, we have directly simulated the contribution of the liquid and solid water to the excess phase. This is accomplished using three-dimensional collocations between the COSMIC/FORMOSAT-3 RO observations and the TRMM PR measurements, which have allowed us to perform realistic simulations of the effects of precipitation in actual RO observables (see Fig. 1). This represents a novel approach to the assessment of the positive refractivity bias with respect to previous studies.

The contribution from precipitation on the phase delay of the signal is due to the scattering of the propagating wave by nonspherical raindrops. The delay induced by raindrops (or frozen hydrometeors) with respect to that of free space can be linked to the scattering term of refractivity in Eq. (1). For the case in this study, the coherent propagation of plane waves is described as the sum of the effects of all the raindrops in a unit volume with various sizes. Formally, the scattered field can be expressed as

where Eⁱ is the incident field, E^s is the scattered field, and T is the “transmission matrix” describing the characteristics of the rain medium (Oguchi, 1983). The propagation through rain can be considered as a propagation through an effective medium with two characteristic axes, characterized by the two eigenvalues of T, λ₁ and λ₂:

where l is the propagated distance.

Raindrops fall following gravity and are flattened due to the air drag, becoming approximately oblate-shaped (e.g., Pruppacher and Beard, 1970; Beard and Chuang, 1987). Here, we do not take into account the canting angle effect (raindrops being tilted by wind), for simplicity and because in this situation its effect is secondary. Therefore, ${\mathit{\lambda }}_{\mathrm{1},\mathrm{2}}=-i{k}_{\mathrm{eff}}^{\mathrm{h},\mathrm{v}}$, where the k_eff is the effective propagation constant of the medium, that is, a complex number, and 1 and 2 indicate the characteristic axes of the medium. For the case in this study, the two characteristic axes correspond to h and v (horizontal and vertical).

The effective propagation constant can be expressed as (e.g., Bringi and Chandrasekar, 2001)

where k₀ is the propagation constant in the homogeneous atmosphere, n_p is the number of particles per unit volume, e_i indicates the unit polarization vector for the linear states, and $\mathit{f}(\widehat{i},\widehat{i})$ is the scattering amplitude vector in the forward-scattering configuration. The real part of the effective propagation constant induces a phase shift, while the imaginary part induces an attenuation. At L band, the attenuation due to the scattering by rain can be neglected. The expression of k_eff is defined for a number of identical particles but can be generalized to a size distribution of particles defined by N(D). Also, the $\mathit{f}\left(\widehat{i},\widehat{i}\right)$ can be expressed as the scattering amplitude matrix, S, using the Jones notation (Jones, 1941). The scattering amplitude matrix (2×2) relates the scattered field components to the incident field components in the far-field approximation. For a right-hand circularly polarized (RHCP) propagating field, as it correspond to GNSS-transmitted signals, a mean effective propagation constant can be defined by

Hence, the specific phase shift induced only by the raindrops to a circularly polarized incident wave is

in units of mm km⁻¹, where λ is the wavelength (mm), S_hh, vv is the co-polar components of the forward-scattering amplitude matrix in a linear base of polarization, N(D) is the particle size distribution (mm⁻¹ m⁻³), and D is the diameter of the particles (mm). The forward-scattering amplitude matrix is computed for each scatterer and depends on the scatterer's size, composition, orientation, and shape (see Bringi and Chandrasekar, 2001, for a detailed explanation). For this study, the T-matrix code is used in order to compute S for raindrops of all sizes between 0.1 and 8 mm of diameter (Mishchenko et al., 1996). For the particle shapes, the Beard and Chuang (1987) model is used, which relates the diameter of the each particle to the relationship between its two characteristic dimensions (i.e., its axis ratio). The complex permittivity for liquid water is obtained from Liebe et al. (1991). The N(D) is obtained at each point from the TRMM products using the same one used to retrieve rain rate from the TRMM PR reflectivity measurements, which is usually approximated with a gamma model (e.g., Kozu et al., 2009).

Using the three-dimensional collocations, we can therefore compute the phase delay that is solely due to precipitation in the following way:

• For each collocated event, we have the precipitation information interpolated into the set of RO ray trajectories. The precipitation information (for example, rain rate, water content) is used, directly or indirectly, to infer the N(D) at each point of these trajectories.

• With the N(D), we can compute the specific ΔΦ^rain along each ray using Eq. (6) and integrate this quantity along each ray path:

  $$\begin{array}{}\text{(7)}& {\mathrm{\Phi }}^{\mathrm{rain}}=\underset{L}{\int }\mathrm{\Delta }{\mathrm{\Phi }}^{\mathrm{rain}}\left(l\right)\mathrm{d}l\end{array}$$

  in units of millimeters, where L is the ray-path length in kilometers.

For each occultation event that has been 3-D-collocated with the TRMM PR, we can have the approximate vertical profiles of precipitation-induced delay along with all the currently provided information, such as the total excess phase delay, the strength of the signal, and the retrieved vertical thermodynamic products. To give the reader context, ray-path length below 15 km can be of the order of hundreds of kilometers, and therefore the amount of liquid water content that is crossed is significant. In big precipitating systems, the total integrated liquid water content along the ray paths can exceed 50 kg m⁻².

4.1 Precipitation-induced phase delay impact

In this section we want to assess the impact that the precipitation-induced phase delay has on RO retrievals. To do so we have designed a study that consists of retrieving the bending angle (phase matching method; Jensen et al., 2003) and refractivity (inverse Abel transform; Fjeldbo et al., 1971) profiles from the total excess phase delay to compare it with the retrieval results when the precipitation-induced delays are removed from the original total excess phase. Therefore, the precipitation-induced delays obtained in the previous section are removed from the actually observed phase delays, obtaining two profiles called rain (original) and rain-free (where the precipitation-induced delay has been removed).

Figure 4(a) Actual SNR (black) corresponding to the RO event C001.2008.345.00.43.G03 (UCAR ID), along with the simulated precipitation-induced phase delay (blue) as a function of time; (b, c) fractional bending angle and refractivity differences between the outputs from the retrieval using the rain-affected profiles and the rain-removed ones, as a function of the impact height (bending angle) and of the geometric height (refractivity). Black lines represent the mean of the 65 cases, while orange shading is the standard deviation.

Download

The bending angle and refractivity retrieval were attempted on both rain and rain-free excess phases from a total of 65 cases collocated with heavy precipitation events. The bending angle profiles calculated by phase matching were smoothed with 200 m windows and compared at the same impact height (corresponding impact parameter minus the collocated radius of Earth).

An example of 1 of the 65 collocated cases is shown in Fig. 4. In panel (a) of Fig. 4 we show an example of the actual occultation SNR (black) together with the precipitation-induced phase delay (blue shaded), as a function of time with respect to the start of the occultation. Note that the maximum precipitation-induced phase delay is of the order of hundreds of millimeters, while the total excess phase at the lower layers of the atmosphere is of the order of kilometers. This case corresponds to a precipitating cell in the Indian Ocean (11^∘ N, 72^∘ E), with an approximate extension of 10 000 km² and rain rate exceeding 20 mm h⁻¹. The combination of extension and intensity makes this case an interesting one, inducing an excess phase larger than 110 mm. The cases selected for this work are those with the largest precipitation-induced excess phase and are all around 100 mm. They are representative of the variety of collocated cases, combining different intensities and extensions.

In panels (b) and (c) of Fig. 4, we show the mean (black line) and standard deviation (orange shading) of the difference between the retrieval using the actual measurements and those after removing the precipitation-induced excess phase, both for the bending angle (panel b) as a function of the impact height and for refractivity (panel c) as a function of geometric height. Because of the integration nature of the inverse Abel transform, the standard deviation (orange shading in panels b and c of Fig. 4) in the retrieved refractivity is much smoother than the one in bending angle profiles. If precipitation had a systematic effect on RO phase delays, a positive bias would be observed in the rain-affected bending angle and refractivity when compared with the rain-free bending and refractivity for the same case. However, this effect is absent in Fig. 4b and c.

The results of the nonexistent mean positive bias shown in Fig. 4b and c suggest that the positive bias found in the retrieved refractivity compared to the weather analyses and reanalyses is not due to neglecting the scattering term in the refractivity. Furthermore, it can be seen how, on average, the impact of taking/not taking into account the precipitation-induced delays when assessing the retrieval increases the variability, implying that the changes in removing precipitation contribution from the signal propagation can be both positive and negative rather than only negative. Since the bending angle and refractivity retrieval process depends mostly on the vertical gradient of the excess phase, the precipitation-induced excess Doppler, which can be positive or negative, will on average lead to unbiased retrieval results. This extra-excess Doppler can be seen as the result of local horizontal inhomogeneity in the refractivity field.

Differently from temperature and pressure, the liquid and ice water content is localized in a small region (compared to the ray travel distance) and might not be contributing along the whole propagation ray path of an occultation. Furthermore, the region where liquid and ice water is present might be far from the tangent point. Yet the refractivity retrieved from an RO observation is located around the RO tangent point and is considered to have horizontal resolution of about 200 km (e.g., Kursinski et al., 1997). Even though the RO observations are integral quantities, most of the contribution from dry and wet air in the bending angle comes from near the tangent point.

In addition, the RO retrievals rely on the spherical symmetric atmosphere approximation. While it has been proven to work properly for the standard RO thermodynamic products, liquid and ice water content contributions to the excess phase cannot be well captured under such an assumption. In consequence, the effect of liquid/ice water content on RO refractivity retrieval results in small errors (e.g., less than 1 % in refractivity standard deviation), which does not introduce obvious biases in either the bending angle or the refractivity (Fig. 4). Thus, the scattering terms in Eq. (1) should not be the cause for the positive N bias in the presence of heavy precipitation.

In the previous section we have shown that the scattering term in Eq. (1) should not be the main cause for the refractivity bias. In this section we test another hypothesis: the bias comes from the problems of large-scale analyses and reanalyses in representing the thermodynamics of precipitation scenarios, especially under high specific-humidity conditions.

Figure 5Fractional refractivity difference between the observations and analyses as a function of the observed specific humidity at a height of 6 km. Panels (a), (c), (e) correspond to the no-rain cases, and panels (b), (d), (f) correspond to the rain cases. Panels (a) and (b) show the results for the comparison of observations and ERA-Interim, panels (c) and (d) show the results for ECMWF high-resolution analysis, and panels (e) and (f) show the results for the GFS.

Download

We have used the data described in Sect. 2.2 to assess the refractivity bias as a function of the RO-retrieved specific humidity for precipitating and non-precipitating scenarios. In Fig. 5 we show the results for the two analyses (ECMWF high-resolution and GFS operational analyses) and one reanalysis (ERA-Interim). We can see how the fractional refractivity difference increases with specific humidity regardless of precipitation. Therefore, the refractivity bias is more correlated with increasing specific humidity than with precipitation itself. However, high specific-humidity conditions are associated with precipitation.

We further classify the collocated COSMIC RO profiles as belonging to four different categories: no rain with low specific-humidity conditions; no rain with high specific humidity; rain with low specific humidity; and rain with high specific humidity. In this case, the threshold for low and high specific humidity is that the RO-retrieved q is lower than 0.5 g kg⁻¹ and higher than 2.7 g kg⁻¹, respectively, in the cases with no rain and that the RO-retrieved q is lower than 0.5 g kg⁻¹ and higher than 2.7 g kg⁻¹, respectively, in the cases with rain. The q and the fractional refractivity difference are evaluated at a height of 6.5 km. These thresholds are based on the lower and higher 20th and 80th percentiles of data with no rain and rain. For these four classifications, the regional fractional refractivity differences are shown in Fig. 6, for the comparison with ERA-Interim, ECMWF high-resolution analysis, and the GFS analysis. Here, the globe is divided in hexagons of a diameter of 30^∘, and only those with a minimum of 15 observations inside are shown. The typical range of observations per bin is between 800 and 10 000 for the no rain with low specific-humidity scenarios, between 15 and 600 for the no rain and high specific-humidity regimes, between 15 and 80 for the rain with low specific humidity, and between 15 and 250 for the rain and high specific-humidity regimes.

The results in Fig. 6 confirm the results anticipated in Fig. 5, i.e., that the fractional refractivity bias can be linked to high specific-humidity conditions rather than to precipitation itself. From the regional dependence of the fractional refractivity bias, some other conclusions can be extracted. The first one is that when there is no rain and the specific humidity is low, the fractional refractivity difference is very small regardless of location and the analyses in use.

The second conclusion one can extract from Fig. 6 is that when specific humidity is high, the fractional refractivity difference is positive and reaches large values (>1 %) for all the analyses in use and regardless of the presence of precipitation. In particular, high specific-humidity observations are concentrated in the tropics, so the largest positive refractivity biases are in this region, in agreement with Fig. 3.

Figure 6Regional averaged differences (color scale) between the observed and the analysis refractivity. Panels (a)–(d) correspond to ERA Interim, panels (e)–(h) correspond to ECMWF high-resolution, and panels (i)–(l) correspond to the GFS. Panels (a), (e), (i), (b), (f), and (j) correspond to the free-of-rain data, where in (a), (e), and (i) the observed specific humidity at a height of 6.5 km is lower than 0.5 g kg⁻¹ and in (b), (f), and (j) is larger than 2.7 g kg⁻¹. Panels (c), (g), (k), (d), (h), and (l) represent rain-affected data, where in (c), (g), and (k) the observed specific humidity at a height of 6.5 km is lower than 0.5 g kg⁻¹ and in (d), (h), and (l) is larger than 2.7 g kg⁻¹. The grid here corresponds to hexagons with a diameter of 30^∘. Only those with a minimum of 15 observations are shown.

Download

The third conclusion is that precipitation under low specific-humidity conditions is rarely observed in the tropics. Under these conditions, the fractional refractivity difference has a more complicated behavior and no clear positive bias is observed, but there is variability depending on the location of the observations.

Finally, in addition to the positive fractional refractivity difference bias linked to high specific-humidity conditions, we show that it is also dependent on the analysis in use. For example, the bias is larger in ERA-Interim than in ECMWF high-resolution analysis and GFS operational analysis, showing the different performance of the analyses and the reanalysis, with a smaller bias for the higher-resolution analyses. On the other hand, for no rain and low specific humidity, the performance of the different analyses and the reanalysis is similar. Overall, the fact that the bias is positive is an indication that models tend to be biased dry. This is in agreement with Hersbach et al. (2015), who noticed a dry bias in ERA-Interim which was attributed to a problem in assimilating microwave radiances affected by rain.

Data sets used for this work are available through third-party repositories. COSMIC RO data are available at the COSMIC Data Analysis and Archive Center at UCAR (https://cdaac-www.cosmic.ucar.edu/cdaac/rest/tarservice/data/cosmic, last access: 13 August 2018). GPM IMERG Final Precipitation L3 Half Hourly 0.1 degree × 0.1 degree V05 files are available at the Goddard Earth Sciences Data and Information Services Center (GES DISC; https://doi.org/10.5067/GPM/IMERG/3B-HH/05, last access: 13 August 2018). Tropical Rainfall Measuring Mission (TRMM) Combined Precipitation Radar and Microwave Imager Rainfall Profile L2 1.5 hours V7 files are available at GES DISC (https://disc.gsfc.nasa.gov/datacollection/TRMM_2B31_7.html, last access: 13 August 2018).

RP and EC led the simulations and data analysis. KNW and COA performed the inversion of the RO profiles. FJT provided the RO-precipitation co-locations. RP, EC, COA, and MTJ performed the RO refractivity and humidity analysis. All the authors contributed to the design of the study, the interpretation of the results, and the writing of the paper.

The authors declare that they have no conflict of interest.