A Numerical Simulation Study of Secondary Ice Productions in a Squall Line Case

: Secondary ice productions (SIPs) can produce ice crystals with a number concentration much higher than that of ice nucleating particles in mixed-phase clouds and therefore inﬂuence cloud glaciation and precipitation. For midlatitude continental mesoscale convective systems (MCSs), how SIPs affect the microphysical properties and precipitation is still not clear. There are few studies of SIPs in midlatitude continental MCSs. This study investigates the roles of three SIPs (rime splintering, freezing drop shattering, and ice-ice collisional breakup) on a squall line case in North China on 18 August 2020 using the WRF model with a modiﬁed Morrison double-moment bulk microphysical scheme. Including SIPs, especially ice-ice collisional breakup, in the model simulations markedly improves the simulated convective area and convective precipitation rate of the squall line, while slightly improving the area and precipitation of the stratiform region. Within the mixed-phase layer in both the convective and stratiform regions of the squall line, ice-ice collisional breakup is the dominant process to generate ice crystals. In contrast, rime splintering generates an order of magnitude fewer ice crystals than ice-ice collisional breakup, while freezing drop shattering plays a negligible role due to the lack of large drops. Ice multiplication through ice-ice collisional breakup and rime splintering produces numerous snowﬂakes and graupel. This leads to enhanced depositional growth and weaker riming, which in turn weakens rime splintering. It is recommended to add SIP parameterization to the model.

Rime splintering occurs when supercooled droplets rime onto graupel or snow particles and subsequently break up to form many small ice crystals in the temperature range of −3 • C to −8 • C [8]. Simultaneous observations of high concentrations of pristine column ice, supercooled droplets, and graupels at the temperature range of −3 • C to −8 • C have been observed in cumulus clouds [26][27][28], weakly convective cells embedded in mixed-phase mid-level stratus and frontal clouds [29][30][31], and the edge of a tropical mesoscale convective system (MCS) with weak updraft in Africa [18].These observations provide strong evidence for rime splintering occurring in relative week convections.In the strong updraft region of deep convection, large amounts of small ice crystals have been observed, but numerical modeling studies [32,33] have shown that the overall contribution of rime splintering is significantly lower than that of freezing drop shattering to the production of secondary ice in tropical MCSs.Numerical simulations by Aleksić et al. [34] showed that rime splintering induced rapid glaciation within a maritime convective cloud but played a less significant role in a continental convective cloud due to the lower liquid water content and narrower droplet spectrum.The modeling work of cumulus clouds in New Mexico found that rime splintering can sufficiently explain the observed rates of glaciation of cumulus clouds [20].Their subsequent research has revealed that increasing cloud condensation nuclei (CCN) number concentration can lead to smaller droplets, weaker rime splintering, and hence slower glaciation [35].But the model simulation of UK summertime cumulus clouds showed that doubled CCN produced more numerous small cloud droplets and fewer graupel, resulting in an increased rime splintering rate [27], as opposed to the result in Phillips, Choularton, Blyth and Latham [35].The different results may be because the two cases have different liquid water content (LWC) and use different parameters to calculate the rate of the rime splintering process.The case in Phillips, Choularton, Blyth and Latham [35], is a continental cumulus and only considers the rime splintering of large droplets (diameters larger than 24 µm), while the case in Huang et al. [27] is closer to the ocean and also takes into account the rime splintering of small droplets (diameters smaller than 13 µm).
Freezing drop shattering is the production of small ice fragments upon the freezing of large drops.At present, there is no clear understanding of the condition for occurrence, but this process tends to occur when there are large drops.Laboratory works have shown that large (80-400 µm diameter) supercooled drops emitted small fragments upon freezing in the temperature range from −5 • C to −30 • C, and larger drops shatter more than the small ones [36][37][38].Freezing drop shattering is usually observed in tropical maritime convective clouds, which have warm bases and large supercooled drops [19,21,39].According to observations by Korolev et al. [19] in four tropical maritime MCSs during their mature stages, large ice particles precipitate and melt in the warm layer to generate large drops.After being recirculated by turbulent updrafts, those large drops collide with the aged ice to freeze and shatter.The observations over the coastal waters of Washington State by Rangno [40] suggest that ice particle concentrations in supercooled stratiform frontal clouds containing drizzle drops (50-500 µm in diameter) can be modestly enhanced (up to 20%) due to the fragmentation of freezing drops where ice is first appearing.Modeling by Sullivan et al. [25] has shown that freezing drop shattering needs a warm cloud base and strong updraft, which tend to help form large supercooled drops.Modeling of tropical maritime convective clouds confirms that the ice number is increased by an order of magnitude and is close to observation when including freezing drop shattering [21,22].
Ice-ice collisional breakup is the generation of multiple ice crystals caused by the collision between ice-phase particles and plays a significant role in clouds with a large number of ice-phase particles.Although recent studies have observed naturally fragmented ice crystals, there is no indication of whether it is caused by ice-ice collisional breakup or by sublimation fragmentation [41].So far, studies often compare the ice crystal number concentrations from observations with those from cloud model simulations that include ice-ice collisional breakup.The results have consistently proved that only the simulations including ice-ice collisional breakup can reproduce the observed number concentration of ice-phase particles in a summer convective line with cold cloud bases over the U.S. high plains [42], in wintertime mixed-phase orographic clouds over the Swiss Alps [43,44], and in the summer clouds over the Antarctic coast [45].These studies also consistently found that ice-ice collisional breakup can produce more than 90% of all nonhomogeneous ice (including heterogeneous nucleation and SIPs) and that the contribution of freezing drop shattering to ice production is negligible.But note that the insignificant effect of freezing drop shattering may be because the parameterizations of drop freezing used in these studies may underestimate the number of freezing drops [46] and therefore may underestimate the secondary ice crystal number concentration from freezing drop shattering.The study of Arctic stratocumulus using a cloud model by Sotiropoulou et al. [47] showed that ice-ice collisional breakup and rime splintering are weak when activated alone because the Arctic is a pristine region (low INP concentration) with a low concentration of ice and low riming rates.Only the combination of rime splintering and ice-ice collisional breakup can explain the observed ice crystal number concentration because the ice crystals produced by rime splintering further enhance ice-ice collisional breakup.
It is generally seen that in different cloud conditions, the dominant SIP may be different.Sullivan et al. [48] investigated the initial conditions and enhancement of the ice crystal number concentration through the three SIPs with a parcel model.The results suggested that the initiation and enhancement of the ice crystal number concentration through rime splintering and freezing drop shattering are determined by a specific set of thermodynamic conditions (such as a warm cloud base and modest updraft) rather than a minimum INP concentration.Too strong updrafts are not conducive to the occurrence of rime splintering and freezing drop shattering because the residence time in the favorable temperature zone is not enough.Freezing drop shattering also needs a warmer cloud base than rime splintering.The results also suggested that a certain concentration of INPs is required to initiate ice-ice collisional breakup.
The effects of SIPs on surface precipitation are studied in some model simulations.Earlier studies focused on the effect of only rime splintering.The model results by Connolly et al. [49] showed that the mean accumulated surface precipitation has no significant change when the rate of rime splintering becomes higher in a deep convective cloud case because higher rime splintering is offset by slower accretion and raindrop freezing rates, and therefore, there is no change in precipitation.When rime splintering occurs in a shallow convective cloud with cloud top temperatures higher than −8 • C, the maximum precipitation intensity is not significantly impacted while the spatial extent of the precipitation increases [50].Recently, studies have considered freezing drop shattering and ice-ice collisional breakup, as well as rime splintering in the models (e.g., [25,42]).For a summertime cold frontal cloud system with a cloud top around 7 km, where rime splintering dominates ice formation, it is found that the precipitation rate is increased in the convective towers and decreased in the non-convective gap regions by SIPs [25].Precipitation rates in summertime and wintertime mixed-phase orographic clouds are impacted when ice-ice collisional breakup dominates ice formation.Summertime accumulated surface precipitation was reduced by 20-40% [42], while wintertime localized regions of high precipitation rates were reduced [43].A sensitivity study showed that a higher rate of ice-ice collisional breakup can decrease the maximum accumulated precipitation in a deep convective event [23].
Freezing drop shattering has a relatively low rate and thus has little effect on precipitation in the cloud types in these modeling studies.For the clouds where freezing drop shattering dominates, the effect of freezing drop shattering on precipitation remains to be studied.
Previous studies mainly focused on the evidence of SIPs and the effects of SIPs on the production of ice crystals, cloud glaciation, and surface precipitation for various cloud types.Notably, there is a scarcity of studies investigating SIPs in midlatitude continental mesoscale convective systems (MCSs), like squall lines, characterized by coexisting convective and stratiform regions.The mechanisms and microphysical feedbacks involving SIPs in these MCSs also need to be studied, especially for the convective region of the MCSs, where both liquid and ice phases have high mass concentrations.The main purpose of this study is to investigate the role of SIPs, including rime splintering, freezing drop shattering, and ice-ice collisional breakup, in a squall line case in North China.We investigate the dominant SIP and the effects of SIPs on the microphysical properties and precipitation of the squall line.Mechanisms or feedbacks are also elucidated.In Section 2, we introduce the case overview and the model setup.Section 3 shows the comparisons between observations and simulations of the squall line.Section 4 examines how the different SIPs affect the squall line system and the microphysical feedback.Sections 5 and 6 discuss the unresolved issues and summarize the main conclusions, respectively.

Observation Data
The observed precipitation rate data are from GPM (Global Precipitation Measurement) IMERG Final Precipitation L3 Half Hourly 0.1 degree × 0.1 degree V06 data (hereafter GPM-Final), a Level 3 NASA product retrieved from Multi-Satellites [51].They are available from https://gpm.nasa.gov/data(accessed on 1 January 2020).These GPM-Final data are climatologically adjusted using the monthly GPCC (Global Precipitation Climatology Centre) gauge analysis product, which is derived from measurements at approximately 6700 gauge stations around the globe [51,52].In comparison to other precipitation datasets such as CMORPH [53], TMPA [51], and others, GPM-Final offers precipitation measurements at a higher spatial (0.1 • ) and temporal (half-hourly) resolution.Furthermore, Jiang et al. [54] showed that precipitation based on GPM-Final data and gauge-based precipitation across Mainland China agree quite well.For the analysis in this study, only precipitation rates higher than 0.01 mm h −1 are analyzed.
The observed radar reflectivity can be accessed from the website of the National Meteorological Center of China (http://www.nmc.cn/publish/radar/huabei.html, accessed on 1 January 2020).We only have access to the composite radar reflectivity images at a time resolution of 6 min.The radar data and the precipitation data for a squall line case in North China are compared to the model simulations with or without SIPs in this study.

Overview of the Case
A squall line case occurred in the North China area on 18 August 2020 and produced heavy precipitation.This case is associated with a cold frontal system.The prefrontal sounding at 14:00 BJT (06:00 UTC) averaged for the region of 39-42 • N and 117-121 • E is shown in Figure 1, based on Final (FNL) reanalysis data (https://rda.ucar.edu/datasets/ds083.2/,accessed on 1 January 2020).The lower troposphere of this region has sufficient water vapor and a modest Convective Available Potential Energy (CAPE) of 850 J kg −1 (parcels from the surface).The lifting condensation level is about 0.5 km at 21 • C.These conditions are favorable for convective triggering.The squall line was then triggered on the leading edge of the cold frontal system and developed into a mature stage in the period around 15:00-17:00 BJT.The squall line then moved in the northeast direction and started to decline.The early stage of this MCS is embedded in the cold frontal cloud.In the mature stage, the MCS moves faster than the cold frontal cloud and then separates.This type of MCS, triggered by cold frontal systems, occurs frequently in North China.

Model Setup
This squall line case is simulated with the Weather Research and Forecast (WRF) model, version 4.0 [55].FNL reanalysis data with a horizontal resolution of 1° × 1° every 6 h are used as the initial field and lateral boundary conditions.The case is simulated from

Model Setup
This squall line case is simulated with the Weather Research and Forecast (WRF) model, version 4.0 [55].FNL reanalysis data with a horizontal resolution of 1 • × 1 • every 6 h are used as the initial field and lateral boundary conditions.The case is simulated from 02:00 to 20:00 BJT on 18 August 2020.Figure 2 shows the satellite cloud image at 15:00 BJT and the two nested domains designed for this study.The horizontal grid spacings of the domains are 9 and 3 km, respectively.The domains have 50 layers in the vertical direction, with resolutions of about 50~300 m in the boundary layer and 400~600 m in the free troposphere.Simulation results in domain 2 output every 10 min for the analysis of the structure and processes.All the simulations use the four-order Runge-Kutta algorithm on the time integration scheme, with a time resolution of 18 s.Morrison double-moment bulk microphysical scheme [56], Yonsei University boundary layer scheme [57], RRTM longwave radiation scheme [58], and Dudhia shortwave radiation scheme [59] are used in three domains.The Kain-Fritsch cumulus parameterization [60] is only used in domain 1.Details of the simulation setup can be seen in Table 1.The Morrison double-moment bulk microphysical scheme includes five hydrometeor species: cloud droplet, rainwater, ice crystal, snow, and graupel [56].The cloud droplet concentration is set to be constant at 250 cm −3 , which is a default value in the Morrison scheme.The number concentrations of the other hydrometeors are predicted.Mass concentrations are all predicted.All hydrometeors are assumed to be spheres with Marshall-Palmer distributions when calculating the microphysical process rates.The maximum number concentration of ice crystal uses the default value of 300 L −1 .A sensitivity test with The Morrison double-moment bulk microphysical scheme includes five hydrometeor species: cloud droplet, rainwater, ice crystal, snow, and graupel [56].The cloud droplet concentration is set to be constant at 250 cm −3 , which is a default value in the Morrison scheme.The number concentrations of the other hydrometeors are predicted.Mass concentrations are all predicted.All hydrometeors are assumed to be spheres with Marshall-Palmer distributions when calculating the microphysical process rates.The maximum number concentration of ice crystal uses the default value of 300 L −1 .A sensitivity test with the maximum number concentration of ice crystal to be 1000 L −1 is also performed.The results show that this value does not affect the analysis of the effect of SIPs on concentrations of ice crystal (Figure S1) and precipitation rate (Figure S2).
In the Morrison scheme, ice crystals are produced by both primary ice processes (PIPs) and SIPs.PIPs, including homogeneous freezing of cloud droplets and heterogeneous nucleation (immersion freezing, contact freezing, and deposition and condensation freezing nucleation), are considered in the scheme.The parameterizations of PIPs are described in the Supplementary Materials (Text S1) in detail.For SIPs, only rime splintering was originally included in the scheme.Freezing drop shattering and ice-ice collisional breakup processes are added to the Morrison scheme in this study.Details of the SIP parameterizations are described in Section 2.4.
Table 2 is a summary of each case.This study mainly conducts a comparison between two simulations: one without SIPs (referred to as the NOSIP case) and another involving all three SIPs (referred to as the ALLSIP case).To dissect the individual effects and explore the relative importance of these three SIPs, three separate cases are simulated, each focusing on only one SIP.Specifically, the RS case only considers the rime splintering process, the DS case solely involves the freezing drop shattering process, and the CB case only includes the ice-ice collisional breakup process.

Parameterizations of SIPs
The rate of change of secondary ice number concentration through rime splintering is positively correlated with the rate of riming in the Morrison scheme and can be expressed as follows: where q rime is the mass concentration of rimed droplets by snow and graupel [8].The fragment number per milligram of rime, N RS , is set to 3.5 × 10 5 .The weight function of ambient temperature, w RS , has a form of piecewise linear distribution, which is set to 0 outside −3 to −8 • C and has a maximum value of 1 at −5 • C (Figure 3a).The secondary ice mass concentration is diagnosed using the number concentration, assuming all the ejected splinters have a diameter of 10 µm.This study adds freezing drop shattering to the Morrison scheme using the parameterization format in Sullivan et al. [25], which assumes all freezing drops produce the same number of ice crystals and is applicable to bulk microphysics schemes.This treatment is not dependent on drop size, as in the bin schemes [21,22,39].The rate of change of secondary ice number concentration through freezing drop shattering is as follows: where n f reeze is the number of freezing drops.In Sullivan et al. [25], only immersion freezing is considered to calculate n f reeze .Korolev et al. [19] found that the number concentration of freezing drops measured above −4 • C cannot be explained by immersion freezing alone but also by ice-rain collisions in tropical maritime MCSs.Therefore, ice-rain collision (the parameterization is described in Supplementary Materials, Text S2) is also considered to calculate n f reeze in this study, similar to Phillips et al. [22].The fragment number by freezing drop shattering, N DS , is set to 10, and the shattering probability, p DS , is set to be a normal distribution of temperature centered at −15 • C with a maximum value of 10% and a standard deviation of 5 • C (Figure 3b), as in Sullivan et al. [25].The mass concentration of secondary ice crystals through freezing drop shattering is calculated using the same method as in rime splintering in this study.Parameterization of ice-ice collisional breakup is also added to the Morrison scheme in this study.The number and mass concentrations of ice crystals produced by ice-ice collisional breakup have the same format as in Phillips et al. [24]: where N CB is the breakup factor and represents the number of ice crystals produced by one collisional breakup, n collision is the total number of collisions between ice phase particles, and q collided is the mass of collided ice phase particles.In this study, there are five ice-ice collisional types: graupel-graupel, snow-graupel, ice crystal-graupel, snow-snow, and ice crystal-snow.The calculations of N CB , n collision , and q collided are summarized in Table 3. are calculated based on the method described in Sotiropoulou et al. [45].

N CB ∂n collision ∂t ∂q collided ∂t
Graupel-Graupel Type I in Phillips et al. [24] Format in Passarelli [61], but with parameters of graupel instead of snow Format in Verlinde et al. [62], but with parameters of graupel instead of snow Snow-Graupel Type II in Phillips et al. [24] Format in Ikawa and Saito [63] Format in Ikawa and Saito [63] Ice crystal-Graupel Type II in Phillips et al. [24] Format in Reisner [64], but with parameters of graupel instead of raindrop \ Snow-Snow Type III in Phillips et al. [24] Format in Passarelli [61], which is already concluded in default Morrison scheme Format in Verlinde et al. [62] Ice crystal-Snow Type III in Phillips et al. [24] Format in Reisner [64], but with parameters of snow instead of raindrop \ N CB is calculated using the method in Phillips et al. [24] and has a format as follows: where A is a measure of the number density (per unit area) of breakable vapor-grown branches or other asperities on the colliding surfaces in the region of contact, γ is a power exponent, C is the asperity-fragility coefficient (including a collection for sublimation effect), a = πD 2 is the equivalent spherical area of the collided particle, 2 is the collisional kinetic energy between two colliding particles.In Phillips et al. [24], the ice-ice collisions are categorized into three broad types, and each type has a set of parameters of A, γ, and C. Type I is the collision of graupel with graupel, where A depends on temperature, and γ and C are two constants.Type II is the collision of ice crystals or snow with graupel, and Type III is the collision of ice crystals or snow with other ice crystals or snow.In Type II and Type III, A and γ are functions depending on the rimed fraction and habits of ice crystals and snow (either dendrite or planar), and C is a constant.
When applying the method in Phillips et al. [24] to calculate the N CB in this study, graupel-graupel collision uses the type I parameters in Phillips et al. [24], snow-graupel and ice crystal-graupel collisions use the type II parameters, while snow-snow and ice crystal-snow collisions use the type III parameters in Phillips et al. [24].Because all icephase particles are assumed to be spherical in the Morrison scheme, the calculations of A and γ for collisions involving ice crystal and snow use the functions for planar habit, which has a closer density and fall speed to the spherical shape of ice particles than the dendrites.a and K 0 are calculated using the default spherical shape, fall velocity, and effective radius calculated from the Morrison scheme.The rimed fraction is not calculated in the Morrison scheme.Sensitivity tests on the rimed fraction indicate that the rimed fraction has little effect on ice crystal number concentration (Figure S3).The rimed fraction of ice crystal and snow is set to be 0.3 in this study, representing a moderately rimed condition.are calculated using the method described in Sotiropoulou et al. [45].Here we briefly describe the methods used for each of the five collision types.∂n collision ∂t by snow-snow is calculated in the original Morrison scheme when calculating the rate of snow aggregation.The calculation is based on Passarelli [61], which is related to the size distribution and fall speed of snow particles.
by snow-snow is added to the Morrison scheme and is calculated by the analytical solution for self-collection derived by Verlinde et al. [62], depending on the size distribution and fall speed parameters of snow particles.∂n collision ∂t and ∂q collided ∂t by graupel-graupel is added to Morrison scheme using a similar method for the snow-snow collision, but with the size distribution and fall speed parameters replaced with those of graupel.∂n collision ∂t and ∂q collided ∂t by snow-graupel are added to the Morrison scheme and are calculated based on Ikawa and Saito [63], which depends on the size distribution and fall speed parameters of snow and graupel.
∂n collision ∂t by ice crystal-snow and ice crystal-graupel is added to the Morrison scheme and is calculated using the "continuous collection" approach based on Reisner [64], but with the size distribution and fall speed parameters replaced with those of snow or graupel.These two collision types do not change the mass of ice crystals, and therefore ∂q collided ∂t does not need to be calculated.

Comparison of Precipitation between Observations and Simulations
Figure 4 shows comparisons of the spatial distribution of the observed and simulated precipitation rates.10 snapshots are shown from 10:00 to 19:00 BJT.The observed precipitation rate (the left column of Figure 4) shows that the organized convective clouds formed at 10:00 BJT and gradually intensified while moving in the northeast direction.These organized convective clouds continue to grow and line up as a very typical squall line from 15:00 to 17:00 BJT.This squall line begins to decline after 17:00 BJT.The simulated precipitation rate (the middle column for the NOSIP case and the right column for the ALLSIP case) generally shows a similar squall line system forming during the time period from 10:00 to 19:00 BJT, but the location of the simulated squall line is approximately 200 km to the west of the observed squall line, and the time is about 2 h late.The simulated convective clouds line up and develop into a typical squall line during 17:00 to 19:00 BJT, which closely resembles the observed squall line during 15:00 to 17:00 BJT.

Comparison of Precipitation between Observations and Simulations
Figure 4 shows comparisons of the spatial distribution of the observed and simulated precipitation rates.10 snapshots are shown from 10:00 to 19:00 BJT.The observed precipitation rate (the left column of Figure 4) shows that the organized convective clouds formed at 10:00 BJT and gradually intensified while moving in the northeast direction.These organized convective clouds continue to grow and line up as a very typical squall line from 15:00 to 17:00 BJT.This squall line begins to decline after 17:00 BJT.The simulated precipitation rate (the middle column for the NOSIP case and the right column for the ALLSIP case) generally shows a similar squall line system forming during the time period from 10:00 to 19:00 BJT, but the location of the simulated squall line is approximately 200 km to the west of the observed squall line, and the time is about 2 h late.The simulated convective clouds line up and develop into a typical squall line during 17:00 to 19:00 BJT, which closely resembles the observed squall line during 15:00 to 17:00 BJT.In Figure 5, we focus on the precipitation area, area-averaged precipitation rate, and area-integrated precipitation rate (precipitation mass in the region per hour, a multiplication of precipitation area, area-averaged precipitation rate, and water density) during 14:00 to 19:00 BJT, which covers the typical squall line period.The precipitating area is partitioned into the convective region and stratiform region of the squall line.The observed and simulated precipitation properties are compared in each region.A precipitation rate threshold in the range of 10-25 mm h −1 has been used in the literature for partitioning convective and stratiform regions of squall lines [65].Here, in this study, we use a precipitation rate threshold of 15 mm h −1 .The convective region is defined as the strong precipitation area (precipitation rate above 15 mm h −1 ), and the stratiform region is defined as the weak precipitation area (precipitation rate below 15 mm h −1 and larger than 0.01 mm h −1 ). 10 mm h −1 and 20 mm h −1 are also tested as the thresholds for partitioning convective and stratiform regions, and the results are shown in Figures S4 and S5.Regardless of the threshold used, the convective area is one order of magnitude smaller than the stratiform area, and the ALLSIP case shows a better simulation of the precipitation in the convective region.Details will be shown below.In Figure 5, we focus on the precipitation area, area-averaged precipitation rate, and area-integrated precipitation rate (precipitation mass in the region per hour, a multiplication of precipitation area, area-averaged precipitation rate, and water density) during 14:00 to 19:00 BJT, which covers the typical squall line period.The precipitating area is partitioned into the convective region and stratiform region of the squall line.The observed and simulated precipitation properties are compared in each region.A precipitation rate threshold in the range of 10-25 mm h −1 has been used in the literature for partitioning convective and stratiform regions of squall lines [65].Here, in this study, we use a precipitation rate threshold of 15 mm h −1 .The convective region is defined as the strong precipitation area (precipitation rate above 15 mm h −1 ), and the stratiform region is defined as the weak precipitation area (precipitation rate below 15 mm h −1 and larger than 0.01 mm h −1 ). 10 mm h −1 and 20 mm h −1 are also tested as the thresholds for partitioning convective and stratiform regions, and the results are shown in Figure S4 and S5.Regardless of the threshold used, the convective area is one order of magnitude smaller than the stratiform area, and the ALLSIP case shows a better simulation of the precipitation in the convective region.Details will be shown below.For the convective region, it is seen in Figure 5a that the observed and simulated convective areas all reach their maximum during 14:00-19:00 BJT.The maximum of the observed convective areas appears at 15:00 BJT, while the maximum of the simulated convective area appears at 17:10 BJT for the NOSIP case and 17:00 BJT for the ALLSIP case, about 2 h later than the observed, as has been seen in Figure 4.The maximum of the convective area is respectively 1.02 × 10 4 , 1.27 × 10 4 , and 1.03 × 10 4 km 2 in the observation, NOSIP case, and ALLSIP case.The time-averaged convective area is shown in Figure 5g and is, respectively, 6.51 × 10 3 , 9.63 × 10 3 , and 8.03 × 10 3 km 2 in the observation, NOSIP case, and ALLSIP case.The root-mean-square error between the ALLSIP case and observation is 2.86 × 10 3 km 2 , and that between the NOSIP case and observation is 4.11 × 10 3 km 2 .The convective area in the ALLSIP case is closer to the observation in terms of the maximum and time-averaged convective area, as well as the root-mean-square error.
The observed and simulated area-averaged precipitation rates in the convective region all show typical pulsating features on a time scale of about 1 h but are out of phase (Figure 5b).The time-averaged precipitation rate during 14:00 to 19:00 BJT is, respectively, 23.01, 25.47, and 24.19 mm h −1 in the observation, NOSIP, and ALLSIP cases (Figure 5h).The precipitation rate simulated by the ALLSIP case is also closer to the observation, and the root-mean-square error between the ALLSIP case and observation is 3.27 mm h −1 , which is also smaller than that between the NOSIP case and observation (3.4 mm h −1 ).For the convective region, it is seen in Figure 5a that the observed and simulated convective areas all reach their maximum during 14:00-19:00 BJT.The maximum of the observed convective areas appears at 15:00 BJT, while the maximum of the simulated convective area appears at 17:10 BJT for the NOSIP case and 17:00 BJT for the ALLSIP case, about 2 h later than the observed, as has been seen in Figure 4.The maximum of the convective area is respectively 1.02 × 10 4 , 1.27 × 10 4 , and 1.03 × 10 4 km 2 in the observation, NOSIP case, and ALLSIP case.The time-averaged convective area is shown in Figure 5g and is, respectively, 6.51 × 10 3 , 9.63 × 10 3 , and 8.03 × 10 3 km 2 in the observation, NOSIP case, and ALLSIP case.The root-mean-square error between the ALLSIP case and observation is 2.86 × 10 3 km 2 , and that between the NOSIP case and observation is 4.11 × 10 3 km 2 .The convective area in the ALLSIP case is closer to the observation in terms of the maximum and time-averaged convective area, as well as the root-mean-square error.
The observed and simulated area-averaged precipitation rates in the convective region all show typical pulsating features on a time scale of about 1 h but are out of phase (Figure 5b).The time-averaged precipitation rate during 14:00 to 19:00 BJT is, respectively, 23.01, 25.47, and 24.19 mm h −1 in the observation, NOSIP, and ALLSIP cases (Figure 5h).The precipitation rate simulated by the ALLSIP case is also closer to the observation, and the root-mean-square error between the ALLSIP case and observation is 3.27 mm h −1 , which is also smaller than that between the NOSIP case and observation (3.4 mm h −1 ).
The area-integrated precipitation rate in the convective region, which is governed by both the convective area and the precipitation rate, is shown in Figure 5c.The area-integrated precipitation rate in the ALLSIP case is closer to the observed, as expected.In general, including SIPs improves the simulated convective area and precipitation rate in the convective region.
For the stratiform region, Figure 5d-f shows the observed and simulated area, areaaveraged precipitation rate, and area-integrated precipitation rate.Both the NOSIP case and the ALLSIP case underestimate the stratiform area (Figure 5g) and overestimate the areaaveraged precipitation rate (Figure 5h) in the stratiform region, which can also be seen in Figure 4.The area-integrated precipitation rate in the stratiform region simulated by NOSIP and ALLSIP cases is both less than the observed rate (Figure 5i).The underestimation of the stratiform area has been a common issue in the model simulation of squall lines [56,[66][67][68][69].
Here, in this study, including SIPs does not improve the simulation of the stratiform area or the precipitation rate in the stratiform region.
Figure 6 shows the histogram of the precipitation in six bins.We have shown in Figure 5 that including SIPs can improve the precipitation averaged in the convective region but does not improve the precipitation averaged in the stratiform region.Here we classify the precipitation rate into six bins and compare the simulations and observations in each bin between 14:00 and 19:00 BJT. Figure 6a shows the precipitation area, Figure 6b shows the difference in the area-averaged precipitation rate between the simulated case and the GPM observation, and Figure 6c shows the area-integrated precipitation rate.For each bin in the convective region (precipitation rate in the range of 15-20 and >20 mm h −1 ), the ALLSIP case improves the simulation of the convective area, area-averaged precipitation rate, and area-integrated precipitation rate.For each bin in the stratiform region (precipitation rate in the range of 0.01-1, 4, 5-10, and 10-15 mm h −1 ), the ALLSIP case nearly does not improve the simulation of stratiform area, area-averaged precipitation rate, and the area-integrated precipitation rate.For the bin with a precipitation rate of 0.01-1 mm h −1 , the precipitating area contributes most to the area of the stratiform region.The simulated area in this bin in the NOSIP case and the ALLSIP case is significantly smaller than that observed, which may explain the large difference in stratiform area between simulations and observations in Figure 5d.In general, considering SIPs in the model indeed improves the simulation of the convective area and precipitation rate in each bin of the convective region and therefore shows better model performance in simulating the convective area.
Three cases, including single SIP (RS, DS, and CB), also improve the simulations of convective area, area-average precipitation rate, and area-integrated precipitation rate in the convective region, but do not improve the simulations in the stratiform region (Figure S6).The CB case has a better simulation of the convective region compared to the RS case and the DS case.But it can be seen that none of these cases can simulate the convective region as well as the ALLSIP case.

Comparison of Radar Reflectivity between Observations and Simulations
Any overestimation of strong precipitation by the model might not necessarily be incorrect, as it could be attributed to observational uncertainty, as highlighted by Shahi [70].
To address this, we conducted a comparison between simulated and observed radar reflectivity.Figure 7 presents a spatial distribution comparison of observed and simulated composite radar reflectivity, featuring 10 snapshots from 10:00 to 19:00 BJT-covering the same time period as the precipitation comparison in Figure 4.The composite radar reflectivity, calculated following Stoelinga [71], reveals that the location and timing of the simulated squall line are posterior to the observed and positioned westward, aligning with the spatial distribution of precipitation rate (Figure 4).A noteworthy observation is the distinct feature where the area with composite radar reflectivity exceeding 50 dBZ (indicating strong convections and strong precipitation) in the ALLSIP case closely aligns with the observation.In contrast, the NOSIP case exhibits a significantly larger area with reflectivity exceeding the observation and the ALLSIP case across all 10 snapshots.ulated squall line are posterior to the observed and positioned westward, aligning with the spatial distribution of precipitation rate (Figure 4).A noteworthy observation is the distinct feature where the area with composite radar reflectivity exceeding 50 dBZ (indicating strong convections and strong precipitation) in the ALLSIP case closely aligns with the observation.In contrast, the NOSIP case exhibits a significantly larger area with reflectivity exceeding the observation and the ALLSIP case across all 10 snapshots.The composite radar reflectivity in observations and simulations between 14:00 and 19:00 BJT is classified into six bins.The ratio of the area of each bin to the total area of the studied region is compared between observations and simulations.For the observation, the area for each bin and the total area of the studied region are obtained using the image recognition technique because we only have images of radar reflectivity.We manually select the area in the image to replicate the simulated area (37°N-43°N and 115°E-123°E) as close as possible, as can be seen in Figure 7.It is found that the fraction of area for each The composite radar reflectivity in observations and simulations between 14:00 and 19:00 BJT is classified into six bins.The ratio of the area of each bin to the total area of the studied region is compared between observations and simulations.For the observation, the area for each bin and the total area of the studied region are obtained using the image recognition technique because we only have images of radar reflectivity.We manually select the area in the image to replicate the simulated area (37 • N-43 • N and 115 • E-123 • E) as close as possible, as can be seen in Figure 7.It is found that the fraction of area for each bin does not significantly alter when shifting the studied region in the simulated radar reflectivity by 0.5 • to the east, west, south, and north of 37 • N-43 • N and 115 • E-123 • E (Figure S7).This suggests that the slight difference between the studied area in the simulations and observation has little impact on the results of the above analysis.
Figure 8 shows the histogram of the composite radar reflectivity in 6 bins from 14:00 to 19:00 BJT.The fraction of area for each bin simulated by the ALLSIP case is closer to the observation, especially in the bins of 50-60 dBZ, as already seen in Figure 7.For the bins of 40-50 dBZ and 50-60 dBZ, the NOSIP case overestimates the fraction of area by a significant amount.For the bins of 10-20 dBZ, 20-30 dBZ, and 30-40 dBZ, the NOSIP case underestimates the fraction of the area.Note that large composite radar reflectivity usually exists in the convective region with relatively strong precipitation.Results shown in Figure 8 are consistent with the results shown in Figure 8 that the ALLSIP case has a better simulation of the convective area (defined as precipitation rate larger than 15 mm h −1 ) and that the NOSIP case overestimates the convective area by a significant amount (Figure 8a).Small composite radar reflectivity usually belongs to the stratiform region.It is seen in Figure 8 that the fraction of area for smaller radar reflectivity is still underestimated even after adding SIPs.The underestimation of the stratiform area has been a common issue in the model simulation of squall lines, which has also been discussed in Figure 5. Processes other than SIPs should be considered to resolve this issue.
Atmosphere 2023, 14, x FOR PEER REVIEW 17 of 36 bin does not significantly alter when shifting the studied region in the simulated radar reflectivity by 0.5° to the east, west, south, and north of 37°N-43°N and 115°E-123°E (Figure S7).This suggests that the slight difference between the studied area in the simulations and observation has little impact on the results of the above analysis.
Figure 8 shows the histogram of the composite radar reflectivity in 6 bins from 14:00 to 19:00 BJT.The fraction of area for each bin simulated by the ALLSIP case is closer to the observation, especially in the bins of 50-60 dBZ, as already seen in Figure 7.For the bins of 40-50 dBZ and 50-60 dBZ, the NOSIP case overestimates the fraction of area by a significant amount.For the bins of 10-20 dBZ, 20-30 dBZ, and 30-40 dBZ, the NOSIP case underestimates the fraction of the area.Note that large composite radar reflectivity usually exists in the convective region with relatively strong precipitation.Results shown in Figure 8 are consistent with the results shown in Figure 8 that the ALLSIP case has a better simulation of the convective area (defined as precipitation rate larger than 15 mm h −1 ) and that the NOSIP case overestimates the convective area by a significant amount (Figure 8a).Small composite radar reflectivity usually belongs to the stratiform region.It is seen in Figure 8 that the fraction of area for smaller radar reflectivity is still underestimated even after adding SIPs.The underestimation of the stratiform area has been a common issue in the model simulation of squall lines, which has also been discussed in Figure 5. Processes other than SIPs should be considered to resolve this issue.We now investigate if the ALLSIP case has a better simulation at all snapshots in the convective region.Figure 9 shows the time series of the fraction of area with radar reflectivity larger than 30 dBZ, 40 dBZ, and 50 dBZ.These dBZ thresholds tend to represent convective regions of the squall line.At each snapshot between 14:00 and 19:00 BJT, the fraction of the convective area simulated by the ALLSIP case is closer to the observation than the NOSIP case.This time-resolving result in Figure 9 is in agreement with the histogram analysis in Figure 8.We now investigate if the ALLSIP case has a better simulation at all snapshots in the convective region.Figure 9 shows the time series of the fraction of area with radar reflectivity larger than 30 dBZ, 40 dBZ, and 50 dBZ.These dBZ thresholds tend to represent convective regions of the squall line.At each snapshot between 14:00 and 19:00 BJT, the fraction of the convective area simulated by the ALLSIP case is closer to the observation than the NOSIP case.This time-resolving result in Figure 9 is in agreement with the histogram analysis in Figure 8. Considering the contrast between precipitation and radar reflectivity, including SIPs indeed can improve the simulation of the convective area and precipitation rate in the convective region.This indicates that future simulations of squall lines should consider SIPs in the models.However, it is noteworthy that the simulation tends to underestimate the stratiform area, and despite including SIPs, there is no observable improvement in simulating the stratiform region.

Model Uncertainty
In order to prove that the better simulation of the convective area and precipitation rate in the convective region in the ALLSIP case is caused by including SIPs in the model rather than the model uncertainty, five additional ensemble simulations are performed, respectively, in the NOSIP case and the ALLSIP case.The initial temperature field at each grid point is perturbed randomly, ranging from −1 K to +1 K, as the model setting in Qu et al. [33].Figure 10 shows the time average of precipitation area, time-and area-averaged precipitation rate, and the time-average of area-integrated precipitation rate of the convective and stratiform regions for the GPM observation, 6 ensembles of the NOSIP case, and 6 ensembles of the ALLSIP case.The 6 ensemble results are very consistent with each other.It can be seen in Figure 10 that the error bars of the 6 ensembles are very small.The difference between ensemble simulations is smaller than the difference between the NOSIP case and the ALLSIP case in the convective region.This indicates that including Considering the contrast between precipitation and radar reflectivity, including SIPs indeed can improve the simulation of the convective area and precipitation rate in the convective region.This indicates that future simulations of squall lines should consider SIPs in the models.However, it is noteworthy that the simulation tends to underestimate the stratiform area, and despite including SIPs, there is no observable improvement in simulating the stratiform region.

Model Uncertainty
In order to prove that the better simulation of the convective area and precipitation rate in the convective region in the ALLSIP case is caused by including SIPs in the model rather than the model uncertainty, five additional ensemble simulations are performed, respectively, in the NOSIP case and the ALLSIP case.The initial temperature field at each grid point is perturbed randomly, ranging from −1 K to +1 K, as the model setting in Qu et al. [33].Figure 10 shows the time average of precipitation area, time-and areaaveraged precipitation rate, and the time-average of area-integrated precipitation rate of the convective and stratiform regions for the GPM observation, 6 ensembles of the NOSIP case, and 6 ensembles of the ALLSIP case.The 6 ensemble results are very consistent with each other.It can be seen in Figure 10 that the error bars of the 6 ensembles are very small.The difference between ensemble simulations is smaller than the difference between the NOSIP case and the ALLSIP case in the convective region.This indicates that including SIPs indeed improves the simulations of convective area and precipitation rate in the convective region in this case.The microphysical feedbacks that lead to improved precipitation and radar reflectivity in the ALLSIP case are analyzed below.
Atmosphere 2023, 14, x FOR PEER REVIEW 19 of 36 SIPs indeed improves the simulations of convective area and precipitation rate in the convective region in this case.The microphysical feedbacks that lead to improved precipitation and radar reflectivity in the ALLSIP case are analyzed below.

General Features of the Simulated Cases
Figure 11 shows the time-averaged cross-sections for the NOSIP and ALLSIP cases.To get a time-averaged cross-section, we first define the leading edge of the storm at each snapshot as the leading line that is oriented in the northeast-southwest direction and also passes a maximum radar reflectivity of 35 dBZ.We then get a time-averaged cross-section based on the relative distance to the leading edge of the storm for all 13 snapshots between 17:00 and 19:00 BJT, as can be seen in Figure 11.Cross sections for single snapshots are shown in Figures S8-S10 for reference.

General Features of the Simulated Cases
Figure 11 shows the time-averaged cross-sections for the NOSIP and ALLSIP cases.To get a time-averaged cross-section, we first define the leading edge of the storm at each snapshot as the leading line that is oriented in the northeast-southwest direction and also passes a maximum radar reflectivity of 35 dBZ.We then get a time-averaged cross-section based on the relative distance to the leading edge of the storm for all 13 snapshots between 17:00 and 19:00 BJT, as can be seen in Figure 11 ), graupel mass (black line, unit: g kg −1 ), and ice crystal mass (green line, unit: g kg −1 ) in the NOSIP case.(c) The time-averaged vertical cross-section of rain water mass (shading, unit: g kg −1 ) and cloud water mass (black line, unit: g kg −1 ) in the NOSIP case.Panels (d-f) are the same as in panels (a-c), but for the ALLSIP case.Note that the time-averaged vertical cross-section is calculated based on the relative distance to the leading edge of the storm for all 13 snapshots between 17:00 and 19:00 BJT.
For both the NOSIP case and the ALLSIP case, the simulated squall line has a cloud top height of about 15 km, the −38 • C layer at about 11 km, and the 0 • C layer at about 5 km (Figure 11a,d).Here, we can define the layer above 11 km as the anvil layer.The layer at 5-11 km can be defined as the mixed-phase layer.The layer below 5 km can be defined as the warm layer.Note that the definition of the mixed-phase layer here does not imply that every grid in this layer is mixed-phase.The definition here mainly considers the height range and temperature range that are possibly mixed-phase for the convenience of analysis.As can be seen in Figure 12a, the fraction of the grids that are mixed-phase (with both liquid and ice phase particles) can be as high as 89% at about 6 km.The WBF process typically only occurs in regions where the air is supersaturated over ice but subsaturated with water in the mixed-phase layer [72].Figure 12b shows the fraction of the grids satisfying both the conditions of subsaturation for liquid water and supersaturation for ice, and Figure 12c shows the fraction of the grids where the WBF process could occur.The WBF process only occurs in less than 50% of the grids in the mixed-phase layer.Figure 11a,d shows that the maximum radar reflectivity in each column is distributed mainly in the warm layer (below 5 km).The anvil layer (above 11 km) mainly contains ice crystals and a little amount of snow and graupel (Figure 11b,e).The mixed-phase layer (5-11 km) mainly contains snow and graupel and also contains some supercooled water.The warm layer mainly contains rainwater and cloud water (Figure 11c,f), as well as melting snow and graupel.
The time-averaged cross sections in the ALLSIP case (Figure 11d-f) have some differences in the mixed-phase layer and warm layer compared with the NOSIP case (Figure 11a-c).In the mixed-phase layer (5-11 km), a very interesting feature is that the ALLSIP case has significantly higher radar reflectivity than the NOSIP case (Figure 11d), due to the significantly increased snow (Figure 11e).The mass of ice crystals is also increased in the mixed-phase layer in the ALLSIP case compared with the NOSIP case.In the warm layer, the radar reflectivity and rainwater mass are decreased in the ALLSIP case (Figure 11d,f).The smaller radar reflectivity in the ALLSIP case in the warm layer is consistent with the smaller composite radar reflectivity in the ALLSIP case as seen in Figure 7.The decreased rainwater and graupel mass account in the ALLSIP case account for the smaller radar reflectivity, compared with the NOSIP case.The difference in the hydrometeor con- Figure 11a,d shows that the maximum radar reflectivity in each column is distributed mainly in the warm layer (below 5 km).The anvil layer (above 11 km) mainly contains ice crystals and a little amount of snow and graupel (Figure 11b,e).The mixed-phase layer (5-11 km) mainly contains snow and graupel and also contains some supercooled water.The warm layer mainly contains rainwater and cloud water (Figure 11c,f), as well as melting snow and graupel.
The time-averaged cross sections in the ALLSIP case (Figure 11d-f) have some differences in the mixed-phase layer and warm layer compared with the NOSIP case (Figure 11a-c).In the mixed-phase layer (5-11 km), a very interesting feature is that the ALLSIP case has significantly higher radar reflectivity than the NOSIP case (Figure 11d), due to the significantly increased snow (Figure 11e).The mass of ice crystals is also increased in the mixed-phase layer in the ALLSIP case compared with the NOSIP case.In the warm layer, the radar reflectivity and rainwater mass are decreased in the ALLSIP case (Figure 11d,f).The smaller radar reflectivity in the ALLSIP case in the warm layer is consistent with the smaller composite radar reflectivity in the ALLSIP case as seen in Figure 7.The decreased rainwater and graupel mass account in the ALLSIP case account for the smaller radar reflectivity, compared with the NOSIP case.The difference in the hydrometeor concentrations between the two simulations and the related microphysical mechanisms will be discussed in detail below.

Ice Crystal Properties
Figure 13 shows the ice crystal number concentration (a,d), mass concentration (b,e), and effective radius (c,f) averaged over the convective region and stratiform region for the five cases.The number and mass concentrations of ice crystals progressively increase in the DS, RS, and CB cases in the mixed-phase layer, for both the convective and stratiform regions.This indicates that ice-ice collisional breakup plays the most significant role in ice crystal production among the three SIPs.The DS case can produce ice crystals with a similar number, mass concentration, and effective radius as the NOSIP case, while the CB case is similar to the ALLSIP case.The RS case produces an ice crystal number and mass concentration that is one order of magnitude lower than that in the ALLSIP case.This also suggests that freezing drop shattering played little role in this squall line case, while ice-ice collisional breakup played the dominant role, almost covering the role of rime splintering.Below, only the microphysical properties in the NOSIP and ALLSIP cases are analyzed in detail.In the convective region, as can be seen from Figure 13a,b, SIPs in the ALLSIP case lead to the most pronounced increase of ice number and mass concentrations by about 5 orders of magnitude at the −8 °C (~6 km) level, compared to the NOSIP case.The number and mass concentration of ice crystals increase by 27% (inset of Figure 7a) and 7% (inset of Figure 13b) in the anvil layer, with the effective radius decreasing by 19% (Figure 7c), which means that there are more but smaller ice crystals in the anvil layer.In the stratiform In the convective region, as can be seen from Figure 13a,b, SIPs in the ALLSIP case lead to the most pronounced increase of ice number and mass concentrations by about 5 orders of magnitude at the −8 • C (~6 km) level, compared to the NOSIP case.The number and mass concentration of ice crystals increase by 27% (inset of Figure 7a) and 7% (inset of Figure 13b) in the anvil layer, with the effective radius decreasing by 19% (Figure 7c), which means that there are more but smaller ice crystals in the anvil layer.In the stratiform region, SIPs also increase the number and mass concentrations of ice crystals most significantly at the −8 • C (~6 km) level, by 5 and 4 orders of magnitude, respectively (Figure 13d,e).The number and mass concentrations in the anvil layer increase by 47% and 33%, respectively, and hence the effective radius of ice crystals decreases by 18% (Figure 7f).

Secondary and Primary Ice Production Rates
Now we focus on the ALLSIP case to study the SIPs and their contributions to ice production.Figure 14 shows the spatial distribution of the vertically averaged rates of three SIPs in the ALLSIP case.Rime splintering and freezing drop shattering occur mainly in the convective region and only in limited areas of the stratiform region (Figure 14a,b), and also have rates at least one order of magnitude lower than ice-ice collisional breakup.This is because supercooled water is required in the rime splintering and freezing drop shattering.In contrast, the ice-phase particles are spread throughout the squall line (Figure 11e) therefore, ice-ice collisional breakup extensively occurs in both the convective and stratiform regions (Figure 14c).However, the rate of ice-ice collisional breakup in the convective region is about an order of magnitude larger than that in the stratiform region.
Atmosphere 2023, 14, x FOR PEER REVIEW 23 of 36 and also have rates at least one order of magnitude lower than ice-ice collisional breakup.This is because supercooled water is required in the rime splintering and freezing drop shattering.In contrast, the ice-phase particles are spread throughout the squall line (Figure 11e) therefore, ice-ice collisional breakup extensively occurs in both the convective and stratiform regions (Figure 14c).However, the rate of ice-ice collisional breakup in the convective region is about an order of magnitude larger than that in the stratiform region.Figure 15 shows the rates of SIPs and PIPs averaged over the convective region and the stratiform region.In the mixed-phase layer of the convective region, rime splintering produces only 3% of ice crystals, freezing drop shattering contributes less than 1%, while ice-ice collisional breakup of snow-graupel produces 90% of ice crystals, which is at least one order of magnitude higher than the other types of collisional breakup.SIPs produce ice crystals much more than PIPs from heterogeneous nucleation (5%) (Figure 15a).In the anvil layer of the convective region, rime splintering and freezing drop shattering do not occur, and collisional breakup contributes less than 1% to ice production.PIPs dominate over SIPs in ice production (Figure 15b).The process that produces the most ice crystals is homogeneous nucleation of cloud droplets (98%), rather than heterogeneous nucleation (1%).
In the mixed-phase layer of the stratiform region, supercooled cloud water and graupel are less than that in the convective region (see Figure 11e, f), so rime splintering produces only 1% of ice crystals, and freezing drop shattering contributes less than 1% (Figure 15c), which is consistent with Figure 14.Ice crystal production is also dominated by iceice collisional breakup of snow-graupel (90%).In the anvil layer of the stratiform region (Figure 15d), the ice crystal production rate is about an order of magnitude smaller than in the convective region.Homogeneous nucleation of cloud droplets also contributes most to the ice crystal number concentration (97%), which is consistent with the convective re- Figure 15 shows the rates of SIPs and PIPs averaged over the convective region and the stratiform region.In the mixed-phase layer of the convective region, rime splintering produces only 3% of ice crystals, freezing drop shattering contributes less than 1%, while ice-ice collisional breakup of snow-graupel produces 90% of ice crystals, which is at least one order of magnitude higher than the other types of collisional breakup.SIPs produce ice crystals much more than PIPs from heterogeneous nucleation (5%) (Figure 15a).In the anvil layer of the convective region, rime splintering and freezing drop shattering do not occur, and collisional breakup contributes less than 1% to ice production.PIPs dominate over SIPs in ice production (Figure 15b).The process that produces the most ice crystals is homogeneous nucleation of cloud droplets (98%), rather than heterogeneous nucleation (1%).In the NOSIP case, only the PIPs occur.For both the convective region and stratiform region, ice crystals in the mixed-phase layer are generated by the heterogeneous nucleation process, but ice crystals in the anvil layer are mainly generated by the homogeneous nucleation process.Compared with the NOSIP case, the PIP rates in the ALLSIP case can be reduced by SIPs in both the mixed-phase layer and the anvil layer (Figure 15).The relevant microphysical feedback will be explained in Section 4.4.

SIP-Related Microphysical Feedbacks in the Convective Region
Figure 16 shows the number concentration, mass concentration, and effective radius of snow, graupel, cloud water, and rainwater averaged in the convective region in the NOSIP and ALLSIP cases.The ALLSIP case generally has more ice-phase content and less liquid-phase content than the NOSIP case, indicating that SIPs accelerate cloud glaciation.In the mixed-phase layer (5 to 11 km), snow and graupel, respectively, have a 358% and 36% increase in the vertically averaged number concentration in the ALLSIP case, compared with the NOSIP case.Snow also has a 51% increase in mass concentration, but graupel has a 13% decrease in mass concentration.Cloud water and rainwater, respectively, decrease by 41% and 40% in mass concentration.In the anvil layer (above 11 km), the number and mass concentrations of snow and graupel change little.In the warm layer (below 5 km), the number concentration of rainwater is increased by 10%, while the number concentration of cloud water is fixed in this study and not shown in the figure.The masses of rainwater and cloud water decreased by 5% and 16%, respectively.In addition, the effective radius of all the particle types has decreased, indicating the important role of SIPs in affecting particle sizes.SIPs mainly occur and have a strong influence on the microphysical properties in the mixed-phase layer.In the following, we mainly focus on the SIP-related microphysical feedbacks in the mixed-phase layer.The subsequent influences on the anvil layer and the warm layer in the convective region will also be discussed.In the mixed-phase layer of the stratiform region, supercooled cloud water and graupel are less than that in the convective region (see Figure 11e, f), so rime splintering produces only 1% of ice crystals, and freezing drop shattering contributes less than 1% (Figure 15c), which is consistent with Figure 14.Ice crystal production is also dominated by ice-ice collisional breakup of snow-graupel (90%).In the anvil layer of the stratiform region (Figure 15d), the ice crystal production rate is about an order of magnitude smaller than in the convective region.Homogeneous nucleation of cloud droplets also contributes most to the ice crystal number concentration (97%), which is consistent with the convective region.
In the NOSIP case, only the PIPs occur.For both the convective region and stratiform region, ice crystals in the mixed-phase layer are generated by the heterogeneous nucleation process, but ice crystals in the anvil layer are mainly generated by the homogeneous nucleation process.Compared with the NOSIP case, the PIP rates in the ALLSIP case can be reduced by SIPs in both the mixed-phase layer and the anvil layer (Figure 15).The relevant microphysical feedback will be explained in Section 4.4.

SIP-Related Microphysical Feedbacks in the Convective Region
Figure 16 shows the number concentration, mass concentration, and effective radius of snow, graupel, cloud water, and rainwater averaged in the convective region in the NOSIP and ALLSIP cases.The ALLSIP case generally has more ice-phase content and less liquid-phase content than the NOSIP case, indicating that SIPs accelerate cloud glaciation.In the mixed-phase layer (5 to 11 km), snow and graupel, respectively, have a 358% and 36% increase in the vertically averaged number concentration in the ALLSIP case, compared with the NOSIP case.Snow also has a 51% increase in mass concentration, but graupel has a 13% decrease in mass concentration.Cloud water and rainwater, respectively, decrease by 41% and 40% in mass concentration.In the anvil layer (above 11 km), the number and mass concentrations of snow and graupel change little.In the warm layer (below 5 km), the number concentration of rainwater is increased by 10%, while the number concentration of cloud water is fixed in this study and not shown in the figure.The masses of rainwater and cloud water decreased by 5% and 16%, respectively.In addition, the effective radius of all the particle types has decreased, indicating the important role of SIPs in affecting particle sizes.SIPs mainly occur and have a strong influence on the microphysical properties in the mixed-phase layer.In the following, we mainly focus on the SIP-related microphysical feedbacks in the mixed-phase layer.The subsequent influences on the anvil layer and the warm layer in the convective region will also be discussed.In the mixed-phase layer of the convective region, the increased ice crystal number from SIPs leads to increased snow number primarily through the enhanced autoconversion of ice crystals, as shown in Figure 17a.SIPs also lead to a significant increase in the rate of ice-rain collisions, which plays a secondary role in snow number production.In the mixed-phase layer of the convective region, the increased ice crystal number from SIPs leads to increased snow number primarily through the enhanced autoconversion of ice crystals, as shown in Figure 17a.SIPs also lead to a significant increase in the rate of ice-rain collisions, which plays a secondary role in snow number production.Meanwhile, enhanced ice-rain collision becomes the dominant process that leads to a graupel number increase in the ALLSIP case, whereas nucleation of rain and snow-rain collision are the dominant processes in the NOSIP case (Figure 17c).The increased number of ice crystals and snow leads to significantly enhanced depositional growth, which is the primary reason for the increased mass of ice crystals and snow in the ALLSIP case.The depositional growth rate of ice crystals increases from 10 −5 g m −2 s −1 in the NOSIP case to 10 −3 g m −2 s −1 in the ALLSIP case, meanwhile, the depositional growth rate of snow increases from 0.7 g m −2 s −1 in the NOSIP case to 2.1 g m −2 s −1 in the ALLSIP case (Figure 17b).The enhanced depositional growth consumes a large amount of liquid water, and thus the riming rates of snow and graupel both decrease (Figure 17b,d).Therefore, the contribution of riming to snow mass is less important, whereas the contribution of depositional growth is more important and becomes the primary reason for snow mass increase in the ALLSIP case.The weaker riming of the graupel is the dominant reason for the graupel mass decrease in the ALLSIP case.
Atmosphere 2023, 14, x FOR PEER REVIEW 26 of 36 of ice crystals and snow leads to significantly enhanced depositional growth, which is the primary reason for the increased mass of ice crystals and snow in the ALLSIP case.The depositional growth rate of ice crystals increases from 10 −5 g m −2 s −1 in the NOSIP case to 10 −3 g m −2 s −1 in the ALLSIP case, meanwhile, the depositional growth rate of snow increases from 0.7 g m −2 s −1 in the NOSIP case to 2.1 g m −2 s −1 in the ALLSIP case (Figure 17b).The enhanced depositional growth consumes a large amount of liquid water, and thus the riming rates of snow and graupel both decrease (Figure 17b,d).Therefore, the contribution of riming to snow mass is less important, whereas the contribution of depositional growth is more important and becomes the primary reason for snow mass increase in the ALLSIP case.The weaker riming of the graupel is the dominant reason for the graupel mass decrease in the ALLSIP case.In the anvil layer, SIPs do not cause direct feedback but can affect PIPs by influencing the liquid water amount transported aloft from the mixed-phase layer.Less liquid water leads to a decreased homogeneous nucleation rate in the anvil layer (Figure 15b), meanwhile, the smaller liquid droplets aloft also lead to smaller ice crystals with lower fall velocity in the anvil layer (Figure 13c).The increase in ice crystals due to lower fall velocity offsets the decrease in ice production from a slower homogeneous nucleation rate, result- In the anvil layer, SIPs do not cause direct feedback but can affect PIPs by influencing the liquid water amount transported aloft from the mixed-phase layer.Less liquid water leads to a decreased homogeneous nucleation rate in the anvil layer (Figure 15b), meanwhile, the smaller liquid droplets aloft also lead to smaller ice crystals with lower fall velocity in the anvil layer (Figure 13c).The increase in ice crystals due to lower fall velocity offsets the decrease in ice production from a slower homogeneous nucleation rate, resulting in a small increase in ice crystal number, as seen in Figure 13a.The deposition of ice crystals is weak, and hence the mass of ice crystals is increased slightly.The slightly increased number of snow and graupel is caused by the slightly enhanced autoconversion of ice crystals and ice-rain collisions, as well as by lower fall velocity due to smaller sizes (see Figure 16).More numerous but smaller snow and graupel particles make their mass change insignificant.
The microphysical properties of the warm layer in the convective region are also affected by the SIP-related microphysical feedbacks in the mixed-phase layer.In the warm layer, the rain number is increased in the ALLSIP case mainly due to the melting of more numerous snow particles (Figure 17e).But the melting of snow does not directly contribute much to the rain mass, because of the smaller snow size in the ALLSIP case.Rain mass is decreased mainly due to the weaker warm rain process (Figure 17f).The warm rain process becomes weaker because of two factors.One is smaller rain from the melting of smaller snow and graupel.The other is the decreased cloud water mass.
The microphysical feedbacks in the mixed-phase layer caused by SIPs could in turn affect the SIPs.It is found that the riming rate of snow and graupel are both decreased in the ALLSIP case.Therefore, rime splintering becomes weaker because the parameterization of rime splintering depends on the mass of rimed liquid produced by snow and graupel.Moreover, the vertically-integrated rate of rime splintering averaged in the convective region is 6.6 × 10 4 m −2 s −1 in the ALLSIP case (with three SIPs), and 8.7 × 10 4 m −2 s −1 in the RS case (only with rime splintering).This indicates that rime splintering is negative feedback.Freezing drop shattering is also a negative feedback because the rain droplets are decreased in cases with SIPs.Moreover, freezing drop shattering has a negligible effect on ice multiplication due to the limited number of large droplets.In addition, SIPs lead to more numerous but smaller snow and graupel.The increased number concentrations of snow and graupel lead to an increased collision rate, whereas smaller snow and graupel can lead to a smaller collision rate.Note that, whether it is positive or negative feedback, all SIPs lead to ice multiplication, but to different degrees.Ice-ice collisional breakup has the highest contribution to ice production in this squall line case.

SIP-Related Microphysical Feedbacks in the Stratiform Region
Figure 18 shows the number concentration, mass concentration, and effective radius of snow, graupel, cloud water, and rainwater averaged in the stratiform region in the NOSIP and ALLSIP cases.In the mixed-phase layer and anvil layer, the main hydrometeor is snow, in addition to ice crystals.Snow number and mass increased by 230% and 26%, respectively, in the ALLSIP case in the mixed-phase layer, while they increased by 26% and 26% in the anvil layer.Graupel and supercooled water both have little amount in the mixed-phase and anvil layers and change little by the SIPs.In the warm layer, rain and cloud water have similar amounts.The rainwater number increases by 44% in the ALLSIP case, while the cloud water number is fixed.Rain mass increases by 28%, but cloud water mass has no significant change.In addition, the effective radius of snow and rain decreases, while that of graupel and cloud water changes little.In the following, we analyze the SIP-related microphysical feedbacks in the stratiform region and mainly focus on the changes in snow and rain.
In the mixed-phase layer, where SIPs mainly occur, the microphysical reasons for the increased snow number and mass are the enhanced local autoconversion of ice crystal process and deposition of snow, respectively (Figure 19a,b), similar to the convective region.The processes leading to the increase in snow number and mass in the anvil layer are also the autoconversion of ice crystals and the deposition of snow, respectively.Previous studies show that the microphysical properties of the stratiform region are determined by a combination of transport from the convective region and local microphysical processes [68,73,74].As can be seen from Figure 11a,d, there are air flows from the upper levels of the convective region to the stratiform region, implying the transport of hydrometeors from the convective region to the stratiform region.Therefore, the increased number and mass concentration of snow in the stratiform region may also be caused by the increased transport from the convective region, where snow number and mass are higher.However, this study does not aim to investigate the relative importance of transport and local microphysics.Our focus is on the effect of SIPs on microphysical properties.In the mixed-phase layer, where SIPs mainly occur, the microphysical reasons for the increased snow number and mass are the enhanced local autoconversion of ice crystal process and deposition of snow, respectively (Figure 19a,b), similar to the convective region.The processes leading to the increase in snow number and mass in the anvil layer levels of the convective region to the stratiform region, implying the transport of hydrometeors from the convective region to the stratiform region.Therefore, the increased number and mass concentration of snow in the stratiform region may also be caused by the increased transport from the convective region, where snow number and mass are higher.However, this study does not aim to investigate the relative importance of transport and local microphysics.Our focus is on the effect of SIPs on microphysical properties.In the warm layer of the stratiform region, the increased rain number and mass are mainly due to the increased melting of snow directly (Figure 19c,d), which is the dominant process among rainwater production.The rain mass is increased also through the enhanced melting of graupel and warm rain processes, but these two processes are not the dominant causes.In addition, the melting of snow occurs in a wide range of stratiform regions, while the melting of graupel and warm rain processes only occur in the limited stratiform region close to the convective region.This is most likely because the formation and growth of graupel is mainly limited to the convective region, whereas the formation and growth of snow can spread over both the convective and stratiform regions.
The stratiform region has limited supercooled water in the mixed-phase layer, and therefore the feedbacks involving deposition and riming as in the convective region have In the warm layer of the stratiform region, the increased rain number and mass are mainly due to the increased melting of snow directly (Figure 19c,d), which is the dominant process among rainwater production.The rain mass is increased also through the enhanced melting of graupel and warm rain processes, but these two processes are not the dominant causes.In addition, the melting of snow occurs in a wide range of stratiform regions, while the melting of graupel and warm rain processes only occur in the limited stratiform region close to the convective region.This is most likely because the formation and growth of graupel is mainly limited to the convective region, whereas the formation and growth of snow can spread over both the convective and stratiform regions.
The stratiform region has limited supercooled water in the mixed-phase layer, and therefore the feedbacks involving deposition and riming as in the convective region have relatively small effects.In addition, the ice multiplication is weaker in the stratiform region than in the convective region (as can be seen in Figure 13).

Discussion
This study is focused on the mature stage of a midlatitude continental MCS.The SIP-related microphysical feedback is very important at this stage.However, a detailed analysis of the SIPs in the early stages is not included in this study.It is possible that different SIPs work at different stages of the same cloud.Qu et al. [32] found that, in a tropical maritime deep convective cloud, freezing drop shattering contributes most to the ice multiplication during the developing stage of cloud evolution, while ice-ice collisional breakup dominates when clouds become nearly glaciated.Rime splintering is not important in this case.Consistently, this study found that ice-ice collisional breakup dominates in the mature stage of a midlatitude continental MCS.Given that this MCS is embedded in a cold frontal cloud during its initial stage, another MCS case or an ideal case may be more suitable to explore which SIP is dominant and whether these microphysical feedbacks exist in the initial stage of a midlatitude continental MCS.
Freezing drop shattering is not important in this case because supercooled liquid water is only sufficient in limited areas of the cloud in the mature stage.This is quite different from cumulus congestus in tropical regions (e.g., [21]) and tropical MCSs [33], where supercooled liquid water is the dominant phase and freezing drop shattering is very important.In contrast to tropical maritime clouds, measurements in midlatitude continental clouds show that updraft cores do not develop large (≥70 µm) supercooled drops [30,39].Less LWC and fewer large drops in this midlatitude continental MCS may result in a small role of freezing drop shattering.Future studies could investigate if a limited number of large drops is the main reason for freezing drop shattering to be unimportant.
The SIP-related microphysical feedback described in the mature stage of this squall line case is different from the feedback in Arctic stratocumulus clouds.Here, the microphysical feedbacks associated with ice-ice collisional breakup weaken rime splintering.This study also found that rime splintering has little effect on ice-ice collisional breakup because the rates of ice-ice collisional breakup in the ALLSIP case and CB case have little difference.However, in the Arctic stratocumulus cloud, rime splintering can enhance ice-ice collisional breakup, and ice-ice collisional breakup is weak when activated alone [47].The different feedbacks may be caused by different water contents in different clouds.The ice water content of MCSs is sufficient to initiate ice-ice collisional breakup at a high rate, and there is also sufficient liquid water content to maintain a high rate of ice-ice collisional breakup by allowing secondary ice crystals to grow into snow and graupel through deposition.However, Arctic clouds have very few ice-phase particles and liquid-phase particles, so rime splintering is needed to generate more ice-phase particles to maintain the high collisional breakup rate.
This study found that including SIPs (mainly ice-ice collisional breakup) results in an increase of about 27% in the ice number concentration and a decrease of 19% in the radius of ice particles in the anvil layer.This kind of change in cloud microphysical properties is similar to the aerosol effect on anvil clouds, which may result in a significant change in cloud radiative properties [75].The simulations without ice-ice collisional breakup process in Fan et al. [75] show that the ice crystal number concentration in the anvil layer increases by a factor of 3-5 and the size of ice particles is reduced by up to 50% when the initial aerosol concentrations increase by a factor of 6.It seems that SIPs and increased aerosol number concentrations have similar effects on the anvil layer, but the aerosol effect is stronger than the SIP effect.The relative importance and joint impact of SIPs and aerosol effects on cloud anvils can be studied in the future.
The study of Han et al. [68] mentioned that the simulated precipitation rate in the stratiform region of squall lines is often lower, while that in the convective region is often higher than the observation, but the reason is still unclear.Several mechanisms, including SIPs, have been proposed, but ice-ice collisional breakup is not considered in their study.
Our study shows that adding ice-ice collisional breakup can lead to less precipitation in the convective region and more precipitation in the stratiform region.Moreover, including three SIPs (especially ice-ice collisional breakup) improves the simulations of the convective area and precipitation rate in the convective region.Note that the parameterizations of three SIPs utilized in the study are based on theoretical results or laboratory results.The ALLSIP case that uses these parameterizations has a better simulation of the surface precipitation and radar reflectivity for this squall line.In the future, more model simulations for more cases are needed to investigate if considering SIPs can improve the simulations in other clouds.
The microphysical feedback in the convective region caused by SIPs could change the latent heats, which could have an impact on the dynamics.As can be seen from Figure 20, the maximum vertical velocity in the convective region increases in the mixed-phase layer in the ALLSIP case.The stronger updrafts may be caused by the SIP-induced phase transition from water vapor and liquid to ice phase, which promotes latent heat release and convection invigoration.However, surface precipitation still decreases as the cloud is invigorated, possibly because the microphysical feedbacks are strong and can dominate over the dynamic effects.
three SIPs (especially ice-ice collisional breakup) improves the simulations of the convective area and precipitation rate in the convective region.Note that the parameterizations of three SIPs utilized in the study are based on theoretical results or laboratory results.The ALLSIP case that uses these parameterizations has a better simulation of the surface precipitation and radar reflectivity for this squall line.In the future, more model simulations for more cases are needed to investigate if considering SIPs can improve the simulations in other clouds.
The microphysical feedback in the convective region caused by SIPs could change the latent heats, which could have an impact on the dynamics.As can be seen from Figure 20, the maximum vertical velocity in the convective region increases in the mixed-phase layer in the ALLSIP case.The stronger updrafts may be caused by the SIP-induced phase transition from water vapor and liquid to ice phase, which promotes latent heat release and convection invigoration.However, surface precipitation still decreases as the cloud is invigorated, possibly because the microphysical feedbacks are strong and can dominate over the dynamic effects.

Conclusions
Three SIPs, including rime splintering, freezing drop shattering, and ice-ice collisional breakup, were investigated in a squall line case in North China on 18 August 2020 using the WRF model with the Morrison double-moment bulk microphysical scheme, which originally includes rime splintering.Freezing drop shattering and ice-ice collisional breakup processes are added to the Morrison scheme here.This study mainly investigates the relative importance of rime splintering, freezing drop shattering, and ice-ice collisional breakup in the mature stage of the squall line case.
The simulations show that including SIPs in the model, especially ice-ice collisional breakup, improves the simulations of convective area and convective precipitation rate, but slightly improves the simulations in the stratiform region.SIPs mainly occur in the mixed-phase layer, exhibiting the most pronounced effect on ice production at about the −8 °C (~6 km) level.Here, ice-ice collisional breakup and rime splintering can respectively increase ice crystal number concentration by up to 5 and 4 orders of magnitude compared to the simulation with no SIPs.Freezing drop shattering has a negligible effect on ice production because there are few large droplets in the mature stage.Ice-ice collisional breakup affects both the convective and stratiform regions of the squall line case, while rime splintering mainly affects the convective region but occurs in a limited area in the stratiform region.

Conclusions
Three SIPs, including rime splintering, freezing drop shattering, and ice-ice collisional breakup, were investigated in a squall line case in North China on 18 August 2020 using the WRF model with the Morrison double-moment bulk microphysical scheme, which originally includes rime splintering.Freezing drop shattering and ice-ice collisional breakup processes are added to the Morrison scheme here.This study mainly investigates the relative importance of rime splintering, freezing drop shattering, and ice-ice collisional breakup in the mature stage of the squall line case.
The simulations show that including SIPs in the model, especially ice-ice collisional breakup, improves the simulations of convective area and convective precipitation rate, but slightly improves the simulations in the stratiform region.SIPs mainly occur in the mixed-phase layer, exhibiting the most pronounced effect on ice production at about the −8 • C (~6 km) level.Here, ice-ice collisional breakup and rime splintering can respectively increase ice crystal number concentration by up to 5 and 4 orders of magnitude compared to the simulation with no SIPs.Freezing drop shattering has a negligible effect on ice production because there are few large droplets in the mature stage.Ice-ice collisional breakup affects both the convective and stratiform regions of the squall line case, while rime splintering mainly affects the convective region but occurs in a limited area in the stratiform region.
A chain of microphysical changes occurs in the convective region and stratiform region.Following the ice multiplication by ice-ice collisional breakup and rime splintering, autoconversion of ice crystals and ice-rain collision are enhanced, leading to more numerous snow and graupel.The increased number of ice-phase particles enhances the deposition process, which increases the mass of ice and snow while decreasing the mass of supercooled water through WBF.Subsequently, this weakens the riming process and decreases the mass of the graupel.Notably, all the ice-phase particles have a smaller size in the ALLSIP case compared to the NOSIP case.
Moreover, this chain of microphysical changes can feedback to rime splintering and ice-ice collisional breakup in different ways.Rime splintering is limited by itself and also limited by ice-ice collisional breakup because the riming due to the two SIPs leads to a lower rime splintering rate, as set in the parameterization of rime splintering.Therefore, in the presence of ice-ice collisional breakup, rime splintering plays a secondary role in secondary ice production in the mature stage of the system.In addition, SIPs lead to more numerous but smaller snow and graupel.The increased number concentrations of snow and graupel lead to an increased collision rate, whereas smaller snow and graupel can lead to a smaller collision rate.Note that, whether it is positive or negative feedback, all SIPs lead to ice multiplication, but to different degrees.The ice-ice collisional breakup has the highest contribution to ice production in this squall line case.
The ice-ice collisional breakup has a slower rate than the primary ice production in the anvil layer.Therefore, it does not directly affect ice production in the anvil layer.However, ice-ice collisional breakup can affect PIPs by transporting fewer and smaller droplets aloft from the mixed-phase layer in the convective region.The fewer droplets and smaller droplets lead to fewer and smaller ice crystals through PIPs, which have a lower fall velocity.The lower fall velocity of smaller ice crystals leads to an increase in ice number, offsetting the decrease in ice production from PIPs, resulting in an increase of 27% and 7% in ice crystal number and mass in the convective region, respectively.The stratiform region has increased in ice number and mass, which is also due to the transport of ice from the convective to the stratiform region.
The rainwater in the warm layer and the surface precipitation rates are also affected by ice-ice collisional breakup.In the convective region, warm rain is the dominant process for rain mass.Ice-ice collisional breakup, on the one hand, can lead to more numerous and smaller snow in the mixed-phase layer, resulting in more and smaller raindrops after melting.On the other hand, ice-ice collisional breakup can lead to decreased cloud water mass through enhanced WBF.Therefore, the warm rain process becomes weaker, and the rain mass decreases.In the stratiform region, the melting of snow is the dominant process for rain mass.Ice-ice collisional breakup increases snow number and mass, and the enhanced melting of snow directly leads to increased number and mass concentrations of rain in the warm layer.As such, there is a smaller surface precipitation rate in the convective region and a larger rate in the stratiform region.Compared with observations, including SIPs (especially ice-ice collisional breakup) in the model improves the simulations of convective area and convective precipitation rate, but slightly improves the simulations in the stratiform region.
This study only simulated one midlatitude continental squall line case.In the future, it is necessary to consider whether the roles of the three SIPs are different among different MCSs.In addition, the relative importance and joint impact of SIPs and aerosol effects on cloud anvils can be studied in the future.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/atmos14121752/s1.Text S1 and Text S2 provide detailed descriptions of parameterizations of primary ice processes and the formation of freezing drops.Figure S1.The number concentration (a,d), mass concentration (b,e), and the effective radius of ice crystals (c,f) averaged in the convective region (a-c) and stratiform region (d-f) during 17:00 to 19:00 BJT for the NOSIP, ALLSIP, NOSIP-Nm1000, ALLSIP-Nm1000 cases.The small boxes in these panels are the profiles of the 11-18 km (anvil layer) in a linear scale.Figure S2: The temporal average of precipitation area (a), temporal and area-average precipitation rate (b), and temporal average of area-integrated precipitation rate (c) of the convective and stratiform region for the GPM observation, NOSIP case, ALLSIP case, NOSIP-Nm1000 case, and ALLSIP-Nm1000 case. Figure S3.The number concentration of ice crystals averaged in the convective and stratiform region during 17:00 to 19:00 BJT for the CB-rf0.2,CB-rf0.3, and CB-rf0.4cases.Figure S4.(a-c): The time series of convective area, area-average precipitation rate, and area-integrated precipitation rate for the convective region.(d-f): The same as in (a-c), but for the stratiform region.(g-i): The temporal average of area, area-average precipitation rate, and the area-integrated precipitation rate of the convective and stratiform region, respectively.The convective and stratiform region are defined by threshold of 10 mm h −1 .Figure S5.The same as Figure S4, but the convective region and stratiform region are defined by of 20 mm h −1 .Figure S6.The temporal average of precipitation area (a), temporal and area-average precipitation rate (b), and the temporal average of area-integrated precipitation rate (c) of the convective and stratiform region for the GPM observation and five simulations.The analysis time range is 14:00-19:00 BJT. ), graupel mass (black line, unit: g kg −1 ), and ice crystal mass (green line, unit: g kg −1 ) averaged along the horizontal axis in (b).(e): The vertical cross section of rain water mass (shading, unit: g kg −1 ) and cloud water mass (black line, unit: g kg −1 ) averaged along the horizontal axis in (b); Panels (f-j) are the same as in panels (a-e), but for the ALLSIP case. Figure S9.The same as Figure S8, but at 18:00 BJT. Figure S10.The same as Figure S8, but at 19:00 BJT.

Figure 2 .
Figure 2. The visible image from Himawari-8 at 15:00 BJT on 18 August 2020 and the domain design of the WRF model.Domains 1 and 2 are abbreviated to d01 and d02, respectively.

Figure 2 .
Figure 2. The visible image from Himawari-8 at 15:00 BJT on 18 August 2020 and the domain design of the WRF model.Domains 1 and 2 are abbreviated to d01 and d02, respectively.

Figure 3 .
Figure 3.The distribution of weight function w RS in the rime splintering (a) and shattering probability p DS in the freezing droplet shattering (b).

∂n collision ∂t and ∂q collided ∂t are predicted in
the bin microphysics scheme in Phillips et al. (2017a).In this study, ∂n collision ∂t and ∂q collided ∂t

Figure 4 .
Figure 4. Spatial distribution of precipitation rate for GPM (left column), NOSIP case (middle column), and ALLSIP case (right column) for the 10 snapshots from 10:00 to 19:00 BJT.Only precipitation rates higher than 0.01 mm h −1 are plotted in this figure.

Figure 4 .
Figure 4. Spatial distribution of precipitation rate for GPM (left column), NOSIP case (middle column), and ALLSIP case (right column) for the 10 snapshots from 10:00 to 19:00 BJT.Only precipitation rates higher than 0.01 mm h −1 are plotted in this figure.

Figure 5 .
Figure 5. (a-c) The time series of convective area (a), area-average precipitation rate (b), and areaintegrated precipitation rate (c) for the convective region.(d-f) The same as in (a-c), but for the stratiform region.(g-i) The temporal average of area, the area-average precipitation rate, and the area-integrated precipitation rate of the convective and stratiform regions, respectively.The convective and stratiform regions are defined by a threshold of 15 mm h −1 .

Figure 5 .
Figure 5. (a-c) The time series of convective area (a), area-average precipitation rate (b), and areaintegrated precipitation rate (c) for the convective region.(d-f) The same as in (a-c), but for the stratiform region.(g-i) The temporal average of area, the area-average precipitation rate, and the areaintegrated precipitation rate of the convective and stratiform regions, respectively.The convective and stratiform regions are defined by a threshold of 15 mm h −1 .

Figure 6 .
Figure 6.(a) Histogram of the precipitation area of each bin.(b) Histogram of the difference in area-average precipitation rate in each bin between the simulated case and the GPM observation.(c) Histogram of the area-integrated precipitation rate in each bin.The analysis time period is 14:00-19:00 BJT.

Figure 7 .
Figure 7. Spatial distribution of composite radar reflectivity (unit: dBZ) for observations (left column), NOSIP case (middle column), and ALLSIP case (right column) for the 10 snapshots from 10:00 to 19:00 BJT.The composite radar reflectivity is calculated as the maximum radar reflectivity in each column.The color bar used for composite radar reflectivity is exactly the same for observations and simulations.

Figure 7 .
Figure 7. Spatial distribution of composite radar reflectivity (unit: dBZ) for observations (left column), NOSIP case (middle column), and ALLSIP case (right column) for the 10 snapshots from 10:00 to 19:00 BJT.The composite radar reflectivity is calculated as the maximum radar reflectivity in each column.The color bar used for composite radar reflectivity is exactly the same for observations and simulations.

Figure 8 .
Figure 8.The histogram of the composite radar reflectivity.The analysis time period is 14:00-19:00 BJT for the observation, NOSIP case, and ALLSIP case.The studied area in the NOSIP and ALLSIP cases is 37° N-43° N and 115° E-123° E. The studied area in the observation is shown in the left column of Figure 7.The fraction is the ratio of the area of each bin to the total area of the studied region.

Figure 8 .
Figure 8.The histogram of the composite radar reflectivity.The analysis time period is 14:00-19:00 BJT for the observation, NOSIP case, and ALLSIP case.The studied area in the NOSIP and ALLSIP cases is 37 • N-43 • N and 115 • E-123 • E. The studied area in the observation is shown in the left column of Figure 7.The fraction is the ratio of the area of each bin to the total area of the studied region.

Figure 9 .
Figure 9.The time series of the fraction of area with the radar reflectivity larger than 30 dBZ (a), 40 dBZ (b), and 50 dBZ (c).The fraction is the ratio of the area for each bin to the total area of the studied region.

Figure 9 .
Figure 9.The time series of the fraction of area with the radar reflectivity larger than 30 dBZ (a), 40 dBZ (b), and 50 dBZ (c).The fraction is the ratio of the area for each bin to the total area of the studied region.

Figure 10 .
Figure 10.The temporal average of precipitation area (a), temporal and area-average precipitation rate (b), and temporal average of area-integrated precipitation rate (c) of the convective and stratiform region for the GPM observation, 6 ensembles of the NOSIP case, and 6 ensembles of the ALLSIP case.The analysis time range is 14:00-19:00 BJT.Color bars are the means of the ensembles, and the error bar represents the maximum and minimum value in the 6 ensemble simulations.

Figure 10 .
Figure 10.The temporal average of precipitation area (a), temporal and area-average precipitation rate (b), and temporal average of area-integrated precipitation rate (c) of the convective and stratiform region for the GPM observation, 6 ensembles of the NOSIP case, and 6 ensembles of the ALLSIP case.The analysis time range is 14:00-19:00 BJT.Color bars are the means of the ensembles, and the error bar represents the maximum and minimum value in the 6 ensemble simulations.

36 Figure 11 .Figure 11 .
Figure11shows the time-averaged cross-sections for the NOSIP and ALLSIP cases.To get a time-averaged cross-section, we first define the leading edge of the storm at each snapshot as the leading line that is oriented in the northeast-southwest direction and also passes a maximum radar reflectivity of 35 dBZ.We then get a time-averaged cross-section based on the relative distance to the leading edge of the storm for all 13 snapshots between 17:00 and 19:00 BJT, as can be seen in Figure11.Cross sections for single snapshots are shown in Figures S8-S10 for reference.Atmosphere 2023, 14, x FOR PEER REVIEW 20 of 36

36 Figure 12 .
Figure 12.(a) The fraction of the grids belonging to the mixed-phase (with both liquid and ice phase particles) to all grids in each precipitation rate bin.(b) The fraction of the grids satisfying both the conditions of subsaturation for liquid water and supersaturation for ice in all grids in each precipitation rate bin.(c) The fraction of the grids where the WBF process could occur to all grids in each precipitation rate bin.Between the two solid gray lines is the mixed-phase layer (5-11 km).The gray dashed line represents the separation line between the convective region and the stratiform region.The analysis time range is 17:00-19:00 BJT.

Figure 12 .
Figure 12.(a) The fraction of the grids belonging to the mixed-phase (with both liquid and ice phase particles) to all grids in each precipitation rate bin.(b) The fraction of the grids satisfying both the conditions of subsaturation for liquid water and supersaturation for ice in all grids in each precipitation rate bin.(c) The fraction of the grids where the WBF process could occur to all grids in each precipitation rate bin.Between the two solid gray lines is the mixed-phase layer (5-11 km).The gray dashed line represents the separation line between the convective region and the stratiform region.The analysis time range is 17:00-19:00 BJT.
Atmosphere 2023, 14, x FOR PEER REVIEW 22 of 36 ice collisional breakup played the dominant role, almost covering the role of rime splintering.Below, only the microphysical properties in the NOSIP and ALLSIP cases are analyzed in detail.

Figure 13 .
Figure 13.Vertical profile of the number concentration (a,d), mass concentration (b,e), and effective radius (c,f) of ice crystals averaged in the convective region (a-c) and stratiform region (d-f) during 17:00 BJT to 19:00 BJT for five cases.The small boxes in these panels are the profiles of the 11-18 km (anvil layer) in a linear scale for the NOSIP case and the ALLSIP case.

Figure 13 .
Figure 13.Vertical profile of the number concentration (a,d), mass concentration (b,e), and effective radius (c,f) of ice crystals averaged in the convective region (a-c) and stratiform region (d-f) during 17:00 BJT to 19:00 BJT for five cases.The small boxes in these panels are the profiles of the 11-18 km (anvil layer) in a linear scale for the NOSIP case and the ALLSIP case.

Figure 14 .
Figure 14.The distribution of vertically integrated rime splintering rate (a), freezing drop shattering rate (b), and total ice-ice collisional breakup rate (c) in the ALLSIP case at 17:00 BJT.The black line in (a-c) represents the boundary of the convective region.

Figure 14 .
Figure 14.The distribution of vertically integrated rime splintering rate (a), freezing drop shattering rate (b), and total ice-ice collisional breakup rate (c) in the ALLSIP case at 17:00 BJT.The black line in (a-c) represents the boundary of the convective region.

Figure 15 .
Figure 15.The vertically integrated production rates of ice crystals averaged of the convective region (a,b) and stratiform region (c,d) during 17:00 BJT and 19:00 BJT for the ALLSIP case.(a,c) are for the mixed-phase layer.(b,d) are for the anvil layer.Note that GGCB, SGCB, IGCB, SSCB, and ISCB represent graupel-graupel, snow-graupel, ice crystal-graupel, snow-snow, and ice crystal-snow collisional breakup.RS represents rime splintering.DS represents freezing drop shattering.HE and HO represent heterogeneous nucleation and homogeneous nucleation.

Figure 15 .
Figure 15.The vertically integrated production rates of ice crystals averaged of the convective region (a,b) and stratiform region (c,d) during 17:00 BJT and 19:00 BJT for the ALLSIP case.(a,c) are for the mixed-phase layer.(b,d) are for the anvil layer.Note that GGCB, SGCB, IGCB, SSCB, and ISCB represent graupel-graupel, snow-graupel, ice crystal-graupel, snow-snow, and ice crystal-snow collisional breakup.RS represents rime splintering.DS represents freezing drop shattering.HE and HO represent heterogeneous nucleation and homogeneous nucleation.

Atmosphere 2023 , 36 Figure 16 .
Figure 16.Vertical profile of the number concentration (a,d,g), mass concentration (b,e,h,j), and effective radius (c,f,i,k) of snow (a-c), graupel (d-f), rain water (g-i), and cloud water (j,k) averaged in the convective region during 17:00 BJT and 19:00 BJT for the NOSIP and ALLSIP cases.

Figure 16 .
Figure 16.Vertical profile of the number concentration (a,d,g), mass concentration (b,e,h,j), and effective radius (c,f,i,k) of snow (a-c), graupel (d-f), rain water (g-i), and cloud water (j,k) averaged in the convective region during 17:00 BJT and 19:00 BJT for the NOSIP and ALLSIP cases.

Figure 17 .
Figure 17.Simulated vertically integrated production rates of the number concentration (a,c,e) and mass concentration (b,d,f) of snow (a,b), graupel (c,d) and rainwater (e,f) averaged in the convective region during 17:00 BJT and 19:00 BJT.Note that the warm rain process includes the autoconversion of cloud droplets to raindrops and the collection of cloud droplets by raindrops.

Figure 17 .
Figure 17.Simulated vertically integrated production rates of the number concentration (a,c,e) and mass concentration (b,d,f) of snow (a,b), graupel (c,d) and rainwater (e,f) averaged in the convective region during 17:00 BJT and 19:00 BJT.Note that the warm rain process includes the autoconversion of cloud droplets to raindrops and the collection of cloud droplets by raindrops.

Figure 18 .
Figure 18.Vertical profile of the number concentration (a,d,g), mass concentration (b,e,h,j), and effective radius (c,f,i,k) of snow (a-c), graupel (d-f), rain water (g-i), and cloud water (j,k) averaged in the stratiform region during 16:00 BJT and 19:00 BJT for the NOSIP and ALLSIP cases.

Figure 19 .
Figure 19.Simulated vertically integrated production rates of the number concentration (a,c,e) and mass concentration (b,d,f) of snow (a,b), graupel (c,d) and rainwater (e,f) averaged in the stratiform region during 17:00 BJT and 19:00 BJT.

Figure 19 .
Figure 19.Simulated vertically integrated production rates of the number concentration (a,c,e) and mass concentration (b,d,f) of snow (a,b), graupel (c,d) and rainwater (e,f) averaged in the stratiform region during 17:00 BJT and 19:00 BJT.

Figure 20 .
Figure 20.Vertical profiles of maximum vertical velocity in the convective region during 17:00 and 19:00 BJT.

Figure 20 .
Figure 20.Vertical profiles of maximum vertical velocity in the convective region during 17:00 and 19:00 BJT.
Figure S7.The histogram of the composite radar reflectivity.But the analysis region is shifting the simulation analysis region (37 • N-43 • N and 115 • E-123 • E) by 0.5 • to the east (a), west (b), south (c), and north (d), respectively.Figure S8.(a,b): Maximum (Composite) radar reflectivity at 17:00 BJT in the NOSIP case.(b) is the part of the red box in (a); (c): The vertical cross section of radar reflectivity (shading, unit: dbz), wind (arrows), and temperature (red line, unit: • C) averaged along the horizontal axis in (b).(d): The vertical cross section of snow mass (shading, unit: g kg −1

Table 1 .
Simulation setup in the Weather Research and Forecasting model (WRF).

Table 3 .
The calculation method of three parameters in the equations of ice-ice collisional breakup rate.