Assessment of GNSS radio occultation refractivity under heavy precipitation

A positive bias at heights between 3 and 8 km has been observed when comparing the radio occultation retrieved refractivity with that of meteorological analyses and re-analyses, in cases where heavy precipitation is present. The effect of precipitation in RO retrievals has been investigated as a potential cause of the bias, using precipitation measurements interpolated into the actual three dimensional RO raypaths to calculate the excess phase induced by precipitation. The study consisted in comparing the retrievals when such extra delay is removed from the actual measurement and when it is not. The results 5 show how precipitation itself is not the cause of the positive bias. Instead, we show that the positive bias is linked to high specific humidity conditions regardless of precipitation. This study also shows a regional dependence of the bias. Furthermore, different analyses and re-analyses show a disagreement under high specific humidity conditions and in consequence, heavy precipitation.


Introduction
Radio Occultation (RO) technique uses opportunistic Global Navigation Satellite System (GNSS) signals to sound the atmosphere.The signal trajectory, travelling from GNSS satellites to Low Earth Orbiters (LEO), is bent due to the index of refraction gradients of the atmosphere.Such bending can be inferred using the phase derivative observable (Doppler shift) obtained by dedicated receivers in the LEOs.Under the assumption of a spherically symmetric atmosphere, the bending angle profile can be integrated to a vertical profile of the refractive index, n(h), through Abel inversion (e.g.Kursinski et al., 1997;Hajj et al., 2002).
Refractivity is defined to account for the deviations of the index of refraction from unity, and is related to geophysical parameters by (e.g.Thayer, 1974;Kursinski et al., 1997): where P is the total pressure (mbar), T is temperature (K), e is the partial water vapour pressure (mbar), n e is the electron density (m −3 ), f is the frequency (Hz), and W w,i are the liquid and ice water contents (g • m −3 ), respectively.These terms are classified as dry, wet, ionospheric and scattering terms.The dry term is dominant below 60-90 km, while the wet term becomes significant in the lower troposphere.The ionospheric term becomes dominant above 60-90 km, and its leading contribution is removed by a combination of two frequencies used by GNSS satellites (L1 = 1.575GHz; L2 = 1.228GHz) (Vorob'ev and Krasil'nikova, 1994).The scattering terms (i.e.W w,i ) are generally much smaller compared to the other refractivity terms in the lower troposphere.Therefore, they are usually neglected in the retrieval of the atmospheric variables, and when RO measurements are assimilated into the Numerical Weather Prediction (NWP) models.
A commonly used method to retrieve temperature, pressure and water vapour from RO observations is the one dimensional variational retrieval (1DVAR).It consists in obtaining the most probable atmospheric variable combining a priori atmospheric information with the observations in a statistically optimal way (Healy and Eyre, 2000).Usually, these a priori values are obtained from global meteorological analyses or reanalyses.On the other hand, bending angle and refractivity profiles are directly assimilated into NWP (e.g.Healy et al., 2005;Cucurull et al., 2007), with a high positive impact in the weather forecasts (Cardinali and Healy, 2014).
An unavoidable link exists between NWP models and RO retrieved temperature, pressure and moisture, due to the fact that RO products use a priori information from the models, and models assimilate RO observations.Yet, differences exists between their products, and its understanding is important in order to detect weaknesses and potentially improve the performance of models.
In this study we compare RO refractivity observations with the global weather analyses and re-analyses, in the presence of precipitation.These analyses have coarse spatial resolution, which has a direct impact in the treatment of heavy precipitation.At these scales, convective processes need to be parametrized.In turn, convective parametrization (CP) has been identified as one of the major source of errors in the modelling of heavy precipitation (e.g.Arakawa, 2004).RO technique offers unique potential to study the interaction between heavy precipitation and vertical thermodynamic processes within the atmosphere, since their signals can penetrate into thick clouds and their products have high vertical resolution.Recent investigations by Cardellach et al. (2014Cardellach et al. ( , 2017) ) and Padullés et al. (2016) have shown potential to retrieve vertical precipitation information adapting RO receivers to collect polarimetric observables (Pol-RO).Therefore, Pol-RO emerge as a technique that could provide relevant simultaneous information of precipitation and thermodynamics (e.g.moisture), to advance in the understanding of the processes linking vertical structure of moisture and heavy precipitation.
While such products are not yet available, in this study we investigate precipitation induced features in standard (nonpolarimetric) RO products using collocations (i.e.space and time coincidences) between the COSMIC/FORMOSAT-3 mission (Anthes et al., 2008) and the Tropical Rainfall Measurement Mission (TRMM) (Kummerow et al., 2000) and Global Precipitation Measurement (GPM) (Hou et al., 2014) missions, and we compare such features with those of analyses and re-analyses.
The refractivity from analyses and re-analyses is derived using the temperature, pressure and moisture that they provide, and Equation 1.
A clear positive bias in the RO refractivity with respect to that of some weather analyses and re-analyses is observed between 3 and 8 km height when precipitation is present in the surroundings of the observation.Previous studies have noted similar biases, for example Lin et al. (2010); Yang and Zou (2012); Zou et al. (2012); Yang and Zou (2016).These studies linked the bias with the liquid and ice water content present in the observation site, suggesting that the scattering term from Equation 1should not be neglected, but used to correct RO refractivity observations instead.However, our approach in this study is different and takes into account the 3-D structure of precipitating medium.Here, the impact of precipitation is assessed directly in the Doppler shift observable, using three dimensional collocations of precipitation structures and realistic RO ray trajectories, together with computational simulations of the effect of the scattering of the propagating signal by liquid and solid water particles.Afterwards, the causes of the observed bias are discussed with focus on the performance of the used analyses and re-analyses, especially under high specific humidity conditions.The reason to proceed this way is because solely comparing the RO observations with data from analyses and re-analyses, one could not make a clear distinction on whether the bias is due to the observation technique limitations or the weather analyses limitations.
This paper is structured in the following way.The details of the data and collocations used for this study are explained in section 2. In section 3 the bias in the comparison between RO observations and analyses and re-analyses is introduced.section 4 presents the results of the assessment of the precipitation induced delay into the RO observables.And in section 5 the specific humidity is assessed as the source of the refractivity bias.Finally, section 6 contains a discussion on the results.

RO, analyses, and precipitation data
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 vapour partial pressure at every 0.1 km of altitude, between 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 re-analysis (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, respectively.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 is the Level-2 orbital 2B31 products, that provide vertical structure information of precipitation and has a limited swath coverage.The used GPM data is the final run of the Integrated Multi-satellitE Retrievals for GPM (IMERG) products (Huffman et al., 2017), that provide 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 is 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 minutes 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 re-analyses.Therefore, the analyses and re-analysis products might or might not be associated to 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 lead to minor effects solely, since the water vapor field is spatio-temporally smoother than the cloud water content field.

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 minutes.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 evidences 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 3 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 Figure 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 into the RO signal propagation and its retrievals, as described in section 4.
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, specially above the melting layer, have larger uncertainty and must be treated carefully.

Collocations with GPM IMERG
In order to improve the statistics of collocated profiles we have performed a larger scale collocation using the GPM IMERG products (every 30 minute with spatial coverage between 60°S and 60°N) and all the COSMIC/FORMOSAT-3 RO products of 2015 and 2016.We can greatly expand the number of collocations by considering only the surface precipitation rate.For each of the COSMIC/FORMOSAT-3 RO events, the corresponding IMERG product has been identified, and the precipitation retrieval has been linked to the RO event.This results in 481,252 RO events from which the surface precipitation in its surroundings has been identified, with a time resolution of ±15min.For each event, the mean rain rate, the maximum rain rate and the number of pixels with non-zero rain rate, in a region of 2 • ×2 • , is stored along with the vertical RO profiles of refractivity, temperature, pressure, water vapour pressure, and the corresponding collocated weather analyses and re-analyses products.

Refractivity bias
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 re-analyses.In Figure 2 the bias is shown, for the comparison between the GPM IMERG collocated RO products and the three different analyses and re-analyses introduced in section 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 While the bias is clearly seen for the 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 re-analyses agree (i.e.no significant bias), as well as among themselves.
In Figure 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 it is where the bias is maximum, as seen in Figure 2. Here, the globe is divided in hexagons with a diameter of approximately 20 degrees, and the events are separated according to their R : R = 0 mm/h; 0 < R < 2.5 mm/h; and shows how the positive bias is present globally under heavy precipitation, although is larger in certain regions, and it depends on the analysis in use.Common features for all three re-analyses are, for example, the positive bias under heavy precipitation that is present in the West Pacific warm pool, the eastern part of the pacific, Indian ocean, the equatorial part of the Atlantic, and over South America and central Africa.These regions are associated to extreme precipitation features (Liu and Zipser, 2015), either to large extension precipitation events or to precipitation systems with a high deep convective cores.
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., 2006Xie et al., , 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 confirms that the RO observation geometry is not a contributing factor to the positive bias.

4 Precipitation induced delay
Once other observational known issues are discarded as plausible sources of the bias, the influence of the scattering term in Equation 1 is assessed.In order to further investigate its importance, we have simulated the contribution of the liquid and solid water directly into the excess phase.This is accomplished using 3-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 100 of precipitation in actual RO observables (see Figure 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 non-spherical 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 Equation 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 i 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 , λ 1 and λ 2 : 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, λ 1,2 = −ik h,v eff , where the k eff is the effective propagation constant of the medium, that is 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 0 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 f ( î, î) 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 f ( î, î) 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 −1 , where λ is the wavelength (mm), S hh,vv are the co-polar components of the forward scattering amplitude matrix in a linear base of polarization, N (D) is the particle size distribution (mm −1 m −3 ), 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 with 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, etc.), directly or indirectly, is used 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 Equation 6, and integrate this quantity along each ray path: in units of mm, where L is the ray-path length in km.
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 a context, the ray-paths 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 −2 .

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 in 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 the rain (original) and the rain-free (where the precipitation induced delay has been removed).
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 200m windows and compared in the same impact height (corresponding impact parameter minus the collocated radius of Earth).An example of one of the 65 collocated cases is shown in Figure 4.In panel (a) of Figure 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 and 72°E), with an approximate extension of 10,000 km 2 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 Figure 4 we show the mean (black line) and standard deviation (orange shade) of the difference between the retrieval using the actual measurements and those after removing the precipitation induced excess phase, both for bending angle (panel b) as a function of the impact height and refractivity (panel c) as a function of geometric height.Because of the integration nature of inverse Abel transform, the standard deviation (orange shade in panels (b) and (c) of Figure 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 the right panels of Figure 4.
The results of nonexistent mean positive bias shown in the right panels of Figure 4 suggest that the positive bias found in the retrieved refractivity compared to the weather analyses and re-analyses is not due to the neglect of 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 of 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 a RO observation is located around the RO tangent point, and considered to have an 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 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 both bending angle and refractivity Figure 4. Thus the scattering terms in Equation 1should not be the cause for the positive N-bias in the presence of heavy precipitation.
5 Specific humidity as a source of refractivity bias In the previous section we have shown that the scattering term in Equation 1should 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 re-analyses in representing the thermodynamics of precipitation scenarios, specially under high specific humidity conditions.
We have used the data described in subsection 2.2 to assess the refractivity bias as a function of the RO retrieved specific humidity, for precipitating and non-precipitating scenarios.In Figure 5 we show the results for the two analyses (ECMWF high resolution and GFS operational analyses) and one re-analysis (ERA-Interim).We can see how the fractional refractivity difference increase 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 to precipitation.
We further classify the collocated COSMIC RO profiles into 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 Figure 6, for the comparison with ERA-interim re-analysis, ECMWF high resolution analysis and the GFS analysis.Here the globe is divided in hexagons of a diameter of 30 deg, 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 Figure 6 confirm the results anticipated in Figure 5, i.e. 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 Figure 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 bias are in this region, in agreement with Figure 3.
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 a 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 re-analysis than in ECMWF high resolution analysis and GFS operational analysis, showing the different performance of the analyses and re-analysis, 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 re-analysis 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.

Summary and discussion
A systematic positive bias in the fractional refractivity difference has been identified when comparing RO retrieved refractivity with that of weather analysis and re-analysis when heavy precipitation was present in the surroundings of the observation.In this paper, the bias has been shown to be linked to the performance of models under high specific humidity conditions rather than with precipitation itself.
This conclusion has been reached after: (1) assessing the impact of precipitation directly into the RO observables (e.g.Doppler shift and bending angle), simulating the contribution of realistic three dimensional precipitation structures into the actual RO ray trajectories, and comparing the retrievals after such a contribution is removed; and (2) evaluating the refractivity bias between RO observations and weather analyses under different humidity and precipitation conditions.This approach is novel with respect to previous studies assessing the same bias.First, precipitation has been shown to have little impact on the positive fractional refracitivity bias between RO observations and analyses and re-analyses.Differences in bending angle and refractivity between rain and rain-removed profiles can be both positive and negative, with no clear bias on average (see right panels in Figure 4).If precipitation, through the scattering term 75 in Equation 1, had a systematic positive impact into the RO retrieved refractivity with respect to when precipitation is not present, the study performed in subsection 4.1 would have shown a positive bias as well.Therefore, precipitation itself does not explain such an impact, but the combination of thermodynamic variables associated with heavy precipitation, might.This statement does not mean that precipitation does not enhance the local refractivity where it occurs (which it is), but that such local enhancements do not necessarily lead to a bias.80 Figure 6.Regional averaged differences (colorscale) between the observed and the analysis refractivity.The first row corresponds to ERA interim, the middle row corresponds to ECMWF high resolution and the bottom one to the GFS.The two left rows corresponds to the free of rain data, where in the first row the observed specific humidity at a height of 6.5 km is lower than 0.5 g/kg, and in the second row it is larger than 2.7 g/kg.The two right columns represent rain affected data, where in the third column the observed specific humidity at a height of 6.5 km is lower than 0.5 g/kg, and in the last row it is larger than 2.7 g/kg.The grid here corresponds to hexagons with a diameter of 30 deg.
Only those with a minimum of 15 observations inside them are shown.
The fractional refractivity bias between RO observations and analyses has been linked to high specific humidity conditions.
The bias appears both in rain and no-rain conditions, and it depends on the analyses and on the geographic region.The spatial resolution of the analyses and re-analyses may also be a factor, since ERA-Interim shows a larger bias than ECMWF high resolution analysis, although they should be based on the same physical model.However, both ECMWF high resolution and GFS operational analyses still exhibit a significant N-bias with increasing humidity.The fact that the bias is not seen in Figure 2 for the no-rain cases is because most of the no-rain cases have low specific humidity conditions, and they weight much more for the mean value of the fractional refractivity difference.On the other hand, the rain cases have a larger contribution in the high specific humidity region (see Figure 5), contributed mostly by tropical precipitation.This is also seen in the right panels in Figure 6, where precipitation with very high specific humidity conditions is mostly observed around the equator.
The bias in fractional refractivity between observations and analyses implies that the retrieved temperature and moisture 10 will be also biased with respect to models.The positive refractivity bias is associated with a combination of colder retrieved temperature with respect to analyses, and a higher retrieved specific humidity than the one in the analyses.This is consistent with Vergados et al. (2015), who showed that ERA-interim is systematically drier than RO in the tropics.Also, the fact that a difference exists between the different analyses used for this study, and that RO thermodynamic retrievals depend on the analyses or re-analysis in use, imply that a difference between the retrievals obtained by different processing centers will exist under such conditions if they use different analyses or re-analyses.
These results stress the need for a better thermodynamic characterization of high specific humidity scenarios, likely to be associated to heavy precipitation.The heights at which the bias is maximum is also consistent with the findings of, e.g.Holloway and Neelin (2009), who argue that heavy precipitation is controlled by the free tropospheric water vapor, and this dependence is not well captured in large scale models.These models are known to have issues with the parameterization of convective processes, hence further investigation in this direction is required.This is the aim of polarimetric radio occultations, which will provide joint products of temperature, pressure, moisture and an indication of the amount of precipitation (mostly sensitive to the heaviest) at each vertical level (Cardellach et al., 2017) with the objective to advance in the understanding of heavy precipitation events, closely linked with high specific humidity conditions.Alternatively, further investigations are being conducted with the aim to make the RO retrievals less dependent on models, which would improve the retrievals itself and provide more independent information of such scenarios.
Competing interests.The authors declare that they have no conflict of interest.

Figure 1 .
Figure 1.Three dimensional collocation of a 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.51450TRMM PR product.Here the precipitation structure is shown in a 3 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 interpolated precipitation information (rain rate) into the RO plane is shown in the 2-Dimensional projection in the latitude-height plane.

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

Figure 3 .
Figure 3. Regional averaged fractional difference between the RO observed refractivity and that from (top row) Era interim re-analysis; (middle row) ECMWF high resolution analysis; and (bottom row) NCEP GFS operational analysis; for a height of 6 km.The compared profiles are classified between no-rain profiles (left column; R = 0 mm/h), low and moderate precipitation (middle column; 0 < R < 2.5 mm/h) and heavy rain profiles (right column; R > 2.5 mm/h).The grid corresponds to hexagons with a diameter of about 20 deg.Only those with a minimum of 25 observations inside them are shown.

Figure 4 .
Figure 4. (left) 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; (right) 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 shade is the standard deviation.

Figure 5 .
Figure 5. Fractional refractivity difference between the observations and analyses as a function of the observed specific humidity at a height of 6 km.The left column corresponds to the no-rain cases and the right column corresponds to the rain cases.The top row shows the results for the comparison of observations and ERA-interim, the middle row show the results for ECMWF high resolution analysis and the bottom row shows the results for the GFS.

Acknowledgements.
Ramon Padulles and K. N. Wang's research was supported by an appointment to the NASA Postdoctoral Program at the Jet Propulsion Laboratory, administered by Universities Space Research Association under contract with NASA.The work conducted at ICE-CSIC/IEEC was supported by the Spanish grant ESP2015-70014-C2-2-R. Part of Cardellach's contribution has been supported by the Radio Occultation Meteorology Satellite Application Facility (ROM SAF) which is a decentralised operational RO processing centre under EUMETSAT.The JPL co-authors acknowledge support from the NASA US Participating Investigator (USPI) program.TRMM and GPM data were obtained courtesy of the NASA Precipitation Processing System (PPS).Part of this work was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration.©2018.All rights reserved.