Microphysical processes producing high ice water contents (HIWCs) in tropical convective clouds during the HAIC-HIWC field campaign: evaluation of simulations using bulk microphysical schemes

Regions with high ice water content (HIWC), composed of mainly small ice crystals, frequently occur over convective clouds in the tropics. Such regions can have median mass diameters (MMDs) <300 µm and equivalent radar reflectivities <20 dBZ. To explore formation mechanisms for these HIWCs, high-resolution simulations of tropical convective clouds observed on 26 May 2015 during the High Altitude Ice Crystals – High Ice Water Content (HAIC-HIWC) international field campaign based out of Cayenne, French Guiana, are conducted using the Weather Research and Forecasting (WRF) model with four different bulk microphysics schemes: the WRF single‐moment 6‐class microphysics scheme (WSM6), the Morrison scheme, and the Predicted Particle Properties (P3) scheme with one- and two-ice options. The simulations are evaluated against data from airborne radar and multiple cloud microphysics probes installed on the French Falcon 20 and Canadian National Research Council (NRC) Convair 580 sampling clouds at different heights. WRF simulations with different microphysics schemes generally reproduce the vertical profiles of temperature, dew-point temperature, and winds during this event compared with radiosonde data, and the coverage and evolution of this tropical convective system compared to satellite retrievals. All of the simulations overestimate the intensity and spatial extent of radar reflectivity by over 30 % above the melting layer compared to the airborne X-band radar reflectivity data. They also miss the peak of the observed ice number distribution function for $\mathrm{0.1}<D_{\max }<\mathrm{1}$ mm. Even though the P3 scheme has a very different approach representing ice, it does not produce greatly different total condensed water content or better comparison to other observations in this tropical convective system. Mixed-phase microphysical processes at −10 ^∘C are associated with the overprediction of liquid water content in the simulations with the Morrison and P3 schemes. The ice water content at −10 ^∘C increases mainly due to the collection of liquid water by ice particles, which does not increase ice particle number but increases the mass/size of ice particles and contributes to greater simulated radar reflectivity.

1 Introduction

High concentrations of small ice particles ingested into jet engines can cause power loss and damaging events (Lawson et al., 1998; Mason et al., 2006). They can also cause air data probe failures (Duviver, 2010). Regions with high ice water content (HIWC), mainly composed of small ice crystals with median mass diameters (MMDs) as low as 300 µm, frequently occur over oceanic convective systems (Mason and Grzych, 2011; Ackerman et al., 2015; Leroy et al., 2016b). Such HIWC regions, with relatively low equivalent radar reflectivities (Z_e) (often less than 20 dBZ; Mason et al., 2006; Fridlind et al., 2015; Protat et al., 2016; Wolde et al., 2016; Leroy et al., 2017), are hard to detect with pilot radars onboard commercial aircraft and are thus potentially hazardous.

In order to explore the processes responsible for the occurrence of HIWC regions and the associated ice crystal properties within tropical convection, the High Altitude Ice Crystals – High Ice Water Content (HAIC-HIWC) international field campaigns (Dezitter et al., 2013; Strapp et al., 2016a) and the HIWC RADAR campaign (Ratvasky et al., 2019) were conducted. The first HAIC-HIWC field campaign took place near Darwin, Australia, during the monsoon season of 2014, with the second out of Cayenne, French Guiana, in May 2015. The HIWC RADAR campaign was conducted out of Florida in August 2015. Data collected during these campaigns have been being analyzed to understand HIWC conditions (Leroy et al., 2015, 2016a, 2017; Protat et al., 2016; Wolde et al., 2016; Korolev et al., 2020), to develop warning products that can identify HIWC regions (Yost et al., 2018; Bedka et al., 2019; Harrah et al., 2019; Haggerty et al., 2020), and to characterize the high-altitude HIWC environment to assess a new ice crystal aircraft certification envelope.

In general, total condensed water content (TWC) values in these campaigns reached as high as 4.1 g m⁻³ averaged over 0.93 km (0.5 nautical mile) distance scales, and even up to about 2 g m⁻³ over 185 km distance scales. Average MMDs in HIWC zones greater than or equal to 1 g m⁻³ increased with temperature, from ∼326 µm at −50 ^∘C to ∼708 µm at −10 ^∘C (Strapp et al., 2021). Leroy et al. (2015, 2016a, 2017) showed that MMDs decrease with increasing TWC and decreasing temperature, indicating small ice crystals are responsible for HIWC regions at high altitudes for both the Darwin and Cayenne datasets. Wolde et al. (2016) found the relationship between the ice water content (IWC) and radar equivalent reflectivity factor followed a power-law fit with coefficients dependent on temperature. However, they found the pilot X-band weather radar on the Canadian National Research Council (NRC) Convair-580 aircraft did not have adequate sensitivity to detect HIWC regions when calibrated using the NRC X-band research radar. Nguyen et al. (2019) proposed a retrieval method for IWC using the specific differential phase (K_dp) and differential reflectivity ratio (Z_dr) data from X-band dual-polarization airborne radar. This method was demonstrated to be superior to the power-law fits between IWC and reflectivity as accounting for Z_dr reduced the dependency of IWC on the variation in the shapes and orientation of ice particles. Ladino et al. (2017) concluded that secondary ice production (SIP) plays a dominant role in the formation of the observed high concentration ice crystals with ice nucleating particles making only a minor contribution. Korolev et al. (2020) proposed that a new “freezing-drop-shattering” mechanism generated small SIP particles above the melting layer at temperatures between 0 and −15 ^∘C in both oceanic tropical mesoscale convective systems (MCSs) and midlatitude frontal clouds. In this SIP mechanism, large liquid drops are transported through the melting layer to a supercooled environment by convective turbulent updrafts and then collide with aged ice, freeze, and shatter. Keinert et al. (2020) used a laboratory study to indicate that bubble bursts dominate the SIP in sea salt drizzle droplets while droplet shattering controls the SIP in pure water droplets.

Several numerical simulation studies on tropical MCSs sampled during the HAIC-HIWC projects have been conducted using different numerical models and different microphysics schemes. Franklin et al. (2016) showed that the Met Office Unified Model (UM) with a single-moment microphysics scheme overestimated the radar reflectivity above the freezing level due to the errors of simulated updraft dynamics and particle sizes, hypothesizing that a double-moment microphysics scheme would improve the model's representation of the observed variability of the ice particle size distribution (PSD). The Stanford et al. (2017) WRF simulations of four tropical deep convection events sampled during the HAIC-HIWC Darwin campaign showed that three microphysics schemes (one bin and two double-moment bulk schemes; Lynn et al., 2005; Thompson et al., 2008; Morrison et al., 2009) produced larger MMDs for TWC >1 g m⁻³ at temperatures between −10 and −40 ^∘C and a high bias in convective radar reflectivity compared to observations. They hypothesized these differences resulted from errors in parameterized hydrometeor PSDs, single-ice-particle properties (e.g., shape and density), and parameterized microphysical processes. The Qu et al. (2018) simulation of a tropical convective system on 16 May 2015 with the Environment and Climate Change Canada's Global Environmental Multiscale (GEM) model and the Milbrandt–Yau double-moment cloud microphysical scheme (MY2) (Milbrandt and Yau, 2005) also produced IWC and ice particle number concentration differing from observations during the HAIC-HIWC Cayenne project, which they hypothesized was due to the poor representation of SIP processes in the microphysics scheme.

Other tropical convective clouds have also been observed and simulated. For example, McFarquhar and Heymsfield (1996) found the numbers of smaller particles (D<100 µm) close to the convection were 1 order of magnitude higher than the numbers found further away according to data obtained during the Central Equatorial Pacific Experiment (CEPEX). Meanwhile, McFarquhar and Heymsfield (1997) indicated the shapes of the ice PSDs in the tropics substantially differ from those in the midlatitudes, especially at temperatures $<-\mathrm{40}$ ^∘C. Lohmann et al. (1995) showed that the average IWC simulated by a coarser resolution ($\sim \mathrm{125}\mathrm{km}\times \mathrm{125}\mathrm{km}$) general circulation model (GCM) agreed well with the observed IWC during CEPEX, especially with respect to the relationship between IWC and temperature, whereas the model underestimated the variability of simulated IWC within each temperature bin. Chen et al. (1997) indicated the main sources of ice particles are frozen cloud droplets and interstitial aerosol particles, and the number concentration of ice particles is influenced strongly by the amount of condensation nuclei in convective inflows according to the simulations and a sensitivity experiment for cases during CEPEX by using a one-dimensional microphysical model. Ackerman et al. (2015) conducted three 3D cloud-resolving model (CRM) simulations of MCSs observed on 23 January 2006 during the Tropical Warm Pool International Cloud Experiment (TWP-ICE) with bulk and bin microphysics schemes but were not able to produce HIWC regions (IWC >2 g m⁻³ and Z_e<30 dBZ). Lang et al. (2011) greatly reduced model bias of excessively large reflectivity values (e.g., 40 dBZ) in the middle and upper troposphere in the simulation of a continental convective case observed during the Tropical Rainfall Measuring Mission (TRMM) Large-Scale Biosphere–Atmosphere Experiment in Amazonia (LBA) through modifying a single-moment bulk microphysics scheme in the Goddard Cumulus Ensemble model; however, there was much less improvement for an oceanic MCS observed during the TRMM Kwajalein Experiment (KWAJEX).

As indicated by the review above, the numerical studies on the HIWC phenomenon to date have not been able to capture the HIWC phenomenon well. This has been attributed to biases in particle properties, parameterized PSDs, and microphysical processes. The lack of knowledge about processes generating HIWC regions suggests that further numerical simulations are needed to explore the microphysical pathways producing HIWCs. Qu et al. (2018) indicated MY2 greatly overestimated the graupel content and hypothesized that HIWC will be better estimated by the next generation of microphysics schemes (e.g., the Predicted Particle Properties (P3) microphysics scheme; Morrison and Milbrandt, 2015). Therefore, this study will be the first test of the P3 scheme in simulating the HIWC phenomenon and in comparing P3 to other bulk microphysics schemes to determine whether P3 is better able to predict HIWCs in a high-resolution numerical weather prediction (NWP) context. In this study, a tropical oceanic MCS on 26 May 2015, which was well sampled during the Cayenne field campaign, is simulated at a high resolution with 1 km horizontal grid spacing. The numerical simulation experiments and their evaluation are described in this paper. In an upcoming companion paper, attention will be focused on sensitivity experiments varying some parameters within the microphysics scheme to enhance understanding of processes leading to the formation of small crystals in HIWC regions.

The next section describes the tropical oceanic MCS sampled on 26 May 2015. Section 3 introduces the collected data and how they were processed. Section 4 shows the simulated fields and their evaluation against observations. A summary and conclusions are presented in Sect. 5.

2 Case description

The tropical MCS observed on 26 May 2015 initiated and developed over the tropical Atlantic Ocean north of Cayenne, French Guiana. The MCS was not associated with obvious synoptic-scale flow features, such as clearly identifiable highs or lows. From the soundings at Cayenne, a deep moist absolutely unstable layer (MAUL) existed over the area where the MCS occurred, and the MAUL was maintained as the MCS developed, consistent with the conceptual model of a MCS proposed by Bryan and Fritsch (2000). There were mainly easterly (westerly) winds below (above) 350 hPa. The first convection initiated and developed over the ocean near the coast in the early morning. The convection moved eastward over the course of the day due to the upper westerly winds. Subsequently, new convective cells continually initiated and developed in a similar location and gradually moved eastward to merge with old convective cells that were present over the ocean to form a large and long-lived MCS.

The convective system was sampled by two research aircraft, the French SAFIRE Falcon 20 and Canadian NRC Convair 580, during the HAIC-HIWC field campaign. Figure 1 shows the observed brightness temperature from GOES-13 geostationary satellite channel 4 (10.8 µm) at 10:45 UTC on 26 May 2015, tracks of the two flights (Fig. 1a), and the height of the aircraft above mean sea level and air temperature at the flight levels (Fig. 1b). Both aircraft sampled close to the convective core of the storm as shown by the tracks of the aircraft through the lowest cloud-top brightness temperature (Fig. 1a). The SAFIRE Falcon 20 sampled at three height levels, i.e., ∼ 7, ∼10, and ∼11.5 km, corresponding to temperatures of about −10, −30, and −45 ^∘C, respectively. The NRC Convair 580 sampled mainly at ∼7 km and a temperature of around −10 ^∘C (Fig. 1b).

Figure 1 (a) Observed brightness temperature (K, shaded) from GOES-13 geostationary satellite channel 4 (10.8 µm) at 10:45 UTC on 26 May 2015 and flight tracks of SAFIRE Falcon 20 (black thick line) and NRC Convair 580 (blue thick line). (b) Height (km, thick solid curves) of aircraft above mean sea level and air temperature (^∘C, thin dashed curves) at the flight level for SAFIRE Falcon 20 (black) and NRC Convair 580 (blue). The white mark “⊗” in (a) represents the location of the release of a radiosonde at Cayenne.

3 Data and method

3.1 Data

The SAFIRE Falcon 20 was equipped with cloud microphysics instrumentation, including a Cloud Droplet Probe (CDP2), Two Dimensional Stereo Imaging Probe (2D-S), Precipitation Imaging Probe (PIP), and Isokinetic Evaporator Probe (IKP-2, Strapp et al., 2016b). The NRC Convair 580 was equipped with an X-band (9.41 GHz) cloud airborne research radar (NAX, Wolde et al., 2016) including three antennae (nadir, zenith, and side-looking) and similar cloud microphysics instrumentation. The two optical array probes, 2D-S and PIP, recorded 2D images of ice crystals nominally in the size range of 10–1280 and 100–6400 µm, respectively. The diode resolutions of 2D-S and PIP are 10 and 100 µm, respectively. The size distribution data with uncertainty of 10 %–100 % (Baumgardner et al., 2017) are processed following the general approach described in McFarquhar et al. (2017), with only center-in particles accepted, and corrections for out-of-focus particles (Korolev, 2007), shattered particles (Field and Heymsfield, 2003; Field et al., 2006; Korolev and Field, 2015), and particles partially within the photodiode array applied (Heymsfield and Parrish, 1978). Due to a poorly defined depth of field for small particles (Baumgardner et al., 2012) and the potential of shattered artifacts only n(D) for D_max>50 µm are considered here. Composite PSDs ranging from 0.05 to 12.845 mm merged from the 2D-S and the PIP were derived at a 5 s time resolution using the 2D-S data below 800 µm, the PIP data above 1200 µm and complementary linear weights for both probes over the range of 800–1200 µm (Fontaine et al., 2017; Leroy et al., 2017). The IKP-2 bulk TWC probe was designed specifically for these campaigns to measure the high-speed, high-TWC environment, up to at least 10 g m⁻³ at 200 m s⁻¹, with a target accuracy of 20 % (Strapp et al., 2016b; Leroy et al., 2017).

Radar reflectivity data from the X-band airborne research radar installed on the NRC Convair 580 (Wolde et al., 2016), TWC measured by the IKP-2, and PSDs measured by the 2D-S and PIP installed on the SAFIRE Falcon 20 and on the NRC Convair 580 are used to statistically evaluate the model simulations. It should be noted that observed radar reflectivity profiles were sampled along the flight track of NRC Convair 580, whose horizontal locations differ from those of the flight track of the SAFIRE Falcon 20 (Fig. 1a). Part of observed TWC/PSD samples at −10 ^∘C and all observed TWC/PSD samples at −45 and −30 ^∘C are from SAFIRE Falcon 20 (Fig. 1b), implying that observed TWC/PSD at −45 and −30 ^∘C may be not fully consistent with observed radar reflectivity in the upper levels in this study. Thermodynamic and wind profiles observed by a radiosonde released at Cayenne and the Cayenne GOES-13 Satellite Cloud Products Data (https://doi.org/10.5065/D6NC5ZX6, Langley Research Center, 2016) are also used.

Figure 2 The model domain configuration (color shaded fields represent terrain elevation, in m). The horizontal grid spacings of d01 and d02 are 3 and 1 km, respectively. The white mark “⊗” represents the location of the release of a radiosonde at Cayenne.

3.2 Model setup

The WRF model version 4.1.3 (Skamarock et al., 2019) is used to simulate the tropical oceanic MCS event on 26 May 2015. Two one‐way nested domains with 3 and 1 km horizontal grid spacing and 51 vertical levels are adopted (Fig. 2). The ERA5 reanalysis data available every 1 h with 0.25^∘ × 0.25^∘ horizontal grid spacing (https://doi.org/10.5065/BH6N-5N20, last access: 13 October 2019) are used for initial and boundary conditions. The model is run from 00:00 to 18:00 UTC on 26 May 2017 for 18 h with a spin-up time of the first 6 h. Physical parameterization schemes include the revised Rapid Radiative Transfer Model (RRTMG) longwave and shortwave radiation scheme (Iacono et al., 2008), the Yonsei University (YSU) planetary boundary layer (PBL) scheme (Hong et al., 2006), the MM5 similarity surface layer scheme (Beljaars, 1995), and the unified Noah land-surface scheme (Tewari et al., 2004). The cumulus parameterization scheme is not activated in this study.

Four bulk microphysics schemes, namely the WRF single-moment 6-class (WSM6) microphysics scheme (Hong and Lim, 2006), the Morrison double-moment scheme (Morrison et al., 2009) and the Predicted Particle Properties (P3) microphysics scheme with one- and two-ice options (Morrison and Milbrandt, 2015; Milbrandt and Morrison, 2016) are used for separate simulations. The simulations using the WSM6 and Morrison microphysics schemes are referred to as the WSM6 and MORR runs hereafter, respectively. These microphysics schemes are used because the PSDs of ice species are parameterized differently in these schemes. Both the WSM6 and Morrison schemes predict the mixing ratios of five cloud hydrometeor species, including cloud water, rainwater, cloud ice, snow, and graupel¹, while the number mixing ratios for all species except cloud water are also predicted in the Morrison scheme. The P3 scheme predicts bulk ice properties (e.g., mean particle density) rather than predicting separate species of ice with fixed properties (e.g., cloud ice, snow, and graupel). P3 uses one or more “free” ice categories to represent all ice-phase hydrometeors, which can eliminate the unphysical “conversion” processes between different traditional ice categories (Morrison and Milbrandt, 2015; Milbrandt and Morrison, 2016). In this study, the options of one- and two-ice categories in the P3 scheme are used, referred to as P3-1ICE and P3-2ICE hereafter, respectively. Technically P3-1ICE and P3-2ICE are two configurations of the same scheme, but they have notably different treatments of ice which is the basis on which all the microphysics schemes were chosen, so these are referred to as different schemes. Output data in the model domain d02 with 1 km horizontal grid spacing are analyzed in this study. It should be noted that cloud ice, snow, and graupel in WSM6 and MORR and both categories of ice in P3-2ICE are treated as a single category of ice particles to compare with the observed ice particles, because the observed ice particles are not separated into different categories.

In WSM6, the main microphysical processes associated with ice particles include ice nucleation, deposition, sublimation, homogeneous and heterogeneous freezing, collection by liquid and other ice categories, autoconversion to other ice categories, melting, and sedimentation (Hong and Lim, 2006). MORR also considers the rime-splintering process in addition to the same microphysical processes as WSM6 (Morrison et al., 2009). Both P3-1ICE and P3-2ICE consider the main microphysical processes associated with ice particles including ice nucleation, deposition, sublimation, homogeneous and heterogeneous freezing, collection by liquid, self-collection, melting, and sedimentation, while collision and collection between ice categories and rime-splintering processes are considered only in P3-2ICE (Morrison and Milbrandt, 2015; Milbrandt and Morrison, 2016).

3.3 Estimation of X-band radar reflectivity

The computations of simulated radar reflectivity are performed using the Rayleigh approximation which is applicable at the X-band given the size of typical ice particles (Ryzhkov et al., 2020). The relations for reflectivity from rain (Z_r), graupel (Z_g), snow (Z_s), and ice (Z_i) are derived in detail in Appendix A. The total equivalent radar reflectivity factor (Z_e) in units of dBZ can thus be attained using

$$\begin{array}{}\text{(1)}& Z_{\mathrm{e}}=\mathrm{10}\times {\log }_{\mathrm{10}}\left({\displaystyle \frac{Z_{\mathrm{r}}+Z_{\mathrm{g}}+Z_{\mathrm{s}}+Z_{\mathrm{i}}}{\mathrm{1}{\mathrm{mm}}^{\mathrm{6}}{\mathrm{m}}^{-\mathrm{3}}}}\right).\end{array}$$

Figure 3 Skew-T plot of observed and simulated Cayenne radiosonde data at 12:00 UTC on 26 May 2015. One full wind barb represents 10 knots (∼5.14 m s⁻¹). The location where Cayenne radiosonde was released is shown in Fig. 2.

4 Results

4.1 Evaluation of simulated sounding, brightness temperature, and radar reflectivity

Figure 3 shows a Skew-T plot of the observed and simulated thermodynamic and wind profiles over Cayenne at 12:00 UTC on 26 May 2015. The observed profiles of air temperature and dew-point temperature show a very moist environment especially between 800 and 350 hPa. Regardless of the choice of microphysical scheme, the moist environment is well simulated although the layer between 500 and 350 hPa is drier with a maximum dew point depression (T−T_d) of ∼5 ^∘C in P3-2ICE and almost zero in the observations. This upper-level drier layer is associated with the initial condition from ERA5 reanalysis data (not shown). The MORR scheme is moister between 500 and 350 hPa, which is more consistent with observations as the maximum T−T_d is less than 2 ^∘C (Fig. 3). As for the wind profiles, all the simulations predict the observed easterly and westerly winds in the lower and upper troposphere respectively, though there are some biases around 300 hPa and up to 250 hPa, where simulated and observed winds are in the opposite directions and wind speeds differ by ∼12 m s⁻¹.

Figure 4 Simulated and observed brightness temperature (BT, K, shaded) from GOES-13 channel 4 (10.8 µm) at 10:45 UTC on 26 May 2015.

Figure 4 shows simulated and observed brightness temperatures (BT) at 10:45 UTC on 26 May 2015. The simulated BT is calculated using the Community Radiative Transfer Model (CRTM, https://www.jcsda.org/jcsda-project-community-radiative-transfer-model, last access: 1 November 2019) using the assumptions consistent with those in the different microphysics schemes, including the characteristics of the cloud species and PSDs. The storm coverage (BT <232 K, yellow to deep red areas) in MORR (∼59 % of domain) is larger than that in the observations (∼47 % of domain), while WSM6 produces a smaller storm coverage (∼33 % of domain) compared to observations. The storm coverages (BT <232 K) in P3-1ICE (∼51 % of domain) and P3-2ICE (∼46 % of domain) are closer to that of the observations (Fig. 4). The lower brightness temperature areas representing deep convection (BT <212 K, red areas in Fig. 4) are larger in MORR (∼25 % of domain), P3-1ICE (∼20 % of domain), and P3-2ICE (∼13 % of domain) and smaller in WSM6 (∼5 % of domain) than in the observations (∼9 % of domain).

Figure 5 Frequency (%, shaded) of simulated and observed BT from GOES-13 channel 4 from 06:15 to 17:45 UTC on 26 May 2015. Contours represent cumulative frequencies 5 %, 25 %, 50 %, 75 %, and 95 %, respectively.

To examine the storm evolution, the frequency distributions of simulated and observed BT from 06:15 to 17:45 UTC on 26 May 2015 are displayed in Fig. 5. By 10:15 UTC, the dominant frequency (>4 %) of simulated and observed BT is around 280 K, indicating that clear regions dominate the domain in the early stage of MCS. The subdominant frequency (2 %–4 %) around 220 K illustrates the deep convection. This subdominant frequency in all of the simulations is consistent with observations, even though the BT ranges are all slightly different (i.e., 216–232 K in WSM6, 202–226 K in MORR, 200–228 K in P3-1ICE, 206–228 K in P3-2ICE, and 208–230 K in observations). There are 25 % of brightness temperatures <226 K in WSM6, <212 K in MORR, <214 K in P3-1ICE, <216 K in P3-2ICE, and <220 K in the observations at 10:15 UTC, indicating MORR, P3-1ICE, and P3-2ICE overpredict strong convective areas, and WSM6 underpredicts strong convective areas compared to the observations. After 10:15 UTC, the frequency around 220 K becomes dominant in the observations and simulations indicating the MCS develops and deep convective areas enlarge, as the maximum frequency exceeds 10 % in the observations and in MORR. The frequency ranges (>4 %) cover 222–246 K in WSM6, 206–234 K in MORR, 208–238 K in P3-1ICE, 210–236 K in P3-2ICE, and 212–232 K in the observations. There are 50 % of BT <238 K in WSM6, <220 K in MORR, <225 K in P3-1ICE, <228 K in P3-2ICE, and <227 K in the observation at 14:15 UTC. Overall this indicates that the storm coverage in MORR, P3-1ICE, and P3-2ICE is more consistent with the observations than that in WSM6. The frequency of brightness temperatures 214–224 K is over 10 % after 12:30 UTC in MORR, and during 15:00–16:00 UTC in the observations. This means deep convective areas as defined by BT metric are larger in MORR than the observations at an early stage in the system. Overall, all the simulations generally reproduce the storm coverage and evolution of this tropical MCS with the average bias (difference between simulation and observation) in storm coverage (BT <232 K) of $\sim -\mathrm{34.3}$ % in WSM6, ∼30.0 % in MORR, ∼12.9 % in P3-1ICE, and ∼2.3 % in P3-2ICE, indicating there is a relatively larger bias in WSM6 in which the storm coverage is smaller. WSM6 underestimates the observed deep convection areas (BT <212 K) by ∼55.5 % ² and the overestimates in MORR, P3-1ICE and P3-2ICE are ∼175.4 %, ∼178.2 %, and ∼76.5 %, respectively (Fig. 5).

Figure 6 Diagram of sampling method for radar reflectivity profiles. R=100 km here. The two magenta points at the flight track and model grid point represent the horizontal locations of observed and simulated radar reflectivity profiles. See text for details.

To obtain a statistical comparison between simulated and observed radar reflectivity, contoured frequency by altitude diagrams (CFADs) (Yuter and Houze, 1995) are used. Because the observed reflectivity is only available along the flight track, equivalent locations for sampling the reflectivities from the modeled fields must be determined. The sampling method here is based on the flight track and observed and simulated BTs. Because there exists a bias between the simulated and observed BTs, they are first spatially normalized respectively within the same region as the model domain d02. Using normalized BTs, here paired samples of observations and simulations are found whose locations are about the same distance from the convective cores. Only the vertical profiles of observed reflectivity in which the observed TWC at the flight locations are larger than 0.1 g m⁻³ are selected as observational samples for the statistical comparison. Across an area with the horizontal location of observed reflectivity profiles as the center and 100 km as the range (Fig. 6), the simulated vertical profiles of reflectivity at model grid points where simulated TWC at the flight level is larger than 0.1 g m⁻³ and the normalized simulated brightness temperature is closest to the normalized observed brightness temperature are selected as simulation samples for statistical comparison. It should be noted that radii from 20 to 200 km in 20 km intervals were tested, and the results were similar. A radius of 100 km was adopted because the standard deviation between the observed and simulated reflectively was the lowest when using this radius threshold. The sampling method used here is similar to the one used by Borderies et al. (2018).

Figure 7 Frequency (%, shaded) of simulated (a: WSM6, b: MORR, c: P3-1ICE, d: P3-2ICE) and observed (e) X-band radar reflectivity. Contours represent cumulative frequencies 25 %, 50 %, 75 %, and 95 %, respectively. The gap at 6–7 km in (e) is due to observed radar reflectivity data not available at the first few range gates (within ∼500 m) from the aircraft.

In general, the simulations overestimate the X-band radar reflectivity above the melting layer (∼4.7 km). Figure 7 shows CFADs and cumulative CFADs of simulated and observed X-band radar reflectivity above 5 km. The CFADs are shown only above 5 km, because the formulae for calculating the simulated radar reflectivity in Sect. 3.3 do not consider the effect of melting, and this study mainly focuses on HIWC regions. From Fig. 7e, 95 % of the cumulative observed reflectivities are <30 dBZ above 6 km. The most frequently (frequency >10 %) observed radar reflectivity is around 25 dBZ at heights of 5–7 km and ∼15 and ∼20 dBZ at a height of ∼8 km (Fig. 7e). The simulated radar reflectivity shows broader distributions with larger values than the distributions of the 95 % cumulative observed reflectivities <30 dBZ above 6 km), and maxima in radar reflectivity can reach 50 dBZ, especially for WSM6, P3-1ICE, and P3-2ICE (Fig. 7a, c, and d). There are 95 % of the cumulative simulated reflectivities <44 dBZ in WSM6, <41 dBZ in MORR, <45 dBZ in P3-1ICE, and <47 dBZ in P3-2ICE above 6 km. The simulated radar reflectivity in MORR has a narrower distribution with 70 % of reflectivities between 34 and 42 dBZ at 5 km (Fig. 7b), which better resembles the observation with 70 % of reflectivities between 24 and 36 dBZ at 5 km (Fig. 7e). The other simulations have broader distributions with 70 % of the reflectivities between 30 and 44 dBZ in WSM6 (Fig. 7a), between 17 and 46 dBZ in P3-1ICE (Fig. 7c), and between 25 and 48 dBZ in P3-2ICE (Fig. 7d) at 5 km. The radar reflectivity in all simulations extends above 14 km, whereas the observed radar reflectivity is mainly below 14 km. Examining the reflectivity and Doppler velocity from zenith-viewing Doppler airborne radar shows a peak-to-peak correlation between them. This suggests there may be stronger updrafts in the simulations associated with the higher extended simulated reflectivity. Overall, all the simulations overestimate the intensity and spatial extent of radar reflectivity. By examining each component of reflectivity, the overestimation of radar reflectivity above the melting layer in WSM6 and MORR mainly results from the overprediction of graupel (not shown), which is similar to the tropical MCS simulations of Lang et al. (2011) and Qu et al. (2018). The P3 scheme, which was expected to yield better estimates of HIWCs, does not reduce the biases in simulated radar reflectivities. It should be noted that the NRC Convair 580 operations avoided the cloud regions with high reflectivity due to safety regulations, and thus it did not approach high reflectivity regions (red zones on the pilot’s radar) within 30 nautical miles (∼55.56 km). However, the BT sampling method has been used to minimize these aircraft sampling biases.

Figure 8 Observed (red) and simulated (blue, a: WSM6, b: MORR, c: P3-1ICE, d: P3-2ICE) ice particle size distributions at the level of −45 ^∘C. The red and blue dashed lines indicate the 25th and 75th percentiles, and the red and blue solid lines represent the median.

Figure 9 Observed (red) and simulated (blue, a: WSM6, b: MORR, c: P3-1ICE, d: P3-2ICE) ice particle size distributions at the level of −30 ^∘C. The red and blue dashed lines indicate the 25th and 75th percentiles, and the red and blue solid lines represent the median.

Figure 10 Observed (red) and simulated (blue, a: WSM6, b: MORR, c: P3-1ICE, d: P3-2ICE) ice particle size distributions at the level of −10 ^∘C. The red and blue dashed lines indicate the 25th and 75th percentiles, and the red and blue solid lines represent the median.

4.2 Cloud microphysical properties

Samples are selected to examine the observed and simulated cloud microphysical properties using the same method as used for sampling radar reflectivity profiles. Compared to the observed ice PSDs, the simulated PSDs have different shapes and variability depending on the microphysics scheme. Figures 8–10 show observed and simulated ice PSDs at the levels of −45, −30, and −10 ^∘C, respectively. These are the three temperature levels at which the in situ observations were focused. The simulated ice PSDs are the composite PSDs of all the ice-phase hydrometeors predicted in each microphysics scheme.

Table 1 Number concentration of median n(D) for 0.1 mm $<D_{\max }<\mathrm{12.845}$ mm (N_{0.1−12.845 mm}) and contributions of small (0.1 mm $<D_{max}<$ 0.3 mm, C_{0.1−0.3 mm}), medium (0.3 mm $<D_{\max }<\mathrm{1}$ mm, C_0.3−1 mm), and large (1 mm $<D_{\max }<\mathrm{12.845}$ mm, C_{1−12.845 mm}) particles at −45, −30, and −10 ^∘C. The bold text indicates the dominant contribution.

At the −45 ^∘C level (Fig. 8), compared to the observations, WSM6 underestimates the median number distribution function n(D) by ∼50 % ($\sim \mathrm{3}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹) near the maximum dimension (D_max) of 0.1 mm (Fig. 8a), while P3-1ICE and P3-2ICE overestimate the median n(D) by ∼282 % ($\sim \mathrm{17}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹) and ∼199 % ($\sim \mathrm{12}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹) respectively at this D_max (Fig. 8c and d). MORR has similar median n(D) magnitude to the observations with both $\sim \mathrm{6}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹ near the D_max of 0.1 mm (Fig. 8b). WSM6, MORR, and P3-2ICE underpredict the observed number concentration for 0.1 mm $<D_{\max }<\mathrm{12.845}$ mm (N_{0.1−12.845 mm}) by ∼17 %, ∼50 %, and ∼16 % respectively, while P3-1ICE overpredicts it by ∼5 times (Table 1). Small particles (0.1 mm $<D_{\max }<\mathrm{0.3}$ mm) consistently make dominant contributions to the total number concentration in observations and simulations, while P3-1ICE produces too many small particles (∼86.6 %) with the overprediction of total number concentration (Table 1). Compared to the observed PSDs, PSDs in WSM6 and MORR have a smaller spread (Fig. 8a and b), and PSDs in P3-1ICE and P3-2ICE have a larger spread (Fig. 8c and d). Based on the shapes, magnitudes, and spreads of PSDs at −45 ^∘C, PSDs in MORR are most consistent with the observations among the simulations.

At the −30 ^∘C level (Fig. 9), all the shapes of simulated PSDs are similar to those at −45 ^∘C, while the magnitudes are smaller by about an order of magnitude (Fig. 8). WSM6 and MORR have similar PSD characteristics in terms of spread and magnitude. The median n(D) in WSM6 and MORR have similar magnitude (∼10⁵ m⁻³ mm⁻¹) at D_max of 0.1 mm compared to the observations (Fig. 9a and b). P3-1ICE and P3-2ICE overestimate the median n(D) by ∼6 times ($\sim \mathrm{5.7}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹) and ∼264 % ($\sim \mathrm{2.5}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹) respectively at D_max of 0.1 mm (Fig. 9c and d). None of the simulations capture the peak of the observed PSD with the median of $\sim \mathrm{2.6}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹ near D_max of 0.3 mm (Fig. 9). There are no obvious peaks in PSDs for 0.1 mm $<D_{\max }<\mathrm{1}$ mm in WSM6, MORR, and P3-2ICE (Fig. 9a, b, and d). There is a PSD peak with a median n(D) of $\sim \mathrm{9.6}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹ near D_max of 0.17 mm in P3-1ICE (Fig. 9c). Medium particles (0.3 mm $<D_{\max }<\mathrm{1}$ mm) are dominant in the observations and WSM6, while small particles make the main contributions in MORR, P3-1ICE, and P3-2ICE (Table 1).

At the −10 ^∘C level (Fig. 10), all simulations underestimate the median PSD for D_max<1 mm, especially P3-2ICE, with an underestimate of ∼94 % ($\sim \mathrm{4.0}\times {\mathrm{10}}^{\mathrm{4}}$ m⁻³ mm⁻¹) at D_max of 0.1 mm. Large particles contribute 58.6 % of N_{0.1−12.845 mm} in P3-2ICE, which is very different from the observations (Table 1). Compared to the observations, MORR has almost the same PSD spread (range between the 25th and 75th percentiles) and median for D_max of 0.1 mm (Fig. 10b). All of the simulations miss the peak of the observed PSD with a median of $\sim \mathrm{1.5}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ mm⁻¹ near D_max of 0.3 mm (Fig. 10). There are enough samples with TWC >1 g m⁻³ at −10 ^∘C to examine PSDs for 0.1 g m${}^{-\mathrm{3}}<$ TWC <1 g m⁻³ and TWC >1 g m⁻³ separately. Both results indicate that all of the simulations miss the peak of the observed PSD near D_max of 0.3 mm (Figs. S1 and S2), which is consistent with that for TWC >0.1 g m⁻³ (Fig. 10).

Overall, the simulated PSDs at the three temperature levels shown in Figs. 8–10 have biases in various degrees with respect to their shapes, magnitudes, and spreads compared to observations. The observations are concentrated around the smaller crystal sizes than are the simulations for most of temperature levels, except P3-1ICE at −45 and −30 ^∘C. These qualitative results do not change through examining the PSDs normalized by TWC at −45, −30, and −10 ^∘C (not shown). It should be noted that the PSD comparison strongly reflects different assumptions built into the schemes regarding PSD shapes, namely inverse exponential PSDs are assumed in WSM6 and MORR, whereas a gamma PSD with a shape parameter (μ) that varies with the slope (λ) following the observations of Heymsfield (2003) that are used to formulate the P3 scheme.

Figure 11 Violin plots of observed (red) and simulated (blue: WSM6, green: MORR, orange: P3-1ICE, magenta: P3-2ICE) (a) total water content (TWC, g m⁻³), (b) number concentration (# m⁻³) within 0.1 mm $<D_{\max }<\mathrm{3}$ mm (N_0.1−3 mm), (c) effective diameter (D_e, mm), and (d) vertical velocity (m s⁻¹) at temperatures of −10, −30, and −45 ^∘C. As for the black box-and-whisker plots, the extremes of the whiskers indicate the 5th and 95th percentiles, the lower and upper limits of the boxes correspond to the 25th and 75th percentiles, the dividing line represents the median value, and the white points represent the average value of samples between the 5th and 95th percentiles. The width of the shaded area represents the proportion of the data located there, and only areas between the 5th and 95th percentiles are shown. The numbers below the upper horizontal axis represent the total number of samples. The red triangles above the bottom horizontal axis indicate that the difference between the experiment mean and the observation mean is not statistically significant (not passing the significant test for p<0.05).

Figure 11 shows statistical distributions of observed and simulated TWC, number concentration for particle sizes of 0.1 mm $<D_{\max }<\mathrm{3}$ mm (N_0.1−3 mm), effective diameter (D_e), and vertical velocity using violin plots (or box-percentile plots, Esty and Banfield, 2003) at temperatures of −10, −30, and −45 ^∘C, respectively. The shaded areas of violin plots represent the proportion of the samples outlining the kernel probability densities. It should be noted that the number concentrations for 0.2 mm $<D_{\max }<\mathrm{3}$ mm were also examined, and the conclusions are consistent. Thus, only the results using N_0.1−3 mm are shown here. There are several different definitions of D_e (McFarquhar and Heymsfield, 1998). Given the main purpose here is to compare the particle sizes between observation and simulation for simplicity, the definition of D_e given by

$$\begin{array}{}\text{(2)}& D_{\mathrm{e}}={\displaystyle \frac{{\int }_{\mathrm{0.1}\mathrm{mm}}^{\mathrm{12.845}\mathrm{mm}}{\sum }_{k=\mathrm{1}}^K{D}^{\mathrm{3}}{n}_k\left(D\right)\mathrm{d}D}{{\int }_{\mathrm{0.1}\mathrm{mm}}^{\mathrm{12.845}\mathrm{mm}}{\sum }_{k=\mathrm{1}}^K{D}^{\mathrm{2}}{n}_k\left(D\right)\mathrm{d}D}}\end{array}$$

is used, where K is the number of ice species. Since D is maximum dimension, only one number distribution function that includes all data is used for the observations, and therefore K=1 for observations.

Generally all the simulations, especially MORR, reproduce the TWC reasonably well within the same order of magnitude as observations at the three temperature levels, with biases within 38 % at −10 ^∘C. In particular at the −10 ^∘C level the 25th and 75th percentiles of TWC in all the simulations cover the same order of magnitude as the observations (Fig. 11a). The differences in N_0.1−3 mm among the simulations are quite large (Fig. 11b). At the three temperature levels, WSM6 and especially MORR underestimate the number concentration (Fig. 11b). This is associated with the underpredicted small particles and overpredicted large particles in WSM6 and MORR compared to the observations (Figs. 8–10). Thus, WSM6 and MORR produce larger D_e compared to the observations at the three temperature levels (Fig. 11c). At −30 ^∘C, the median N_0.1−3 mm values in WSM6 ($\sim \mathrm{0.4}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³) and MORR ($\sim \mathrm{0.2}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³) are underestimated by ∼50 % and ∼75 %, respectively, consistent with the underestimate of particle number near the peak of the observed PSD at D_max of ∼0.3 mm (Fig. 9a and b). The N_0.1−3 mm at −45 ^∘C in P3-1ICE is about 1 order of magnitude larger than observed (Fig. 11b) mainly due to many more small particles for 0.1 mm $<D_{\max }<\mathrm{0.4}$ mm (Fig. 8c). Similarly, compared to the observations, P3-1ICE overestimates N_0.1−3 mm at −30 ^∘C, with the median overestimated by ∼129 % (∼10⁵ m⁻³) (Fig. 11b). This is explained by an overestimate of particle number at D_max of ∼0.1 mm (Fig. 9c). Accordingly, P3-1ICE produces smaller D_e than the observations at −45 and −30 ^∘C with an underestimate of median D_e by ∼38 % and ∼46 % respectively (Fig. 11c). Due to the larger PSD spreads in P3-2ICE at −45 and −30 ^∘C (Figs. 8d and 9d), N_0.1−3 mm and D_e in P3-2ICE accordingly have a larger spread than the observations (Fig. 11b and c). This occurs even though the spread of TWC between P3-2ICE and observations is similar (Fig. 11a). At −10 ^∘C, values of N_0.1−3 mm from all of the simulations, in particular P3-2ICE, are about 1 order of magnitude smaller than observed (Fig. 11b), implying larger mean particle size than observed (Fig. 11c). This is mainly attributed to the underestimate of small particle number for D_max<1 mm, especially near the peak of the observed PSD at D_max of ∼0.3 mm, and the overestimate of large particles in all of the simulations (Fig. 10). The simulated vertical velocity is in general stronger than in the observations, especially at −45 and −10 ^∘C (Fig. 11c), corresponding to the higher extent of simulated radar reflectivity (Fig. 7).

4.3 Cloud microphysical processes

As discussed above, all the four microphysics schemes underpredict the number concentration by about 1 order of magnitude at −10 ^∘C compared to the observations (Fig. 11b), although they predict similar TWC to the observed TWC (Fig. 11a). In this section, the processes producing HIWCs are determined. WSM6 is a single-moment scheme in which the number concentration of ice particles is not predicted directly. Through examining the number concentration derived diagnostically from the water content of each hydrometeor species in WSM6, it is found that their distributions are similar to those in MORR, especially the distributions of ice particles (not shown). Thus, only the double-moment MORR and P3 schemes are examined in detail here. Regions with IWC >1 g m⁻³ are defined as HIWC regions in this study. There are very limited observed and simulated HIWC samples at temperatures of −45 and −30 ^∘C (Fig. 11a). Through comparing between simulated sampling profiles with TWC >0.1 g m⁻³ at −45 and −30 ^∘C and simulated sampling profiles with TWC >1 g m⁻³ at −10 ^∘C, it is found that their main microphysical processes are the same at the same vertical levels. Hence, the profiles of water content, number concentration, and microphysical processes with HIWC regions at −10 ^∘C are used as examples to discuss here. The subsamples whose observed and corresponding simulated TWCs at −10 ^∘C are larger than 1 g m⁻³ are selected from the total samples (1778) at −10 ^∘C (Fig. 10) to conduct a composite analysis. There are 509, 488, and 427 paired samples selected for MORR, P3-1ICE, and P3-2ICE, respectively (Table 2).

Table 2 Averaged TWC, N_0.1−3 mm and air vertical velocity (w) of HIWC samples (TWC >1 g m⁻³) at −10 ^∘C. The differences in all averages between the simulations and observations pass the significant test for p<0.05.

^a Observation. ^b Model.

From Table 2, for the observation–MORR paired HIWC samples the average observed and simulated TWCs at −10 ^∘C are similar, ∼1.6 and ∼1.8 g m⁻³, respectively, while the simulated N_0.1−3 mm ($\sim \mathrm{1.78}\times {\mathrm{10}}^{\mathrm{4}}$ m⁻³) is about 1 order of magnitude less than the observed N_0.1−3 mm ($\sim \mathrm{1.72}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³). The average simulated air vertical velocity for the HIWC points (∼0.32 m s⁻¹) is about twice as large as the observations (∼0.16 m s⁻¹). For the observation–P3-1ICE paired samples, the average observed and simulated TWCs at −10 ^∘C are ∼1.7 and ∼1.8 g m⁻³ respectively, while the simulated N_0.1−3mm ($\sim \mathrm{1.25}\times {\mathrm{10}}^{\mathrm{4}}$ m⁻³) is underpredicted compared to the observed N_0.1−3 mm ($\sim \mathrm{1.98}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³) by ∼94 %. The simulated air vertical velocity (∼1.12 m s⁻¹) is 3.67 times greater than observed (∼0.24 m s⁻¹). Similarly, for the observation–P3-2ICE paired samples, the averaged observed and simulated TWCs at −10 ^∘C are ∼1.7 and ∼1.9 g m⁻³ respectively, while the simulated N_0.1−3 mm ($\sim \mathrm{4.63}\times {\mathrm{10}}^{\mathrm{3}}$ m⁻³) is about 2 orders of magnitude less than the observed N_0.1−3 mm ($\sim \mathrm{2.02}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³), and the simulated air vertical velocity (∼2.05 m s⁻¹) is about 15 times larger than the observation (∼0.13 m s⁻¹). Therefore, both MORR and P3, in particular P3-2ICE, substantially underpredict the ice particle number for 0.1 mm $<D_{\max }<\mathrm{3}$ mm and overpredict the vertical motion in the HIWC regions, which is associated with stronger and higher-extended simulated radar reflectivity (Fig. 7).

Figure 12 Vertical profiles of averaged water content (a, c, e, solid curves, g m⁻³) and total number concentration (b, d, f, solid curves, # m⁻³) of each hydrometeor category in MORR (a, b – cloud water: blue; rain water: red; cloud ice: magenta; snow: orange; graupel: green; and sum of ice species: black), P3-1ICE (c, d – cloud water: blue; rain water: red; and cloud ice of first category: black), and P3-2ICE (e, f – cloud water: blue; rain water: red; cloud ice of first category: orange; cloud ice of second category: green; and sum of two ice categories: black). The thick black dashed lines represent vertical profiles of vertical velocity (m s⁻¹). The vertical thin black dashed lines indicate vertical velocity of 0 m s⁻¹. The horizontal thin black and gray dashed lines indicate the heights of temperatures at 0, −10, −30, and −45 ^∘C from bottom to top, respectively. The numbers of samples used to calculate the averages are shown in Table 2.

From the vertical profiles of water content for each hydrometeor class in MORR (Fig. 12a), graupel is dominant above the 0 ^∘C layer up to ∼8 km height, with snow being the second most important class. Snow is the dominant hydrometeor category above 8 km in terms of mass content. The aforementioned underestimate of ice particle number concentration at −10 ^∘C is associated with the total number concentration ($<\mathrm{4}\times {\mathrm{10}}^{\mathrm{4}}$ m⁻³) of cloud ice, snow, and graupel (Fig. 12b). According to Eq. (A17), under the condition of similar simulated and observed TWCs, an underestimate of ice particle number concentration, especially graupel, leads to mass concentrating at the large end of PSD resulting in large reflectivities. Moreover, not only is the vertical motion overpredicted in the HIWC regions as mentioned above, but also the height with air vertical velocity >0 m s⁻¹ in MORR can be up to 15 km with a maximum of ∼1 m s⁻¹ at 13 km (Fig. 12a and b). These characteristics are consistent with the stronger intensity and higher distribution of simulated radar reflectivity shown in Fig. 7b.

The total ice number concentration in P3-1ICE is ∼10⁴ m⁻³ at −10 ^∘C, and the maximum is $\sim \mathrm{4}\times {\mathrm{10}}^{\mathrm{5}}$ m⁻³ at −45 ^∘C (Fig. 12d). The air vertical velocity >0 m s⁻¹ in P3-1ICE extends from heights of 1.5 to 15.5 km with a maximum of ∼1.4 m s⁻¹ at 12.5 km (Fig. 12c and d), which is stronger than MORR. Compared to P3-1ICE, there is one more “free” ice category considered in P3-2ICE. The difference in mean mass-weighted diameters between the two ice categories is over 500 µm. The vertical motion in P3-2ICE is stronger with a maximum of ∼2.1 m s⁻¹ at 6 km (Fig. 12e). Although the IWCs are similar between P3-1ICE and P3-2ICE, the total number concentration of the two ice categories in P3-2ICE is less than 5×10³ m⁻³ at −10 ^∘C, which is less than that in P3-1ICE (Fig. 12d and f). The average liquid water contents (LWCs) in MORR, P3-1ICE, and P3-2ICE at −10 ^∘C are ∼0.029, ∼0.094, and ∼0.187 g m⁻³, respectively, while the observed LWC from the Cloud Droplet Probe (CDP) is less than 0.008 g m⁻³. Thus, all of the simulations, especially P3-2ICE, overpredict LWC at −10 ^∘C and air vertical velocity above the 0 ^∘C layer. The distributions of vertical velocity (Fig. 11d) and Doppler velocity from zenith-viewing Doppler airborne radar (not shown) confirm that the model produces stronger updrafts than observations. It should be noted that both the NRC Convair 580 and SAFIRE Falcon 20 avoided cloud regions with strong updrafts where the presence of liquid phase is expected. However, the BT sampling method has been used to minimize these sampling biases.

Figure 13 Vertical profiles of (a–c) averaged mass conversion rates (solid lines, 10⁻³ g m⁻³ s⁻¹) and (d–f) averaged number concentration tendencies (solid lines, # m⁻³ s⁻¹) of (a, d) cloud ice, (b, e) snow, and (c, f) graupel due to different microphysical processes in MORR. Only the profiles whose column-maximum conversion rates are larger than 10⁻⁶ g m⁻³ s⁻¹ and column-maximum number concentration tendencies are larger than 1 m⁻³ s⁻¹ are shown in (a–c) and (d–f), respectively. The total microphysics tendencies of cloud ice (a: Qi_TEND, d: Ni_TEND), snow (b: Qs_TEND, e: Ns_TEND), and graupel (c: Qg_TEND, f: Ng_TEND) are shown by black solid curves. The horizontal thin black and gray dashed lines indicate the heights of 0, −10, −30, and −45 ^∘C from bottom to top, respectively. The numbers of samples used to calculate the averages are shown in Table 2. C: cloud ice; R: rain water; I: cloud ice; S: snow; and G: graupel.

From vertical profiles of microphysical conversion rates in MORR (Fig. 13), the main source terms of cloud ice content are ice nucleation from homogeneous and heterogeneous freezing on aerosol at −45 ^∘C and vapor deposition at −30 and −10 ^∘C, and the main sink terms are collection by snow and autoconversion to snow. The net conversion rate of cloud ice (Qi_TEND, sum of all microphysical conversion rates including sedimentation) at −10 ^∘C is negative (Fig. 13a). The net number concentration tendency of cloud ice (Ni_TEND) at −10 ^∘C in MORR is $\sim -\mathrm{20}$ m⁻³ s⁻¹, mainly due to the accretion of cloud ice, autoconversion to snow and sublimation (Fig. 13d). The IWC at −10 ^∘C mainly consists of snow and graupel (Fig. 12a), and their main source terms are vapor deposition and collection of cloud water (Fig. 13b and c). The net snow conversion rate (Qs_TEND) is $\sim \mathrm{0.1}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³ s⁻¹ at −10 ^∘C (Fig. 13b), while there is a negative net conversion rate of graupel (Qg_TEND) ($\sim -\mathrm{0.6}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³ s⁻¹) mainly due to its large sedimentation tendency at −10 ^∘C with a rate of $\sim -\mathrm{0.8}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³ s⁻¹ (Fig. 13c). Autoconversion of cloud ice to snow (∼5 m⁻³ s⁻¹) and the sedimentation term (∼1 m⁻³ s⁻¹) increase the snow particle number; however, they are offset by the self-collection of snow (aggregation), collection of cloud water by snow to form graupel, and sublimation (Fig. 13e). The collection of cloud water by snow to form graupel is the main production term of graupel particle number, while it is offset by the sedimentation and sublimation terms. Finally, the net number concentration tendencies of both snow and graupel (Ns_TEND and Ng_TEND) are near 0 m⁻³ s⁻¹ at −10 ^∘C (Fig. 13e and f). Therefore, the collection of cloud water by graupel is the key source term of total IWC at −10 ^∘C in MORR, which increases the mean mass/size of graupel and does not directly influence its number leading to mass concentrating at the large end of PSD. This is associated with the strong simulated reflectivity above the melting layer.

Figure 14 Vertical profiles of (a, c, e) averaged mass conversion rates (solid lines, 10⁻³ g m⁻³ s⁻¹) and (b, d, f) averaged number concentration tendency (solid lines, # m⁻³ s⁻¹) of ice categories due to different microphysical processes in (a, b) P3-1ICE and (c–f) P3-2ICE. Only the profiles whose column-maximum conversion rates are larger than 10⁻⁶ g m⁻³ s⁻¹ and column-maximum number concentration tendencies are larger than 1 m⁻³ s⁻¹ are shown in (a, c, e) and (b, d, f), respectively. The total microphysics tendencies of first ice category (a, c Qi_TEND; b, d: Ni_TEND) and second ice category (e: Qi2_TEND; f: Ni2_TEND) are shown by black solid curves. The horizontal thin black and gray dashed lines indicate the heights of 0, −10, −30, and −45 ^∘C from bottom to top, respectively. The numbers of samples used to calculate the averages are shown in Table 1. C: ice; R: rain water; I: cloud ice.

From vertical profiles of microphysical conversion rates in P3-1ICE (Fig. 14a), the main production terms of ice content are vapor deposition at −45 and −30 ^∘C, collection of cloud water by ice, vapor deposition, and collection of rain water by ice at −10 ^∘C. High ice number concentrations at the upper levels in P3-1ICE (Fig. 11b) should be associated with homogeneous freezing of cloud droplets (Fig. 14b). As for the first ice category in P3-2ICE, in addition to the same main production terms as in P3-1ICE, there is another source term merging with the second ice category due to similar mean mass-weighted diameters between the two ice categories (Fig. 14c). The net tendency of the two ice categories in P3-2ICE (i.e., Qi_TEND + Qi2_TEND) is $\sim \mathrm{1.05}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³ s⁻¹ at −10 ^∘C, which is much larger than that in P3-1ICE ($\sim \mathrm{0.02}\times {\mathrm{10}}^{-\mathrm{3}}$ g m⁻³). It is mainly due to the stronger collection of cloud water and rain water by ice in P3-2ICE, which may be associated with the greater cloud water and rain water content at −10 ^∘C in P3-2ICE than P3-1ICE (Fig. 12c and e). The aforementioned collection of cloud water and rain water by ice does not increase the ice particle number. Although the deposition nucleation can increase the ice particle number in P3-1ICE and P3-2ICE, it is small (less than 0.5 m⁻³ s⁻¹) at −10 ^∘C. Merging ice categories does not increase the total ice particle number in P3-2ICE. The sedimentation of ice number is ∼3.3 and ∼3.4 m⁻³ s⁻¹ at −10 ^∘C in P3-1ICE and P3-2ICE respectively, which dominates the net ice number concentration tendencies at −10 ^∘C in both P3-1ICE (∼2.8 m⁻³ s⁻¹) and P3-2ICE (∼3.1 m⁻³ s⁻¹). The much lower number concentration in P3-2ICE than P3-1ICE (Table 1) is likely due to aggregation associated with collection between the two ice categories in P3-2ICE (Fig. 14).

To summarize briefly, due to the overprediction of LWC in MORR, P3-1ICE, and P3-2ICE above the melting layer, there exist obvious mixed-phase processes at −10 ^∘C. The IWC at −10 ^∘C increases mainly due to the collection of liquid water by ice particles, which does not increase ice particle number but increases the mass and size of ice particles. The lower ice particle numbers in the simulations could also be associated with excessive aggregation and/or missing SIPs, such as collision-induced breakup and “freezing-drop-shattering” proposed by Korolev et al. (2020). The large ice particles and lower ice particle numbers lead to mass concentrating at the large end of PSD contributing to strong simulated radar reflectivity. Introduction of parameterizations for the missing SIPs may be able to overcome some of these model limitations, as will be examined in the future.

Figure 15 As Fig. 7 but for 3 km (d01) simulations.

5 Summary and conclusions

A tropical oceanic convective system observed on 26 May 2015 during the HAIC-HIWC international field campaign based out of Cayenne, French Guiana, was simulated using the WRF model. Observation data from radiosondes, GOES-13 geostationary satellite, airborne radar, and cloud microphysics instrumentation were used to assess the simulated convective system in terms of the thermodynamic and dynamic environment, storm coverage, evolution and structure, and microphysical properties. The major results are summarized as follows:

1. By comparing simulated and observed soundings, all of the simulations using different microphysics schemes replicate temperature with average bias within 1.6 %, dew-point temperature with average bias within 6 %, wind speed with average bias within 14 %, and wind direction with average bias within 36^∘, with the MORR scheme giving closest agreement with observations.

2. WRF basically reproduces the coverage and evolution of this tropical MCS based on a comparison between simulated and observed brightness temperature with the average bias in storm coverage (brightness temperature <232 K) by $\sim -\mathrm{34.3}$ % in WSM6, ∼30.0 % in MORR, ∼12.9 % in P3-1ICE, and ∼2.3 % in P3-2ICE. Thus, WSM6 underestimates the storm coverage, and P3-2ICE produces the closest storm coverage to the observation.

3. In general, all of the simulations overestimate the intensity and spatial extent of radar reflectivity above the melting layer compared to the observed X-band radar reflectivity. There are 95 % of the simulated reflectivities <44 dBZ in WSM6, <41 dBZ in MORR, <45 dBZ in P3-1ICE, and <47 dBZ in P3-2ICE above 6 km, while 95 % of the observed reflectivities are <30 dBZ above 6 km. The radar reflectivity >0 dBZ in all simulations extends above 14 km, whereas the observed radar reflectivity >0 dBZ is mainly below 14 km.

4. Different microphysics schemes have different shapes, magnitudes, and spreads in the simulation of ice PSDs at the different temperature levels. All of the simulations miss the peak of the observed ice number distribution function for $\mathrm{0.1}<D_{\max }<\mathrm{1}$ mm.

5. Both the WSM6 and MORR schemes underestimate the number concentration of ice particles at the temperature levels of −45, −30, and especially −10 ^∘C with a maximum bias up to 1 order of magnitude, though the simulated total water contents are similar to the observations. This indicates that WSM6 and MORR simulate fewer small particles and more large particles compared to the observations. P3-1ICE overestimates the number concentration with a maximum bias up to 1 order of magnitude, and P3-2ICE generates a larger spread of number concentrations covering about 2 orders of magnitude at temperatures of −45 and −30 ^∘C. Both schemes, and especially P3-2ICE, underestimate the number concentration by about 1 order of magnitude at −10 ^∘C. This indicates the P3 scheme produces more large particles at this level.

6. Mixed-phase processes play an important role at −10 ^∘C due to the overprediction of LWC in MORR, P3-1ICE, and P3-2ICE above the melting layer. Stronger simulated radar reflectivity in MORR above the melting layer results from large graupel particles associated with greater graupel water content and fewer graupel particles compared with in situ observations. Rapid growth of graupel mass mainly through collecting cloud water but with a limited increase in graupel number mainly by conversion of snow to graupel through the collection of cloud water above the melting layer leads to large mean graupel sizes in MORR. Similarly, in P3-1ICE and P3-2ICE the IWC at −10 ^∘C increases mainly due to the collection of cloud water and rain water while the net ice number concentration tendencies are near 0, which generates large mean ice particle sizes. The large ice particles generate strong radar reflectivity, partially explaining the bias of simulated radar reflectivity with P3-1ICE and P3-2ICE.

It should be noted that simulations of deep convection at different model resolutions can be much different. To examine the sensitivity of model resolution, CFADs and cumulative CFADs of radar reflectivity are also calculated using the simulation data from the 3 km domain (Fig. 15). Although there are some differences in specific values of reflectivity, the intensity and distribution of reflectivity from the 3 km simulations (Fig. 15) are basically consistent with those from the 1 km simulations (Fig. 7). Although this does not prove that the conclusions with respect to the differences in behavior among the microphysics schemes will be the same at all resolutions, it does at least indicate that the results in this study have some generality. It provides something useful about the microphysics schemes to numerical forecast guidance for HIWCs in current high-resolution NWP models, which are now routinely run at O(3 km) and now more and more often at O(1 km). However, a caveat is that 1 km is not cloud-resolving O(100 m), and thus horizontal entrainment is still not being resolved at this resolution, which affects the amount of liquid available for riming growth in updrafts. Therefore, there still exist uncertainties in the 1 km simulations.

In conclusion, no one microphysics scheme outperforms the other scheme in simulating this tropical oceanic MCS, as is evident from examining the simulated soundings, brightness temperature, radar reflectivity, ice particle size distributions, total water content, and number concentration. The Morrison scheme underestimates the number concentration at different temperature levels compared to the observations. This indicates that large ice particles, especially graupel, are overpredicted in this scheme, which is similar to the Qu et al. (2018) simulation of a different tropical MCS using a different model and microphysics scheme. Even though the P3 scheme has a much different approach for representing ice, it does not produce greatly different TWC or better comparison to the observations using either one- or two-ice categories. This suggests that other aspects need to be considered, such as microphysical process rate formulations or parameters. To enhance understanding of processes leading to the formation of small crystals in HIWC regions, sensitivity experiments varying parameters within the P3 microphysics scheme (e.g., mass-dimensional relations, size distribution parameters, microphysical conversion rates, or representation of different processes like secondary ice production) will be examined in a future paper.

Appendix A: Formulae for X-band radar reflectivity factor

Most cloud models with bulk microphysics parameterizations predict either mass content (one-moment schemes) or mass content and number concentration (two-moment schemes) for a number of hydrometeor categories. In one-moment schemes, number concentration can be obtained diagnostically. They also commonly assume inverse exponential or gamma size distributions of hydrometeors with respect to particle maximum dimension. This allows the use of simple analytical formulae for converting mass content and number concentration to radar reflectivity factor if the scatterers are small compared to the radar wavelength so that the Rayleigh approximation can be used. A gamma size distribution is represented by

$$\begin{array}{}\text{(A1)}& N\left(D\right)=N_{\mathrm{0}}{D}^{\mathit{\mu }}{e}^{-\mathrm{\Lambda }D},\end{array}$$

where N₀, μ, and λ are the intercept, shape, and slope factors respectively. The total number concentration N_t is thus given by

$$\begin{array}{}\text{(A2)}& N_{\mathrm{t}}=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}N\left(D\right)\mathrm{d}D={\displaystyle \frac{N_{\mathrm{0}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}{{\mathrm{\Lambda }}^{\mathit{\mu }+\mathrm{1}}}}.\end{array}$$

A1 Calculation of radar reflectivity factor for rain

For rain, the radar reflectivity factor Z in the Rayleigh approximation is given by the well known formula

$$\begin{array}{}\text{(A3)}& Z=\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}_{\mathrm{eq}}^{\mathrm{6}}N\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}},\end{array}$$

where D_eq is the equivolume diameter of raindrop, and D_eq=D for liquid water species in bulk schemes. The liquid water content LWC is defined as

$$\begin{array}{}\text{(A4)}& \text{LWC}={\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{\mathit{\rho }}_{\mathrm{w}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{D}^{\mathrm{3}}N\left(D\right)\mathrm{d}D,\end{array}$$

where ρ_w is the density of water. Taking into account Eq. (A2), integration of Eqs. (A3) and (A4) yields

$$\begin{array}{}\text{(A5)}& Z={\displaystyle \frac{N_{\mathrm{t}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{7})}{{\mathrm{\Lambda }}^{\mathrm{6}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}},\end{array}$$

and

$$\begin{array}{}\text{(A6)}& \text{LWC}={\displaystyle \frac{\mathit{\pi }{\mathit{\rho }}_{\mathrm{w}}{N}_{\mathrm{t}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{4})}{\mathrm{6}{\mathrm{\Lambda }}^{\mathrm{3}}\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}},\end{array}$$

which results in the following formula for estimating Z from LWC and N_t with an inverse exponential size distribution assumption (μ=0) in WSM6, MORR, and P3:

$$\begin{array}{}\text{(A7)}& Z={\displaystyle \frac{\mathrm{720}{\text{LWC}}^{\mathrm{2}}}{{\mathit{\pi }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{w}}^{\mathrm{2}}{N}_{\mathrm{t}}}}.\end{array}$$

A2 Calculation of radar reflectivity factor for ice particles

The radar reflectivity factor for ice is defined as

$$\begin{array}{}\text{(A8)}& Z={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{2}}}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{\mathit{\rho }}_{\mathrm{s}}^{\mathrm{2}}\left(D_{\mathrm{eq}}\right)D_{\mathrm{eq}}^{\mathrm{6}}N\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}}\end{array}$$

(Ryzhkov and Zrnić, 2019, Eq. 5.14), where ρ_i is the density of solid ice sphere, assumed here to be 0.917 g cm⁻³, and

$$\begin{array}{}\text{(A9)}& K_{\mathrm{w}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{w}}-\mathrm{1}}{{\mathit{\epsilon }}_{\mathrm{w}}+\mathrm{2}}},\text{ and }{K}_{\mathrm{i}}={\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{i}}-\mathrm{1}}{{\mathit{\epsilon }}_{\mathrm{i}}+\mathrm{2}}},\end{array}$$

where ε_w and ε_i are the dielectric constants of water and solid ice, and $|K_{\mathrm{w}}{|}^{\mathrm{2}}=\mathrm{0.930}$ and $|K_{\mathrm{i}}{|}^{\mathrm{2}}=\mathrm{0.176}$ in this study. Taking into account that the mass of ice particle m_s can be expressed as

$$\begin{array}{}\text{(A10)}& m_{\mathrm{s}}\left(D_{\mathrm{eq}}\right)={\displaystyle \frac{\mathit{\pi }}{\mathrm{6}}}{\mathit{\rho }}_{\mathrm{s}}\left(D_{\mathrm{eq}}\right)D_{\mathrm{eq}}^{\mathrm{3}}.\end{array}$$

Equation (A8) can be written in a different form often used in cloud models (e.g., Hogan et al., 2006):

$$\begin{array}{}\text{(A11)}& Z={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{i}}}}\right)}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{m}_{\mathrm{s}}^{\mathrm{2}}\left(D_{\mathrm{eq}}\right)N\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}},\end{array}$$

where D_eq is the equivalent volume diameter of an ice particle (Ryzhkov and Zrnić, 2019).

The m_s–D_max relations are commonly used (e.g., Locatelli and Hobbs, 1974; Mitchell, 1996; Finlon et al., 2019; Ding et al., 2020), where D_max is the maximal dimension of ice particles. These relations are often represented as power-law dependencies

$$\begin{array}{}\text{(A12)}& m_{\mathrm{s}}\left(D_{\max }\right)=a{D}_{\max }^b.\end{array}$$

Then Eq. (A11) can be rewritten as

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}Z& ={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{i}}}}\right)}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{a}^{\mathrm{2}}{D}_{\max }^{\mathrm{2}b}N\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}}\\ & ={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{i}}}}\right)}^{\mathrm{2}}\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{a}_{\mathrm{0}}^{\mathrm{2}}{D}_{\mathrm{eq}}^{\mathrm{2}b}N\left(D_{eq}\right)d{D}_{eq},\end{array}\end{array}$$

where a₀=aη^b. In Eq. (A13), the difference between maximal dimension D_max and equivolume diameter D_eq is taken into account with a scaling factor $\mathit{\eta }=D_{\max }/D_{\mathrm{eq}}$ assuming a constant aspect ratio of ice particles across the size spectrum.

A2.1 Constant m–D relation across the size spectrum

In WSM6 and MORR, a and b in the m_s–D_max relations are assumed to be constant across the size spectrum for snow and graupel particles, implying densities of snow and graupel are constant. This leads to the following expressions for Z and ice water content IWC:

$$\begin{array}{}\text{(A14)}& Z={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_{\mathrm{i}}}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Gamma }(\mathrm{2}b+\mathit{\mu }+\mathrm{1})N_{\mathrm{t}}}{\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}}{\displaystyle \frac{a_{\mathrm{0}}^{\mathrm{2}}}{{\mathrm{\Lambda }}^{\mathrm{2}b}}}\end{array}$$

and

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}\text{IWC}& =\underset{\mathrm{0}}{\overset{\mathrm{\infty }}{\int }}{a}_{\mathrm{0}}{D}_{\mathrm{eq}}^{b}N\left(D_{\mathrm{eq}}\right)\mathrm{d}{D}_{\mathrm{eq}}\\ & ={\displaystyle \frac{\mathrm{\Gamma }(b+\mathit{\mu }+\mathrm{1})N_{\mathrm{t}}}{\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})}}{\displaystyle \frac{a_{\mathrm{0}}}{{\mathrm{\Lambda }}^b}},\end{array}\end{array}$$

so that

$$\begin{array}{}\text{(A16)}& Z={\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\left({\displaystyle \frac{\mathrm{6}}{\mathit{\pi }{\mathit{\rho }}_i}}\right)}^{\mathrm{2}}{\displaystyle \frac{\mathrm{\Gamma }(\mathit{\mu }+\mathrm{1})\mathrm{\Gamma }(\mathrm{2}b+\mathit{\mu }+\mathrm{1})}{\left[\mathrm{\Gamma }\right(b+\mathit{\mu }+\mathrm{1}){]}^{\mathrm{2}}}}{\displaystyle \frac{{\text{IWC}}^{\mathrm{2}}}{N_{\mathrm{t}}}}.\end{array}$$

It is important that Eq. (A16) is not sensitive to the variability of the prefactor a in the m_s(D_max) power-law relation (A12) and is minimally affected by the variability of the exponent b in Eq. (A12). In WSM6 and MORR, b=3 and μ=0, and thus Eq. (A16) can be simplified as

$$\begin{array}{}\text{(A17)}& Z={\displaystyle \frac{\mathrm{720}}{{\mathit{\pi }}^{\mathrm{2}}{\mathit{\rho }}_{\mathrm{i}}^{\mathrm{2}}}}{\displaystyle \frac{|K_{\mathrm{i}}{|}^{\mathrm{2}}}{|K_{\mathrm{w}}{|}^{\mathrm{2}}}}{\displaystyle \frac{{\text{IWC}}^{\mathrm{2}}}{N_{\mathrm{t}}}}.\end{array}$$

A2.2 Variable m–D relation across the size spectrum

P3 represents each ice category much differently than WSM6 and MORR. In P3, a and b in the m_s–D_max relations (Eq. A12) are variable across the size spectrum (Morrison and Milbrandt, 2015). Therefore, radar reflectivity factor is calculated numerically by using Eq. (A13).