the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

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

### Cuong M. Nguyen

### Mengistu Wolde

### Alexei Korolev

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.

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^{−3}), 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.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^{−1}) 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),

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)

where *n* is the number concentration in reciprocal litres, *k* is wavenumber in
reciprocal metres, *R**e*[] 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),

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

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,

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

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.

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.

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^{−1}, this results in a 100 m
radar sampling volume.

## 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.

*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^{−1}. The NRC *K*_{dp}
estimation algorithm is summarized in the flow chart below.

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.257*Z*^{0.391}) and at $T=-\mathrm{10}$ ^{∘}C (IWC =0.253*Z*^{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^{−3} (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^{−3}) (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^{−1} (and *Z*_{DR}∼1.12) where *Z*_{dr}
started increasing and *ρ*_{hv} decreased with respect to *K*_{dp}. At
*K*_{dp}<1 ^{∘} km^{−1}, *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^{−1}, 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^{−1}) 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*^{2} 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.

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^{−1} and in small ice particles
(MMD in the range 0.3–1 mm) *K*_{dp} is between 0.6 and 1 ^{∘} km^{−1} (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^{−3}, and the root-mean-square differences (hereinafter referred to as the rms differences) are
0.49 and 0.48 g m^{−3} 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.

## 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.

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^{−3} (∼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^{−3}.

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^{−3}, 0.41) and (0.79 g m^{−3}, 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.

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.

*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^{−3}) 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^{−3}, 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.

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^{−1}. At larger *K*_{dp}, IWC saturates at 2.5 g m^{−3} 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.

^{*} 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.

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^{−3} and correlation coefficient increases from 0.69 to 0.72.
The standard deviations of the fitting coefficients (*a*_{1},*b*_{1}) for the
*K*_{dp}-only method and for (*a*_{2},*b*_{2}) 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^{−1} 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^{−3}.

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^{−3}. 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.

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.

The data used in this analysis are stored at the National Research Council database. Please contact the lead author for access to the data.

CN developed the algorithms, processed radar data, and performed the data analysis. MW was the flight mission scientist for all the flights and operated the radar. MW, AK, and Ivan Heckman (ECCC) processed the in situ data. The paper was written by CN. MW and AK have provided content for sections of the paper and also reviewed the draft and made some changes to the original draft.

The authors declare that they have no conflict of interest.

This work is supported by the FAA (NAT-I-8417) and the NRC RAIR programme. We would like to acknowledge Jim Riley, Tom Bond, and Chris Dumont of FAA and Steve Harrah of NASA for their support for this work. We thank the engineering, scientific, and managerial staff from NRC and ECCC, who made the project possible by working long hours during instrument integration and field operations. Special thanks are due to Ivan Heckman (ECCC) for support of processing cloud microphysical data. We would also like to acknowledge the support we received from the HAIC-HIWC community during the flight campaign.

This paper was edited by Brian Kahn and reviewed by Alexander Ryzhkov and two anonymous referees.

Aydin, K. and Tang, C.: Relationships between IWC and polarimetric measurements at 94 and 220 GHz for hexagonal columns and plates, J. Atmos. Ocean. Tech., 14, 1055–1063, 1997.

Bringi, V. N. and Chandrasekar, V.: Polarimetric Doppler Weather Radar: Principles and Applications, Cambridge University Press, 636 pp., 2001.

Davison, C. R., MacLeod, J. D., Strapp, J. W., and Buttsworth, D. R.: Isokinetic Total Water Content Probe in a Naturally Aspirating Configuration: Initial Aerodynamic Design and Testing, 46th AIAA Aerospace Sciences Meeting and Exhibit, Reno, Nevada, AIAA-2008-0435, 2008.

Fontaine, E., Schwarzenboeck, A., Delanoë, J., Wobrock, W., Leroy, D., Dupuy, R., Gourbeyre, C., and Protat, A.: Constraining mass–diameter relations from hydrometeor images and cloud radar reflectivities in tropical continental and oceanic convective anvils, Atmos. Chem. Phys., 14, 11367–11392, https://doi.org/10.5194/acp-14-11367-2014, 2014.

Grzych, M. and Mason, J.: Weather conditions associated with jet engine power loss and damage due to ingestion of ice particles: What we've learned through 2009, 14th Conf. on Aviation, Range and Aerospace Meteorology, Atlanta, GA, Amer. Meteor. Soc., 6.8, 2010.

Hemsfield, A. J.: Precipitation development in stratiform ice clouds: A microphysical and dynamical study, J. Atmos. Sci., 34, 367–381, 1977.

Hogan, R. J., Mittermaier M. P., and Illingworth, A. J.: The retrieval of ice water content from radar reflectivity factor and temperature and its use in evaluating a mesoscale model, J. Appl. Meteorol. Clim., 45, 301–317, 2006.

Korolev, A., Heckman, I., and Wolde, M.: Observation of Phase Composition and Humidity in Oceanic Mesoscale Convective Systems, AMS Cloud Physics Conference, Vancouver, BC, 9–13 July 2018.

Lawson, R. P., Angus, L. J., and Heymsfield, A. J.: Cloud particle measurements in thunderstorm anvils and possible weather threat to aviation, J. Aircraft, 35, 113–121, 1998.

Leroy, D., Fontaine, E., Schwarzenboeck, A., and Strapp, J. W.: Ice crystal sizes in high ice water content clouds. Part I: Mass–size relationships derived from particle images and TWC for various crystal diameter definitions and impact on median mass diameter, J. Atmos. Ocean. Tech., 34, 117–136, 2016.

Lu, Y., Aydin, K., Clothiaux, E. E., and Verlinde, J.: Retrieving Cloud Ice Water Content Using Millimeter- and Centimeter-Wavelength Radar Polarimetric Observables, J. Appl. Meteorol., 54, 596–604, 2015.

Mason, J., Strapp, W., and Chow, P.: The ice particle threat to engines in flight, Proc. 44th AIAA Aerospace Sciences Meeting and Exhibit, Reno, NV, American Institute of Aeronautics and Astronautics, AIAA-2006-206, https://https://doi.org/10.2514/6.2006-206, 2006.

Matrosov, S. Y., Reinking, R. F., Kropfli, R. A., and Bartram, B. W.: Estimation of ice hydrometeor types and shapes from radar polarization measurements, J. Atmos. Ocean. Tech., 13, 85–96, 1996.

Protat, A., Delanoë, J., Strapp, J. W., Fontaine, E., Leroy, D., Schwarzenboeck, A., Lilie, L., Davison, C., Dezitter, F., Grandin, A., and Weber, M.: The measured relationship between ice water content and cloud radar reflectivity in tropical convective clouds, J. Appl. Meteorol. Clim., 55, 1707–1729, 2016.

Rotemberg, J.: A Heuristic Method for Extracting Smooth Trends from Economic Time Series, NBER Working Paper No. 7439, 1999.

Ryzhkov, A., Bukovcic, P., Murphy, A., Zhang, P., and McFarquhar, G.: Ice microphysical retrievals using polarimetric radar data, 10th European Conference on Radar in Meteorology and Hydrology, the Netherlands, 40, 1–6 July 2018.

Ryzhkov, A. V., Zrnić, D. S., and Gordon, B. A.: Polarimetric method for ice water content determination, J. Appl. Meteorol., 37, 125–134, 1998.

Stocker, T. F., Qin, D., Plattner, G.-K., Tignor, M., Allen, S. K., Boschung, J., Nauels, A., Xia, Y., Bex, V., and Midgley, P. M. (Eds.): Climate Change 2013: The Physical Science Basis. Cambridge University Press, 1535 pp., 2013.

Strapp, J.W., Korolev A., Ratvasky T., Potts R., Protat A., May P., Ackerman A., Fridlind A., Minnis P., Haggerty J., Riley J. T., Lilie L. E., and Isaac G.A.: The High Ice Water Content Study of Deep Convective Clouds:Report on Science and Technical Plan, DOT/FAA/TC-14/31, 2016.

Vivekanandan, J., Bringi, V. N., Hagen, M., and Meischner, P.: Polarimetric radar studies of atmospheric ice particles, IEEE T. Geosci. Remote, 32, 1–10, 1994.

Wang, Y. and Chandrasekar, V.: Algorithm for estimation of the specific differential phase, J. Atmos. Ocean. Tech., 26, 2569–2582, 2009.

Wolde, M. and Pazmany, A.: NRC dual-frequency airborne radar for atmospheric research, 32nd Conf. on Radar Meteorology, Albuquerque, NM, Amer. Meteor. Soc., P1R.9, 2005.

Wolde, M. and Vali G.: Polarimetric signatures from ice crystals observed at 95 GHz in winter clouds. Part I: Dependence on crystal form, J. Atmos. Sci., 58, 828–841, 2001.

Wolde, M., Nguyen, C., Korolev, A., and Bastian, M.: Characterization of the Pilot X-band radar responses to the HIWC environment during the Cayenne HAIC-HIWC 2015 Campaign, vol. 8th AIAA Atmospheric and Space Environments Conference, AIAA Aviation, 2016, AIAA 2016-4201, 2016.

^{−3}) retrieval from a dual-polarization side-looking X-band airborne radar.

*Z*

_{dr}and

*K*

_{dp}are used to mitigate the effects of ice crystal shape and orientation on the variation in IWC – specific differential phase (

*K*

_{dp}) joint distribution. Empirical analysis shows that the proposed method improves the estimation bias by 35 % and increases the correlation by 4 % on average, compared to the method using

*K*

_{dp}alone.

^{−3}) retrieval...