Articles | Volume 17, issue 3
Research article
14 Feb 2024
Research article |  | 14 Feb 2024

Forward operator for polarimetric radio occultation measurements

Daisuke Hotta, Katrin Lonitz, and Sean Healy

Global Navigation Satellite System (GNSS) polarimetric radio occultation (PRO) observations sense the presence of hydrometeor particles along the ray path by measuring the difference of excess phases in horizontally and vertically polarised carrier waves. As a first step towards using these observations in data assimilation and model diagnostics, a forward operator for the GNSS-PRO observable ΦDP (polarimetric differential phase shift) has been implemented by extending the existing two-dimensional forward operator for radio occultation bending-angle observations. Evaluation of heavy-precipitation cases showed that the implemented forward operator can simulate the observed ΦDP in synoptic-scale atmospheric river (AR) cases very accurately. For tropical cyclone cases it is more challenging to produce reasonable ΦDP simulations, due to the high sensitivity of ΦDP with respect to displacement of the position of the tropical cyclones. It was also found that snow is the dominant contributor to the simulated ΦDP and that the ability to compute the ray paths in two dimensions is essential to accurately simulate ΦDP.

1 Introduction

The speed of light is slowed down when radio waves pass through the air, and this “retardation” is larger when the air is heavier and more humid. Because of this, as radio waves travel through stratified atmosphere from an emitter on a Global Navigation Satellite System (GNSS) satellite to a receiver on board a low-Earth-orbit (LEO) satellite, they undergo bending (or refraction) to minimise the travel time along the ray. In radio occultation (GNSS-RO) observations, this bending is retrieved from continuous measurement of the phase of the radio waves. As the refraction depends on the density of dry air and the amount of water vapour, measurements of bending can inform us about the thermodynamic properties of the atmosphere along the ray paths. GNSS-RO observations are routinely assimilated at most major numerical weather prediction (NWP) centres and are recognised as an indispensable component of modern NWP systems (e.g. Bonavita2014).

The carrier waves employed in GNSS are circularly polarised to minimise the impact of receivers' antenna alignment on the accuracy and stability of positioning. Because the carrier waves are polarised, it should be possible, in principle, to obtain information on properties of hydrometeors along the rays, just like polarimetric phase-shift measurement from dual-polarised weather radars (e.g. Kumjian2013). Polarimetric radar observations exploit the fact that, when polarised radio waves travel through a medium containing non-spherical objects like hydrometeor particles, the phase is delayed differently in the horizontally and vertically polarised waves due to the directionally differential cross-section of such objects, so that polarimetric differential phase delay contains information about the presence of those particles along the ray. For example, if a large difference between horizontal and vertical phase shifts is observed, that indicates the presence of more and/or larger hydrometeors (and, thus, heavier precipitation) along the ray since large rain droplets are usually oblate. This same principle should also hold for GNSS carrier waves to allow for inference of properties of hydrometeors from GNSS-RO observations if the horizontally and vertically polarised components of the radio waves can be processed separately, provided that depolarisation during the wave propagation through the ionosphere has been accounted for. An additional benefit of using GNSS carrier waves that are usually in the L-band (India's regional navigation system NavIC also uses S-band in addition to L-band), in comparison to X-band to S-band as in most weather radars, is that the relatively low frequencies in the L-band range may make the measurement more sensitive to larger hydrometeors while being insensitive to smaller particles like aerosols and non-precipitating cloud droplets.

Such polarimetric measurement of GNSS-RO observations, which we shall call PRO hereafter, has not been explored until recently but was enabled by the sensor deployed for Radio Occultation and Heavy Precipitation (ROHP) mission on board the Spanish PAZ satellite (Cardellach et al.2019). The PAZ satellite was successfully launched in May 2018 and has already been producing PRO measurements for more than 4 years (as of April 2023), with the observed cases including many heavy-precipitation events.

PRO observation complements the standard RO measurements of a bending-angle profile, which is sensitive to the thermodynamic variables (temperature, pressure and humidity) of the atmosphere, with additional information on the vertical profile of heavy precipitation. Such additional information is provided by measurement of the differential phase shift at each vertical level of the ascending/descending rays, which in turn is enabled by measuring the phase delay at two orthogonal (horizontal and vertical) polarisations.

The promise of PRO measurements is already established by recent studies. Cardellach et al. (2019) confirmed, with real data, that PRO measurements of differential phase shift exhibit stronger signals in the presence of heavier precipitation. Turk et al. (2021) and Padullés et al. (2021) simulated PRO measurements using hydrometeor retrieval products from collocated Global Precipitation Measurement (GPM) measurements, and their results suggested that PRO measurements do contain useful information about vertical structure of deep convective clouds. Murphy et al. (2019) simulated PRO measurements of an airborne instrument using output data from an atmospheric mesoscale model and showed that such measurements can provide useful guidance on validation of microphysics representation in the model.

An important benefit that is unique to PRO observations is that because the regular RO (or bending) measurement and the newly available PRO measurement are carried out simultaneously, profiles of thermodynamic and cloud-related properties can be observed at the same time. Hence, if an accurate observation operator is available that can simulate PRO measurements from state variables of a NWP model, PRO observations can potentially be of great diagnostic value to modelling of physical processes.

As a first step towards PRO assimilation and model validation with PRO measurements, we develop an offline forward operator of PRO measurement for the Integrated Forecasting System (IFS) model of the European Centre for Medium-Range Weather Forecasts (ECMWF).

The paper is structured as follows. Section 2 describes the specification and the main components of the forward operator, clearly presenting the assumptions we made and their potential limitations. Section 3 describes the data and the model used in this study, along with the cases examined. Section 4 presents the results including those from several sensitivity experiments, followed by discussion and conclusions in Sect. 5.

2 Description of the forward operator

The main observable of GNSS-PRO is the differential phase shift ΦDP=ΦH-ΦV, which is the additional excess of the phase delay of the horizontal wave ΦH in comparison to that of the vertical wave ΦV. This can be computed as the integration along the ray path, s, of the specific differential phase shift KDP:

(1) Φ DP = GNSS LEO K DP ( s ) d s ,

where GNSS and LEO symbolically represent, respectively, the position of the transmitter and receiver of GNSS radio signals. KDP indicates how much the phase of the horizontally polarised wave is delayed in comparison to that of the vertically polarised wave as they travel a unit distance. Since RO measurements represent path-integrated quantities, a positive value of ΦDP is an indication of the presence of horizontally oriented hydrometeors somewhere along the ray path.

The main components to computing Eq. (1) are determination of the ray path and estimation of KDP from hydrometeors represented in the model. In IFS, hydrometeors from parametrised convection are represented as their vertical mass fluxes, and thus we need to convert convective mass fluxes to mixing ratios in order to relate them to KDP (see Sect. 2.3).

2.1 Ray tracing

We develop the PRO forward operator by extending the operational two-dimensional (2D) forward operator for RO bending angles (Healy et al.2007). The bending angle, α, computed by this 2D forward operator, can be symbolically written as

(2) α = GNSS LEO d α d s d s ,

with the ray path being identical to the one in Eq. (1). To compute the PRO observable ΦDP, we exploit the analogy between Eq. (1) and Eq. (2) and use the existing code from the operational 2D bending-angle operator to compute the ray path and to integrate the integrands.

The ray tracing follows “Approach 2” of Healy et al. (2007), which is described in detail in their Sect. 3.2. For each ΦDP measurement, the latitude, longitude, height of the tangent point and azimuth angle of the ray are taken from the observed data. The forward operator then makes a 2D slice of the 3D model field in the direction of the azimuth angle centred around the tangent point. The slice comprises vertical columns, each at the model's native vertical full level, of equally spaced locations along the occultation plane. Angular distance of two adjacent columns is set so that their horizontal distance is approximately equal to the typical horizontal grid spacing of the input grid, and the number of vertical columns in the slice is chosen so that the slice horizontally spans  1200 km. When the input model field is on a 0.25× 0.25 regular lat–long grid, for example, the forward operator first computes the latitudes and longitudes of 31 points equally spaced with Δθ=40/6371 rad (corresponding to 40 km physical distance at the Earth's surface), along the great circle centred around the tangent point's horizontal position and with the azimuth angle specified by the observed data, and horizontally interpolates the model fields to these horizontal points to construct the 2D slice. Similarly, when the input model field is on a 0.125× 0.125 regular lat–long grid, the angular spacing is set as Δθ=20/6371 rad, and the slice will contain 61 columns. In this study, for ease of implementation, the horizontal interpolation is done by nearest-grid search.

Once the 2D slice is set up, the forward operator computes the ray path, starting from the tangent point by integrating the ray equations in both directions (towards the receiver on board the LEO satellite and towards the transmitter on board the GNSS satellite). The numerical integration of the ray equations is based on the second-order Runge–Kutta method (Healy et al.2007, used the fourth-order Runge–Kutta method, but the operator was later simplified to adopt the second-order method). The KDP contributions from each type of hydrometeors (see next section) are integrated as each section of the ray is traced, and the accumulated KDP from each section is finally summed up to obtain the total ΦDP. The detail of ray tracing is described in Healy et al. (2007). Vertical interpolation for KDP is performed column-wise. While KDP can be negative in nature in situations, for example, where particles in prolate orientation dominate over oblate ones, in our modelled formulation of KDP (see the next subsection, in particular Eq. 3), it always assumes non-negative values. It is hence desirable to ensure monotonicity in the interpolation, and for this reason we employ a simple linear interpolation in the vertical. Unlike refractivity, exponential decay with height is not assumed for KDP.

During an occultation event, the horizontal position of the tangent point drifts as the ray ascends or descends. While the operational 2D operator for bending angles has accounted for such “tangent-point drift” (see Sect. 4.3) since 2011, Healy et al. (2007) found that assimilating regular RO observations of bending angles can be beneficial, even when the effect of tangent-point drift is neglected. We found, from PAZ data, that in a single occultation event, the tangent point typically drifts  100 km. This can be tolerated for bending angles, but it is not clear if the same can be said for KDP because it is sensitive to hydrometeors, and their horizontal variability is much larger than that of thermodynamic fields. We examine this aspect in Sect. 4.

We remark that the ray tracing implemented in our 2D operator relies on the position of tangent points and the impact parameter provided from the RO data processing centres, but the tangent-point position can only be determined after ray tracing has been done. In this sense, our ray tracing is dependent on externally performed estimation of the ray path. The accuracy of this tangent-point estimation may impact the performance of our ray tracing.

2.2 Relating hydrometeor water content to KDP

Specific differential phase shift, KDP, is induced by the difference in scattering properties of hydrometeor particles for the horizontally and vertically polarised waves. A “first-principle” scattering calculation of KDP from NWP variables would require assumptions about details of cloud microphysics and precipitation that are not currently represented in the forecast model. In IFS, the only prognostic variables related to hydrometeors are their water content (or vertical mass flux; see next subsection), so KDP needs to be estimated from the water content variables from each type of hydrometeors.

In IFS, hydrometeors represented with the resolved (or large-scale) microphysics scheme can be categorised into the following four different kinds: non-precipitating liquid water, non-precipitating ice water, precipitating liquid water (or rain) and precipitating ice water (or snow). We denote the water content of these categories, respectively, by LWC, IWC, RWC and SWC. In addition to these resolved-scale variables, the deep-convection parametrisation scheme also represents precipitating rain and snow separately. We denote the rainwater content and snow water content that are attributable to the deep-convection scheme, respectively, by RWCconv and SWCconv.

To determine KDP from the water content of each category of hydrometeor, here we adopt a simple linear relation between the water content and KDP contribution. For ice (i.e., IWC, SWC and SWCconv), following Padullés et al. (2021), we adopt the following water content-to-KDP formula (which, in turn, is taken from Bringi and Chandrasekar2001):

(3) K DP ( WC ) = 1 2 C ρ × WC × 1 - ar ,

where WC means any one of IWC, SWC or SWCconv (in units of g m−3); KDP(WC) (in units of mm km−1) denotes the KDP contribution from the hydrometeor whose water content is denoted WC; ρ is the particle density in units of g cm−3; and ar (non-dimensional) is the assumed axis ratio of the particles. C is a proportionality constant which is derived from the theory of scattering by spheroid objects which is set as C=1.6 (g cm−3)−2 in this study. For ar and ρ, we assume them to be constant and arbitrarily chose their values as ar =0.5 and ρ=0.2 g cm−3.

The main approximating assumptions behind the formula relating KDP to water content, Eq. (3), are (1) particles are spheroid, (2) particles are small enough in comparison to the wavelength so that Rayleigh scattering dominates, (3) ar and ρ are independent of the particle size (see the derivation in Chap. 7 of Bringi and Chandrasekar2001), and (4) the natural variability of the particle size distribution (PSD) can be ignored (i.e., can be regarded fixed regardless of the atmospheric conditions). Compared with findings from the weather radar community, assumptions (1) and (2) may seem too crude for ice particles, but as we show later, the simulated ΦDP values are quite consistent with observations. We revisit this point in Sect. 5.1. Assumption (3) is admittedly difficult to justify, and we consider this to be an important limitation of our approach (see Sect. 5.3).

It is relevant to note that KDP, at L-band, in which the Rayleigh limit is appropriate, is nearly proportional to the third moment of the PSD (Bringi and Chandrasekar2001, Sect. 7.1), which, in turn, is proportional to the water mass, unlike reflectivity, which is related to the sixth moment of the PSD and thus is sensitive to the contributions from larger but fewer particles. Therefore, in comparison to radar reflectivity, KDP is less sensitive to the precise shape of the PSD. In fact, as shown in Fig. 11 in Turk et al. (2021), KDP depends virtually linearly on the mass content. This figure also suggests that, if the particles are flat enough (i.e., if ar is small enough), KDP for ice particles can be larger than for rain particles. This, along with the fact that the rays of GNSS-(P)RO tend to travel longer in the upper atmosphere above freezing levels than below, may explain the greater contributions to ΦDP from ice and snow particles than from liquid and rain particles (see discussion in Sect. 5.1).

We also note here that the orientation of the ice/snow particles are situation dependent, and hence the axis ratio (ar) would better be allowed to vary. For example, Padullés et al. (2021), when simulating ΦDP using the ice water content retrieved from Global Precipitation Mission (GPM) Microwave Imager (GMI) observations, allowed ar to vertically vary from 1 near the cloud top to 0 below 10 C level and further modified it depending on the polarisation difference (PD) of the observed brightness temperature measured by GMI. Similarly, the particle density ρ should be different depending on the particle size and shape. In our study, however, for simplicity, and due to lack of knowledge about particle orientation and details of particle shape, these effects are not accounted for.

The formula for liquid water (LWC, RWC and RWCconv) should be different from Eq. (3), but here, we use Eq. (3) for LWC, RWC and RWCconv as well. We are aware of the imperfection of this approach and aware that at least the coefficients should be refined. As we discuss in Sect. 4, however, the ΦDP contributions from liquid water are found not to be so dominant. Hence, small changes in tuning these parameters for liquid water will not change the main signals.

2.3 Converting mass flux to water content

In IFS, the amount of hydrometeor is represented (and archived) differently for the resolved (large-scale) scheme and for the parametrised convection scheme. In the resolved microphysics, LWC, IWC, RWC and SWC are directly available as specific water content (in units of kg kg−1), which can be readily converted to mass per volume (in units of g m−3). In the convective scheme, however, the amounts of hydrometeor (RWCconv and SWCconv) are represented only as their vertical mass fluxes (in units of kg m−2 s−1). To convert them into water content mass per volume, some additional assumptions have to be made. Here we follow the approach adopted in RTTOV-SCATT, described in Appendix B of Geer et al. (2007). In this approach, we assume that the particle density ρ is constant, the fall speed of a particle of diameter D is proportional to Dβ for some β and the particle size distribution follows an exponential decay. With these assumptions, one can show (Alan Geer, personal communication, 2022), with calculus involving the Gamma function, that the water content WC, in units of g m−3, can be derived from the vertical mass flux FL as FL/(aρ)1/b, with the parameters a, b and ρ given in Table 8 of Geer et al. (2007). This approach, adopted in RTTOV-SCATT, should be also applicable to models other than the IFS as long as vertical mass fluxes of water substances are available in their forecast or analysis product.

3 Model and data

3.1 Forecast model fields

The forecast data used to simulate ΦDP are produced by running the version Cy47R3 of the IFS model at Tco1279 horizontal resolution (approximately 9 km of grid spacing) with 137 vertical model levels initialised from the operational analysis fields, available twice daily at 00:00 and 12:00 UTC, that are valid at the time closest to the start time of the occultation event simulated. Since the IFS Cy47R3 was operational from 12 October 2021 until 26 June 2023, the forecast fields used in this study are essentially identical to the operational deterministic forecasts (apart from technical differences such as use of different versions of the compiler and some of the software stack) for the event that falls inside this period (tropical cyclone Kompasu, to be precise). For the other events that occurred before this period (see Tables 1 and 2), the forecast fields used in this study are generated from a newer version of the IFS than the then-operational suite. The model fields defined on the model's native Tco1279 Octahedral Gaussian grid are interpolated to a regular 0.125× 0.125 lat–long grid by ECMWF's MARS system before being ingested into the forward operator.

Table 1List of the examined atmospheric river (AR) cases (provided by Ramon Padullés). RO ID is an identification code given following the UCAR convention for each occultation event; time (UTC) is the time at which the occultation begins; and latitude and longitude are those of the “occultation point”, where the excess phase exceeds 500 m for the first time during the occultation event.

Download Print Version | Download XLSX

Table 2As in Table 1 but for tropical cyclone (TC) cases, with an additional column for the TC names.

Download Print Version | Download XLSX

The model fields are available at an hourly interval in time. The forecast fields at two adjacent time steps are ingested to the forward operator, which are then linearly interpolated to the start time of the occultation event. The mass flux variables are archived as time accumulation from the beginning of the forecast integration; for these variables, we take the difference of forecast fields at two adjacent time steps that contain the time of occultation event and divide the difference by 3600 s as if it is an instantaneous value, assuming that their values are constant over the 1 h time interval.

3.2 Observations

We use the version V06 of PAZ Level-1B data processed and calibrated at the Spanish Institute of Space Sciences ICE-CSIC/IEEC (Cardellach et al.2019). The data are available for download from the ROHP-PAZ website (accessible at, last access: 2 February 2024) upon registration. In this dataset, data for each occultation event are provided in a separate netCDF file containing time series of the calibrated ΦDP and the tangent-point height along with other relevant data and metadata that include the azimuth angle at the “occultation point” (i.e., the point at which the excess phase first becomes greater than 500 m). To simulate ΦDP with our forward operator we additionally need information on the latitude and longitude of the tangent point. These data are not included in the PAZ dataset, so we extracted them from the corresponding UCAR-processed level-2 data.

The ΦDP data measured from PAZ are known to undergo height-dependent systematic errors, and the PAZ dataset provides ΦDP data calibrated with two different approaches, one based on antenna pattern and the other based on linear regression. In this study we used the ΦDP profile with antenna pattern calibration (denoted with the variable name dphase_cal_ant in the netCDF files) which was shown in Padullés et al. (2020) to be more accurate than the linear-regression-based calibration.

3.3 Examined cases

Using the forecast model and the observed data described above, ΦDP profiles are simulated from the model and compared with the observations for five atmospheric river (AR) cases and six tropical cyclone (TC) cases, which all exhibit large ΦDP signals in the PAZ observations and are accompanied with heavy precipitation. The AR cases were selected by Michael Murphy, Jennifer Haase and Ramon Padullés, while the TC cases were selected by Ramon Padullés. These cases are proposed for a multi-centre model intercomparison project of ΦDP simulated from NWP output fields. A summary of the cases is given in Tables 1 and 2.

4 Results

4.1 Baseline results

We first examine the overall agreement of the simulated and observed ΦDP profiles and the relative contributions from different categories of hydrometeors. The simulated and observed ΦDP for AR and TC cases, plotted as vertical profiles against the tangent height, are shown in Fig. 1. To use the best possible simulation as the baseline, in these plots the tangent-point drift is fully accounted for (i.e., we use the correct tangent-point position at each tangent height). To indicate the heights below which one can expect ΦDP contributions from rain and liquid water, in each panel, the freezing level that is taken from the PAZ metadata is shown with a thin horizontal black line.

Figure 1Comparison of the observed (blue) and simulated total (purple) ΦDP profiles for the AR cases (top two rows) and the TC cases (bottom two rows). ΦDP contributions from resolved-scale non-precipitating ice (IWC) and liquid (LWC), resolved-scale precipitating rain (RWC) and snow (SWC), and convective scheme rain (RWCconv) and snow (SWCconv) are also shown with different colours depicted in the legend. Thin horizontal black lines in each panel show the freezing level at the tangent point inferred from the retrieved temperature profile obtained from the PAZ RO data. Freezing levels roughly coincide with the levels at which RWC / LWC and RWCconv start to show non-zero contributions.


For the AR cases (Fig. 1, top two rows), simulated ΦDP profiles fit very well to the observed profiles albeit with some overestimation in two of the cases that occurred during January 2021. Considering all the uncertain assumptions that are made in linking hydrometeor water content to KDP (Sect. 2.2), this level of agreement is quite surprising. From Fig. 1 we can also observe that simulated ΦDP is dominated by contributions from resolved-scale snow (SWC; yellow solid lines). Because of the uncertainty in how we estimate KDP from hydrometeor water content (Sect. 2.2), we cannot assert that SWC contribution dominates solely by judging from their dominance in magnitude. However, the shape of the profile of SWC contribution closely resembles that of the observed ΦDP profile for any of the five cases, which should mean that ΦDP is predominantly determined by resolved-scale snow. In contrast, for the TC cases (Fig. 1, bottom two rows), the simulated and observed ΦDP values do not agree well, with the former significantly overestimating the latter. We investigate why the simulation results are so drastically different in AR and TC cases in the rest of this section.

4.2 Sensitivity to model-field displacement

We have seen that the simulated ΦDP agrees well for AR cases but not so much for TC cases. One factor that may explain this sharp contrast is the spatial scale of the phenomena: the horizontal scales of TCs are typically much smaller than those of synoptic disturbances like AR, so that even a small positional error in the model fields may have significant impact on simulated ΦDP in the TC cases, while these may be tolerated in the AR cases. Such a “displacement effect” was suggested earlier in Murphy et al. (2019).

Estimates of the uncertainty in simulated ΦDP that can be attributed to displacement of model fields are plotted in Fig. 2. Here, we took a simple approach and estimated the uncertainty by shifting the latitude and longitude of the tangent points by ±0.1 (which correspond to shifts in position by  10 km, namely, roughly by one grid point) or 0, resulting in nine profiles computed for each event in total. We assume that the spread among such profiles would represent the range of uncertainty that we would have if the forecast had displacement error on the order of one grid point. The IFS version Cy47R3, which we used to generate forecast fields assessed in this study, has tropical cyclone position error of O(∼30) km on average in the 0–12 h forecast range (Forbes et al.2021, their Fig. 6a). Thus, the displacement of O(10) km that we examine here should yield reasonable estimates of ΦDP's sensitivity to the background field uncertainties of particularly successful TC forecasts.

Figure 2 shows that the simulated ΦDP values are insensitive to the forecast displacement in the AR cases but are more sensitive in the TC cases. This high sensitivity to the displacement error can explain the poorer fit between the simulated and observed ΦDP, albeit not the systematic overestimation in the TC cases.

Figure 2Uncertainty of simulated total ΦDP (purple) and their contributions from resolved-scale snow (SWC; teal) and parametrised convective snow (SWCconv; pale blue), computed for each of the AR cases (top two rows) and each of the TC cases (bottom two rows). Here, the uncertainty is estimated by shifting or not shifting the latitude and longitude of the tangent points by ±0.1. This results in computing 3×3=9 profiles in total, and the range between the minimum and maximum of these nine profiles is shown by the shading. The unperturbed profiles are shown by solid lines.


4.3 Impact of tangent-point drift

The results shown so far have been computed by fully taking into account the effect of tangent-point drift (i.e., by changing the horizontal position of the tangent point for each tangent-point height). In practice, this can be prohibitively expensive because, each time the tangent-point position changes, the 2D slice has to be re-generated. It is thus desirable to reduce the frequency of tangent-point position updates to minimise the number of 2D slices to be created as long as the accuracy is not too degraded.

Here, in addition to the “full drift” approach shown above, we explored two more approaches: “no drift”, in which the drift of tangent point is not accounted for, and “11-batch”, in which 11 neighbouring tangent-point heights are grouped into a batch which shares the same 2D slice. In the 11-batch approach, rays in each batch are assumed to share the same tangent-point horizontal position, which is the sixth point of the 11 tangent points within the batch. The ECMWF's operational system uses the 11-batch approach to assimilate bending angles.

Profiles of ΦDP simulated with the three approaches handling the tangent-point drift are shown in Fig. 3. As we can expect from the insensitivity of simulated ΦDP with respect to the horizontal displacement in AR cases, the different approaches yield ΦDP values that are equally consistent with the observations. In contrast, in TC cases, simulated ΦDP profiles are highly dependent on how the tangent-point drift is handled. Contrary to a naive expectation, however, the “full-drift” approach, which is the most expensive but should be the most accurate, does not necessarily result in simulations most consistent with the observations. This is likely because the overall error is dominated by the errors that result from displacement error, and thus the impact from refining tangent-point drift is obscured.

Figure 3Impact of different approaches to account for tangent-point drifts. Profiles of ΦDP simulated with the three different ways to handle the tangent-point drift (see text for detail) are plotted for each of the AR cases (top two rows) and each of the TC cases (bottom two rows), along with the observed ΦDP.


4.4 Limitations of 1D operator

Most global NWP centres rely on a one-dimensional (1D) forward operator in simulating and assimilating RO bending angles. A 1D forward operator computes bending-angle observables by only using the atmospheric profile at the tangent point, assuming that atmospheric thermodynamic fields within the occultation slice can be regarded horizontally uniform, and hence an identical atmospheric profile can be used for all the columns in the 2D slice. It is thus interesting to see how a 1D operator would perform in simulating the PRO observable ΦDP. Such an assessment will also clarify how essential 2D ray tracing is in simulating realistic ΦDP.

Our 2D operator can easily operate in “1D mode” to emulate a 1D operator. To do this, we just set the derivative of the horizontal ray position to zero when integrating the ray equation (which is equivalent to assuming zero horizontal gradient of refractivity within the 2D slice). For simplicity, the tangent-point drift is ignored in our 1D computation.

The results of 1D computation are summarised in Fig. 4. Unlike the 2D results (Fig. 1), simulated ΦDP values are highly inconsistent with the observed ΦDP, even for the AR cases. The extreme case is Hurricane Larry (centre panel in the bottom-most row), in which the simulated ΦDP is almost zero except very near to the surface.

Figure 4As in Fig.  1 but with the profiles simulated with 1D computation.


To understand why, cross-sections of resolved-scale snow water content (SWC), which is the dominant contributor to KDP, are informative (Fig. 5). In any of the cases, the distribution of SWC is far from being horizontally uniform, violating the assumption of the 1D computation. In the case of Hurricane Larry, the tangent point just happens to be inside the eye where there is no cloud and hydrometeors at all, so that 1D computation using only the profile at the tangent point completely misses out the hydrometeors in its vicinity.

Figure 5Resolved-scale snow water content (in g m−3) along the rays.


The poor fit of the simulated and observed ΦDP highlights the importance of the capacity to perform 2D computation in accurately simulating ΦDP observable. The importance of 2D computation is not limited to PRO observations but also to regular RO observations of bending angles because the assumption of horizontal homogeneity can be questionable inside and in the vicinity of weather disturbances associated with heavy precipitation where the variability of thermodynamic parameters, especially humidity fields, tend to be large. This point will become increasingly important as the model's horizontal resolution gets even higher to better resolve precipitation systems.

5 Summary and discussion

A forward operator for the GNSS-PRO observable ΦDP has been implemented by extending the existing 2D forward operator for GNSS-RO bending-angle observations assuming a linear relation between hydrometeor water content and KDP. The implemented forward operator has been tested with five atmospheric river (AR) cases and six tropical cyclone (TC) cases which all accompanied heavy precipitation and strong ΦDP signals in the actual observations. Despite all the simplifications, the implemented forward operator is found to be able to produce simulated ΦDP profiles that are remarkably consistent with the corresponding observed profiles in most of the AR cases. Quantitatively speaking, however, the simulated ΦDP profiles still tend to somewhat overestimate the observed profiles, and we speculate that such disagreement is likely attributable to imperfect assumptions on hydrometeor properties like the oblateness and the tilting angle (represented with the parameter ar), as we acknowledged in Sect. 2.2, and the uncertainty of hydrometeor variables in the forecast fields. A more rigorous validation with more samples from different cases would be in order to draw a quantitative conclusion.

In contrast to the success with AR cases, we found TC cases to be much more challenging, with the simulated ΦDP values systematically overestimating the observed ΦDP. Several additional ΦDP simulations have been conducted with varying configurations of the implemented operator to understand its capacity and limitation. In this section we highlight and discuss the following main findings from this study.

5.1 Why does snow dominate?

From the results in Fig. 1 it was found that simulated ΦDP is dominated by the contribution from snow particles. This does not immediately mean that ΦDP (or KDP) is dominated by snow in reality because our modelled KDP–water content relation involves multiple uncertain assumptions, and the hydrometeor representation in the IFS model is also subject to forecast uncertainties. Nevertheless, the remarkable agreement between the observed and simulated ΦDP in AR cases suggests that the dominance of snow contributions to ΦDP is likely realistic.

This finding is in line with previous findings by Turk et al. (2021) and Padullés et al. (2021), who simulated ΦDP profiles observed by PAZ using liquid and solid water content retrieval from collocated cloud-sensitive measurements from other satellites. They found that the liquid phase hydrometeor alone cannot explain the observed ΦDP values, especially at high tangent-point heights above the freezing level. It appears, nevertheless, that liquid hydrometeor contributions to ΦDP being negligible in comparison to snow has not been reported before.

In the weather radar community, it is widely accepted that KDP per mass of snow particles is an order of magnitude smaller than that of liquid rain (Doviak and Zrnić1993); therefore snow is largely undetectable from KDP measurements (e.g. Kumjian2013). Because of this, most forward operators that have been developed to simulate or assimilate radar KDP observations only consider warm rain conditions (e.g. Li and Mecikalski2013; Yokota et al.2016; Kawabata et al.2018a, b). From this perspective, it is surprising that PRO measurement likely predominantly senses the presence of snow rather than rain.

The geometry of GNSS-PRO may be a factor contributing to the high sensitivity of ΦDP with respect to snow: the rays travel long distances at high altitudes except in the very vicinity of the tangent points as they go through the atmosphere from an emitter aboard a GNSS satellite toward the receiver aboard the PAZ satellite, so that, even if KDP values for snow/ice particles are relatively small compared to those for rain droplets (in fact, KDP for snowflakes that are flat enough can be larger than for rain droplets; see Fig. 11 of Turk et al.2021), the integration of such smaller KDP over a much longer path than for rain droplets that are confined to lower altitudes below the freezing level can dominate in their contributions to the total ΦDP.

Another factor that may explain this apparent contradiction would be the difference in carrier wave frequencies of GNSS and weather radars. In GNSS, the L-band is chosen as the carrier frequency since radio waves at these frequencies are less prone to attenuation by hydrometeors, thus allowing for signals to propagate in all sky conditions. In the L-band, the frequency is  1.4 GHz, corresponding to the wavelength of λ 20 cm. In contrast, in weather radars, the carrier frequency of the C- to X-band is typically chosen to maximise backscatter from hydrometeors, with a much shorter wavelength of λ= 3–5 cm. The longer wavelength of GNSS carrier waves makes the phase measurement more sensitive to the overall bulk properties of the hydrometeor particles than to their detailed shapes (Turk et al.2021), whereby making KDP at these frequencies more sensitive to snow particles than at lower frequencies, which may justify the simple linear KDP/ WC relation.

5.2 Sensitivity to displacement

In Sect. 4.2 we saw that the poor OB fit for tropical cyclone (TC) cases is partly due to the high sensitivity of ΦDP to the displacement of clouds. Even a small shift in the latitude and longitude of only 0.1, which is equivalent to around 10 km (just one grid point of the deterministic high-resolution model), can lead to completely different simulations for TC cases. While this is helpful as it may inform the model about its incorrect TC positions through observations, it poses a methodological challenge for the data assimilation system.

Consider, for example, a scenario where the observed ΦDP is greater than the background ΦDP due to the misplaced TC position. In such a situation, positive O-B departures can be “corrected” in many different ways, such as locally increasing snow along the ray (which would be a wrong correction), changing the refractivity or temperature so that the ray passes through areas of intense snow (also a wrong correction), or shifting the position of the TC in the background fields (which would be the right correction). Out of these many possibilities, the data assimilation method needs to correct the background fields in the right way, but given the sparsity of GNSS-PRO observations, it is not obvious whether the information content provided by such observations is rich enough to constrain the correction in the right direction.

Apart from the sparsity of observations, correcting displacement of the background fields is also difficult because it is known to induce non-Gaussianity in the probability distribution of the background errors (e.g. Chen and Snyder2007; Aonashi and Eito2011).

5.3 Overestimation of ΦDP in TC cases

In this study, we have assumed a linear relation between KDP and hydrometeor water content variables as we discussed in Sect. 2.2. Despite such a simple assumption, our forward operator achieves remarkably good simulations for AR cases. Yet, the simulated ΦDP values are systematically overestimated for TC cases, which deserves to be explored.

In our forward operator we assumed that the axis ratio (ar) is constant at the arbitrarily chosen value of 0.5. While this choice resulted in ΦDP simulations that are in remarkably good agreement with the observations in AR cases, its validity may need to be reconsidered for TC cases. There are several pieces of observational evidence that ar should be larger (i.e., snow particles should be less horizontally oriented) in deep convective clouds than in stratiform clouds because strong turbulent mixing inside deep convection randomises particle orientation (e.g. Gong and Wu2017). It may thus be worth allowing ar to vary in our formulation depending on the strength of mixing or vertical velocity of the background fields (Ramón Padullés, personal communication, 2023).

5.4 Importance of 2D ray tracing

We found that, unlike the successful ΦDP simulations with the 2D operator, the 1D operator fails to reproduce the observed ΦDP, even for AR cases where the 2D operator performs very well, which underlines the importance of incorporating the 2D aspect in ray tracing calculation. This is in contrast to the case of regular GNSS-RO bending-angle assimilation where a 1D operator is considered to be accurate enough to allow for extraction of meaningful information content from observations, although additional benefit is demonstrated with a 2D operator.

At the moment ECMWF is the only operational NWP centre to perform 2D ray tracing in assimilating GNSS-RO observations operationally. Our results suggest that, when other centres start investigation on GNSS-PRO assimilation, they would need to start by first extending their RO forward operator to adopt 2D ray tracing. Depending on how the code is parallelised, this alone can be non-trivial work.

5.5 Future directions

This study investigated characteristics of our implementation of a forward operator for the GNSS-PRO observable ΦDP. To the authors' knowledge, this is the first ΦDP forward operator ever implemented for an NWP model. While our first implementation demonstrated promising results, especially for the synoptic-scale atmospheric river (AR) cases, the results also identified several challenges that warrant further investigations.

The key challenge in assimilating PRO observations would be to account for displacement error of the background, and this will be particularly important for smaller-scale phenomena such as tropical cyclones. While the currently operational 4DVar is known to be able to correct position errors in the background by assimilating dense observations like all-sky microwave radiances (e.g. Duncan et al.2022), it is not clear if such a correction is possible with horizontally sparse observations like GNSS-PRO, and further methodological development along this line might be needed. We speculate that such correction will likely be possible if PRO observations are made available from a constellation of receiver satellites (e.g. Turk et al.2019) to allow for dense sampling with multiple different azimuth angles.

The linear relation between KDP and hydrometeor contents that we adopted is found to be quite successful despite its simplicity, but its limitations have also been identified. To better account for a wider range of weather situations, it would be worth exploring a more elaborate KDP–WC relation. To this end, integration with RTTOV-SCATT would be beneficial because that allows the assumptions on microphysical properties like particle size distribution to be more consistent across different components of the NWP system.

In this study we focused on simulating the polarimetric differential phase shift ΦDP as the observable of GNSS-PRO, but ΦDP is not the only GNSS-PRO observable. Wang et al. (2021) introduced polarimetric bending angles as an alternative observable quantity and showed that polarimetric bending angles can be less prone to issues with the multipath effect, which may be beneficial especially for measurements at low altitude. It would be thus be worthwhile to also explore building forward operator for polarimetric bending angles.

Code and data availability

The PAZ data are available for download from (Cardellach et al.2022) upon registration. Some additional metadata for PAZ data were retrieved from UCAR-processed level-2 data, which are publicly available from (UCAR COSMIC Program2018). The forecast data were produced with ECMWF's IFS Cy47R3 suite (, ECMWF2019). We have made the IFS forecast data used in this study available for download at (Hotta et al.2023). The forward operator for PRO observations developed in this study is based on the ROPP code (Culverwell et al.2015,, which is available free of charge from (last access: 5 February 2024) after agreeing to license conditions and completing user registration. We plan to put the PRO code in the ROPP package in a later release, but the code is still under development. Every detail necessary to reproduce the PRO code is documented in Sect. 2. To convert hydrometeor vertical mass fluxes to water content variables, we used code from RTTOV-SCATT (Saunders et al.2020), which can be obtained at (last access: 5 February 2024) after agreeing to license conditions and completing user registration.

Author contributions

DH contributed to the conceptualisation, investigation, methodology, software, visualisation and writing. KL contributed to the conceptualisation, methodology, writing and software. SH contributed to conceptualisation, software, interpretation and writing.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The authors thank Estel Cardellach and Ramon Padullés of ICE-CSIC/IEEC in Spain and Joe Turk of NASA/JPL for kindly providing PAZ data and for guiding us on the use and interpretation of their data. Michael Murphy and Jennier Haase are also acknowledged for selecting the AR cases examined. The authors also thank Peter Bechtold, Alan Geer, Robin Hogan and a number of other colleagues at ECMWF for their support. The authors thank Josep M. Aparicio and Joe Turk for their careful review of the manuscript and for constructive comments. This study is an outcome from a collaboration between ECMWF and Japan Meteorological Agency (JMA). The authors acknowledge the financial support of the Space-related Overseas Fellowship Program offered by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of the government of Japan.

Financial support

This research has been supported by the Ministry of Education, Culture, Sports, Science and Technology (MEXT) of the government of Japan.

Review statement

This paper was edited by Peter Alexander and reviewed by Joe Turk and Josep M. Aparicio.


Aonashi, K. and Eito, H.: Displaced ensemble variational assimilation method to incorporate microwave imager brightness temperatures into a cloud-resolving model, J. Meteorol. Soc. Jpn., 89, 175–194, 2011. a

Bonavita, M.: On some aspects of the impact of GPSRO observations in global numerical weather prediction, Q. J. Roy. Meteor. Soc., 140, 2546–2562, 2014. a

Bringi, V. N. and Chandrasekar, V.: Polarimetric Doppler weather radar: principles and applications, Cambridge University Press,, 2001. a, b, c

Cardellach, E., Oliveras, S., Rius, A., Tomás, S., Ao, C., Franklin, G., Iijima, B., Kuang, D., Meehan, T., Padullés, R., de la Torre Juárez, M., Turk, F. J., Hunt, D. C., Schreiner, W. S., Sokolovskiy, S. V., Van Hove, T., Weiss, J. P., Yoon, Y., Zeng, Z., Clapp, J., Xia-Serafino, W., and Cerezo, F..: Sensing heavy precipitation with GNSS polarimetric radio occultations, Geophys. Res. Lett., 46, 1024–1031, 2019. a, b, c

Cardellach, E., Padullés, R., and Oliveras, S.: Radio Occultation and Heavy Precipitation with PAZ (ROHP-PAZ), ICE-CSIC/IEEC [data set],​ (last access: 5 February 2024), 2022. a

Chen, Y. and Snyder, C.: Assimilating Vortex Position with an Ensemble Kalman Filter, Mon. Weather Rev., 135, 1828–1845, 2007. a

Culverwell, I. D., Lewis, H. W., Offiler, D., Marquardt, C., and Burrows, C. P.: The Radio Occultation Processing Package, ROPP, Atmos. Meas. Tech., 8, 1887–1899,, 2015 (code available at:, last access: 5 February 2024). a

Doviak, R. J. and Zrnić, D. S.: Doppler Radar and Weather Observations, Academic Press, San Diego, CA, 2nd edn., ISBN 978-0-12-221422-6, 1993. a

Duncan, D. I., Bormann, N., Geer, A. J., and Weston, P.: Assimilation of AMSU-A in All-Sky Conditions, Mon. Weather Rev., 150, 1023–1041,, 2022. a

ECMWF: IFS Documentation CY47R3, (last access: 5 February 2024), 2019. a

Forbes, R., Laloyaux, P., and Rodwell, M.: IFS Upgrade Improves Moist Physics and Use of Satellite Observations, ECMWF Newsletter, 169, 17–24, 2021. a

Geer, A. J., Bauer, P., and Lopez, P.: Lessons learnt from the 1D+ 4D-Var assimilation of rain-and cloud-affected SSM/I observations at ECMWF, ECMWF Tech. Memo., 535, 49 pp., 2007. a, b

Gong, J. and Wu, D. L.: Microphysical properties of frozen particles inferred from Global Precipitation Measurement (GPM) Microwave Imager (GMI) polarimetric measurements, Atmos. Chem. Phys., 17, 2741–2757,, 2017. a

Healy, S. B., Eyre, J. R., Hamrud, M., and Thépaut, J. N.: Assimilating GPS radio occultation measurements with two‐dimensional bending angle observation operators, Q. J. Roy. Meteor. Soc., 133, 1213–1227, 2007. a, b, c, d, e

Hotta, D., Lonitz, K., and Healy, S.: Deterministic ECMWF forecast data with additional precipitation flux parameters, European Centre for Medium-Range Weather Forecasts [data set],, 2023. a

Kawabata, T., Bauer, H.-S., Schwitalla, T., Wulfmeyer, V., and Adachi, A.: Evaluation of Forward Operators for Polarimetric Radars Aiming for Data Assimilation, J. Meteorol. Soc. Jpn., Ser. II., 96A, 157–174, 2018a. a

Kawabata, T., Schwitalla, T., Adachi, A., Bauer, H.-S., Wulfmeyer, V., Nagumo, N., and Yamauchi, H.: Observational operators for dual polarimetric radars in variational data assimilation systems (PolRad VAR v1.0), Geosci. Model Dev., 11, 2493–2501,, 2018b. a

Kumjian, M. R.: Principles and applications of dual-polarization weather radar. Part I: Description of the polarimetric radar variables, Journal of Operational Meteorology, 1, 226–242,, 2013. a, b

Li, X. and Mecikalski, J. R.: Evaluation of The Sensitivity of The Dual-Polarization Doppler Warm-Rain Radar Data Assimilation to Radar Forward Operators for a Convective Storm, J. Meteorol. Soc. Jpn., Ser. II, 91, 287–304,, 2013.  a

Murphy, M. J., Haase, J. S., Padullés, R., Chen, S. H., and Morris, M. A.: The potential for discriminating microphysical processes in numerical weather forecasts using airborne polarimetric radio occultations, Remote Sens.-Basel, 11, 2268,, 2019. a, b

Padullés, R., Ao, C. O., Turk, F. J., de la Torre Juárez, M., Iijima, B., Wang, K. N., and Cardellach, E.: Calibration and validation of the Polarimetric Radio Occultation and Heavy Precipitation experiment aboard the PAZ satellite, Atmos. Meas. Tech., 13, 1299–1313,, 2020. a

Padullés, R., Cardellach, E., Turk, F. J., Ao, C. O., de la Torre Juárez, M., Gong, J., and Wu, D. L.: Sensing Horizontally Oriented Frozen Particles With Polarimetric Radio Occultations Aboard PAZ: Validation Using GMI Coincident Observations and Cloudsat a Priori Information, IEEE T. Geosci. Remote, 60, 1–13, 2021. a, b, c, d

Saunders, R., Hocking, J., Turner, E., Havemann, S., Geer, A., Lupu, C., Vidot, J., Chambon, P., Köpken-Watts, C., Scheck, L., Stiller, O., Stumpf, C., Borbas, E., and Brunel, P.: RTTOV-13 science and validation report, NWP-SAF report NWPSAF-MO-TV-046, EUMETSAT NWP-SAF, Met Office, Exeter, UK, (last access: 5 February 2024), 2020. a

Turk, F. J., Padullés, R., Ao, C. O., Juárez, M. d. l. T., Wang, K.-N., Franklin, G. W., Lowe, S. T., Hristova-Veleva, S. M., Fetzer, E. J., Cardellach, E., Kuo, Y.-H., and Neelin, J. D.: Benefits of a closely-spaced satellite constellation of atmospheric polarimetric radio occultation measurements, Remote Sens.-Basel, 11, 2399,, 2019. a

Turk, F. J., Padullés, R., Cardellach, E., Ao, C. O., Wang, K.-N., Morabito, D. D., Juarez, M. d. l. T., Oyola, M., Hristova-Veleva, S., and Neelin, J. D.: Interpretation of the precipitation structure contained in polarimetric radio occultation profiles using passive microwave satellite observations, J. Atmos. Ocean. Tech., 38, 1727–1745, 2021. a, b, c, d, e

UCAR COSMIC Program: PAZ Data Products, UCAR/NCAR - COSMIC [data set],, 2018. a

Wang, K.-N., Ao, C. O., Padullés, R., Turk, F. J., Juárez, M. d. l. T., and Cardellach, E.: The effects of heavy precipitation on polarimetric radio occultation (PRO) bending angle observations, J. Atmos. Ocean Tech., 39, 161–173,, 2021. a

Yokota, S., Seko, H., Kunii, M., Yamauchi, H., and Niino, H.: The tornadic supercell on the Kanto Plain on 6 May 2012: Polarimetric radar and surface data assimilation with EnKF and ensemble-based sensitivity analysis, Mon. Weather Rev., 144, 3133–3157, 2016. a

Short summary
Global Navigation Satellite System (GNSS) polarimetric radio occultation (PRO) is a new type of GNSS observations that can detect heavy precipitation along the ray path between the emitter and receiver satellites. As a first step towards using these observations in numerical weather prediction (NWP), we developed a computer code that simulates GNSS-PRO observations from forecast fields produced by an NWP model. The quality of the developed simulator is evaluated with a number of case studies.