Forward operator for polarimetric radio occultation measurements

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. Bonavita, 2014).

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. Kumjian, 2013). 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 ${\mathrm{\Phi }}_{\mathrm{DP}}={\mathrm{\Phi }}_{\mathrm{H}}-{\mathrm{\Phi }}_{\mathrm{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 K_DP:

$$\begin{array}{}\text{(1)}& {\mathrm{\Phi }}_{\mathrm{DP}}=\underset{\mathrm{GNSS}}{\overset{\mathrm{LEO}}{\int }}{K}_{\mathrm{DP}}\left(s\right)\mathrm{d}s,\end{array}$$

where GNSS and LEO symbolically represent, respectively, the position of the transmitter and receiver of GNSS radio signals. K_DP 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 K_DP 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 K_DP (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

$$\begin{array}{}\text{(2)}& \mathit{\alpha }=\underset{\mathrm{GNSS}}{\overset{\mathrm{LEO}}{\int }}\left({\displaystyle \frac{\mathrm{d}\mathit{\alpha }}{\mathrm{d}s}}\right)\mathrm{d}s,\end{array}$$

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 $\mathrm{\Delta }\mathit{\theta }=\mathrm{40}/\mathrm{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 $\mathrm{\Delta }\mathit{\theta }=\mathrm{20}/\mathrm{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 K_DP contributions from each type of hydrometeors (see next section) are integrated as each section of the ray is traced, and the accumulated K_DP 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 K_DP is performed column-wise. While K_DP can be negative in nature in situations, for example, where particles in prolate orientation dominate over oblate ones, in our modelled formulation of K_DP (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 K_DP.

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

Specific differential phase shift, K_DP, is induced by the difference in scattering properties of hydrometeor particles for the horizontally and vertically polarised waves. A “first-principle” scattering calculation of K_DP 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 K_DP 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 K_DP from the water content of each category of hydrometeor, here we adopt a simple linear relation between the water content and K_DP contribution. For ice (i.e., IWC, SWC and SWCconv), following Padullés et al. (2021), we adopt the following water content-to-K_DP formula (which, in turn, is taken from Bringi and Chandrasekar, 2001):

$$\begin{array}{}\text{(3)}& K_{\mathrm{DP}}\left(\text{WC}\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}C\mathit{\rho }\times \text{WC}\times \left(\mathrm{1}-\text{ar}\right),\end{array}$$

where WC means any one of IWC, SWC or SWCconv (in units of g m⁻³); K_DP(WC) (in units of mm km⁻¹) denotes the K_DP contribution from the hydrometeor whose water content is denoted WC; ρ is the particle density in units of g cm⁻³; 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⁻³)⁻² 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⁻³.

The main approximating assumptions behind the formula relating K_DP 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 Chandrasekar, 2001), 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 K_DP, at L-band, in which the Rayleigh limit is appropriate, is nearly proportional to the third moment of the PSD (Bringi and Chandrasekar, 2001, 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, K_DP is less sensitive to the precise shape of the PSD. In fact, as shown in Fig. 11 in Turk et al. (2021), K_DP 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), K_DP 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⁻¹), which can be readily converted to mass per volume (in units of g m⁻³). 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⁻² s⁻¹). 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⁻³, can be derived from the vertical mass flux FL as ${\left\{\text{FL}/\left(a\mathit{\rho }\right)\right\}}^{\mathrm{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 1 List 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.

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

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 https://paz.ice.csic.es/, 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 1 Comparison 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 K_DP (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 K_DP 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 2 Uncertainty 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 $\mathrm{3}\times \mathrm{3}=\mathrm{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 3 Impact 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 4 As 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 K_DP, 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 5 Resolved-scale snow water content (in g m⁻³) 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 K_DP. 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 K_DP) is dominated by snow in reality because our modelled K_DP–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 K_DP 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 K_DP measurements (e.g. Kumjian, 2013). Because of this, most forward operators that have been developed to simulate or assimilate radar K_DP observations only consider warm rain conditions (e.g. Li and Mecikalski, 2013; 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 K_DP values for snow/ice particles are relatively small compared to those for rain droplets (in fact, K_DP 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 K_DP 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 K_DP at these frequencies more sensitive to snow particles than at lower frequencies, which may justify the simple linear $K_{\mathrm{DP}}/$ WC relation.

5.2 Sensitivity to displacement

In Sect. 4.2 we saw that the poor O−B 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 Snyder, 2007; Aonashi and Eito, 2011).

5.3 Overestimation of Φ_DP in TC cases

In this study, we have assumed a linear relation between K_DP 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 Wu, 2017). 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 K_DP 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 K_DP–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.