Microphysical processes producing high ice water contents (HIWCs) in tropical convective clouds during the HAIC-HIWC field campaign: dominant role of secondary ice production

Abstract. High ice water content (HIWC) regions in tropical deep convective clouds, composed of high concentrations of small ice crystals, were not reproduced by Weather Research and Forecasting (WRF) model simulations at 1 km horizontal grid spacing using four different bulk microphysics schemes (i.e., 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) for conditions encountered during the High Altitude Ice Crystals (HAIC) and HIWC experiment. Instead, overestimates of radar reflectivity and underestimates of ice number concentrations were realized. To explore formation mechanisms for large numbers of small ice crystals in tropical convection, a series of quasi-idealized WRF simulations varying the model resolution, aerosol profile, and representation of secondary ice production (SIP) processes are conducted based on an observed radiosonde released at Cayenne during the HAIC-HIWC field campaign. The P3 two-ice category configuration, which has two “free” ice categories to represent all ice-phase hydrometeors, is used. Regardless of the horizontal grid spacing or aerosol profile used, without including SIP processes the model produces total ice number concentrations about 2 orders of magnitude less than observed at −10 ∘C and about an order of magnitude less than observed at −30 ∘C but slightly overestimates the total ice number concentrations at −45 ∘C. Three simulations including one of three SIP mechanisms separately (i.e., the Hallett–Mossop mechanism, fragmentation during ice–ice collisions, and shattering of freezing droplets) also do not replicate observed HIWCs, with the results of the simulation including shattering of freezing droplets most closely resembling the observations. The simulation including all three SIP processes produces HIWC regions at all temperature levels, remarkably consistent with the observations in terms of ice number concentrations and radar reflectivity, which is not replicated using the original P3 two-ice category configuration. This simulation shows that primary ice production plays a key role in generating HIWC regions at temperatures <-40 ∘C, shattering of freezing droplets dominates ice particle production in HIWC regions at temperatures between −15 and 0 ∘C during the early stage of convection, and fragmentation during ice–ice collisions dominates at temperatures between −15 and 0 ∘C during the later stage of convection and at temperatures between −40 and −20 ∘C over the whole convection period. This study confirms the dominant role of SIP processes in the formation of numerous small crystals in HIWC regions.


been proposed to represent the fragmentation due to collision of ice particles. In these methods, the number of ice fragments 50 produced per collision is simply fit to the observations of Takahashi et al. (1995) (Sullivan et al., 2017(Sullivan et al., , 2018b, set to a constant value (Hoarau et al., 2018) or dependent upon the initial kinetic energy and colliding particles' size and rimed fraction (Phillips et al., 2017a, b). The shattering of freezing droplets, the first proposed SIP mechanism, generates small ice splinters following droplet freezing, in which a closed ice shell is formed, freezes inward, and subsequently shatters (Langham and Mason, 1958;Mason and Maybank, 1960;Lauber et al., 2018). Lauber et al. (2018) summarized previous laboratory studies to 55 show that overall fragmentation frequency and the number of ejected splinters per fragmenting droplet increase with increasing droplet diameter. Although the splinter production rates vary with temperature among the laboratory experiments conducted by different research groups, most laboratory studies showed that the freezing droplet fragmentation may be most efficient at temperatures below −10 • C, outside the temperature range of H-M process (Mason and Maybank, 1960;Brownscombe and Thorndike, 1968;Takahashi andYamashita, 1969, 1970;Kolomeychuk et al., 1975;Pruppacher and Schlamp, 1975;Takahashi, 60 1975;Lauber et al., 2018;Keinert et al., 2020).
Several numerical studies investigated the roles of different SIP mechanisms in ice particle production in different types of clouds in different regions (Phillips et al., 2017a(Phillips et al., , 2018Sullivan et al., 2018a;Fu et al., 2019;Sotiropoulou et al., 2020Sotiropoulou et al., , 2021. Phillips et al. (2017a) indicated that the average ice number concentration at temperatures between ∼0 and −30 • C increased by one to two orders of magnitude with the inclusion of fragmentation in ice-ice collisions in a cloud-resolving model simu-have indicated that HIWC phenomenon cannot be captured well by numerical models (Franklin et al., 2016;Stanford et al., 85 2017;Qu et al., 2018;Huang et al., 2021). Huang et al. (2021) (Morrison and Milbrandt, 2015;Milbrandt and Morrison, 2016). All of their simulations overestimated the radar reflectivity and underestimated the number concentration of 90 ice particles in HIWC regions compared to the observations. They hypothesized that these biases could be attributed to the poor representation of SIP processes in the microphysics schemes. As a companion paper of Huang et al. (2021), the roles of different SIP mechanisms, namely the H-M mechanism, shattering of freezing droplets, and fragmentation of ice-ice collisions, in the formation of numerous small crystals in HIWC regions are investigated in the current study through a series of sensitivity experiments with the P3 microphysics scheme. 95 The next section describes the implementation of SIP mechanisms in the P3 scheme and sensitivity experiments conducted in this study. The results of the sensitivity experiments are discussed in Section 3, followed by a summary and conclusions presented in section 4.

Methodology
2.1 Implementation of SIP parameterizations 100 Milbrandt and Morrison (2016) expanded the P3 scheme to include multiple "free" ice categories, in which particle populations with different sets of bulk properties are allowed to coexist and the detrimental effects of bulk property dilution, where information from particles' different growth paths is lost due to a single set of bulk properties, can be reduced. Their simulation results indicated that at least two ice categories are required to correctly represent the rime splintering process and reduce the bulk property dilution effects. Therefore, the P3 scheme with two ice categories (P3-2ICE) is adopted in this study. The three The parameterization of freezing-droplet shattering is implemented using the numerical formulation proposed by Phillips et al. (2018), which combined observations from previous laboratory studies and considered the physics of collisions. Two modes of the scheme, fragmentation during heterogeneous drop freezing (mode 1) and accretion of raindrops (mode 2), are 120 considered in this study. The number of fragments per frozen drop in mode 1 is dependent on raindrop size and freezing temperature, and it is dependent on collision kinetic energy in addition to raindrop size and freezing temperature in mode 2 (Phillips et al., 2018). The raindrop-freezing fragmentation scheme is implemented in P3 by adopting a bin-emulating approach, in which bulk particle size distributions are discretized into bins for the calculations of microphysical process rates (Saleeby and Cotton, 2008;Morrison, 2012). A more detailed description of this parameterization is found in Phillips et al. (2018).

125
The ice-ice collection or aggregation (collision and coalescence) process was considered in the original P3 scheme, but fragmentation during ice-ice collision was not. The physically-based parameterization of ice multiplication by breakup during ice-ice collision proposed by Phillips et al. (2017b) is adopted and implemented in the P3 scheme using a bin-emulating approach. The scheme is based on an energy conservation principle, in which the number of new fragments per collision is dependent on collision kinetic energy, temperature, and colliding particles' size and rimed fraction. Parameters in the scheme 130 are estimated based on previous laboratory and field experiments, with more details found in Phillips et al. (2017b). The collection (aggregation) efficiency (E agg ) between ice particles follows the laboratory study of Connolly et al. (2012), in which E agg are 0.09, 0.21, 0.6, 0.1, 0.08, 0.02 at temperatures of −5, −10, −15, −20, −25, −30 • C, respectively. When temperature > −5 • C and < −30 • C, E agg is set to 0.09 and 0.02 respectively, and otherwise E agg is linearly interpolated between temperatures. As with most bulk schemes, the collision efficiency between ice particles is assumed to be 1, implying 135 E agg is equal to the coalescence efficiency (E coal ). Therefore, ice-ice collision breakup efficiency is equal to 1 − E coal . Field et al. (2006) indicated that a constant aggregation efficiency of 0.09 (E agg = 0.09) produced good agreement with aircraft observations, and Morrison and Grabowski (2010) assumed E agg = 0.1 in their study. The main results and conclusions do not change in sensitivity experiments using a constant E agg = 0.1. Therefore, only results using E agg following Connolly et al. (2012) are shown in this paper.

Numerical experiments
Idealized experiments that consume less computing resources are conducted first, and then the optimal configurations from these idealized studies are used to rerun a real-case experiment to examine whether changes to the default P3 scheme can improve the simulation of HIWC phenomenon.

145
The WRF version 4.1.3, as used by Huang et al. (2021), is employed in a three-dimensional quasi-idealized framework to simulate the tropical oceanic convection observed on 26 May 2015. The input sounding used for the initial horizontally-uniform thermodynamic environment is from a radiosonde released at Cayenne at 00:00 UTC 26 May 2015 (Fig. 1a). The sounding has a deep moist absolutely unstable layer and mainly easterly (westerly) winds below (above) 350 hPa. The surface-based convective available potential energy of the sounding is 2378 J kg −1 .

150
The model domain is 200 × 100 km 2 with horizontal grid spacings between 125 and 1000 m and the model top is 18 km with 71 vertical levels. The model time step is 1 s. Three dimensional subgrid-scale mixing is calculated using a 1.5-order turbulent kinetic energy scheme (Skamarock et al., 2019) instead of a planetary boundary layer (PBL) scheme. An ocean surface is assumed and the surface moisture and sensible/latent heat fluxes are estimated using the MM5 similarity surface layer scheme (Jiménez et al., 2012). The other physical parameterization schemes, including longwave and shortwave radiation scheme, 155 land-surface scheme and cumulus parameterization scheme, are not activated in the idealized simulations. The P3 two-ice microphysics scheme is adopted and the detailed setups are described in Section 2.2.2.
Updraft nudging (Naylor and Gilmore, 2012) is adopted to initiate convection within the horizontally-uniform thermodynamic environment. The updraft (w t ) within a spheroid with 10-km horizontal radius (x r = y r = 10 km) and 1.5-km vertical radius (z r = 1.5 km) centered at z c = 1.5 km is determined by with where ∆t is the small model time step, α = 0.5 s −1 , w max = 10 m s −1 , β =  . At the −10 • C level, Ni/IWC covers three orders of magnitude between 10 3 and 10 6 g −1 and ∼53.9% and ∼45.8% of samples are between 10 4 and 10 5 g −1 and between 10 5 and 10 6 g −1 , respectively. With an increase of either upward or downward vertical velocity, Ni/IWC increases, passing the significance test for p < 0.05 ( Fig. 3a). At the −30 • C level, Ni/IWC 210 covers two orders of magnitude between 10 4 and 10 6 g −1 , and ∼85.4% of samples are between 10 5 and 10 6 g −1 (Fig. 3b).
With an increase of either upward or downward vertical velocity, Ni/IWC increases at −30 • C (Fig. 3b, passing the significance test for p < 0.05), but the slope is less than that at −10 • C (Figs. 3a and b). At the −45 • C level, ∼98.4% of Ni/IWC samples are between 10 5 and 10 6 g −1 (Fig. 3c). Ni/IWC does not appear to increase with vertical velocity at −45 • C (Fig. 3c, not passing the significance test for p < 0.05) in contrast to results at temperatures of −10 and −30 • C. At the −10 • C level, the simulations excluding SIP processes produce Ni/IWC mainly covering two orders of magnitude between 10 2 and 10 4 g −1 (Figs. 4a1-d1), which is about two orders of magnitude less than observed (Fig. 3a). With an increase of upward vertical velocity or decrease of downward vertical velocity, Ni/IWC has a decreasing trend (Figs. 4a1-225 d1), which is also different from the observations (Fig. 3a). In addition, the radar reflectivities in these simulations are mainly greater than 35 dBZ at −10 • C (Figs. 4a1-d1), which is obviously overestimated compared to the observations in HIWC regions where 95% of the cumulative observed reflectivities are < 30 dBZ (Huang et al., 2021). Therefore, at −10 • C the experiments without SIP processes using any horizontal grid spacing (Figs. 4a1-c1) or any aerosol profile (Figs. 4b1 and d1) cannot produce HIWC regions consistent with observations. There are no obvious differences among these simulations, although the number 230 concentration of cloud droplets is reduced by ∼74.5% in the NoSIP250m experiment using the aerosol profile based on UHSAS observations (not shown), which is closer to the observed.
Overall, the simulations without SIP processes underestimate Ni/IWC and overestimate radar reflectivity at temperatures of −10 and −30 • C, that is, they cannot produce HIWC regions at these temperature levels. These simulations can produce HIWC not passing the significance test for p < 0.05), which is inconsistent with the observations (Fig. 3a). About 72% of radar reflectivities in HM250m are greater than 30 dBZ, which are overestimated compared to observed in HIWC regions at −10 • C, where 95% of the cumulative observed reflectivities are < 30 dBZ (Huang et al., 2021). In RFZB250m, Ni/IWC covers three orders of magnitude between 10 3 and 10 6 g −1 (Fig. 5b1), which is consistent with the observations (Fig. 3a). However, RFZB250m does not produce the observed relationship between Ni/IWC and downward vertical velocity, which is similar to 260 HM250m. About 63% of radar reflectivities in RFZB250m are less than 30 dBZ (Fig. 5b1). Thus, RFZB250m can produce HIWC regions at −10 • C to a certain extent. In IICB250m (Fig. 5c1), Ni/IWC covers four orders of magnitude between 10 2 and 10 6 g −1 , which is underestimated compared to observations especially at larger vertical velocities. Meanwhile, IICB250m fails to capture the observed relationship between Ni/IWC and upward vertical velocity. However, IICB250m produces ∼30% of samples with HIWC characteristics, that is, high Ni/IWC > 10 5 g −1 and radar reflectivities < 20 dBZ (Fig. 5c1). In SIPs250m 265 that includes all three SIP mechanisms (Fig. 5d1), Ni/IWC covers the same range as observed (i.e., between 10 3 and 10 6 g −1 ), and the increase of Ni/IWC with greater upward or downward vertical velocity is also captured well. Around 96% of radar reflectivities in SIPs250m are less than 30 dBZ, indicating that SIPs250m produces HIWC regions at −10 • C remarkably consistent with the observations.

285
To further examine the role of SIP mechanisms in different locations of the convective storm, we analyze a vertical cross section through the convective core (based on the maximum composite reflectivity) of microphysical process rates relevant to ice particle production including: the H-M mechanism, shattering of freezing droplets, fragmentation of ice-ice collisions, and other microphysical processes (i.e., primary ice nucleation, homogeneous and heterogeneous freezing of cloud droplets and rain). Results from SIPs250m are shown in Fig. 6 for regions with IWC > 1 g m −3 . The H-M process (mainly at −5 • C) and 290 shattering of freezing droplets (mainly at temperatures between −5 and −20 • C) dominate ice particle production (> 58%) in the strong updraft core regions where there is plentiful LWC. Fragmentation during ice-ice collisions is dominant (∼100%) in the other HIWC regions (Fig. 6). In general, total ice particle production rates are about 4 times larger in the strong updraft regions (w >10 m s −1 ) than those in other HIWC regions. The importance of freezing fragmentation enhanced by updrafts is consistent with an observational study on mixed-phase clouds at temperatures > −10 • C in the Arctic (Luke et al., 2021).  Because SIP processes need to be triggered by preexisting ice and the ice-ice collision process is strongly dependent on Ni, the relative contribution of SIP processes to ice particle production should be different at different stages of convection. Figure 7 shows the time evolution of the microphysical process rates relevant for ice particle production including the H-M mechanism, shattering of freezing droplets, fragmentation of ice-ice collisions, and other microphysical processes (i.e., primary ice nucleation, homogeneous and heterogeneous freezing of cloud droplets and rain) in regions with IWC > 1 g 305 m −3 at different temperatures in SIPs250m. It indicates that the roles of SIP processes in ice particle production in HIWC regions vary during the evolution of convection. At the early stage of convection (t < 40 min), primary ice production (mainly homogeneous freezing of cloud droplets) dominates ice particle production (> 50% of total ice particle production rate) at temperatures less than −40 • C, fragmentation of ice-ice collisions is dominant (> 50% of total ice particle production rate) at temperatures between −40 and −20 • C and at 0 • C, and shattering of freezing droplets plays the key role (> 50% of total ice 310 particle production rate) at temperatures between −15 and −5 • C (Fig. 7). With the development of convection, Ni increases, and the fragmentation of ice-ice collisions becomes dominant (> 50% and maximum close to 100% of total ice particle production rate) at temperatures less than 0 • C. The H-M process also plays a role (∼5% of total ice particle production rate) https://doi.org/10.5194/acp-2021-781 Preprint. Discussion started: 27 October 2021 c Author(s) 2021. CC BY 4.0 License.
in the ice particle production around −5 • C. Therefore, primary ice production is dominant in HIWC regions at the very early stage of convection at temperatures less than −40 • C, shattering of freezing droplets dominates ice particle production in HIWC regions at temperatures greater than −15 • C during the early stage of convection, and fragmentation during ice-ice collisions is dominant at temperatures greater than −15 • C during the later stage of convection and at temperatures less than −20 • C over the whole convection period.

Improvement in real-case simulation
To examine whether the new P3 two-ice scheme including all three SIP mechanisms can improve the simulation of HIWC g −1 , it covers three orders of magnitude between 10 3 and 10 6 g −1 with ∼75.5% of Ni/IWC values < 10 5 g −1 (Fig. 8a3). From the observed radar reflectivity shown in Fig 10 5 and 10 6 g −1 in P3-2ICE_SIP and observations, respectively (Fig. 8b3). There are ∼85.4%, ∼93.0%, and ∼99.1% of the simulated radar reflectivities in P3-2ICE_SIP < 30 dBZ at −10 • C, < 20 dBZ at −30 • C, and < 15 dBZ at −45 • C, respectively, (Figs. 8b1-b3), which is very consistent with the observed. These results are also consistent with those in the quasi-idealized simulation SIPs250m (Fig. 5d1-d3). Therefore, the real-case simulation using the new P3 two-ice scheme including all three SIP mechanisms can successfully reproduce the HIWC regions of the observed tropical oceanic convective system. It also 345 confirms the dominant role of SIP processes in HIWC regions with high concentration of small ice crystals.
A previous study (Huang et al., 2021) used the WRF model at 1-km horizontal grid spacing with four different bulk microphysics schemes to simulate tropical deep convective clouds observed during the HAIC-HIWC field campaign. The simulations overestimated the intensity and spatial extent of radar reflectivity above the melting layer and failed to reproduce the observed 350 high concentrations of small ice crystals in HIWC regions, in which there are numerous small ice crystals with MMDs of 200-300 µm, Z e often < 20 dBZ, and IWCs often > 1.5 g m −3 . To explore formation mechanisms for HIWC regions and biases in the WRF simulations, a series of quasi-idealized sensitivity experiments on the model resolution, aerosol profile, and SIP processes are conducted based on an observed sounding from a radiosonde released at Cayenne during the HAIC-HIWC field campaign. The P3 two-ice scheme, which has two "free" ice categories to represent all ice-phase hydrometeors, is used.

355
The main results are summarized as follows: (1) By comparing simulations to observations, regardless of the horizontal grid spacing (1 km, 250 m and 125 m) or aerosol profile used (default constant profile in original P3 scheme or aerosol profile based on UHSAS measurements from HAIC-HIWC field campaign), without SIP processes the model produces total ice number concentrations about two orders of magnitude less than observed at −10 • C and about an order of magnitude less than observed at −30 • C. These simulations also 360 overestimate the radar reflectivity at −10 and −30 • C. Although the simulations can produce HIWC regions at −45 • C, they overestimate the ice number concentration compared to observations.
(2) Three simulations turning on one of three SIP mechanisms separately (i.e., the Hallett-Mossop mechanism, fragmentation during ice-ice collisions, and shattering of freezing droplets) can produce higher ice number concentrations at −10 and −30 • C, but they do not fully replicate observations of HIWCs, with the results of the simulation with shattering of freezing droplets 365 most closely resembling the observations.
(3) The simulation including all three SIP processes successfully produces HIWC regions at all temperature levels in terms of ice number concentration and radar reflectivity. Based on a vertical cross section of ice particle production rates through the mature convection, the H-M mechanism (mainly at −5 • C) and shattering of freezing droplets (mainly at temperature between −5 and −20 • C) dominate ice particle production in strong updraft core regions where there is plentiful LWC, and   I  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I ,  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I : l • • : : : : . . . , : : : . .  I  I  I  I  I  I  I  I  I  I  I  I  I  I  I Fig. 2c of IWC (gray contours: 0.05, 1, and 3 g m −3 from thin to the thick), LWC (green contours: 0.05, and 1 g m −3 from thin to the thick), vertical velocity (vertical vectors), and the microphysical process rates relevant for ice particle production processes including (red) H-M mechanism, (blue) shattering of freezing droplets, (orange) fragmentation of ice-ice collision, and (magenta) other microphysical processes (i.e., ice nucleation, homogeneous and heterogeneous freezing of cloud droplets and rain) in regions with IWC > 1 g m −3 at different temperatures in SIPs250m at t = 60 min. The pie charts denote ice particle production rates summed over all the source terms with the area of each color proportional to the ice particle production rate. . Time evolution of the averaged microphysical process rates relevant for ice particle production processes including (red) H-M mechanism, (blue) shattering of freezing droplets, (orange) fragmentation of ice-ice collision, and (magenta) other microphysical processes (i.e., ice nucleation, homogeneous and heterogeneous freezing of cloud droplets and rain) in regions with IWC > 1 g m −3 at different temperatures. The pie charts denote ice particle production rates summed over all the source terms with the area of each color proportional to the ice particle production rate.  Fig. 1 of Huang et al. (2021)), during the Cayenne field campaign, and the observed temperature ranges of samples at the three levels are −12.6 to −7.9 • C, −30.4 to −29.7 • C, −44.7 to −43.6 • C, respectively. The simulation using the original P3 two-ice scheme is shown in a1-a3, and the simulation using the P3 two-ice scheme including the three SIP processes is shown in b1-b3. The simulations at the three temperature levels are interpolated from the 1-km model outputs. The scatters of simulations are color-coded according to the magnitude of radar equivalent reflectivity factor (dBZ).