The impacts of secondary ice production on microphysics and dynamics in tropical convection

. Secondary ice production (SIP) is an important physical phenomenon that results in an increase of ice particle concentration and can therefore have a significant impact on the evolution of clouds. In this study, idealized simulations of a 15 mesoscale convective systems (MCS) was conducted using a high-resolution (250-m horizontal grid spacing) mesoscale model and a detailed bulk microphysics scheme in order to examine the impacts of SIP on the microphysics and dynamics of a simulated tropical MCS. The simulations were compared to airborne in situ and remote sensing observations collected during the High Altitude Ice Crystals - High Ice Water Content (HAIC-HIWC) field campaign in 2015. It was found that the observed high ice number concentration can only be simulated by the models which include SIP processes. Inclusion of SIP processes 20 in the microphysics scheme is crucial for the production and maintenance of high ice water content observed in tropical convection. It was shown that SIP can enhance the strength of the existing convective updrafts and result in the initiation of new updrafts above the melting layer. Agreement between the simulations and observations highlights the impacts of SIP on the maintenance of tropical MCSs in nature and the importance of including SIP parameterizations in models.

The primary effect of SIP is the enhancement of ice particle concentration which, depending on environmental conditions, may exceed the concentration of PIP ice particles by several orders of magnitude (e.g., Hobbs and Rangno, 1985;Ladino et 35 al., 2017). Such an enhancement of ice particle concentration may have a significant effect on the phase composition, cloud dynamics, precipitation rate, and cloud radiative properties, impacting the energy balance and hydrological cycle on regional and global scales.
At present, seven mechanisms are recognized as sources of secondary ice in clouds. These include the fragmentation of freezing droplets (hereafter FFD) (e.g., Kleinheins et al. 2021), rime splintering (i.e., the Hallett-Mossop process, hereafter 40 HM) (e.g., Hallett and Mossop, 1974), fragmentation due to ice-ice collisions (e.g., Vardiman 1978;Takahashi et al. 1995), ice fragmentation due to thermal shock (e.g., Dye and Hobbs, 1968), fragmentation of sublimating ice (Oraltay and Hallett, 1989), activation of INPs in transient supersaturation around freezing drops (e.g., Prabhakaran et al., 2020), and break-up of freezing water drops on impact with ice particles (James et al. 2021). A detailed description of the first six SIP mechanisms and the status of associated laboratory studies are discussed in the review of Korolev and Leisner (2020). It was found that 45 HM and FFD are the most experimentally studied SIP mechanisms, and in which production rates of secondary ice have been quantified. However, a detailed analysis of previous experiments by Korolev and Leisner (2020) revealed a large diversity of the ice production rates, which led to the conclusion that these SIP processes need to be studied further. The other four mechanisms have a limited number of laboratory experiments, and cover only a fraction of environmental conditions (e.g., fragmentation during ice collisions, fragmentation of sublimating ice), or only demonstrated the general feasibility of SIP 50 mechanisms (e.g., fragmentation due to thermal shock, activation of INPs in transient supersaturation around freezing drops).
All these led Korolev and Leisner (2020) to the conclusion that the relative contributions of each of the six SIP mechanisms in the enhancement of ice concentrations remain uncertain.
For the last few years, there were many new efforts on systematic studies of the effect of SIP on cloud microphysics with the help of cloud simulations (e.g., Phillips et al. 2017aPhillips et al. , 2018Sullivan et al. 2018;Hoarau et al. 2018;Fu et al. 2019;55 Sotiropoulou et al. 202055 Sotiropoulou et al. , 2021Dedekind et al. 2021;Hawker et al. 2021;Huang et al., 2021, 2022 andothers). Most of these modeling efforts were focused on matching simulated moments of particle size distributions (PSDs) with those observed in situ. In many ways, the implementation of SIP in numerical models was hindered by the lack of consensus on parameterizations of SIP mechanisms.
One of the main objectives of this work is to identify and simulate the occurrence of high ice water content (IWC) associated 60 with the enhancement of ice particle concentrations from SIP processes. Cloud environments with high IWC (> 1 g m -3 ) pose a hazard for civil aviation and may result in engine power loss, stall, or damages (e.g., Lawson et al. 1998;Mason et al. 2006, Mason andGrzych, 2011). The phenomenon of high IWC is well documented from in situ observations in tropical mesoscale convective systems (MCSs) (e.g., Heymsfield and Palmer 1986, Lawson et al. 1998, Gayet et al. 2012, Fridlind et al. 2015Leroy et al. 2017, Strapp et al. 2021. Several previous modeling studies using different cloud microphysical parameterizations 3 attempted to reproduce high IWCs. Ackerman et al. (2015) used a 1D model to explore microphysics in tropical MCSs.
Simulations performed with 3D models (Franklin et al. 2016;Stanford et al. 2017;Qu et al. 2018) pointed to the inaccuracies in the estimation of cloud PSD, IWC, and ice category comparing to the observations. Huang et al. (2021) conducted highresolution simulations of tropical convection and found significant overestimates of radar reflectivity and underestimates of total ice crystal concentration (Ni). Adding SIP in high-resolution simulations, Huang et al. (2022) found significant 70 improvement of simulated Ni compared to the in situ observations.
In most of the previous numerical studies investigating SIP, the microphysics schemes used were based on the traditional approach of representing ice-phase hydrometeors whereby they are partitioned into various predefined categories (e.g., pristine ice, snow, graupel, etc.) with prescribed physical properties. This approach has several inherent limitations and problems, including a limited range of ice properties (e.g., bulk density) that can be represented, inconsistent physical processes applied 75 to the categories, and the need to parameterize conversion between categories, an artificial process which can not be constrained from observations and is purely ad hoc. To address this problem, Morrison and Milbrandt (2015) proposed a new approach and developed a new microphysics parameterizationthe Predicted Particle Properties (P3) scheme -whereby all ice-phase hydrometeors are represented by a single "free" ice category whose physical properties evolve continuously. While flexible in this regard, one limitation of the original P3 scheme was that it could not represent more than one population of ice 80 particles (with different bulk properties) at a given time and grid location. The scheme was thus generalized to allow for a user-specified number of "free" ice categories, each of which have properties that evolve continuously and can represent any ice type (Milbrandt and Morrison 2016).
The P3 scheme was used in the tropical convection simulations of Huang et al. (2022). While they found significant improvement of simulated Ni compared to the in situ observations by adding three SIP mechanisms, their simulations were 85 limited to two ice categories. This is the minimum number of categories required for including SIP processes, since at least two categories are needed to represent the co-existence of newly formed small ice splinters and pre-existing large ice particles.
However, as will be shown below, the use of more than two ice categories may be beneficial or even necessary to model the impacts of SIP in deep convection.
This study is focused on the examination of the effects of SIP on the microphysics and dynamics of a simulated tropical 90 MCS. Quasi-idealized simulations of a MCS were conducted using a near cloud-resolving configuration (250-m horizontal grid spacing) of a 3D dynamical model with the P3 microphysics scheme. Model configurations with up to four free ice categories were tested. The enhancement of ice particle concentration by SIP is represented by HM and FFD mechanisms. In the absence of a consensus on SIP parameterizations, these two processes were described by two specific parameterizations proposed in the literature, which provide a sufficient enhancement of Ni above the melting layer consistent with in situ 95 observations in the MCSs (Korolev et al. 2020). Simulated ice PSD, Ni, IWC, radar reflectivity, and Doppler velocity were compared against in situ and remote sensing observations collected during the HAIC-HIWC field campaign in 2015 (Strapp et al. 2021). Without looking for exact match between model simulations and observations, this study aimed to show whether the simulation with SIP produces better estimation of the observed microphysics compared to the simulation without SIP.
The remainder of the paper is structured as follows. The next section describes the observation data used for evaluation. In 100 section 3, the setup of the model, the microphysics scheme, and the parameterizations of SIP are detailed. Section 4 describes the choice of control simulation with regard to the number of ice categories in P3. Section 5 assesses the impact of SIP on the formation of ice clouds based on the control simulation. The role of SIP in strengthening and sustaining tropical convection is discussed in section 6. This is followed by an assessment of the impact of ice-ice collection efficiency on the simulation. The final section offers a summary of the study and conclusions. 105

Observation data
In situ data employed in this study were collected from the National Research Council of Canada (NRC) Convair-580 and The measurements of PSDs were performed by three particle probes, which covered different particle size ranges. The software employed a retrieval algorithm of partially viewed particle images (Heymsfield and Parrish 1979;Korolev and Sussman 2000), which allowed the enhancement of particle statistics and extended the maximum size of the composite PSD up to 12.8 mm.
All particle probes were equipped with anti-shattering tips to mitigate the effect of ice shattering on the measurements of ice particle concentration. Residual shattering artifacts were identified and filtered out with the help of the inter-arrival time 120 algorithm (Field et al. 2006;Korolev and Field 2015).
The bulk IWC was measured by an Isokinetic Probe (IKP: Davison et al. 2008). The IKP allowed measurements of IWC up to 10 g m 3 at the aircraft speed 200 m s 1 . Such a high upper limit of IWC well exceeded maximum IWC (~5 g m 3 ) measured during the HAIC-HIWC campaign and ensured that the measured IWC never exceeded the IKP saturation level.
Both aircraft were equipped with the same instruments for measurements of PSDs and bulk IWC and used the same 125 processing algorithms applied to these measurements. Such arrangement minimized differences in systematic errors specific to different types of instruments and synchronized data processing. Therefore, if any potential biases in data acquisition and data processing existed, they would be the same for both data sets collected from the NRC Convair-580 and SAFIRE Falcon-20.
Besides comparisons with IWC and Ni, model results are also compared with reflectivity and Doppler velocity measured 130 by the NRC aircraft X-band radar (NAX) installed on the NRC Convair-580 (Wolde and Pazmany 2005). Statistics of the NAX data included the MCS cloud segments that precipitated down to the ground surface level. Cloud segments with outflow cirrus having radar returns disconnected from the ground were excluded from the statistics.

Atmospheric model and initialization 135
The model used in this study is the Global Environmental Multiscale (GEM) model (Côté et al. 1998;Girard et al. 2014).
GEM is used for operational numerical weather prediction (NWP) in Environment and Climate Change Canada (ECCC) as well as research in ECCC and Canadian universities. The dynamical core of GEM is formulated based on the non-hydrostatic fully-compressible primitive equations with a terrain-following hybrid vertical grid. As such, it can be run at cloud-resolving (sub-km grid spacing) scales. It can be run on global or limited-area domains and is capable of one-way nesting. In this study, 140 an idealized model configuration was used to simulate tropical deep convection, with a horizontal grid spacing of 250 m in a simulation domain of 160 km  160 km, with 83 vertical levels over a tropical ocean surface. The horizontal grid spacing of 250 m is nearly at the cloud-resolving scale (Bryan et al. 2003;Lebo and Morrison, 2015). It is also close to the corresponding distances (110 to 180 m) of the 1 Hz in situ observation data from the aircraft which flew mostly at 110-120 m s -1 for Convair-580 and 150-180 m s -1 for Falcon-20. To resolve the vertical profiles near the tropical melting layer, vertical grid spacings of 145 approximately 100 m were used between the altitudes of 4 and 7.5 km. The readers are referred to Table 1 for more details of dynamics/numerics and physics configurations used in the model.
The atmospheric initial conditions were horizontally homogeneous, based on an initial sounding taken from the operational global GEM analysis at 12 UTC on 15 May 2015 at 6.769 N and 49.551 W. The initial profile ( Fig. 1) had 1697 J kg 1 of convective available potential energy (CAPE). The GEM analysis also provided the initial sea surface temperature. The 150 location and time were chosen based on the occurrence of an extensive mesoscale system that formed in this region and was observed ( Fig. 2) during the HAIC -HIWC field campaign (Strapp et al. 2021).
To initiate the model storm, the updraft nudging method of Naylor and Gilmore (2012) was used to force convection during the first 15 min of the simulation. Three distinct updrafts 15 km apart from each other in the western part of the simulation domain were initialized. Each updraft was forced by perturbing the vertical air velocity ( ) in a spheroid with a horizontal 155 radius of 10 km and vertical radius of 1.5 km, centered at 1.5 km altitude: where β is the distance from the center of the spheroid normalized by its radius, α is an inverse nudging timescale (0.5 s 1 ), dts is the model time step, is the maximum updraft speed (10 m s 1 ) for nudging.

Cloud microphysics scheme
All cloud microphysical processes in the GEM simulations were represented by the P3 two-moment bulk microphysics scheme, where up to four (free) ice categories were used. For each ice category, there are four prognostic (i.e. advected) mixing ratio variables: the total ice mass, the rime mass, the bulk volume, and the total number. From the prognostic variable fields, various bulk physical properties can be computed. The size distribution of each category is represented by a complete gamma 165 function. The liquid-phase component of P3 consists of two-moment categories for cloud droplets and rain. For details on the representation of each hydrometeor category and the parameterized microphysical processes, readers are referred to Morrison and Milbrandt (2015); for details on the multi-category configuration, see Milbrandt and Morrison (2016).
It should be noted that there have been several further key developments to P3 since the first version of the multi-icecategory scheme, along with various minor modifications. Major developments include a triple-moment treatment of rain 170 (Paukert et al. 2019), introduction of a prognostic liquid fraction for wet ice (Cholette et al. 2019), and a triple-moment treatment of ice (Milbrandt et al. 2021). The version of P3 used in this study does not include these major modifications. The impacts of these components of P3 on SIP may be examined in future work. It is possible, for example, that triple-moment ice, which in principle results in a better representation of the PSD dispersion, may be important for some aspects of modeling SIP and its impacts. However, such work is outside of the scope of this study. 175

Parameterization of SIP
The following two SIP mechanisms were examined in this study. Other mechanisms will be considered in the future.

Rime splintering/Hallet-Mossop (HM)
The pristine version of P3 used in this study (i.e., prior to SIP-related modifications examined here) includes a parameterization of the HM mechanism if two or more ice categories are used. The requirement for at least two categories is 180 to prevent dilution of ice particle properties when two populations of ice are forced to be represented by a single size distribution and set of physical properties (see Milbrandt and Morrison 2016). The parameterized HM process produces a maximum of 350 ice crystals per mg of collected liquid water, with crystal sizes of 10 µm, during riming of rain within a temperature range of -3°C > T > -8°C, with the peak value at -5°C and varying linearly to 0 at the extreme temperature ranges.
This ice multiplication parameterization has been used in several traditional (fixed ice category) microphysics schemes (e.g., 185 Reisner et al. 1998). This is similar to the original HM parameterization in P3 as used in Milbrandt and Morrison (2016).
One modification for HM parameterization in this study is to exclude the use of collected liquid water from the situation where raindrops (100 µm < D < 3500 µm) were collected/nucleated by small ice particles (D < 100 µm). From the point of view of SIP mechanism, it is more appropriate to apply this part of collected liquid water to the FFD mechanism.
The parameterization of the FFD mechanism was implemented following Lawson et al. (2015): where Nf is the average number of ice fragments per drop, and D is the drop diameter in micrometers. The parameterization of the FFD process was applied for raindrops (100 µm < D < 3500 µm) which were nucleated by ice particles (D < 100 µm).
Following Keinert et al. (2020) the activity of the FFD process was limited to the temperature range -25°C < T < -2°C. Within 195 the temperature range, the current parameterization is not temperature dependent. Further studies are undergoing for exploring the impacts of variation due to temperature.

Ice collection efficiency
Although not directly part of SIP, ice aggregation is another process that impacts Ni and may therefore be important in affecting high IWC. There are several different approaches to parameterize the ice collection efficiency, which is a key 200 parameter in the aggregation parameterization (Hallgren and Holster 1960;Lin et al. 1983;Cotton et al. 1986;Ferrier et al. 1994Ferrier et al. , 1995Milbrandt and Yau 2005a,b;Seifert and Beheng 2006). Khain and Pinsky (2018) showed that the collection efficiencies in these parameterizations vary by more than two orders of magnitude for any given temperature between -40°C and 0°C. Sensitivity simulations conducted in the frame of this current work showed a high sensitivity of the modeled IWC and Ni to the collection efficiency of ice. The collection efficiency used in this study follows Cotton et al. (1986). Further 205 discussion of the sensitivity of the modeling results to the ice aggregation parameterization is presented in section 7.

Establishment of the control configuration
The number of ice categories used in the P3 scheme can impact the overall simulation results (Milbrandt and Morrison, 2016). SIP results in large quantities of small ice splinters, which can be co-located with pre-existing larger ice particles, and thus bimodal or multi-modal ice size distributions may occur. This cannot be represented with the single-ice-category 210 configuration since P3 uses complete gamma size distributions for each hydrometeor category. Thus, SIPor any other ice initiation process -could result in the "dilution" of the bulk particle properties of existing ice where the microphysics scheme tries to represent two or more populations of ice particles within a single size distribution and with a single set of bulk physical properties (e.g., mean size). Thus, before examining the impacts of SIP in the simulation, a control configuration must be established based on the minimum number of ice categories needed to represent ice particle evolution in P3 with sufficient 215 detail. This is determined by the number of categories beyond which adding more does not change the simulated average profiles for more than 15% for the fields of interest (e.g., IWC and Ni). To address this a series of sensitivity tests was conducted using the baseline version of P3 (with no SIP), varying the number of ice categories from 1 to 4, and using P3 with SIP included, varying the number of categories from 2 to 4. Table 2 summarizes the complete list of experiments conducted in this study. For the experiments starting with "BASE" 220 (denoting the baseline P3 configuration), no SIP is activated. The experiments starting with "SIP" use the new parameterization for both FFD and HM process. In the following analysis, we focus on the simulation results between 90 and 150 min, when the convective systems are more complex and more closely resemble the observed MCS.
As described in section 3, three convective cells were initiated in the eastern part of the simulation domain. Figure  The initial formation of the three convective updrafts can still be seen at 30 min (Fig. 3a). By 60 min, these updrafts started to merge, forming a larger system (Fig. 3b). This system then moved westward (towards the right of the domain) and developed into a sustained system (Fig. 3c). By 180 min the convection began to weaken (Figs. 3d, h). For the altitude range between 5 and 12 km, adding one more ice category to one, two and three-category baseline simulation produce maximal changes of 12%, 23% and 7% for IWC respectively (Fig. 4a). Similarly, the maximal changes of 25%, 31% and 14% are found for Ni (Fig. 4b). The radar reflectivity (Fig. 4f) for the one-ice category run (BASE-1ICE) is about 4 to 6 dBZ lower than the three or four-category simulations in the same altitude range. This is likely caused by the fact that a single ice category is not sufficient to represent the co-existence of large and small ice particles and results in reduction 240 of the concentration of large ice particles, leading to lower radar reflectivity. With regard to the number of ice categories for SIP simulations, a similar conclusion to the baseline simulations was found. Adding one more ice category to two and threecategory SIP simulations produce maximal changes of 31% and 9% for IWC respectively (Fig 5a). For Ni, the maximal changes of 70% and 15% are found (Fig 5b). Therefore, at least three ice categories in P3 appear to be necessary and sufficient to examine the impacts of including SIP processes. Note that adding more ice categories for baseline and SIP simulations will 245 not change significantly the morphology of the storm. It is of passing interest to note that the similarity of the three and fourice category results is consistent with the 1D kinematic simulations in Milbrandt and Morrison (2016). However, given that the use of more ice categories in P3 is generally preferable in principle (though there is added computational cost), and that the four-ice-category simulations were already performed, BASE-4ICE is taken as the control run for the sensitivity studies to follow; this simulation is referred to as CTR. Correspondingly, we focus on the four ice category simulation including SIP 250 (SIP-4ICE) for direct comparison to CTR. Figure 6 shows simulated average profiles for the four-ice category simulation including SIP processes, SIP-4ICE (see Table 2) as well as for CTR. As seen in Fig. 6a, the SIP-4ICE simulation has at least 100% higher IWC compared to the control 255 run above 6 km. SIP-4ICE has significantly higher Ni than CTR (Fig. 6b), with differences reaching two to three orders of magnitude near 6 km and about one order of magnitude above 11 km. Figure 6c shows the maximum vertical air velocity in the domain. The simulations are very similar below the melting layer (~4.5 km), however, above the melting layer, wmax of SIP is 2 to 5 m s -1 higher. This suggests that SIP enhances convection due to the sudden production of a large number of small ice particles resulting in the rapid freezing of rainwater and depletion 260 of water vapor by diffusional growth. Both effects result in latent heating, which invigorates the convection. This can be inferred from Fig. 6e showing that RWC in CTR is reduced form 0.12 g m -3 at 4.2 km to 0.05 g m -3 at 5 km, whereas RWC form SIP-4ICE is changed from 0.12 g m -3 at 4.2 km to 0.01 g m -3 at 5 km. Another potential mechanism of convection enhancement above the melting layer will be discussed in section 6.

Domain-averaged profiles
The medians of radar reflectivity of the SIP simulation is 5 to 10 dBZ lower compared to those of the CTR between the 265 altitudes of 5 and 12 km (Fig. 6f). This is because the SIP simulation has smaller ice particles, despite the higher IWC values, due to the higher Ni.
In order to confirm that the simulation differences illustrated in Fig. 6 are indeed the result of the inclusion of parameterized SIP processes, and not simply a chance set of changes that could result from perturbations in the model due to some minor code change, a set of 10 ensemble simulations were run for each configuration (CTR and SIP-4ICE) but with perturbed initial 270 conditions. Figure 7 shows the ensemble results for the two configurations at 120 min. Each configuration includes 11 members (1 unperturbed + 10 perturbed members). The initial temperature profile is randomly perturbed with the maximum range of ±1°C at all model levels. The ensemble results show high consistency. The SIP simulation consistently produces much higher IWC and Ni. The vertical velocity wmax of an ensemble member is the single maximal value of a give model level. Therefore, the profiles of wmax can be noisy. The ensemble profiles of wmax of CTR and SIP-4ICE overlap, especially below the altitude 275 of ~6 km. However, the averaged maximal vertical velocity wmax (solid lines) diverges above the altitude of ~6 km, which indicates that the SIP simulations generate stronger updrafts in general. The RWC and the radar reflectivity are both strongly reduced by 5 to 10 dBZ in the SIP simulations between 5 and 10 km. Similar remarks can be made for all other times between 90 and 150 min (not shown). The results summarized in Fig. 7 lend support to the idea that the differences between CTR and SIP-4ICE are indeed due to the effects of SIP on the microphysical and thermodynamic fields, and not merely another model 280 realization of a chaotic weather system. Figure 8 shows comparisons of the probability density function of Ni, F(Ni), calculated from the simulations and measured during airborne in situ observations at two different altitude (H) ranges. The Ni measured in situ is calculated for size range between 40 μm and 12.5 mm. The ice particles smaller than 40 μm is not counted due to large uncertainty of the instrument 285 for the size range (Baumgardner et al. 2017). The F(Ni) measured in situ were averaged over all HAIC-HIWC flights for the clouds with IWC > 0.01 g m -3 . Figure 8a shows a comparison of F(Ni) at altitudes 6 < H < 7 km. The CTR simulation (blue line) shows a significant underestimation of Ni compared to the measured values (black line). The concentrations corresponding to the modal value of F(Ni) in CTR are nearly two orders of magnitude lower than those obtained from in situ observations. The SIP-4ICE simulation (red line) shows good agreement with the measured values. 290

Ice number concentration
The general behaviour of the functions F(Ni) for 11 < H < 12 km (Fig. 8b) is similar to that obtained for 6 < H < 7 km in

Ice water content
Similar to the results discussed above, Fig. 9 shows comparison of probability density functions of IWC, F(IWC), obtained from model simulations and aircraft observations at altitudes 6 < H < 7 km (Fig. 9a) and 11 < H < 12 km (Fig. 9b). As seen in 300 Fig. 9a, F(IWC) from SIP-4ICE is in good agreement with the observations for IWC smaller than 3.25 g m 3 . In contrast, the simulated frequency of encountering high IWC in CTR is about ½ to 1/500 of the observed frequency between IWC of 1 and 2 g m -3 . There is no data with IWC higher than 2.5 g m -3 in CTR.
SIP-4ICE produces some points with IWC > 3.25 g m 3 , which were not observed by the instruments. The F(IWC) of these high IWC conditions from SIP-4ICE are below 1.7 10 -5 . For the Convair-580 aircraft, the number of observed 1-s average 305 data points with IWC > 0.01 g m 3 is 59,893. This sets the limit of F(IWC) of the observation data at 1.7  10 -5 . If the campaign lasted much longer, it is conceivable that these high IWC conditions might have eventually been observed. Figure 9b shows similar results but for higher altitudes (between 11 and 12 km). The frequency of IWC between 0.3 and 1.3 g m -3 in CTR is about 2 to 3 orders of magnitude lower than the observed frequency. There is no data with IWC higher than 1.30 g m -3 from CTR. SIP-4ICE produces closer estimates compared to the observation data as they both have IWC up to 310 ~3 g m -3 . Between 11 and 12 km, Ni of SIP-4ICE is considerably improved compared to CTR as shown in Fig. 8b. However, the IWC of SIP-4ICE at these altitudes is still underestimated by 1 to 2 orders of magnitude beyond IWC of 0.7 g m -3 compared to the observation. One possible reason for this underestimate is the differences in the sampling of data. At higher altitudes between 11 and 12 km, there are often extensive areas with thin ice clouds near the convection. The data from SIP-4ICE used in the statistics include these areas if the IWC of the grid cell larger than 0.01 g m -3 . In contrast, the HAIC-HIWC campaign 315 targeted conditions with high IWC. The thinner ice clouds with IWC between 0.01 and 0.3 g m -3 might not have been sufficiently sampled as they are less relevant to the extreme conditions causing safety issues for aviation. This difference might partly explain the higher F(IWC) below 0.3 g m -3 and lower F(IWC) above 0.3 g m -3 for SIP-4ICE compared to the observations. Another possible explanation of the underestimation is the uncertainty of the strength of simulated convections.
In this study, the maximal updraft nudging speed =10 m s -1 is used as the default value. Simulations with different 320 are also tested. Using =15 m s -1 in SIP-4ICE simulation will produce 63% higher averaged IWC between 10 and 11 km than that produced by default SIP-4ICE simulation with =10 m s -1 . For SIP-4ICE simulation, the impact of for lower altitudes between 6 and 7 km is negligible. The CTR simulation with =15 m s -1 produces higher IWC compared to the default CTR simulation with =10 m s -1 (20% and 50 % higher for altitudes range of 6-7 km and 10-11 km respectively). However, these values are still significantly smaller than those of the default SIP-4ICE simulation. 325 The corresponding simulated radar reflectivity of the cross-section indicated by black lines in Fig. 10 is shown in Fig. 11. One significant difference between the CTR (Figs. 10a-c) and SIP-4ICE simulation (Figs. 10d-f) is that the reflectivity from the simulation with SIP is significantly lower than that of the control run between altitudes of approximately 5 and 10 km. This is due to the higher Ni values and thus smaller ice particle sizes in SIP-4ICE. 335 Figure 12 shows a comparison of the frequency distribution of radar reflectivity for CTR and SIP-4ICE (Figs. 12a, b) and for the NRC Convair-580 X-band radar (Fig. 12c). For the results of the two simulations, only the atmospheric columns with both ice water path and rain water path larger than 1 g m -2 are selected. The X-band radar data in Fig. 12c was averaged over all research flights during the HAIC-HIWC campaign. Figure 12d shows the simulated and observed median values of the reflectivity. At altitudes higher than 10.3 km reflectivity in both the CTR and SIP-4ICE runs is lower than the measured 340 reflectivity. However, between 5 and 10 km, SIP-4ICE has values closer to the observations with a maximum overestimation of 4 dBZ at 5 km. CTR clearly overestimates the reflectivity by 5 to 15 dBZ between 5 and 10 km.

Longwave radiation and radar reflectivity
For SIP-4ICE between 5 and 8 km, there is a large amount of cases with reflectivity larger than 25 dBZ which is not found in the observation. Figure 13 shows the comparison of Ni and IWC against reflectivity for CTR (panel a, b), d) and observation (e, f). Although the distributions of Ni and IWC of SIP-4ICE are closer to the observation than those of CTR, a significant amount of cases (35.8%) from SIP-4ICE has still reflectivity larger than 25 dBZ, whereas only 0.4% from the observation is above 25 dBZ. Figure 14 shows the 2D histograms for IWC and reflectivity (panel a-d), Ni and reflectivity (e-h) for each of the 4 ice categories in SIP-4ICE. The ice category 1 and 3 have in general lower reflectivity (1% and 3% higher than 25 dBZ respectively) and higher Ni (peak value at 10 4.4 and 10 4 m -3 respectively). By contrast, the ice category 2 and 4 have higher 350 reflectivity (23% and 24% higher than 25 dBZ respectively) and lower Ni (peak value at 10 3.4 and 10 2.8 m -3 respectively). This means that the ice particles in category 2 and 4 have much larger mean sizes than those of category 1 and 3. Most of the ice particles in category 2 and 4 with reflectivity larger than 25 dBZ are lightly rimed aggregated ice with rime fraction of ~20%.
The mean-mass diameter of these ice particles could reach several mm. It is probably that the high reflectivity above 25 dBZ are the results of these large ice sizes. 355 Several possible reasons could explain the high reflectivity in the model simulation. Firstly, the model might overproduced the large ice particles due to excessively large ice-ice collection efficiency (which is a highly uncertain parameter). Despite the important role of ice-ice collection process in clouds, there is no consensus in the scientific literature on how to quantify it. A further discussion on the impact of ice-ice collection efficiency on the simulation results is presented in the section 7.
Secondly, current P3 version uses a diagnostic shape parameter for the gamma size distribution of ice. The shape parameter, 360 which is a measure of the relative spectral dispersion, might give a wide distribution that implies too many large ice particles, hence higher reflectivity. Note, higher moments such as reflectivity are more sensitive to the tail of the size distribution. The triple-moment version of P3 (Milbrandt et al. 2021, JAS) uses a prognostic shape parameter which in principle should give better size distributions for those large ice particles. However, this comparison is beyond the scope of the current study, since the triple-moment version of P3 was not yet available when this current study began. Finally, P3 scheme uses Rayleigh 365 scattering approach to calculate the radar reflectivity. It is not an advanced instrument simulator that can take into account the attenuation, multiple and Mie scattering, etc. For the X-band radar on Convair-580, the ice particles of several mm are in a transition region between Rayleigh and Mie scattering which is not symmetric but with a stronger forward scattering lobe. To simulate the radar reflectivity with Rayleigh scattering approach might overestimate the reflectivity for the large ice particles.
The radar reflectivity is a useful indicator to understand the model performance. However, considering that a rigorous 370 instrument simulator is not used in this study, it is better to be cautious with regard to the comparison with the reflectivity.Note that at the altitude of the melting layer (~4.5 km) none of the simulations reproduce the distinct bright band that is clearly apparent in the observation data. This is due to the fact that the version of P3 used in this study does not properly represent the transition state of melting ice with a wet surface and ice core (nor does P3 artificially boost the reflectivity contribution from ice during melting in order to mimic a bright-band effect). As mentioned in section 3.3, a newer version of P3 includes a 375 prognostic variable for the liquid mass content for each ice-phase category, which allows for mixed-phase particles and a corresponding improvement in the calculation of reflectivity in the melting zone (to be shown in a forthcoming publication).

Vertical Doppler velocity
The Doppler velocity from the simulation is calculated using mass weighted fall speed of all hydrometeors subtracted from the vertical wind speed. Figure 15 shows the simulated Doppler velocity for CTR and SIP-4ICE for the vertical cross section 380 shown in Fig. 11. The Doppler velocity below the melting layer is mostly negative due to the high fall speed of the rain. Above the melting layer the average Doppler velocity gradually increases with altitude from -2 to 0 m s -1 between 5 and 14 km. The gradual increase of the Doppler velocity is primarily linked to the size of ice particles. Positive values of the Doppler velocity are associated with convective cloud regions where the updraft velocity exceeds the falls speed of the ice particles. One difference between the control and SIP simulations is that the SIP-4ICE simulation has higher Doppler velocity between 5 and 385 8 km for all three times analyzed. Figure 16 shows similar results to Fig. 12 but for the Doppler velocity. To mitigate the effect of anvils the diagrams in Figs. 16a and b used the same mask as that for radar reflectivity in Figs. 12a and b. Both CTL (Fig. 16a) and SIP-4ICE (Fig.   16b) show distribution patterns very similar to those obtained from the observations (Fig. 16c). Figure 16d shows comparisons of the simulated and observed median values of the Doppler velocity versus altitude. The SIP-4ICE simulation produces very 390 close results to the measurements between the altitude of 5 and 9 km with maximum difference of ±0.3 m s -1 . The CTL simulation overestimates the negative Doppler velocity by approximately 0.9 m s -1 in the same range of altitudes compared to the measured values. This result is consistent with the systematic underestimation of Ni in CTR, and therefore, overestimation of the mean particle sizes and fall speeds.

Role of SIP in tropical convection 395
In the previous section it was shown that wmax is higher in the SIP-4ICE simulation than in the control simulation above the melting layer (Fig. 6c), and that the inclusion of SIP results in enhanced formation of high IWC regions above the melting layer. Altogether these results suggest that SIP plays an important role in the microphysics and thermodynamics of tropical convection, at least for the MCS examined. In order to explore SIP impacts in more detail, in this section, the rates and locations of SIP within the model storm and initiation of secondary convection are analyzed. 400 Figure 17 shows the active SIP areas and rates at different altitudes for the simulated MCS. Similar to Hu et al. (2021), the values of w are used to distinguish three different situations: 1) SIP within updrafts (marked by red lines, w > 3 m s -1 ), 2) SIP within downdrafts (marked by yellow lines, w < 3 m s -1 ), and 3) SIP with moderate vertical wind velocities (marked by blue lines, -3 m s -1  w  3 m s -1 ). 405

SIP productions at different altitudes
For the employed SIP parameterizations, FFD is active in the range of altitudes 5 < H < 8.8 km (Fig. 17d), whereas HM occurs in a narrower range of 5.2 < H < 6.8 km (Fig. 17a). The range of the SIP activation altitudes are primarily determined by the temperature ranges of the FFD and HM processes (section 2) and the temperature of the cloud base. The vertical velocity has a lesser effect on the SIP activation altitudes, and it may displace the upper and lower boundaries of the SIP regions within approximately 200 m. Within these altitude ranges, SIP is mostly found in the area with moderate vertical wind velocities, 410 followed by the area within updrafts, and occurs less in downdrafts for both HM and FFD processes (Fig. 17a, d). The active SIP areas of FFD are usually smaller than those of HM process in the range of altitudes 5.5 < H < 6.2 km for all three situations.
To activate the FFD process, the employed parameterization requires the presence of large size raindrops (100 µm < D < 3500 µm) together with the presence of small ice particles (D < 100 µm). This condition likely restricts the FFD process to a smaller area than that of HM process. 415 The vertical wind speed has a significant impact on the rate of SIP. Figure 17e shows the average SIP rates (m -3 s -1 ) within the active FFD areas shown in Fig. 17d. For most altitudes, the rate of FFD increases with an increase of w and it is at least one order of magnitude higher in updrafts than in the downdrafts. This is related to the lower number of precipitation-sized drops in downdrafts compared to updrafts.
Another factor is related to the effect of w on the residence time of drops within the range of SIP activation altitude ∆ . 420 For a drop with a terminal fall velocity ( ) the residence time can be assessed as = abs(∆ /( ( ) − )) .
Depending on the sign of ( ( ) − ) the drop will be moved upward through ∆ or downward. For the extreme situation when ( ) = the drop is suspended in an updraft indefinitely, it can freeze and generate secondary ice or mechanically interact with other cloud particles and thereby change its fall velocity.
The rate of the HM process (Fig. 17b) is higher in updrafts above 6.0 km where the riming process is active. At altitudes 425 below 6.0 km, the rates are similar in updrafts, downdrafts and areas with moderate vertical velocities (-3 m s -1  w  3 m s -1 ). Figure 17c shows the total SIP rate (m -1 s -1 ) from the HM process which is the product of the area of active HM and mean SIP rate from HM shown in Figs. 17a,b. The total SIP rate shows how many ice particles are produced by SIP horizontally across the domain per 1 m of vertical layer per second. The total SIP rate from the FD process is shown in Fig. 17f. Below the altitude of 6.3 km, both the HM and FFD processes in the areas with moderate vertical wind velocities show the highest total SIP rate, 430 followed by the area within updrafts. The lowest total rates are found in downdrafts. The larger active SIP areas associated with moderate vertical wind velocities (Fig. 17a, d) contribute significantly to the high total SIP rates across the domain. At altitudes above 6.3 km, updraft regions contribute more to the total SIP rates. This is due to the high average rates in updrafts (Fig. 17b, e), since the corresponding SIP areas are smaller than those with moderate vertical velocities.
The results obtained show that the vertical extent ∆ of FFD is deeper and its rate is higher than those of HM. This finding 435 leads to the conclusion that in the simulations of this study the overall contribution of FFD in the production of the secondary ice in tropical MCSs is significantly higher than HM.

Role of SIP in initiating of secondary convection
The freezing of rain into ice generates latent heating which should enhance the existing convection. This may explain the increase of wmax above 6 km in the Figs. 6c and 7c in SIP-4ICE compared to CTR. As explained in the previous subsection, below 6.3 km altitude, both FFD and HM are more active outside of the major updrafts originating below the melting layer.
High activity of the FFD and HM processes might eventually initiate new updrafts in stratiform regions inside MCSs above the melting areas.
A Lagrangian trajectory analysis was used to trace the cloud parcels affected by SIP. For this analysis, 1152 air parcels were selected with active SIP at an altitude of 5.6 km from the SIP-4ICE simulation between 90 and 150 min. Each selected 445 parcel is traced backward (t < 0 min) and forward (t > 0 min) for 15 min (t = 30 min). These parcels are then classified into two different groups. The first group included the parcels which at t = -15 min had altitudes within 5 < H < 6 km and at t = 15 min ending their trajectories at H > 6.5 km. The second group included parcels with the same initial altitudes as the first group at t = 0 min. However, the altitude of their trajectories remained at H < 6 km at the end of forward tracing (t = 15 min). The total number of parcels of the first category was 47 and that of the second category was 1105. 450 Figure 18 shows the time history of mean values of environmental and microphysical parameters of the parcels for the two categories. Most parcels, which underwent SIP processes between -5 and 5 min, were located at the same altitude of 5.6 km at t = 0 min. These parcels started to rise from t = -3 min, and eventually reached 7.5 km at t = 15 min (Fig. 18c). On the other hand, the other group of parcels did not rise throughout the 30 min analysis period (hereafter named non-rising parcels). The mean altitude of these parcels decreased by about 700 m (Fig. 18c). 455 Figure 18e shows that the rising parcels had an initial positive vertical speed from t = -3 to -1 min. However, their potential temperature differences with respect to the environmental values at the same altitudes (Δθ) were generally negative for t < -1 min, and decreased slightly from t = -3 to -1 min (Fig. 18g). Thus, these parcels were not gaining positive buoyancy at this stage. However, the Δθ of the rising parcels started to increase quickly between t = 0 min and 2 min, becoming positive and reaching a difference of nearly +1 K between t = 3 and 8 min (Fig. 18g). Figures 18i, j also show very high SIP rates of the 460 rising parcels from the FFD process during this time period. For the non-rising parcels, there was an increase in Δθ during the same period but at a much slower rate due to smaller SIP tendency (Figs. 18i, and j). The Δθ remained negative during the whole analysis period, and these parcels were therefore convectively stable.
The main reason for the high SIP rate between t = 0 and 2 min for the rising parcels is that there was a substantial amount of large rain drops available for activating the FFD process (Figs. 18f, and h). The rain mean-mass diameter at t = 0 min was 465 large (560 μm). The RWC was also large (0.36 g m -3 ). The rain mean-mass diameter for the non-rising parcels (502 μm) was slightly smaller than that of rising parcels at t = 0 min. However, the corresponding RWC was quite low (0.02 g m -3 ). With large RWC and rain mean-mass diameter, the rising parcels had a high SIP potential, which eventually led to greater Ni, vapor growth, increased latent heating and buoyancy, thereby enhancing secondary convection.
As mentioned in subsection 6.1, in the simulations of this study, the FFD process plays a dominant role in SIP compared 470 to the HM process. This agrees with what we found in Figs. 18i, j. The sudden increase of the Δθ with respect to the environment is more likely linked to the high SIP rate from FFD. SIP, in particular FFD, may therefore play a role in the initiation of new updrafts above the melting layer. Figure 19 shows an example of a rising air parcel. The black line represents the parcel trajectory from t = -15 to 15 min.
From t = -15 to 0 min, the air parcel had no significant change of altitude. Near t = 0 min, the air parcel was close to an existing 475 updraft (Fig. 19a, red surface) and was in an area where there was rainwater (Fig. 19b, red surface). Shortly after t = 0 min, the air parcel started to rise. The supply of rainwater resulted in an enhanced SIP and led to a higher rate of latent heating and rapid increase the buoyancy.

Impact of ice aggregation
In addition to SIP, which clearly has significant effects on Ni, formation of high IWC, radar reflectivity, and vertical wind 480 velocity, the aggregation of ice particles also plays an important role in determining the microstructure. Aggregation results in a decrease of Ni and increase of radar reflectivity and particle fall velocity. The rate of aggregation is characterized by the iceice collection efficiency, eii. In the frame of this study, the following parameterization of the ice-ice collection efficiency has been employed (Cotton et al. 1986): = min (10 0.035( −237.16)−0.7 , 0.2) (4) 485 where T is the temperature in K.
There is a diversity of parameterizations of employed in models (Khain and Pinsky, 2017), which may vary by up to three orders of magnitude. In this study, the mid-range in (4) is used as the default parameterization. However, the uncertainty in raises a question about the impact of the ice aggregation parameterization on the Ni and high IWC formation.
To explore the effect of the ice aggregation rate on the high IWC formation, a sensitivity test (SIP-COL) was performed 490 with another parameterization, i.e.
As seen from Eq. (5) the new varied linearly from 0.1 to 1.0 within the temperature range -20° < T < -5°C. For T < -20°C, = 0.1, and for -5°C < T < 0°C = 1.0. For all temperatures, in (5) is much higher than that in (4). We want to use this high parameterization to show the significant impact of on the simulation results. 495 Comparisons of distributions of Ni for two SIP simulations with varying and the observations are shown in Fig. 20. For the altitudes between 6 and 7 km, using the new parameterization of aggregation [given by Eq. (5)] (SIP-COL) results in a decrease of modal value of F(Ni) by two orders of magnitude compared to using that described by Eq. (4) (SIP-4ICE). At higher altitudes between 10 and 11 km, the SIP-COL produces ~5 times higher F(Ni) at Ni of 10 5 m -3 and lower F(Ni) by one to two orders of magnitude for Ni > 4×10 5 m -3 compared to those of SIP-4ICE. 500 Figure 21 shows similar results to Fig. 20 but for the distribution of ice mass. Between 6 and 7 km, applying the linear approach for from Eq. (5) to the SIP-4ICE simulation (SIP-COL), the F(IWC) between 0.4 and 2.4 g m -3 is slightly reduced by up to 50%. However, the F(IWC) for extreme situation with IWC larger than 2.5 g m -3 is somehow slightly enhanced.
For altitudes between 11 and 12 km, the SIP-COL experiment produces a lower F(IWC) up to ~1 order of magnitude for most of IWC (> 0.16 g m -3 ) than the SIP-4ICE simulation. Although the estimation of F(IWC) by SIP-COL between 6 and 7 505 km is relatively close to that of SIP-4ICE, SIP-COL produces much lower F(IWC) in higher altitudes between 10 and 11 km than the SIP-4ICE simulation. This indicates that the Ni at the lower altitudes play an important role in determining the IWC at upper altitudes.

Conclusion
The impacts of SIP on the microphysics and dynamics of deep convection have been examined using quasi-idealized near 510 cloud-resolving simulations of a tropical MCS based on storm observations during the HAIC-HIWC field campaign. GEM model simulations using the P3 microphysics scheme were conducted using 250-m horizontal grid spacing and horizontally homogeneous atmospheric initial conditions, with updraft nudging to initiate convection. It was established through sensitivity tests that a minimum of three ice categories in P3 are necessary to examine SIP in detail; four categories were used for most of the simulations. P3 was modified to include rime splintering (HM) and fragmentation of freezing drops (FFD), which have 515 been the most closely examined SIP mechanisms in laboratory studies. The parameterizations of the HM and FFD processes used were based on the information available from previously published results.
In the control configuration with no SIP processes at altitudes of 6 to 7 km, the simulated ice number concentrations were about two orders of magnitude lower than the values obtained from in situ measurements. The simulated frequency of encountering high IWC condition is about ½ to 1/500 of the observed frequency between IWC of 1 and 2 g m -3 . With the SIP 520 mechanisms activated, the model results for these fields were dramatically improved compared to the observations. The Doppler velocities above the melting layer were also notably closer to the measured values, indicating improved ice fall speeds in the simulations with SIP active. SIP was responsible for an increase in ice concentrations of up to three orders of magnitude at altitudes of 6 to 7 km. As a result, the total ice mass was distributed over a much larger number of particles and thus mean particle size was smaller with a lower fall speed. Consequently, ice was more easily transported to higher altitudes, ultimately 525 resulting in sustained cloud regions with high IWC.
Analysis conducted in this study lead to the following general conclusions for the high resolution NWP simulations:

1.
SIP processes play a critical role in the formation and maintenance of high IWC with low reflectivity at upper levels in MCSs.
2. SIP enhances secondary convection above the melting layer due to an increase in buoyancy caused by greater 530 latent heat releasing during vapor deposition on numerous secondary ice particles. Enhanced secondary convection may in turn extend the longevity of MCSs and regions with high IWC.

3.
Aggregation of ice particles results in a decrease of ice number concentration and IWC at upper levels but is very sensitive to details of the parameterization of this process, in particular the collection efficiency, which remains uncertain. 535 In order to minimize errors in interpretation of the results due to unresolved convective updrafts, the simulations conducted in this study were all done with a horizontal grid spacing of 250 m. This is a much higher resolution than current operational numerical weather prediction (NWP) models. However, tests with 1 km grid spacing (Fig. 22) indicated that impacts of including SIP are very similar to those at 250-m grid spacing, where 1 km is close to the grid spacing of several current operational and experimental NWP systems. Further, the P3 microphysics scheme is already used operationally in the 540 Canadian 2.5-km system (Milbrandt et al. 2015). Therefore, the conclusions regarding the importance of including SIP processes in models are not limited to numerical modelling in research mode, but also have important implications for current and/or upcoming operational NWP, in particular for systems that provide numerical guidance for civil aviation operating at cruising altitudes between 10 and 14 km.
Finally, although the simulations conducted with the activated SIP process clearly resulted in improved results compared 545 to the observations, this is not a basis for concluding that the HM and FFD parameterizations used are accurate representations of these physical processes. While the formulations were based on either laboratory experiments or combined modelling and in situ observational study, they are still largely ad hoc. This study further highlights the importance of these processes in deep convection and the need to include them in some fashion in numerical models. However, accurate parameterizations that capture the underlying physics of these mechanisms, not just their bulk effects, continue to be topics of research. 550

Authors Contribution
ZQ, AK and JAM conceptualized the research goals and aims. ZQ, JAM and AK designed the experiments with the support 555 from YH, GMM and HM. ZQ performed the simulations and analysis with the help from AK, JAM, IH, MW and CN. ZQ prepared the manuscript with contributions from all co-authors.