Evaluation of radar reflectivity factor simulations of ice crystal populations from in situ observations for the retrieval of condensed water content in tropical mesoscale convective systems

This study presents the evaluation of a technique to estimate cloud condensed water content (CWC) in tropical convection from airborne cloud radar reflectivity factors at 94 GHz and in situ measurements of particle size distributions (PSDs) and aspect ratios of ice crystal populations. The approach is to calculate from each 5 s mean PSD and flight-level reflectivity the variability of all possible solutions of m(D) relationships fulfilling the condition that the simulated radar reflectivity factor (T-matrix method) matches the measured radar reflectivity factor. For the reflectivity simulations, ice crystals were approximated as oblate spheroids, without using a priori assumptions on the mass–size relationship of ice crystals. The CWC calculations demonstrate that individual CWC values are in the range ±32 % of the retrieved average CWC value over all CWC solutions for the chosen 5 s time intervals. In addition, during the airborne field campaign performed out of Darwin in 2014, as part of the international High Altitude Ice Crystals/High Ice Water Content (HAIC/HIWC) projects, CWCs were measured independently with the new IKP-2 (isokinetic evaporator probe) instrument along with simultaneous particle imagery and radar reflectivity. Retrieved CWCs from the T-matrix radar reflectivity simulations are on average 16 % higher than the direct CWCIKP measurements. The differences between the CWCIKP and averaged retrieved CWCs are found to be primarily a function of the total number concentration of ice crystals. Consequently, a correction term is applied (as a function of total number concentration) that significantly improves the retrieved CWC. After correction, the retrieved CWCs have a median relative error with respect to measured values of only−1 %. Uncertainties in the measurements of total concentration of hydrometeors are investigated in order to calculate their contribution to the relative error of calculated CWC with respect to measured CWCIKP. It is shown that an overestimation of the concentration by about +50 % increases the relative errors of retrieved CWCs by only +29 %, while possible shattering, which impacts only the concentration of small hydrometeors, increases the relative error by about +4 %. Moreover, all cloud events with encountered graupel particles were studied and compared to events without observed graupel particles. Overall, graupel particles seem to have the largest impact on high crystal number-concentration conditions and show relative errors in retrieved CWCs that are higher than for events without graupel particles. Published by Copernicus Publications on behalf of the European Geosciences Union. 2240 E. Fontaine et al.: Evaluation of radar reflectivity factor simulations

Abstract. This study presents the evaluation of a technique to estimate cloud condensed water content (CWC) in tropical convection from airborne cloud radar reflectivity factors at 94 GHz and in situ measurements of particle size distributions (PSDs) and aspect ratios of ice crystal populations. The approach is to calculate from each 5 s mean PSD and flight-level reflectivity the variability of all possible solutions of m(D) relationships fulfilling the condition that the simulated radar reflectivity factor (T-matrix method) matches the measured radar reflectivity factor. For the reflectivity simulations, ice crystals were approximated as oblate spheroids, without using a priori assumptions on the mass-size relationship of ice crystals. The CWC calculations demonstrate that individual CWC values are in the range ±32 % of the retrieved average CWC value over all CWC solutions for the chosen 5 s time intervals. In addition, during the airborne field campaign performed out of Darwin in 2014, as part of the international High Altitude Ice Crystals/High Ice Water Content (HAIC/HIWC) projects, CWCs were measured independently with the new IKP-2 (isokinetic evaporator probe) instrument along with simultaneous particle imagery and radar reflectivity. Retrieved CWCs from the T-matrix radar reflectivity simulations are on average 16 % higher than the direct CWC IKP measurements. The differ-ences between the CWC IKP and averaged retrieved CWCs are found to be primarily a function of the total number concentration of ice crystals. Consequently, a correction term is applied (as a function of total number concentration) that significantly improves the retrieved CWC. After correction, the retrieved CWCs have a median relative error with respect to measured values of only −1 %. Uncertainties in the measurements of total concentration of hydrometeors are investigated in order to calculate their contribution to the relative error of calculated CWC with respect to measured CWC IKP . It is shown that an overestimation of the concentration by about +50 % increases the relative errors of retrieved CWCs by only +29 %, while possible shattering, which impacts only the concentration of small hydrometeors, increases the relative error by about +4 %. Moreover, all cloud events with encountered graupel particles were studied and compared to events without observed graupel particles. Overall, graupel particles seem to have the largest impact on high crystal number-concentration conditions and show relative errors in retrieved CWCs that are higher than for events without graupel particles.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Introduction
Clouds play an important role within the hydrological cycle, radiative transfer, and heat balance of the Earth. Thus, improving knowledge of ice hydrometeor properties and understanding of related processes is important for improving numerical weather forecast and global climate models, as such models use simple schemes to describe the ice hydrometeors' properties. As a consequence, significant differences in the representation of ice properties in clouds (Li et al., 2005(Li et al., , 2007 lead to large variations in the quantification of ice cloud effects on climate evolution (Intergovernmental Panel on Climate Change Fourth Assessment Report). Among different ice properties, the spatiotemporal distribution of cloud condensed water content (CWC) is a key parameter for evaluating and improving numerical weather prediction (Stephens et al., 2002). Increasingly, remote sensing tools are used to study cloud properties such as hydrometeor size distributions (ice or water), liquid and/or ice water content, ice particle shape (spherical, hexagonal, etc.), and precipitation rate. In particular, the radar (at different frequencies of 5. 5, 9.4, 35, 94 GHz) is the most common measurement technique used to measure clouds properties. Radar reflectivity factors are an integral value of all backscattering cross sections from all hydrometeors within the radar sampling volume, which makes the radar a complex measurement device for estimating cloud properties.
The methodology applied in this study to simulate 94 GHz radar reflectivity factors is based on assumptions about individual cloud particle properties. If hydrometeors are droplets, then Mie solutions can be applied to the Maxwell equations; however, this is not the case for non-spherical ice crystals. Discrete dipole approximations (DDA; Draine and Flatau, 1994;Liu, 2008) can be used to calculate backscatter cross sections for complex shapes and thus tackle this question for ice crystals. However, in order to apply DDA simulations to ice crystal radar reflectivity factors, a classification of ice hydrometeor habits is essential. Unfortunately, more than 50 % of ice crystal images classified using automated or manual shape recognition techniques from ground-based and airborne measurements are typically identified as irregular (i.e. they are not identified as a specific habit). This statement still holds when using very-high-resolution imaging such as from the Cloud Particle Imager (2.3 µm resolution; e.g. Mioche, 2010). Hogan et al. (2011) used the oblate spheroid approximation to simulate radar reflectivity factors at 3 and 94 GHz. For their calculation of the ice fraction in horizontally oriented oblate spheroids, these authors used a constant axis ratio of 0.6 and used the mass-size relationship from Brown and Francis (1995), originally from Locatelli and Hobbs (1974), to derive CWCs and radar reflectivity factors. Even though a mass-size relationship for ice crystals with constant coefficients has been utilized, the study of Hogan et al. (2011) claimed minimal errors between measured and simulated radar reflectivity factors, even smaller than the calibration uncertainty of the cloud radar. Fontaine et al. (2014) also used the oblate spheroid approximation for ice crystals to calculate CWC in mesoscale convective systems (MCS), and Drigeard et al. (2015) used the Fontaine et al. (2014) results to simulate radar reflectivity factors at 5.5 GHz. Although simulations of radar reflectivity factors agreed with reflectivity observations at 94 and 5.5 GHz, no direct bulk measurements of CWC were available to evaluate the retrieved CWCs.
A simpler way to calculate CWC from radar reflectivity factors is based on empirical Z-CWC or Z-CWC-T relationships. Such relationships have been established in earlier studies, either with or without direct measurements of CWC (Protat et al., 2016(Protat et al., , 2007Hogan et al., 2006), for different types of clouds and at different geographical locations (tropics, continental, mid-latitude). When no direct simultaneous measurements of CWC are available, Z-CWC (and Z-CWC-T ) relationships are established by using constant mass-size relationships (for CWC calculations from particle size distribution (PSD) measurements) and most of the time these studies use m(D) coefficients suggested by Brown and Francis (1995). As a matter of fact this would mean that mass-size coefficients are constant and are linked neither to temperature nor to the variability of PSDs or the type of clouds. In Fontaine et al. (2014), these simplistic assumptions were not used. This more sophisticated method allows for a more dynamic solution with varying mass-size coefficients related to the variability of microphysics (PSD and aspect ratio of ice crystals). Fontaine et al. (2014) end up with the retrieval of an average CWC calculated as a function of time from multiple possible m(D) solutions for the T-matrix simulations of the radar reflectivity factor (Ze) simulations compared to measured Z. The method has been established and published without validation from simultaneous direct measurements of CWC.
During the High Altitude Ice Crystals (HAIC; Dezitter et al., 2013)/High Ice Water Content (HIWC; Strapp et al., 2016) airborne field campaign performed out of Darwin, the Falcon-20 (F-20) from SAFIRE (Service des Avions Français Instrumentés pour la Recherche en Environnement) measured on the same aircraft 94 GHz Z, PSD, aspect ratios of hydrometeors, and independently CWC. The objective of our study is to use this comprehensive HAIC/HIWC dataset to evaluate the method presented by Fontaine et al. (2014).
The next section of this paper presents the dataset of the first HAIC/HIWC airborne campaign and associated data processing. The principle of the cloud radar reflectivity simulation method  is then briefly recalled. The third section is dedicated to the evaluation of the Fontaine et al. (2014) method by comparing averaged in retrieved CWCs from the T-matrix simulations with direct measurements of CWC IKP . Moreover, the study suggests correction functions for calculated CWC (based on T-matrix simulations of the reflectivity factors) as a function of tem-Atmos. Meas. Tech., 10, 2239-2252, 2017 www.atmos-meas-tech.net/10/2239/2017/ perature, largest size of hydrometeors in PSDs, and ice crystal concentrations. Then, two further sections are dedicated to the estimation of possible uncertainties in the measurements leading to the corrected CWC retrievals and the impact of graupel particles on those retrievals. Finally, the study ends with a discussion and conclusion section.

Data processing
The HAIC/HIWC projects were designed to investigate the microphysical processes responsible for engine damage observed when commercial aircraft divert around convective cores. First HIWC studies (Mason et al., 2006) indicated that this new form of icing (now referred to as ice crystal icing) was due to the production of high concentrations of small ice crystals by deep convection. In this context 23 flights were performed over the Darwin area, mainly over the ocean north of Australia, during the monsoon season. Three of the 23 flights were dedicated to the calibration of the instruments and are not included in this study. The dataset includes more than 17 000 data points where all parameters are synchronized and averaged over 5 s. As one of the priority of the HAIC/HIWC projects was to measure high ice water content and the variability of CWC as a function of distance from the convective cores at typical altitudes flown by commercial aircraft, the flight strategy was to fly long legs at constant altitude at −50, −40, −30, and −10 • C and to get as close as possible to the most convective zone of the MCS (see Fig. 2 in Leroy et al., 2017), thereby avoiding red aircraft radar echoes at normal gain as commercial aircraft would also do. More details on HAIC/HIWC projects and implemented flight strategy can be found in Dezitter et al. (2013), Leroy et al. (2016Leroy et al. ( , 2017, Protat et al. (2016), and . In situ microphysical measurements and radar reflectivity data used in this study were provided by three types of instruments, which were mounted on board the F-20 for the HAIC/HIWC Darwin campaign.
(1) The 94 GHz multi-beam Doppler cloud radar RASTA measured both cloud radar reflectivity factors and 3-D cloud dynamics below and above the flight level. The uncertainty on the measured reflectivity is about ±1 dBZ (e.g. Protat et al., 2009). We only use radar vertical antennas (zenith, nadir) producing vertical profiles of the radar reflectivity along the flight trajectory. Radar reflectivity factors are interpolated at the flight altitude using validated gates (typically 180 m above and below the aircraft) to retrieve the most likely radar reflectivity at flight altitude. RASTA has a vertical resolution of 60 m and 0.7 • beam width.
(3) An Isokinetic Evaporator Probe (IKP-2: ) that provides direct CWC measurements was also used. The IKP-2 is a second-generation version of the prototype IKP (Davison et al., 2008), which was downsized for the F-20. The IKP-2 was developed to provide reliable measurements of CWC in deep convective clouds at temperatures colder than −10 • C, up to at least 10 g m −3 at aircraft speeds of 200 ms −1 , and with a target accuracy of 20 %. The IKP-2 samples the cloud particles isokinetically, evaporates them, and measures the resulting humidity of the evaporated particles and background air. CWC (hereafter CWC IKP given in g m −3 ) is then obtained by subtracting the water vapour background measurement from the IKP-2 total hygrometer signal. System accuracy estimates are better than 20 % for CWC greater than 0.25 g m −3 for temperatures lower than −10 • C . Accuracy increases with decreasing temperature due to the exponential decrease in background humidity, which drives much of the IKP-2 error. For example, at −56 • C, system accuracy was estimated at better than 4 % for CWC greater than about 0.1 g m −3 .
The 2-D-S and the PIP record monochromatic 2-D shadow images of cloud hydrometeors (ice and/or water) along the flight trajectories. The 2-D-S records images of hydrometeors in the size range 10-1280 µm at a 10 µm pixel resolution, whereas the PIP records images in the size range 100-6400 µm and beyond (perpendicular to photodiode array and reconstruction of truncated images parallel to the array) at a 100 µm pixel resolution. For both probes, PSDs were produced by image analysis into number concentrations per unit volume of sampled air as a function of their size.
In this paper, the size of ice hydrometeors is given in terms of the maximum diameter D max (e.g. see Leroy et al., 2016, for definition). The size of truncated images and sampling volume are corrected using the method presented in Korolev and Sussman (2000). This reconstruction method allows extrapolating hydrometeor sizes to a maximum size of 2.56 mm for 2-D-S and to 12.8 mm for PIP.
In addition, many artefacts can bias PSD estimates from 2-D image analysis. Therefore, supplementary post-processing is needed to retain only the natural ice particles. One of the most important causes of artefacts is the shattering of hydrometeors on the tips of OAPs. During the first HAIC/HIWC field campaign in Darwin, the newest antishattering tips were used for 2-D-S and PIP to reduce the shattering from large ice crystals. In addition, analysis of the time-dependent (along the flight trajectory) interarrival time spectra was performed to determine the cut-off time, which separates natural hydrometeors images from artefact particles (Korolev and Isaac, 2005;Field et al., 2006;Korolev and Field, 2015). It has been shown that both mitigation techniques are needed to maximize the removal of shattering artefacts (Jackson et al., 2014). Furthermore, OAP images (for both 2-D-S and PIP) of splashed hydrometeors were removed using the ratio between their projected surface area and the surface defined by the box Lx × Ly (e.g. see Dx and Dy in Leroy et al., 2016), where Lx is the projection of the size of each hydrometeor along the flight trajectory, and Ly is the projection along the array of diodes. Images with an area ratio less than 0.25 were considered as splashed particles and removed. This threshold has been calculated statistically and allows the removal of larger splashed ice crystals as well.
Another important correction is related to the sizing of outof-focus hydrometeors. In our study, the size of out-of-focus particles was corrected using the method presented in Korolev (2007). In addition, noisy pixels (satellite pixels) which may affect hydrometeor images were eliminated thereby applying the method described in Lawson (2011). These singlepixel noisy pixels are not firmly attached to the hydrometeors' images and therefore must be removed to get the best estimation of the true diameter (e.g. D max ).
Finally, high number concentrations of ice particles lead to gaps in the sampling times of OAPs due to insufficient time to record all hydrometeors images. This probe effect (also called OAP overload or dead time) was taken into account and would otherwise lead to an underestimation of the number concentration of hydrometeors. While the 2-D-S probe overload times are directly registered, the PIP overload is estimated by comparing the number of images in the PIP files to the separately registered total particle counts of particles that passed through the laser beam. The ratio of counted particles (1-D information) to recorded particle images (2-D information) is used to correct for the concentration. During an OAP overload, 1-D counted particles may reach 15000 while only about 10 000 have a recorded image, which would result in an uncertainty of 50 % on the concentration of hydrometeors, without overload correction (Fontaine, 2014). Further details on post-processing of 2-D-S and PIP data are given in Leroy et al. (2016). The individual 2-D-S and PIP PSDs were merged into a composite PSD using the algorithm described in Eq. (1). The resolution of the composite PSD is 10 µm (by interpolating the PIP raw PSDs at the 2-D-S resolution), and PSDs are averaged over 5 s time intervals for improved large particle statistics. The transition zone for changing from 2-D-S to PIP data in the composite spectrum (see equation below) is from a D max of 805 µm (median diameter for a size bin) to 1205 µm. N (D max ) is given per litre. where 3 Retrievals of CWC from radar reflectivity simulations

Simulations of radar reflectivity factors: Ze
This section reviews the method used in Fontaine et al. (2014). The principle of the technique is to retrieve CWC from simulations of radar reflectivity factors (Ze), calculated from OAP image information and compared to measured RASTA radar data. The data were collected in tropical MCSs that formed over the African continent and over the Indian Ocean during two aircraft campaigns dedicated to the Megha-Tropiques project (Roca et al., 2015). The main drawback of these experiments is that they do not include direct measurement of CWC. The Fontaine et al. (2014) method uses oblate spheroids (Hogan et al., 2011) to approximate the backscatter cross section (Q back ) of natural hydrometeors. In Eq.
(2) below, Ze is defined in mm 6 m −3 : with where and where Note that in Eq.
(2) λ is the emitted wavelength of the radar and K w-ref is the dielectric constant of liquid water at the same frequency. Q back is a function of the ice fraction (f ice ; Atmos. Meas. Tech., 10, 2239-2252, 2017 www.atmos-meas-tech.net/10/2239/2017/ Eq. 3) in the spheroid and the axis ratio of the oblate spheroid (here denoted As, see Eq. 5). f ice is thereby a function of the mass-size relationship (see Eq. 3). Equation (3) also limits the mass of an ice hydrometeor to the mass of an ice sphere of diameter D max . f ice also allows the calculation of the dielectric properties of the ice spheroids (Maxwell Garnet, 1904;Drigeard et al., 2015). As (Eq. 5) is the average aspect ratio of all hydrometeors and is calculated every 5 s according to N (D max ) and Z. Pi (Eq. 6) is a weighting function that has been defined to account for the volume occupied by the hydrometeors in the sampled volume. As(D max ) for a particle is defined as the ratio of the width (radius perpendicular to D max ) divided by D max . In Fontaine et al. (2014) the calculation of As only considers all hydrometeors with sizes D max ≤ 2005 µ which contribute on average to 95 % of Ze (Drigeard et al., 2015). Since the processing of 2-D-S and PIP has been further improved by Leroy et al. (2016), we decided to consider all hydrometeors from 15 µm to 1.2845 cm for As calculation. Note that a comparison of the As calculation utilized in Fontaine et al. (2014) with the new As calculation results in a decrease of As of less than 3 %, which has a negligible impact on the retrieval results.
Note that the constant aspect ratio of 0.6 used in Hogan et al. (2011) is rather close to the peak of the As frequency distribution presented for African monsoon MCS and also oceanic MCS over the Indian Ocean . The average aspect ratio calculated for the HAIC/HIWC dataset is 0.55, which is then very similar to respective average values for the various datasets sampled over the African continent, United Kingdom, Indian Ocean, and north of Australia.
In Fontaine et al. (2014) the exponent β of the m(D) power law relationship has been constrained as a function of time from the ice particle images. For this study, no a-priori assumptions on the mass-size relationship of hydrometeors have been chosen and therefore a variational approach has been applied to calculate CWC from Ze reflectivity factor simulations. For a given but variable exponent β i the corresponding pre-factor α i is calculated to match the measured reflectivity Z with the simulated Ze. β i is varying stepwise between 1 to 3 by increments of 0.01. Thereby, the range of potential solutions are explored using the oblate spheroids approximation with the T-matrix method (Mishchenko et al., 1996) calculating Q back of each spheroid. Hence, for a given 5 s data point, 201 calculations of α i and 201 calculations of corresponding CWC(β i ) are performed (see Eq. 10): For each 5 s data point, from the 201 possible CWC(β i ) values, an average value CWC is deduced (Eq. 11): where N tot ≤ 201, since the minimum value allowed for α i is the mass of an empty sphere (air density). Figure 1 shows two examples (flights 9 and 12) with all possible CWC(β i ) retrievals (colour band), average CWC (red line), and CWC IKP measured by the IKP (overlaid black line). The example in Fig. 1a shows results from flight 9, in which the F-20 research aircraft flew in the more stratiform part of the cloud system (w ∼ 0 m s −1 ), whereas results from flight 12 shown in Fig. 1b represent a case with more signatures of convective updrafts. Overall, Fig. 1 demonstrates that the variational retrieval method produces a large variability of possible CWC(β i ) for each 5 s data point. In general, the average CWC (red line) is close to CWC IKP . The bandwidth of all possible solutions CWC(β i ) as a function of time is calculated from the difference CWC = max(CWC(β i ))-min(CWC(β i )) between the maximum and the minimum values of CWC(β i ). On average, it is found that CWC accounts for 61 % of CWC (with 64 % for the median relative error) for the entire dataset (20 flights performed over Darwin area). Finally, the calculations reveal that 80 % of the HAIC dataset satisfy the condition CWC IKP = CWC ± 32 %, where no a priori assumptions on mass-size relationships were applied and β i linearly varies between 1 and 3, thereby producing equally eligible solutions CWC(β i ) that are finally averaged to produce a 5 s data point for CWC.
In general, and for each given 5 s data point, maximum CWC is obtained for β = 1 and minimum CWC for β =3. For β = 1, ice hydrometeors below D max = 200 µm (sometimes even below 300 µm) may reach the maximum density of 0.917 g cm −3 , while for β = 3 the density of oblate spheroids is constant as a function of D max (see Eqs. 3-4), where the density of icy spheroids is equal to 0.917 × f ice ).
The impact of β on Ze and on the retrieved CWCs is illustrated in Fig. 2. One can notice that for different values of β (β = 1, 2, 3) the corresponding value of α can be found such that the cumulative sum of Ze as a function of D max normalized by the measured radar reflectivity factor Z (in mm 6 m −3 ; Fig. 2a.) is equal to 1. The respective cumulative mass of ice crystals (as a function of D max ) then is normalized by CWC (Fig. 2b). This CWC ratio may deviate from 1, whereas the normalized cumulative sum of Ze has been equal to 1, independently of chosen β. For β =1, Ze is reached sooner (D max ≈ 1700 µm) than for β = 3 (D max ≈ 3000 µm). The likely explanation is that with increasing β, the backscattered energy is increased for large hydrometeors and the mass contribution of smaller hydrometeors is considerably reduced since the contribution of numerous smaller hydrometeors (compared to larger hydrometeors) on retrieved CWC is decreasing with increasing β.

CWC deviations from T-matrix simulations of reflectivity with respect to IKP direct measurements
This section focuses on the potential error in CWC retrievals from T-matrix simulations of radar reflectivity factors (at frequency of 94 GHz) for populations of ice hydrometeors approximated with oblate spheroids. Therefore, the relative errors of retrieved CWC with respect to reference CWC IKP (measured by the IKP-2 probe) are calculated and then analysed as a function of microphysical properties of ice hydrometeors such as total concentration (N T ; Fig. 3), temperature T (Fig. 4), maximum size of hydrometeors in PSDs (max(D max ); Fig. 5), total cloud water content CWC IKP (Fig. 6), and radar reflectivity factors Z (Fig. 7). Blue lines in Figs Fig. 3 (blue line) it appears that CWC resulting from T-matrix simulations approximating ice hydrometeors with oblate spheroids is poorer at the lower (N T < 100 L −1 ) and higher (N T > 5000 L −1 ) ranges of number concentrations. In particular, the reference CWC IKP is underestimated for small N T and overestimated for larger N T . Furthermore, Fig. 4 (blue line) seems to illustrate that this method increasingly overestimates with decreasing temperature (blue line in Fig. 4). For example, CWC exceeds CWC IKP at 220 K (±5 K) by about 25 %. Finally, the relative errors of CWC with respect to CWC IKP slightly but continuously increase with the maximum size of hydrometeors within the respective data point (blue line in Fig. 5), where the relative error of CWC reaches +25 % when max(D max ) = 1 cm versus ∼ 0 % for max(D max ) = 800 µm.

Correction functions for CWC retrievals from T-matrix simulations of reflectivity
Because of the above findings, three different types of corrections are studied in order to (i) quantify the limitations of the oblate spheroid approximation and (ii) suggest suitable correction functions that use in situ measured quantities over the entire dataset with CWC IKP larger than 0.1 g m −3 . These corrections are performed using the inverse of the original relative errors ( Without the above correction functions, retrieved initial CWCs are larger than CWC IKP by about 19 % on average (with a median value of +16 %; Table 1, first row). Therefore, all three correction functions (Eqs. 9-11) have a median factor of 0.84 in common that reduces the initial CWC such that CWC × f (X; X = N T , T , max(D max )) better matches CWC IKP . The expressions in parentheses of Eqs. (9)-(11) try to redistribute the relative error in CWC from T-matrix simulations over the entire range of observed N T , T , and max(D max ) values, but they have negligible impact on the median relative error itself. Even though no correction functions for CWC have been proposed as a function of CWC IKP and Z, Figs. 6 and 7 illustrate the impact of N T , T , and max(D max ) correction functions (Eqs. 9-11) on the redistribution of the relative error also as a function of CWC IKP and Z. Figure 3 reveals that f (N T ) (Eq. 9) decreases biases of retrieved CWC × f (N T ) over the entire N T bandwidth, while f (T ) (Eq. 10) and f (max(D max )) (Eq. 11) do not change the shapes of the relative error lines as compared to uncorrected CWC relative errors (Fig. 3). Also, the function  f (N T ) (Eq. 9) also generally decreases relative errors of the CWC × f (N T ) retrievals when plotted as a function of T (Fig. 4) and as a function of max(D max ) (Fig. 6). Furthermore, the f (T ) correction function (Eq. 10) reduces the differences between CWC×f (T ) and CWC IKP as a function of T (Fig. 4). However, CWCs corrected as a function of T still show bias when presented as a function of N T (Fig. 3, grey line) or as a function of Z (Fig. 7, grey line). Finally, the f (max(D max )) correction function (Eq. 11) reduces the relative error of retrieved CWC×f (max(D max )) as a function of max(D max ) in Fig. 5 but does not have much impact on the shape of the relative error distributions as a function of N T , T , CWC IKP , and Z (Figs. 3, 4, 5, 7) relative to uncorrected CWC relative errors. In addition, rows 2-4 of Table 1 present mean (average), median, 10th, 25th, 75th, and 90th percentiles of the relative error after applying the correction functions, with an obvious decrease of mean and median values and corresponding shift of relative error distribution percentiles.
Overall, the f (N T ) correction seems most efficient to remove the CWC bias. Heymsfield et al. (2013) showed that to- tal concentrations of ice hydrometeors tend to increase with decreasing temperature in tropical MCS for temperatures −60 • C < T < 0 • C. This evolution of the increasing total concentration of hydrometeors related to decreasing temperatures is therefore suggested as the key to explaining trends in relative CWC errors as a function of N T and T . Figure 8 summarizes the above results by showing probability distribution functions of CWC × f (X; X = N T , T , max(D max )) versus CWC IKP . Imperfections of the corrections described by Eqs. (10) and (11) are clearly visible in Fig. 8c and d, where high CWC IKP values are still overestimated (as in Fig. 8a) by CWC × f (T ) and CWC × f (max(D max )), respectively. The correction function f (N T ) (Eq. 9) produces the best results (Fig. 8b), where the maximum of the probability distribution function in CWC × f (N T ) versus CWC IKP representation follows the line y = x.

Uncertainties in ice particle concentrations and impact on results
This section investigates the impact of uncertainties in crystal concentrations on the CWC retrieval, with a particular focus on shattering of ice crystals. As discussed in Sect. 3.2 the relative errors increase with total number concentration, with overestimations by about 50 % of CWCs for very large concentrations of hydrometeors, which can reach 10 4 hydrometeors per litre in most convective parts of sampled MCS (Fig. 3a). In order to investigate the impact of uncertainties of number concentrations on the retrieved CWC we apply two different types of functions on the measured PSDs, in which both functions increase the number concentrations of measured PSDs.
First, a function f shatt is applied to the PSD in order to increase concentrations of hydrometeors in the first PSD size bin (5-15 µm) by about 50 %, while concentrations of hydrometeors larger than 500 µ remain unchanged. The function f shatt decreases in a logarithmic way with D max from first bin to 500 µm. f shatt is applied to PSD such that N (D max ) = f shatt (D max ) × N (D max ) and aims to produce new PSDs where the optimized probe tips still would have produced shattered crystal fragments and/or removal processing would have failed to remove numerous shattered ice particles. Then, the retrieval method (see Sect. 3.1) is applied to these new PSDs in order to calculate new values for α i and subsequently CWC i (hereafter α i,fshatt and CWC i,fshatt ). For the purpose of this Sect. 4, the method was only applied for β = [1, 2, 3] in order to get a good idea of the maximum impact of possible shattering artefacts. Results are presented in terms of relative errors in Table 2 for α i,fshatt and Table 3 for CWC i,fshatt , respectively. Relative errors in (%) are calculated with respect to coefficients α i and CWC i (for β = [1, 2, 3]) calculated in Sect. 3 for non-modified original N (D max ) size distributions and without correction functions discussed in Sect. 3.3. (Eqs. 9 to 11). The relative errors illustrate that the chosen concentration increase of solely small hydrometeor sizes has very limited impact on retrieved α i and CWC i . Indeed, we observe that the median relative error of the prefactor α i,fshatt with respect to α i is roughly −3 % for β = 1 and 0 % for β = 2 and β = 3. The 1st and 99th percentiles are shown in order to demonstrate that the relative errors in α i,fshatt are small over the entire dataset. Consecutively, median relative errors of CWC i,fshatt with respect to CWC i are of the order of 4, 3, and 1 % for β = 1, 2, and 3, respectively. The 99th percentile of the relative error does not exceed 10 % in retrieved CWC.
Second, a simple concentration uncertainty factor of 1.5 is applied over the entire measured size range, which increases the number concentration by 50 % such that N (D max ) = (1.5 × N (D max )). Note that 50 % is approximatively the missed number of ice crystal images by the PIP due to the probe overload in high concentrations of ice crystals, though data have been corrected for overload times (see Sect. 2 and Fontaine, 2014). Simulations of the reflectivity factor with modified N (D max ) were performed with resulting prefactor α i,50 % and derived CWC i,50 % . Results of the comparison of α i,50 % and CWC i,50 % with α i and CWC i , respectively, are also presented in Table 2 and Table 3. Globally, this second scenario of concentration enhancement of original PSD has a larger impact on retrieved prefactor (α i ) and calculated CWC i than was the case for the first scenario. Indeed, adding 50 % to the concentrations of hydrometeors results in a median decrease of prefactor α i,50 % with respect to α i of −18 % for β = 1 and −20 % for both β = 2 and 3. The forced decrease in α (from α i,50 % to α i, ) goes along with an increase of CWC i,50 % with respect to CWC i by about +29 % (for β =1) and +27 % for β = 2 and 3. Indeed, for a given radar reflectivity factor we simulate (and measure) the same Ze (and Z) for N (D max ) (and N (D max )) with decreasing α (α i,50 % compared to α i, ), while the CWC i,50 % increases by almost 30 % (CWC i,50 % case compared to CWC i ). Hence, if two different size distributions produce an identical Ze, CWC can be significantly different. In other words, Ze may differ significantly for the same CWC associated with two underlying different size distributions even without considering uncertainties in concentration measurements.

Impact of graupel on retrieved CWCs
Graupel are more spherical than pristine ice crystal shapes and aggregates. When approximating graupel particles as oblate spheroids in the T-matrix calculations, their density is close to solid ice density and their aspect ratio is close to 1. Hence, observations of graupel particles dominating the crystal populations are investigated within the entire dataset to estimate their impact on retrieved CWCs. A detection algorithm of virtually spherical particles has been applied to PIP image data to select events where graupel are observed. Most of these events are observed for N T > 2000 L −1 , CWC > 1 g m −3 . For those events, retrieved CWCs also have relative errors larger than 100 %. Some events are also detected for N T < 2000 L −1 , but with relative errors less than 30 %. Table 1 gives relative errors for graupel events only (bracketed bold numbers). It shows mean and median relative errors of retrieved CWCs without corrections and with correction applied from Eqs. (9) to (11). The correction as a function of N T (Eq. 9) has the largest impact on the relative errors (mean and median) for graupel events. Indeed, mean relative errors of retrieved CWCs for all graupel events are 83 % (without correction) and 16 % (after f (N T ) correction). Likewise, the percentages for the median relative errors are 75 and 7 %, respectively, without and with correction. Note that corrections as functions of max(D max ) and T are less efficient in reducing mean and median relative errors of retrieved CWCs for graupel events. In the following, three different microphysical categories are distinguished in order to compare PSD and As(D max ) of observed graupel events with respective non-graupel data: 1. CWC IKP > 2 g m −3 and N T > 2000 L −1 (Fig. 9a and b), 2. 1 < CWC IKP < 2 g m −3 and N T > 2000 L −1 (Fig. 9c and d), and 3. 1 < CWC IKP < 2 g m −3 and N T < 2000 L −1 (Fig. 9e and f) Figure 9 shows median, 25th, and 75th percentiles of PSDs (Fig. 9a, c, and e) and size-dependent As(D max ) (Fig. 9b, d, and f). Graupel and non-graupel events are shown as red and blue lines, respectively. From this figure, there is no evidence of significant differences between median PSDs with and without detection of graupel particles (three left figures). However, the median aspect ratio substantially changes for hydrometeors larger than 1mm, with increased As(D max ) for the graupel particle events compared to non-graupel events. For all three graupel categories defined above, the number of observations with graupel particles are smaller than those without graupel particles. Indeed, for the event category where CWC IKP > 2 g m −3 and N T > 2000 L −1 a total of 451 events without graupel particles and solely 130 events with graupel particles are detected, with median relative errors in uncorrected retrieved CWCs of +40 and +74 %, respec-tively. For the event category with 1 < CWC IKP ≤ 2 g m −3 and N T > 2000 L −1 , there are 1687 and 80 graupel and nongraupel events, respectively, with median relative error of uncorrected retrieved CWCs of +25 and +63 %, respectively. Finally, for the third category (1 < CWC IKP < 2 g m −3 and N T < 2000 L −1 ) only 16 events are detected with graupel particles compared to 1701 events without graupel particles, with corresponding median relative errors of +10 and +9 %, respectively. Graupel particles seem to have smaller impact on this third category as compared to two other categories. Most of the graupel events occur for high crystal number concentrations (N T > 2000 L −1 ), thereby leading to higher median relative errors in uncorrected retrieved CWCs as compared to the non-graupel events.

Discussion and conclusions
The objective of the present study was to evaluate with direct CWC observations the CWC retrieval technique outlined in Fontaine et al. (2014), based on matching simulated radar reflectivity factors from PSD measurements with measured radar reflectivity factors at 94 GHz. Since mass-size relationships are considered to be unknown, the reflectivity simulations explore a wide range of possible solutions of (α i , β i ), varying β in the range of 1 to 3. This produces a series of possible CWCs for a given data point, one for each value of β. From this series, the average value CWC over all values of β is calculated. On average it is found that the difference CWC(β = 1) − CWC(β = 3) is approximately 64 % of the average CWC. Of the Darwin data points, 77 % meet the condition of CWC(β = 3) ≤ CWC IKP ≤ CWC(β = 1), which goes along with the relation CWC IKP = CWC ±32 %. However, the retrieved CWC values are generally larger than CWC IKP by about 16 % (median value), which illustrates that the approximation of ice oblate spheroids tends to underestimate the backscattered energy of real reflectivity measurements of ice hydrometeors at 94 GHz, and a constant factor of 0.84 could be applied as a first-order correction of retrieved CWC.
One of the possible explanations is that the calculation of the average aspect ratio from 2-D images (Eqs. 5 and 6) which has been adopted as the flattening parameter utilized for the approximation of ice oblate spheroids might be somewhat too large. A way to investigate impact and calculation of axis ratio for ice oblate spheroids approximations would be to use As(D max ) instead of As (Eq. 5) in simulations of radar reflectivity factors.
This study also demonstrated that the total concentration of ice hydrometeors could be used as input for a correction algorithm that minimizes differences between CWC and CWC IKP and that this parameter was the best of several parameters evaluated for this purpose. These differences, before correction, were found to increase with increasing ice concentration, with CWC underestimating CWC IKP at low ice concentrations and overestimating CWC IKP at high concentrations.
Another attempt of explanation is the uncertainty of measured total crystal concentrations, which could partly explain large relative errors at high concentrations. Indeed, for very high concentrations we find a higher relative error of about +50 % as compared to IKP measurements. However, according to the results of Sect. 4 of this study, an uncertainty of 50 % in ice crystals concentrations can solely explain 30 % of the relative errors. Likewise, concentration errors related to large crystal shattering could not explain more than 35 % of the relative errors at very high concentrations of ice hydrometeors. However, this cannot explain the negative relative errors for low ice crystal concentrations. Also, at very low ice crystals concentrations oblate spheroids approximation of crystals could be not sufficiently adapted, since for low concentrations real shapes of ice hydrometeors might be even more important due to the lack of any averaging process over all possible shapes and possible orientations as should be more likely the case for higher concentrations.
Moreover, graupel particles are found to substantially affect the relative errors in retrieved CWCs in high total concentration situations as compared to moderate total crystal concentration situations with graupel. In the latter observed conditions, graupel may have been less numerous and their impact on the aspect ratio of ice hydrometeors larger than 1 mm becomes less important than at higher concentrations (N T > 2000). Despite the fact that most graupel events are found during high crystal concentration conditions, the correction function f (N T ) appears to be very efficient in reducing the relative error in retrieved CWCs when graupels are detected. Since the suggested f (N T ) correction has been established for the entire dataset and is not exclusively designed for the graupel event correction of retrieved CWCs, median relative errors of f (N T ) corrected CWCs of all graupel data are of the order of 7 % compared to −1 % for the entire dataset (see Table 1).
This study does not focus on the mass-size relationship directly, but it is clear that the coefficients of the m(D) relationship, and particularly its exponent, considerably impact the simulation of radar reflectivity factors on the one hand and CWC calculation on the other hand. In this context, we recall that Fontaine et al. (2014) obtain different CWCs for one and the same radar reflectivity factor, which is also valid reciprocally, where one and the same CWC is related to a range of radar reflectivity factors, thereby varying the β exponent and corresponding prefactors α in the T-matrix simulations of the radar reflectivity factors. Hence, getting large differences between calculated CWC (related to radar reflectivity simulations) and measured CWC from IKP-2 does not mean that the oblate spheroid approximation for Ze is wrong, but the difference stems primarily from the choice of β (and respective α) of the mass-size relationship. We are aware that a single m(D) power law is not sufficient to assimilate all the complexity of the mass of ice hydrometeors. The main goal Atmos. Meas. Tech., 10, 2239-2252, 2017 www.atmos-meas-tech.net/10/2239/2017/ of this paper is to evaluate the Fontaine et al. (2014) method of calculating an average CWC from all possible solutions for (α, β) without a priori assumption. The method has been already used in further studies: Drigeard et al. (2015) and Alcoba et al. (2015). Finally, the method allows constraining the effective density of ice hydrometeors and also the simulation of radar reflectivity factors at different frequencies: 94, 9.4 and 5.5 GHz. The method will be used in the near future to derive mass-size relationships for different size ranges of the ice hydrometeors' size spectrum.
Data availability. No data sets were used in this article.