Determination of ice water content (IWC) in tropical convective clouds from X-band dual-polarization airborne radar

This paper presents a methodology for ice water content (IWC) retrieval from a dual-polarization side-looking X-band airborne radar. Measured IWC from aircraft in situ probes is weighted by a function of the radar differential reflectivity (Z_dr) to reduce the effects of ice crystal shape and orientation on the variation in IWC – specific differential phase (K_dp) joint distribution. A theoretical study indicates that the proposed method, which does not require a knowledge of the particle size distribution (PSD) and number density of ice crystals, is suitable for high-ice-water-content (HIWC) regions in tropical convective clouds. Using datasets collected during the High Altitude Ice Crystals – High Ice Water Content (HAIC-HIWC) international field campaign in Cayenne, French Guiana (2015), it is shown that the proposed method improves the estimation bias by 35 % and increases the correlation by 4 % on average, compared to the method using specific differential phase (K_dp) alone.

1 Introduction

Ice water content (IWC) and its spatial distribution inside clouds are known for the significant effects they exert on the Earth's energy budget and hydrological cycle (e.g. Stocker et al., 2013). Aside from its significant effect on the atmospheric processes, high ice water content (IWC >1 g m⁻³), which is resultant from a high concentration of small ice crystals in tropical mesoscale convective systems, has been linked to aircraft incidents and accidents (Lawson et al., 1998; Mason et al., 2006; Grzych and Mason, 2010; Strapp et al., 2016). Since early 1990s, over 150 engine rollback and power-loss events have been attributed to the ingestion of ice particles produced in convective clouds (Grzych and Mason, 2010). Many studies have been undertaken to understand the details of the meteorological processes responsible for producing areas of HIWC. Equally important, methods using multiplatform observations from ground, airborne, and space supplemented by weather models are being developed for improving detection and avoidance of high-IWC regions that would be potentially hazardous for aviation (Strapp et al., 2016).

Conventional methods of deducing IWC from radar measurements assume a statistical relationship between the radar reflectivity factor (Z) and IWC. Such relationships are usually obtained based on IWC and Z calculated from in situ measurements of particle size distributions (PSDs) and a size-to-mass parameterization (m(D)) (e.g. Heymsfield, 1977; Hogan et al., 2006). In recent studies (Protat et al., 2016), IWC was measured directly by bulk microphysical probes and Z was measured from either an airborne or ground-based radar. However, all of these studies show large uncertainties in the IWC–Z relationship despite the introduction of additional constraints such as air temperature (T) or the inclusion of refined m(D) in the IWC calculations (Fontaine et al., 2014; Protat et al., 2016).

Lu et al. (2015) conducted an extensive simulation on both millimetre- and centimetre-wavelength radar and concluded that the IWC–Z relationship is very sensitive to ice crystal PSDs (from 1 to 2 orders of magnitude in variability) and as such is not recommended for IWC retrievals. Another approach employs polarimetric observations. The non-spherical geometry of ice crystals provides information on the types and habits of ice crystals (Matrosov et al., 1996; Wolde and Vali, 2001). It has been shown that the radar-specific differential phase (K_dp) is less dependent on PSD; hence, it is potentially useful for IWC retrieval (Vivekanandan et al., 1994; Lu et al., 2015). Aydin and Tang (1997) suggested the possibility of combining K_dp and differential reflectivity ratio (Z_dr) for IWC estimation for clouds composed of pristine ice crystals. However, even for the polarimetric approach, knowledge about ice crystal mass density (ρ) and axis ratio is still needed to obtain accurate estimates of IWC. Simulation results (Lu et al., 2015) show that if only the general type of ice crystals is known, errors in IWC retrieval based on K_dp are within 30 % of their true values. Unfortunately, the aforementioned parameters (ρ and particle axis ratio) are, in general, unknown and additional assumptions are often invoked. Ryzhkov et al. (1998), for instance, took into consideration ice crystal shapes and size–density parameterization of scatterers to reduce the uncertainty in IWC estimates. Modelling work (Ryzhkov et al., 1998) shows that for average-sized pristine and moderately aggregated ice crystals, the ratio between the reflectivity difference $Z_{\mathrm{DP}}=Z_H-Z_V$ and K_dp is practically insensitive to the shape and density of the ice particles and is a good estimator of their mass.

In this paper we present a new method for assessment of IWC based on the K_dp and Z_dr measurements from a side-looking X-band airborne radar in tropical mesoscale convective systems (MCSs). The IWC will be weighted with a function of Z_dr to minimize the dependency of the IWC–K_dp relationship on the particle shape and orientation, hence improving the IWC estimation errors without knowledge of the PSD or density of the ice particles. The proposed method is examined using datasets collected during the High Altitude Ice Crystal – High Ice Water Content (HAIC-HIWC) international field campaign in Cayenne, French Guiana, in May 2015. The campaign was carried out to enhance the knowledge of microphysical properties of high-altitude ice crystals and mechanisms of their formation in deep tropical convective systems in order to address aviation safety issues related to engine icing (Strapp et al., 2016).

2 Background

2.1 Polarimetric parameters characterizing ice crystals

In conventional single-polarization Doppler radar, measured radar reflectivity and radial velocity are used to assess cloud and precipitation spatial variability, precipitation rate, and characteristic hydrometeor types. In dual-polarization radar systems, measurements are made at more than one polarization state (Bringi and Chandrasekar, 2001). Such systems can be configured in several ways depending on the measurement goals and the choice of orthogonal polarization states. In this study, the results and discussions will be limited to the consideration of the linear horizontal and vertical (H∕V) polarization basis. The intrinsic backscattering properties of the hydrometeors for the two polarization states enable the characterization of microphysical properties such as size, shape, and spatial orientation of the cloud/precipitation particles in the radar resolution volume. Hence, using polarization, it is generally possible to achieve more accurate classification of hydrometeor types and estimate hydrometeor amounts such as rainfall rate. Polarimetric backscattering properties of hydrometeors depend on many factors such as radar wavelength, radar elevation angle, particle size, shape, orientation, etc. In this section, we summarize how the differential reflectivity (Z_dr, dB) and the specific differential phase (K_dp, ^∘ km⁻¹) are measured by a polarimetric Doppler radar in the Rayleigh scattering regime and at low radar elevation angles.

In general, the differential reflectivity of an ensemble of n particles of size D and axis ratio r is given by Eq. (1) (Eq. 7.4 in Bringi and Chandrasekar, 2001),

$$\begin{array}{}\text{(1)}& Z_{\mathrm{dr}}=\mathrm{10}{\log }_{\mathrm{10}}\left[{\displaystyle \frac{{\left|S_{hh}\left(r,D\right)\right|}^{\mathrm{2}}}{{\left|S_{vv}\left(r,D\right)\right|}^{\mathrm{2}}}}\right],\end{array}$$

where S_hh and S_vv are the diagonal elements of the backscattering matrices in the forward-scatter alignment (FSA) convention.

The specific differential phase is defined as (Eq. 7.6 in Bringi and Chandrasekar, 2001)

$$\begin{array}{}\text{(2)}& K_{\mathrm{dp}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }n}{k}}Re\left[f_{hh}\left(r,D\right)-f_{vv}\left(r,D\right)\right],\end{array}$$

where n is the number concentration in reciprocal litres, k is wavenumber in reciprocal metres, Re[] stands for the real part of a complex number, and f_hh and f_vv are the forward-scattering amplitudes in m at horizontal and vertical polarization, respectively. Equation (2) shows that Z_dr does not change with increasing number of ice particles while K_dp is proportional to n. Consequently, for a large number of small particles with the axis ratio close to unity (r≈1), Z_dr→0 dB and the second term in Eq. (2) becomes small but K_dp can still be large.

In a simple form of the calculations of Z_dr and K_dp of ice crystals, it is customary to approximate columns as homogeneous prolate spheroids and plates as homogeneous oblate spheroids. In the case of side incidence, the elevation angle is assumed to be close to zero and there is no (or very small) canting in the vertical plane. In the absence of wind shear and turbulence, and assuming a perfectly aligned spheroid model, Z_dr and K_dp can be expressed as functions of ice particle size, axis ratio, and the relative permittivity of the particle (ε). For example, for oblate spheroid ice particles with a particle size distribution, N(D) (Eq. 7.5–7.8 in Bringi and Chandrasekar, 2001),

$$\begin{array}{}\text{(3)}& {\displaystyle }\left|S_{hh}\left(r,D\right)\right|\approx {\displaystyle \frac{k^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{V\left(D\right)|\mathit{\epsilon }-\mathrm{1}|}{\left[\mathrm{1}+\frac{\mathrm{1}}{\mathrm{2}}\left(\mathrm{1}-{\mathit{\lambda }}_o\right)|\mathit{\epsilon }-\mathrm{1}|\right]}},\text{(4)}& {\displaystyle }\left|S_{vv}\left(r,D\right)\right|\approx {\displaystyle \frac{k^{\mathrm{2}}}{\mathrm{4}\mathit{\pi }}}{\displaystyle \frac{V\left(D\right)|\mathit{\epsilon }-\mathrm{1}|}{\left[\mathrm{1}+{\mathit{\lambda }}_o|\mathit{\epsilon }-\mathrm{1}|\right]}},\text{(5)}& {\displaystyle }{Z}_{hh,vv}={\displaystyle \frac{{\mathit{\lambda }}^{\mathrm{4}}}{{\mathit{\pi }}^{\mathrm{5}}{K}_p^{\mathrm{2}}}}\int \mathrm{4}\mathit{\pi }|S_{hh,vv}{|}^{\mathrm{2}}N\left(D\right)\mathrm{d}D,\text{(6)}& {\displaystyle }{K}_{\mathrm{dp}}={\displaystyle \frac{k}{\mathrm{2}}}\int V\left(D\right)N\left(D\right)\mathrm{d}D,\end{array}$$

where K_p is dielectric factor of water at 0 ^∘C ($K_p^{\mathrm{2}}=\mathrm{0.93}$); V(D) is the particle volume; λ_o is the depolarizing factor, which is only a function of the axis ratio $r=b/a$ (for oblate particles; a is the semi-major axis length; and b is the semi-minor axis length (a>b)). The depolarizing factor is defined as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\lambda }}_o=& \mathit{\lambda }\left(\text{oblate}\right)={\displaystyle \frac{\mathrm{1}+f^{\mathrm{2}}}{f^{\mathrm{2}}}}\left(\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{f}}{\tan }^{-\mathrm{1}}f\right);\\ & f^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}-\mathrm{1}\end{array}.\end{array}$$

A similar equation for K_dp can also be derived for prolate spheroid ice particles with symmetry axis parallel to the horizontal plane (Hogan et al., 2006).

On the other hand, the IWC can be defined in terms of the size distribution,

$$\begin{array}{}\text{(8)}& \mathrm{IWC}=\int \mathit{\rho }\left(D\right)V\left(D\right)N\left(D\right)\mathrm{d}D,\end{array}$$

where ρ(D) is the mass density of ice crystals with size D.

2.2 Polarimetric methods for IWC retrieval

An inspection of Eqs. (3) and (4) suggests that for small ice crystal particles, the radar cross section (${\mathit{\sigma }}_{hh,vv}=\mathrm{4}\mathit{\pi }{\left|S_{hh,vv}\right|}^{\mathrm{2}}$) is roughly proportional to the square of the ice particle mass (${\mathit{\rho }}_i^{\mathrm{2}}\left(D\right)V(D{)}^{\mathrm{2}}$), a conclusion also confirmed by results from simulated data (Lu et al., 2015). In addition, according to Lu et al. (2015), for particles with sizes comparable to or larger than the radar wavelength, there is no clear relationship between the radar cross section and ice particle mass due to the Mie resonance effects. In either case, σ_hh,vv is not directly proportional to the particle mass. Hence, the Z–IWC relationship depends strongly on the particle size distribution and the radar frequency. Consequently, using Z only to estimate IWC without knowledge of the PSD can lead to errors as large as 1 order of magnitude. On the other hand, Eq. (6) indicates that if the terms in square brackets (α) are proportional to the ice density (ρ(D)), then the K_dp–IWC relationship is independent of PSD. The proportionality constant depends on several factors such as the ice crystal type, orientation, and the measurement elevation angle. It is shown that the variability of this proportionality constant significantly increases at large elevation angles (Lu et al., 2015). Furthermore, when the exact ice crystal type is known, averaged relative error in the estimated IWC using K_dp can be as small as 10 %, regardless of whether PSD is known or not. If the ice crystal types are unknown but can be generally categorized, the errors can be higher, but mostly less than 30 %. These numbers were averaged from elevations in the interval [0–70^∘]. If IWC is estimated using K_dp at small elevation angles (less than 10^∘) such as from a side-looking antenna, we would expect better results.

For a given radar volume, if the orientation of the ice crystal changes, the K_dp value changes (Eq. 7) while the IWC of the radar volume does not. Consequently, in the case of spatial variability of ice crystal shapes and orientations, the IWC estimation based solely on K_dp may be biased. To mitigate this problem, the measured IWC needs to be modified to include the information of the ice particles' orientation. One way to do this is to weight the measured IWC by a function of ice crystal shapes and orientations before applying a linear regression model to the K_dp–IWC relationship. In a simple approach, the weighting function can be in a form of $Z_{\mathrm{DR}}^a$ (Z_DR is the linear version of Z_dr and a is a constant coefficient) as suggested in Aydin and Tang (1997) (derived from their approximation IWC $\approx K_{\mathrm{dp}}^a{Z}_{\mathrm{DR}}^b$, where a and b are constant coefficients). Proceeding more rigorously, Ryzhkov et al. (1998, 2018) demonstrated that both K_dp and difference reflectivity Z_DP ($Z_{\mathrm{DP}}=Z_H-Z_V$) are dependent on the particles' aspect ratios and orientation, whereas their ratio is very robust with respect to those factors. Indeed, in simulation and modelling work considering 12 different crystal habits, Ryzhkov et al. (2018) showed that the ratio Z_DP∕K_dp in combination with reflectivity can be used to estimate IWC. In detail, for exponential size distribution and with an assumption of $\mathit{\rho }\left(D\right)=\mathit{\alpha }{D}^{-\mathit{\beta }}$, β≈1, $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)\mathrm{IWC}$ is proportional to K_dp. Also, according to Ryzhkov et al. (2018), this approximation is almost insensitive to the ice habit, aspect ratio, and orientation of the ice particles, but is affected by the degree of riming. Hence, it works better for clouds with a low degree of riming. This condition might not be true for all types of ice clouds, but might be suitable for HIWC regions, which are often composed of a high concentration of small ice particles (Leroy et al., 2016).

At Z_DR≈1 (or Z_dr≈0 dB), the weighting function $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$ is close to zero, and hence it can introduce large errors in the estimates. Therefore, there should be a certain threshold for Z_dr to determine how the weighting function would be calculated. In detail, if Z_dr is less than a threshold, the weighting function $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$ is replaced by $\left(\mathrm{1}-Z_{\text{DR-threshold}}^{-\mathrm{1}}\right)$. In this paper, we use Z_DR-threshold=1.12 (see Sect. 6 for more detailed derivation of this threshold value). This threshold is very close to a 1.15 threshold proposed by Ryzhkov et al. (1998) for “cold” storms for temperatures below −5 ^∘C.

In summary, there are two polarimetric methods for IWC retrieval, which will be investigated and compared in this paper. They are expressed as

$$\begin{array}{}\text{(9)}& {\displaystyle }\mathrm{IWC}=a_{\mathrm{1}}{K}_{\mathrm{dp}}+b_{\mathrm{1}},\text{(10)}& {\displaystyle }\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)\mathrm{IWC}=a_{\mathrm{2}}{K}_{\mathrm{dp}}+b_{\mathrm{2}},\end{array}$$

where model parameters (a_i,b_i) will be estimated from measured data. A flow chart of IWC retrieval using K_dp and Z_DR is shown in Fig. 1.

Figure 1 Flow chart of IWC retrieval using (K_dp,Z_DR) and in situ data. In this paper, Z_DR-threshold=1.12.

3 Airborne measurements

During the Cayenne HAIC-HIWC project, the NRC Convair CV-580 conducted 14 research flights in both continental and oceanic mesoscale convective systems (MCSs) with high IWC. All the flights were conducted during daytime only due to flight restrictions. Analysis of all the IR satellite imageries show that the Convair flights were performed close to the peak intensities of the targeted MCSs, which had a typical lifetime of 7.5 and 4.5 h for oceanic and continental MCSs, respectively (Walter Strapp, personal communication, 2019). So the fact that the flights were conducted only during the daytime did not prevent the flights from sampling of the storms during their peak intensities. For this campaign, the Convair aircraft was instrumented by the NRC and Environment and Climate Change Canada with an array of in situ cloud microphysics probes, atmospheric sensors, and the NRC airborne W- and X-band (NAWX) Doppler dual-polarization radars (Wolde and Pazmany, 2005). The unique quasi-collocated in situ and radar data collected during the HAIC-HIWC mission provided a means for developing techniques for detection and estimation of high IWC that could be adopted in operational airborne weather radars.

3.1 Airborne radar data

In this study, dual-polarization radar data from the NRC airborne X-band radar (NAX) (Fig. 2) side-looking antenna are used. Some important radar parameters are given in Table 1. More detailed information on the radar system can be found in Wolde and Pazmany (2005). In the Cayenne project, the radar complex I and Q samples are processed to powers and complex pulse pair products according to the radar parameter specification table, and the products are recorded in binary format. Due to the size of the aircraft radar radome, the NAX dual-polarization parabolic side antenna is relatively small (0.66 m), hence exhibiting some limitations in terms of side lobe performance. The antenna orthomode transducer (OMT)–feedhorn combination is relatively large compared to the parabolic dish. The large feed structure creates some significant side lobes at ±90^∘ planes. As a result, when the side lobes intercept targets with strong returns below the aircraft, such as the earth surface or a storm melting layer, significant returns from the side lobes will contaminate signals coming via the antenna's main lobe. In most situations, the effect is more prominent at a range equal to or greater than the distance where the antenna side lobes hit the ground. At regions where signals are contaminated by ground clutter via the side lobes, the data are intermittent and exhibit large biases. Unfortunately, with the pulse pair data from the Cayenne campaign, methods to separate clutter from the precipitation signals are limited. To overcome this issue, a method is developed to detect regions with strong clutter contamination based on signal correlations between the nadir and zenith returns. If the correlation coefficient exceeds a predefined threshold, the corresponding side data in those regions are discarded. If the width of the discarded data region is relatively small (less than 300 m in radar range) it will be filled through interpolation. In addition, due to the limitation of the radar hardware, the measurements of dual-polarization parameters are not useable below a range of 1000 m from the aircraft, but reflectivity can be measured accurately from 450 m. Hence, in this work, radar profiles were extracted at a horizontal distance of 1000 m from the aircraft. This is not an ideal condition when the in situ data and the radar data are not spatially coincident. However, in most scenarios the advantage of having fine radar sampling volumes with a high order of accuracy in time synchronization between in situ probes overcomes the location offset. At large distances from cloud boundaries and convective cores, the microphysics properties of glaciated clouds can be considered spatially quasi-uniform at scales of the order of a few hundred metres. This is specifically relevant to the measurements in MCSs during the HAIC-HIWC project. Moreover, there was no attenuation correction applied to reflectivity and Z_dr because in ice precipitation regions and at close range attenuation at the X-band is negligible.

Figure 2 The NRC Convair CV-580 and the dual-polarization side-looking X-band radar.

For the Cayenne project, the in situ microphysical data are processed at 1 Hz, which is lower than that of the radar data. Hence, the radar data were decimated to match with temporal resolution of the in situ data. At the Convair CV-580 average ground speed of 100 ms⁻¹, this results in a 100 m radar sampling volume.

Table 1 X-band radar parameters for the Cayenne campaign.

3.2 In situ data

For the project, the NRC Convair CV-580 was equipped with state-of-the-art in situ sensors for measurements of aircraft and atmospheric state parameters and cloud microphysics. There were multiple sensors to measure bulk liquid water content (LWC) and total water content (TWC), with hydrometeor size distribution ranging from small cloud drops to large precipitation particles. A detailed list of the Convair in situ sensors used during the Cayenne HIWC-HIWC project is provided in Wolde et al. (2016). Here we will briefly describe the in situ microphysical sensors used in correlating the airborne radar measurements with regions of HIWC. TWC was measured by an Isokinetic probe (IKP2) that was specifically designed to measure very high TWC (Davison et al., 2008). Alternatively, IWC was estimated from the measured PSDs with the D–M parameterization tuned using IKP2 measurements. In the Cayenne Convair datasets, IWCs calculated from PSDs and measured by IKP2 agreed quite well and the difference between them in the HIWC regions on average did not exceed 15 %. Because the IKP2 data were not available in all flights, estimated IWC from PSDs (IWC_PSD) has been used in this work. Additionally, mean mass diameter (MMD) was also used to characterize the microphysical properties of the high-IWC regions and interpret X-band radar measurements. MMD was calculated from composite particle size distributions measured by SPEC 2D-S and DMT PIP 2-D imaging probes.

As shown in Korolev et al. (2018) in the MCSs studied during the Cayenne HAIC-HIWC project, the fraction of mixed-phase clouds at −15 ^∘C $<T<-\mathrm{5}$ ^∘C did not exceed 4.6 %, and in most mixed-phase cloud regions LWC ≪ IWC. Hence, we did not filter out the very small fraction of liquid observed in our analysis; i.e. we assumed TWC = IWC. This finding significantly simplifies the processing and interpretation of cloud microphysical measurements.

Figure 3 The NRC K_dp estimation algorithm for X-band radar flow chart.

4 K_dp estimation algorithm for X-band airborne weather radar

The radar specific differential phase (K_dp) is defined as the slope of the range profile of the differential propagation phase shift Φ_dp between horizontal and vertical polarization states (Bringi and Chandrasekar, 2001). The measured differential phase shift between the two signals at the H and V polarizations, Ψ_dp, contains both Φ_dp and differential backscatter phase shift δ_dp. If δ_dp is relatively constant or negligible, the profile of Ψ_dp can be used to estimate K_dp. The estimated phase Ψ_dp usually exhibits discontinuities due to phase wrapping, statistical fluctuations in estimation, and the gate-to-gate variation in δ_dp. Because the statistical fluctuations in the estimates of Ψ_dp will be magnified during the differentiation, resulting in a large variance of the K_dp estimates, the following considerations need to be addressed in the K_dp estimation algorithm.

Phase unfolding. Phase wrapping occurs when the total Φ_dp accumulation exceeds the unambiguous ranges. This depends on the system differential phase Φ_dp(0) and the cumulative phase due to the medium. The NAX radar operates in the simultaneous transmission mode (VHS) and the unambiguous range is 360^∘. The system differential phase Φ_dp(0) of NAX is about 64^∘. For the Cayenne dataset, no observations were made when the phase was folded.

δ_dp “bump”. It seems that δ_dp was negligible in the HIWC environment in the Cayenne campaign. We did not observe the presence of significant changes in δ_dp over a short range.

Range filtering. In this work, the range scale was set at 500 m; thus the fluctuations at scales smaller than 500 m will be suppressed.

Once the phase data are quality controlled, filtered, and decimated to match the temporal resolution of the in situ data, a heuristic algorithm similar to one reported in Rotemberg (1999) is applied to the data to extract the Ψ_dp smooth trend and then K_dp is computed from it. This approach does not require an assumption of Φ_dp being a monotonically increasing function as it is in some other existing K_dp retrieval algorithms (Wang and Chandrasekar, 2009); therefore, it would also work well with negative K_dp, which possibly appears in ice clouds. Our preliminary analysis shows that the algorithm can provide estimates with standard deviation no greater than 1 ^∘ km⁻¹. The NRC K_dp estimation algorithm is summarized in the flow chart below.

Figure 4 Panel (a) shows IR GOES-13 image with the overlaid segments of the Convair CV-580 flight track on 26 May 2015. Different time segments of the flight track are shown by different colours. Panel (b) shows X-band side reflectivity from a period of 11:07–11:30 UTC corresponding to white, yellow, and purple segments in (a). A break line at around 7.1 km is the location of contaminating ground clutter via the side antenna's side lobe (Sect. 3.1), which was filtered out.

Figure 5 Time series of (a) IWC_PSD, MMD, (b) K_dp, Z_dr, (c) ρ_hv, and Z_H for 26 May Convair CV-580 flight.

5 Results

In this section, results illustrating the performance of the proposed polarimetric algorithms are presented. In addition to the polarimetric method, we also include results from the conventional IWC–Z relations for comparison. Because the histogram of static temperature (not shown) indicated a bimodal distribution with two centres at around −5 and −10 ^∘C, two IWC–Z relations at $T=-\mathrm{5}$ ^∘C (IWC =0.257Z^0.391) and at $T=-\mathrm{10}$ ^∘C (IWC =0.253Z^0.596) were obtained by fitting power-law curves to scatter plots of all the data points at those two temperature levels (Wolde et al., 2016).

5.1 Case study I: 26 May flight

In this case, a 20 min segment of the Convair flight inside an Oceanic MCS on 26 May 2105 is selected for analysis. The MCS was sampled north-northwest of Cayenne, French Guiana, during early morning hours. Figure 4a shows IR satellite imagery obtained during the flight where the aircraft's flight track is shown in different colours, which represent the aircraft's location at different time segments. The reflectivity field from the NAX side antenna is shown in Fig. 4b. The selected period begins at a point when the aircraft started to sample at the proximity of the convective core of the storm with the lowest cloud-top brightness temperature (white segment in the IR image). The brightness temperature was increasing toward the end of the segment (magenta segment). The aircraft flew between 6.9 and 7.2 km altitude and the static air temperature (T_s) varied from −12.8 to −8.2 ^∘C. In addition to the radar data, IWC_PSD and MMD time series from particle probes are shown in Fig. 5. The IWC estimates from radar data have been decimated to match with the temporal resolution of the in situ data.

The aircraft sampled two regions: a convective region before 11:23 UTC and a stratiform region after 11:25 UTC (Fig. 4b), with IWC in both regions mostly higher than 1.5 g m⁻³ (Fig. 5a). It is worth noticing that the reflectivity measurements along the flight path were fairly constant at ∼20 dBZ and the MMD values were relatively small at HIWC regions (IWC >1.5 g m⁻³) (Fig. 5a). From Fig. 5, it follows that (1) K_dp, in general, is highly correlated with IWC; (2) K_dp increases at the regions dominated by small ice crystals (between 11:08 to 11:16 and 11:24 to 11:27 UTC); and (3) regions with larger MMD exhibit deceasing ρ_hv and increasing Z_dr. In Fig. 6, Z_dr, ρ_hv, and IWC are expressed as functions of K_dp. In this case, there is a break point at K_dp≈1.5 ^∘ km⁻¹ (and Z_DR∼1.12) where Z_dr started increasing and ρ_hv decreased with respect to K_dp. At K_dp<1 ^∘ km⁻¹, Z_dr was mainly flat and IWC linearly increased with respect to K_dp (Fig. 6b). This suggests the pristine ice crystals' axis ratio might be fairly constant, but the particle number density increased, resulting in an enhancement in both K_dp and IWC (shown by a linear IWC–K_dp relationship). From K_dp>1 ^∘ km⁻¹, the Z_dr increment with respect to K_dp was greater, but the IWC increase does not follow the same degree as in the previous segment. If a linear IWC–K_dp relationship derived from the first segment (K_dp<1 ^∘ km⁻¹) is applied to the second portion, IWC will be overestimated. It is not easy to identify the exact reasons of this observation. Many factors could contribute to this circumstance such as changes in ice crystals' size, shape, orientation (e.g. particles with higher axis ratio that are aligned in the horizontal plane), or particle density. In this work we used K_dp and Z_dr to mitigate this dependency and improve estimation of IWC. In Fig. 6b and c, measured IWC and $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$IWC are shown as solid black lines and their linear fitting curves (red lines) are superimposed. The R² goodness-of-fit parameter indicates that a linear regression fits $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$IWC better in comparison to IWC.

Figure 6 Z_dr and ρ_hv (a), IWC_PSD (b), and $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)\mathrm{IWC}$ (c) as functions of K_dp. In panels (b) and (c), mean values and frequency distributions are computed from data points in each K_dp bin of 0.1 and 0.05 ^∘ km⁻¹, respectively. Regression parameters (a₁,b₁) for the K_dp-only method and (a₂,b₂) for the (K_dp,Z_DR) method are estimated from the mean values using a simple linear fitting algorithm.

Figure 7 Panel (a) shows IWC_PSD (black line), estimated IWC using reflectivity (dashed green line), K_dp alone (blue line), and the (K_dp, Z_DR) combination (red line) for the 26 May case. Panel (b) shows estimation biases for the three estimators. Average biases for IWC(Z), IWC(K_dp), and IWC(K_dp, Z_DR) are 0.60, 0.082, and 0.018 g m⁻³ and rms differences are 0.98, 0.49, and 0.48 g m⁻³ for the three algorithms, respectively.

To gauge the performance of the polarimetric methods, results from the conventional IWC–Z estimator are also included in Fig. 7a. This figure shows the measured IWC along the Convair's flight path depicted in black, IWC–Z results in green, and IWC estimates using polarimetric methods in blue and red for K_dp-only and (K_dp, Z_DR) algorithms, respectively. One can observe that the two polarimetric methods agree well with measured IWC while the IWC estimates from just using radar reflectivity exhibit biases as large as 1 order of magnitude. The large errors in the IWC–Z estimator are due to the presence of mixtures of large aggregates and small ice crystal regions as indicated in the PIP images (not shown) in clouds. Large aggregates have a dominant contribution into the radar reflectivity, which explains the positive biases of the IWC–Z estimates. On the other hand, K_dp is not biased toward large aggregates. The magnitude of K_dp in aggregates with MMD >2 mm is usually smaller than 0.4 ^∘ m⁻¹ and in small ice particles (MMD in the range 0.3–1 mm) K_dp is between 0.6 and 1 ^∘ km⁻¹ (Fig. 8a). It follows that estimators utilizing K_dp information would overcome the large aggregate effects in radar volumes. It is worth noting that the two algorithms capture the IWC variation at the end of the segment well. If the in situ measurements are considered as the ground truth, the estimation biases are computed and shown in Fig. 7b. On average, biases are 0.082 and 0.018 g m⁻³, and the root-mean-square differences (hereinafter referred to as the rms differences) are 0.49 and 0.48 g m⁻³ for the K_dp alone and (K_dp, Z_DR) methods, respectively. The correlation coefficients between IWC_PSD and estimated IWCs are 0.66 and 0.70 for the two methods. In this case study, the inclusion of Z_dr improves the accuracy of the IWC estimates.

Figure 8 K_dp (a) and Z_dr (b) as functions of median mass diameter (MMD). Over 17 000 data points from seven selected flights (Sect. 6) during the Cayenne campaign are used.

5.2 Case study II: 23 May flight

For this case, a segment of the Convair flight on 23 May 2015 inside an MCS north of the Surinam coast and over French Guiana (Fig. 9a) is selected. The flight segment consists of a HIWC region of a very high concentration of small ice particles and a region of mixed moderately large aggregates and pristine ice crystals. This affords an excellent example to gauge the performance of the algorithms. In Fig. 9a, the selected segment is displayed in purple. The radar reflectivity field from the side antenna is shown in Fig. 9b. In this segment, the aircraft's altitude was between 6.74 and 6.78 km and T_s ranged from −11 to −8 ^∘C.

Figure 9 Similar to Fig. 4 but for the 23 May case.

In addition to the radar data, IWC_PSD and MMD time series from particle-imaging probes are shown in Fig. 10. The aircraft sampled two small cores where IWC was higher than 1 g m⁻³ (∼18:30 and ∼18:34 UTC). At these high IWC cores, the clouds were dominated by small ice particles (Fig. 11a) and MMD was in the 400 µm range. In contrast, for the flight segment between 18:36 and 18:44 UTC, when the temperature was higher, the aircraft sampled a mixture of large aggregates with sizes exceeding 6 mm and small ice particles (Fig. 11b), where the IWC was less than 0.5 g m⁻³.

The IWC estimates from the methods are depicted in Fig. 12a. In the regions around the two HIWC peaks, results from the three estimators agree quite well with IWC_PSD. There are small biases in the outcomes of the two polarimetric algorithms that can be attributed to the errors of fitting linear regression models to the data and/or the difference in the sampling locations of the radar and the in situ data (Sect. 3.1). In the region after 18:38 UTC, the IWC–Z results show very large errors due to the presence of aggregates in the clouds. The large aggregates dominate the measurements of radar reflectivity, resulting in positive biases of the IWC–Z estimates (Fig. 12b). The errors for this case are as large as 300 % in most estimates. In contrast, both the polarimetric methods provide much better results compared to the conventional IWC–Z method. They capture the variation in IWC at smaller scales (around 18:33:58 UTC) and larger scales (around 18:41:54 UTC) well. This again confirms that these algorithms are robust to the variation in ice crystal type, shape, and distribution. The rms differences and correlation coefficients for the K_dp-only and (K_dp, Z_DR) methods are (0.84 g m⁻³, 0.41) and (0.79 g m⁻³, 0.55), respectively. The combination of K_dp and Z_DR provides better results, which can be seen at the edges of the second IWC peak (indicated by ellipses) in Fig. 12a. At those regions, MMD (Fig. 10) and Z_dr (not shown) values are large. This may be an indication of ice crystals with a high axis ratio aligned in the horizontal plane. When this happens, the algorithm based on K_dp alone will overestimate IWC. On the other hand, the product $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$IWC already includes the particles' shape and orientation effects; thus, estimates based on it should yield better results. When large particles dominated the volume (after 18:36:14 UTC) Z_dr become small (Fig. 8b), and then the (K_dp,Z_DR) estimator provides no advantage over the K_dp-only estimator.

6 Experimental evaluation

In the previous section, two case studies were analysed in detail. In both cases, results from the polarimetric methods show a much better agreement with in situ measurements compared to the IWC estimates from the radar reflectivity factor, especially when larger particles dominate the radar volume. In addition, applying a function of Z_DR to IWC before fitting a linear regression model to the data improves the estimation accuracy and correlation. In this section, more data from different flights collected during the mission were analysed and summarized. Out of total 14 campaign flights, there were seven flights with good data quality (radar and in situ) and with an applicable number of high-IWC data points, and data from those flights were used this analysis.

Figure 10 IWC_PSD and MMD time series for the 23 May case.

Figure 11 Samples of 2-D imagery from the SPEC 2D-S (a) and DMT PIP (b) probes at two time stamps as in Fig. 10. The width of the 2D-S image strip is 1.28 mm and that of the PIP is 6.4 mm. The aircraft's altitude was 6.75 km at A and 6.74 km at B.

Z_DR threshold (Sect. 2.2) is determined from all selected data (17 699 points in total). In order to find an optimal Z_DR threshold from the available data, average bias and rms of IWC estimates are expressed as a function of the Z_DR threshold (Fig. 13). Z_DR threshold was changed within [1.01, 1.2] with 0.1 increments, and bias and rms were computed for each value of Z_DR threshold. In Fig. 13, average bias and rms of IWC estimates from the K_dp-only algorithm, which are independent of Z_DR threshold, are also displayed (blue lines). It follows that the average biases for the two methods are very small (within ±0.08 g m⁻³) and the (K_dp, Z_DR) method provides unbiased estimates at a Z_DR threshold of 1.06. However, rms of the (K_dp,Z_DR) method is quite large at a small Z_DR threshold and reduces with an increasing Z_DR threshold. It becomes saturated at 0.498 g m⁻³, which is slightly below the rms of the K_dp-only algorithm. Considering all the factors, we selected an optimal Z_DR threshold of 1.12, where rms values of the two methods are equal but the average bias is smaller with the (K_dp,Z_DR) method.

Figure 12 Similar to Fig. 7 but for the 23 May case. Average biases for IWC(Z), IWC(K_dp), and IWC(K_dp, Z_DR) are 1, −0.189, and −0.222 g m⁻³ and rms differences are 1.34, 0.84, and 0.79 g m⁻³ for the algorithms, respectively.

Figure 13 Average bias and rms as a function of Z_DR threshold for all data from the seven selected flights.

In Fig. 14, IWC and $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$IWC are expressed as functions of K_dp for selected flights. Linear fits to all the data from selected flights are also shown. The respective y axes are scaled to the maximum value of IWC and $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$IWC for comparison. For most cases, the linear relationships are well approximated up to K_dp=2 ^∘ km⁻¹. At larger K_dp, IWC saturates at 2.5 g m⁻³ and the IWC–K_dp relationship departs from the linear trend. Due to the limited amount of data of large measured K_dp and IWC, identifying the major reasons for this saturation is not attempted. In these scenarios, applying a more sophisticated method (such as a parametric model) will likely reduce errors at high K_dp, but this is beyond the scope of this paper. Here, a simple linear regression model (based on the approximation in Eq. 13) is used and errors are computed from all data points.

Figure 14 IWC_PSD (a) and (1-$Z_{\mathrm{DR}}^{-\mathrm{1}}$)IWC (b) as functions of K_dp for the seven selected flights. Linear fits (dash lines) are also plotted using coefficients computed from all data point (Table 2).

Table 2 Polarimetric method performance for selected flights during the Cayenne 2015 campaign. The MCS types and scales from Strapp et al. (2016) are also listed.

^* For all data points, optimal fitting parameters (0.88, 0.45) were used for the K_dp-only algorithm and (0.13, 0.04) were used for the (K_dp, Z_dr) algorithm.

Figure 15 (a) Combined IWC time series data from the selected flights: measured IWC (black line), estimated IWC using K_dp alone (blue line), and estimated IWC using K_dp and Z_DR (red line). (b) Estimation errors for the two estimators. For all study cases, the aircraft flew in the range of [5.6, 7.5] km and most of the data points were within a temperature range of (−10 ^∘C ±25 ^∘C).

Figure 16 Bias and rms difference as a function of IWC derived from the seven selected flights. Mean values and standard deviation are computed from data points in each IWC bin of 0.2 g m⁻³.

It is also worth noting that the deviation of the $\left(\mathrm{1}-Z_{\mathrm{DR}}^{-\mathrm{1}}\right)$ IWC–K_dp curves from the linear fit is smaller compared to that of the original IWC–K_dp curves. This spread in the IWC–K_dp relationship can be attributed to the properties of ice crystals and the medium's state. In other words, when the dependency of the IWC–K_dp relationship on ice crystal shape and orientation was removed (or partially removed), the spread in the IWC–K_dp relationship around the linear fit should be smaller. This is a very important outcome which helps to reduce estimation errors when a single estimator is used for all the cases. Results for IWC estimates are shown in Table 2 for the two polarimetric methods only. In each row, statistical error analysis is shown for each flight with the optimal fitting model derived from data of that flight. Information about MCS of selected flights is also provided. The last row displays results computed from all selected data of 17 699 points. In all cases, improvement in IWC estimation when Z_dr information is utilized in the algorithm is clear. For all data, the bias changes from −0.07 to −0.045 g m⁻³ and correlation coefficient increases from 0.69 to 0.72. The standard deviations of the fitting coefficients (a₁,b₁) for the K_dp-only method and for (a₂,b₂) for the (K_dp, Z_DR) method are (0.12, 0.33) and (0.032, 0.033), respectively. The uncertainty of the retrieval depends on the uncertainty in the fitting parameters as well as the values of K_dp and Z_DR and their measurement accuracy. Typical values of K_dp and Z_DR for HIWC regions (MMD between 0.3 and 1 mm) are about 1 ^∘ km⁻¹ and 1.12 (Fig. 8). At those typical values, standard deviation of IWC estimates using the (K_dp,Z_DR) algorithm is 0.6 g m⁻³.

Figure 15 shows time series of IWC_PSD from the seven flights and estimated IWC from the two algorithms. As mentioned before, for each algorithm, a single set of fitting parameters is used for the combined data. Evidently, the method utilizing Z_dr yields better results in terms of estimation bias and correlation (Table 2). In Fig. 16, estimation bias and standard deviation are expressed as a function of IWC. It can be seen that inclusion of Z_dr improves estimation bias at all IWC points. On average, an improvement of 35 % in average bias was achieved. As observed in Fig. 16, larger biases happen at IWC greater than 2 g m⁻³. This is attributed to strong departures from the linear model in the joint IWC–K_dp distribution. The inclusion of Z_dr has been proven to be able to mitigate these large errors but not completely fix the problems. To improve the radar-derived IWC estimates further, more additional data processing (such as hydrometeorology classification) and more sophisticated regression models are needed.

7 Conclusions

Accurate detection and estimation of HIWC in tropical mesoscale convective systems are critical for reducing hazards caused by the ingestion of ice particles into the engines of commercial aircraft. The objective of this paper is to find a method to improve IWC retrieval from a side-looking X-band dual-polarization airborne radar. It is shown that the use of the specific differential phase (K_dp) and differential reflectivity ratio (Z_dr) significantly reduces errors in IWC retrieval over the conventional IWC–Z method. In general, the IWC–K_dp relationship can be approximated by a linear model, and IWC retrieval using K_dp captures the IWC variation very well, regardless of the information of PSD. One major drawback of the K_dp algorithm is that it provides large estimation biases when the ice particle's aspect ratio and orientation are changing. To mitigate this effect, Z_dr is used to reduce the dependency of IWC on the variation in ice particles' shapes and orientation. We proposed a method in which IWCs are weighted by a function of Z_dr before applying a linear model to the IWC–K_dp joint distribution. This approach uses an assumption of a low degree of riming within the radar volume. This is suitable for HIWC regions, which are often composed of a very high density of small ice particles. Z_dr at regions of mixtures of small pristine ice crystals and larger particles such as aggregates is generally low and will not be used in the weighting function. Results from selected Convair CV-580 flights from the Cayenne campaign show that the proposed method is able to improve estimation biases by 35 % and correlation by 4 %, on average. In our analysis, a single set of fitting parameters is applied for all the data points. The results can be improved further by including advanced data processing techniques such as ice crystal type classification or using a more sophisticated regression model for the modified IWC–K_dp joint distribution.

Most of HIWC data points used in this are measured at a narrow window of the temperature range (−10 ^∘C ±2.5 ^∘C). More data are needed to study the temperature variability of the proposed method.