Numerical Sensitivity Studies on Effects of Ice Nucleating Processes on Electrification in Thunderstorms

Aerosol particles serve as cloud condensation nuclei (CCN) and ice nuclei (IN) for the formation of both cloud droplets and ice crystals, and thus, they perhaps exert a pronounced influence on the microphysical properties of clouds. As lightning flashes are produced in thunderclouds, considerable progress has been made in understanding the interaction between aerosol and lightning (Mansell and Ziegler ^[1]; Shi et al. ^[2]; Tan et al. ^[3]; Tan et al. ^[4]). However, in contrast to the interaction between CCN and the electrification process that has been well understood, the influence of IN on electrification in thunderstorms remains highly uncertain. It is commonly admitted that ice crystals are formed via either homogeneous or heterogeneous nucleation; therefore, the effects of two nucleation processes on the electrification in thunderstorms need to be studied.

The charge centers reside in regions where ice particle production may be actively involved in electrification (Reynolds and Neill ^[5]; Saunders et al. ^[6]). The main mechanism responsible for thunderstorm electrification is the non-inductive charging separation between ice crystal and graupel/hail (Takahashi ^[7]; Jayaratne and Saunders ^[8]; Saunders et al. ^[9]; Ziegler et al. ^[10]; Saunders and Peck ^[11]). Since the charge distributions are prominently related to the feature of lightning discharges (Carey and Rutledge ^[12]; Coleman et al. ^[13]; Qie et al. ^[14]; Tan et al. ^[15]; Tan et al. ^[16]), the formation of ice crystals is vital to the electrification process and the production of lightning. The size of ice crystals plays an important role in charge transfer processes in thunderclouds (Jayaratne and Saunders ^[8]; Mansell et al. ^[17]). The transferred charge obtained by small ice crystals is very sensitive to the size of ice crystals, while the larger ice crystal stably charges when it experiences further growth. Brooks et al. ^[18] suggested the effect of the rime accretion rate on the separated charge is more appropriate to determine the charge transfer than only the effective water content. As falling velocities are primarily related to the size of ice particles, the diameter plays an important role in the electrification process. Moreover, a sensitive modeling study shows that the non-inductive charge separation is very sensitive to the concentration of ice crystals (Mansell et al. ^[17]). Several observational studies have revealed a strong correlation between lightning flash rate and ice masses (Blyth et al. ^[19]; Latham et al. ^[20]; Deierling et al. ^[21]). Therefore, the microphysical characteristics of ice crystals have a great influence on electrification.

Ice crystal formation in thunderstorms may be caused by freezing droplets at temperatures lower than –12℃(Pruppacher and Klett ^[22]). A large fraction of initial ice crystals are homogeneously produced at temperature below –37℃. Heterogeneous freezing can be available and efficient at lower altitudes and may have important effects on ice initiation and precipitation (Gilmore et al. ^[23]; Van Den Heever et al. ^[24]; Ekman et al. ^[25]; Phillips et al. ^[26]; Pruppacher and Klett ^[27]). Compared to homogeneous nucleation, which is relatively well understood, the heterogeneous nucleation process is much more complicated (Heymsfield et al. ^[28]). IN particles can be obtained via several distinct pathways (Vali et al. ^[29]), for example, deposition freezing (deposition of water vapor as ice directly on a dry aerosol particle), immersion freezing (INs are immersed in supercooled cloud droplets, and then freezing occurs), contact freezing (an aerosol particle contacts the surface of a liquid droplet), and condensation freezing (water vapor condenses on IN, which acts similarly to cloud condensation nuclei). However, despite the small number of heterogeneously formed ice crystals in the cloud (occurring at relatively low levels and at an early stage compared with homogeneous freezing), it is still suggested that heterogeneously formed ice crystals may significantly alter the homogeneous freezing process. Furthermore, there is a water vapor competition between homogeneous freezing of droplets and heterogeneous freezing of IN, leading to a strong reduction in ice crystal concentration (Demott et al. ^[30]). An analysis of winter storms in California showed that immersion freezing is the dominant nucleation process in continental storms, but contact freezing is more significant over the ocean (Lee et al.^[31]).

A modest number of studies have been conducted to investigate the relationship between IN and electrification in thunderstorms. Gonçalves et al. ^[32] developed the numerical simulation using the Brazilian Regional Atmospheric Model System to investigate the effect of bacterial ice nuclei on the frequency and intensity of lightning activity. Their results show that IN bacteria clouds directly affect the thunderstorm microphysical development and lightning formation. However, we still know little about the role of homogeneous and heterogeneous freezing processes on ice crystal production and electrification in thunderstorms.

As the electrification in thunderstorms is intrinsically linked to microphysics and hydrometeor particles, the ice nucleation processes may have an important potential impact on thunderstorm electrification. This study aims to present sensitivity studies of the effect of two different ice nucleating processes on the microphysics and electrification within thunderstorm clouds. For this purpose, a two-dimensional cumulus model with detailed cloud microphysics and an electrification scheme is used. Numerical experiments are mainly tested for the relationship between ice nucleating and electrification in thunderstorms.

The numerical model used in this study was developed by Hu and He ^[33]. It is a non-hydrostatic cumulus model (refer to the Appendix for more details). As described in Shi et al. ^{[2, 34]}, a resolution of 250 m and time steps of 2 s were used to calculate the microphysical and electrification processes in the 76 km×20 km domain. The microphysics scheme had five hydrometeor categories with a gamma function distribution. The five categories were cloud droplets, rain, ice crystal, graupel, and hail. The model predicted the mixing ratio and number concentration of each category. The main cloud physical processes were activation, condensation, evaporation, collision, auto conversion, nucleation, multiplication, melting, and freezing.

For electrification, the Gardiner/Ziegler non-inductive charging parameterization scheme involving charge separation between ice crystal and graupel was used in this simulation (Ziegler et al. ^[10]). The sign of charge acquired by the graupel depends on the cloud water content and the ambient temperature. According to Mansell et al. ^[17]:

where D_g and D_i are the diameters of the colliding particles (graupel and ice crystal). E_r is the rebound probability, and E_gi is the graupel-ice crystal collision efficiency; in this study, E_r = 0.01 and E_gi = 0.7. N is the number concentration.V_i and V_g are the mass-weighted mean terminal speeds for ice crystal and graupel. β is given similar to Mansell et al. ^[17] by:

The charge per collision δ_q is approximated as:

where D_i is the diameter for ice, δL is a parameter related to cloud water content (CWC) and it is given by as:

where τ = (−21/T_r)(T − 273.16) is the scaled temperature used by Ziegler et al. ^[10] to allow the reversal temperature T_r to be varied. The reversal temperature for CWC above 0.1 g m^–3 is set at T_r = –15℃. At temperatures below T_r, graupel (ice) charges negatively (positively) and at higher temperatures, the charging sign is reversed.

Electrification via induction in the model occurs when graupel particles collide with cloud droplets. Inductive collision charging parameterization is based on Ziegler et al. ^[10]:

where Q_eg is the individual charge from graupel, and D_c and D_g are the diameters of cloud droplets and graupel, respectively. V_g is the falling speed of graupel, and N and N_g are the cloud droplet and graupel concentrations, respectively. N_0g is the number concentration intercept for graupel. Γ(x) is the complete gamma function, and E_z is the vertical electric field. The symbols E_gc and E_r denote graupel-cloud droplet collision efficiency and rebound probability, respectively. θ is the polar collision angle. According to Mansell et al.^[17], the coefficients for inductive graupel cloud droplet charging in this study ( E_r = 0.01 and cosθ = 0.4) fall within the moderate to strong range, with E_r spanning from 0.007 to 0.015 and cosθ from 0.2 to 0.5.

It must be stated that lightning discharge processes are considered to restrain the charge accumulation in this study. Lightning discharges are parameterized based on Tan et al. ^{[15, 16]}. Lightning initiation uses the runaway electron threshold for the break-even field, and thereafter, bidirectional channels are propagated in a stochastic step-by-step fashion. The leaders of IC lightning do not reach the ground, and a height threshold (1.5 km or 6 grid points above ground) is used to define a flash to be CG lightning (including +CG and –CG).

It has been suggested that micron cloud particles can be frozen homogeneously at temperatures between 35℃ and 40℃, and a lower temperature even initiates the freezing of cloud particles with a size of approximately 1 μm (Sassen and Dodd ^[35]; Heymsfield and Sabin ^[36]; Heymsfield and Miloshevich ^[37]). The size means the diameter. Homogeneous freezing is described in the model according to the approach of Koop and Murray ^[38]. The coefficient of homogeneous nucleation rate is calculated according to:

where Q_cif is the mixing ratio of nucleated ice crystals. Q_c is the mixing ratio of supercooled liquid drops that can be frozen, including both cloud droplets and rain drops. The probability P(T, t) that drops are frozen at time t can be written as:

where T is the ambient temperature in Kelvin, V is the volume of drops, and J_v(T) is the volume-dependent homogeneous ice nucleation rate coefficient, which is fitted to a polynomial for simpler computation (Koop and Murray ^[38]):

where T_m is 273.15 K, and c_i can be expressed as six parameters. For more details, see Koop and Murray ^[38].

The most common type of IN is mineral dust particles transported into the atmosphere. It has been suggested that the dust particles (serve as IN) have a significant influence on cloud microphysics and dynamics (Van Den Heever et al. ^[24]; Cziczo et al. ^[39]; Fan et al. ^[40]). Soot particles and biological aerosols can also act as IN, but their number concentrations are generally low, and the ice nucleating efficiency is much weaker. Therefore, the heterogeneous freezing parameterization is designed and implemented as immersion freezing. DeMott et al. ^[41] assumed that the concentrations of IN active in mixed-phase cloud conditions can be related to aerosols (dust particles) larger than 0.5 μm. The approach is based on a "global" type of IN collected from multiple locations and may well describe drop freezing in the immersion mode for representative insoluble particles. The concentration of ice formed from supercooled liquid drops via immersion is described in the parameterized form mentioned by DeMott et al. ^[42]:

where N_IN is the number concentration of ice nucleation particles at T, T is the cloud temperature in Kelvin, N_aer is the number concentration of aerosol particles with diameters larger than 0.5 μm, and cf is the calibration factor, which by default is set to 1. a = 0, b = 1.25, c = 0.46, and d = –11.6.

During the model simulations, the number of ice particles formed in a time step ∆t depends on the ice nucleation rate. The equation to calculate the ice nucleation rate Q_i with respect to the time step ∆t(s) has the form:

where N_{i, new} is the number concentration of ice nucleation particles after ∆t, and N_{i, old} is the number concentration of ice nucleation particles before ∆t. This equation means that at the present (new) time step, if a higher ice number is predicted at a given model grid box than that at the last time (old) step, the ice number difference (divided by model time step ∆t) will be added to the prognostic equation. Otherwise, there will be no new ice formation.

Mountain thundercloud studies of electrical evolution in thundercloud cases are simulated. This case occurred between July and August of 1999 at the Langmuir Laboratory for Atmosphere Research in the mountain of central New Mexico (Coleman et al. ^[13]), and the National Centers for Environmental Prediction grids data provided atmospheric sounding profiles on July 31, 1999 around the observing site, as shown in Fig. 1 (Shi et al. ^[43]). Heat and a bubble over a flat terrain were initialized by a temperature disturbance of 3.5 K, and the humidity disturbance was 60%, which resulted in a normal convective cloud. The bubble was located in the domain center at a height of 1 km with a horizontal radius of 5 km and a vertical radius of 1 km. The initial setting for aerosol was 100 cm^–3. Hereafter, the experiments were classified into three cases, which were the homogeneous case (H), immersion case (I), and coupled case (C). In the homogeneous case, only homogeneous freezing was considered; similarly, only immersion freezing was considered in the I case. Moreover, a coupled case (C case) was performed under the condition that both freezing processes were considered.

Figure 1.  Initial sounding data for the case mentioned by Coleman et al. ^[13]. The red solid line represents the environment temperature, the green solid line represents the dew point temperature, and the black solid line represents the state curve. The wind vector is displayed on the right frame. The blue dashed lines are pseudo adiabats, and the red dashed lines are dry adiabats. The lifting condensation level is 5.59℃, and 652.51 hPa.

To investigate the effect of different ice nucleating processes on thundercloud properties, three cases were calculated for thunderstorms. Fig. 2 shows the temporal evolution of maximum and minimum vertical velocities in the three cases. In the H case, the maximum updraft was 9.36 m s^–1, occurring at 6.5 km in the 38th minute. The maximum updraft of the I case was 8.77 m s^–1, occurring at 6.25 km in the 35th minute. The maximum updraft of the C case was 11.66 m s^–1, occurring at 5.5 km in the 32nd minute. The vigorous development of thunderstorms was mainly in the time range from 15 to 60 minutes. The maximum updraft was an indicator of the largest local latent heat release. In the developing stage, the quick increase in updraft velocity mainly resulted from the release of latent. It also can be seen from Fig. 2 that the heterogeneous freezing process considered in the I case and the C case led to a fast convective development due to the faster latent heat release in the I and C cases compared with the H case. Therefore, one can conclude that the heterogeneous freezing process led to a faster development of convection. Moreover, due to the enhanced ice crystal particle production, the maximum updraft in the H case was larger than that in the I case. For the C case simulation, more latent heat release consisted of two nucleating processes, and thus, the maximum updraft was the strongest among the three cases. After that, the decrease in maximum updraft velocity was due to the falling precipitation particles. A slight difference in the maximum updraft velocity in the three cases was in the time interval of 50–80 minutes. This might be caused by the difference in precipitation behavior. From the evolution of the modeled maximum and minimum vertical velocities within the domain, the variations in the downdraft were roughly similar to that in the updraft within the first 40 minutes. The downdraft was closely connected with precipitation particles; the maximum values are shown in Fig. 2, and the space-time positions were (40 min, 6.75 km), (34 min, 6.25 km), and (24 min, 3.75 km), respectively. The maximum downdraft in the C case was lower than those in the cases of H and I, and this was due to the smaller ice particle production in the C case (discussed in the next section). After the thunderstorm cloud development for 40 minutes in the cases of I and C, the downdraft peaked again when the precipitation enhanced. Compared with the I and C cases, there was an absence of a downdraft center in the H case. This was because the larger the graupel particles, the lower the airflow would be. However, the H case produced a large number of small ice crystals at high altitudes, which were difficult to be converted into graupel particles because of their small scale. Therefore, the small number concentration of graupel particles led to the absence of a downdraft center in the H case after 50 minutes.

Figure 2.  Evolution of maximum and minimum vertical velocities in three cases. (a, b, c) The maximum and (d, e, f) minimum values of the represented quantities are given at the top of each plot. The black solid lines represent the isotherm (–40℃ and –15℃, and to 0℃).

The temporal development of the ice particles in the three freezing processes was rather different. Fig. 3 shows the time-height ice crystal mixing ratio in the three cases. In the H case, the ice particles mixing ratios over 5 g kg^–1 resided at altitudes from 5 to 10 km, where the temperature was below –20℃. In the dissipation stage of the thunderstorm, the mixing ratio of ice crystal exceeding 1 g kg^–1 was still available above 5 km. In the I case, because of the immersion freezing process arising, the cloud developed rapidly, and the ice crystal particles appeared after approximately 7 minutes (Fig. 3b). The ice crystal production from the immersion freezing process became active in the region between 2 and 5 km with temperatures from –20℃ to 0℃. Therefore, the immersion freezing process played an important role in ice crystal production at warmer temperatures. In the C case, the earlier appearance of ice crystals in the low region (where the temperatures were high) was associated with immersion freezing. After 32 minutes, the ice crystal resided between 6–8 km, probably attributing to the homogeneous freezing process. Furthermore, Fig. 4 shows the vertical distribution of the maximum ice crystal number concentration in the three cases. The maximum ice crystal number concentrations for the H, the C, and the I cases were 7.97×10⁸ kg ^–1, 9.17×10⁶ kg^–1, and 1.29×10⁶ kg^–1, respectively. The most noticeable difference in Fig. 4 is that the homogeneous freezing process predicted a significantly higher number concentration of ice crystals, while the ice crystal number concentration produced from the immersion freezing process was relatively weak. The peaks of the number concentration of ice crystal in the C case were 1.23×10⁶ kg^–1 (homogeneous freezing) and 9.17×10⁶ kg^–1 (immersion freezing), respectively, which was less than that in the H case and the I case. Therefore, it can be inferred from the above results that homogeneous freezing was dominant at colder temperatures, and immersion freezing produced frozen hydrometeors in advance between the isotherms of 0℃ and –20℃. When the temperature was lower than –12℃, heterogeneous nucleation can begin to occur with the participation of ice nucleation. Therefore, immersion freezing produced high values in the lower layers. In case C, which contained two processes, the competition between homogeneous and heterogeneous freezing appeared in about 32 minutes with the increase of the number concentration of ice crystals. Heterogeneous freezing began to consume cloud water in the early stage of development. When vertical airflow prevailed, the updraft brought the remaining droplets into the low-temperature layer and froze into the ice by homogenization, which not only reduced the efficiency of heterogeneous nucleation but also inhibited the sublimation growth of ice crystals at the bottom. At the same time, small ice crystals with lightweight formed by heterogeneous nucleation were brought into the low-temperature region dominated by homogeneous nucleation and compete with homogeneous nucleation for cloud water in the middle and high levels. Therefore, homogeneous freezing is the major source of ice crystals in the cloud anvil, so the cloud anvil is composed of a large number of small size of ice crystals, which is consistent with the results of previous studies (Ekman et al. ^[25]; Phillips et al. ^[26]; Seifert et al. ^[44]).

Figure 3.  Spatial and temporal distribution of the ice crystal mixing ratio in (a) homogeneous freezing case, (b) immersion freezing case, and (c) couple freezing case. The black solid lines represent the isotherm (at 0℃, –15℃, and –40℃).

Figure 4.  Vertical distribution of the maximum ice crystal number concentration in the three cases.

The vertical variation of the maximum ice nucleation rate in three cases is shown in Fig. 5. The homogeneous freezing process in the H case occured between 7–8 km, while the ice crystal produced from immersion freezing (I case) resided in the height of 2–3 km. The maximum homogeneous nucleation rate (H case) was about 0.35 g kg^–1 s^–1, which was significantly greater than the maximum immersion nucleation rate (2.8×10^–3 g kg^–1 s^–1) in the I case. In the C case, two peaks in height can be found in Fig. 5c. This can be consistent with two nucleation processes. The peak of nucleation rate (the homogeneous freezing process) at about 7 km was 0.01 g kg^–1 s^–1, while the peak of nucleation rate (the immersion freezing process) at about 2 km was 7.8×10^–4 g kg^–1 s^–1. Therefore, both the two nucleation freezing rates were lower than those in the H and the C cases, and thus it can be inferred that water vapor competition can be found between the two nucleation freezing processes. The simulation result is consistent with the relevant studies (Ekman et al. ^[25]; Phillips et al. ^[26]; Chen et al. ^[45]; Li et al. ^[46]; Jensen et al. ^[47]).

Figure 5.  Vertical distribution of the maximum ice nucleation rate in the three cases. (a) Homogeneous nucleation case (case H), (b) immersion nucleation case (case I), and (c) two nucleation processes are considered (case C).

The time evolution of the mixing ratio and number concentration of cloud droplets, rain, and graupel particles for the three cases are presented in Fig. 6. After about 8 minutes of simulation, the cloud droplets were activated from aerosol particles, and after that, the cloud droplets were lifted by the strong updraft. The largest difference in the cloud droplet content among the three cases occured at the 20th minute. The higher cloud droplet content in the H case extended up to the 7 km level (Fig. 5a), indicating less consumption compared with the I case and the C case, in which the cloud droplet content decreased more rapidly with height because of a more efficient collision between cloud droplet and ice particles (Figs. 6b and 6c). The ice crystal in the H case was produced after 30 minutes (Fig. 3a), yielding less water vapor consumption, and therefore, the stronger supersaturation in the H case was responsible for the higher cloud droplet production. Hence, the maximum mixing ratio of cloud droplets in the three cases was 7.4 g kg^–1, 3.2 g kg^–1, and 2.6 g kg^–1, respectively.

Figure 6.  Time-height plots showing the mixing ratio (shaded region) and number concentration of cloud droplet (a, b, and c), rain (d, e, and f) and graupel (g, h, and i) evolution in the three cases. Isotherms (thin horizontal black lines) at 0℃, –15℃, and –40℃ are the same in all panels. Cloud droplet number concentration with contour intervals of 10⁵, 10⁶, 10⁷, 10⁸ kg^–1, rain number concentration with contour intervals of 10 ², 10³, 10⁴, 10⁵ kg^–1, and graupel number concentration with contour intervals of 10 ², 10³, 10⁴, 10⁵ kg^–1.

Two microphysical processes (including auto-conversion of cloud droplet rain and melting of graupel) had a great impact on the growth of rain. Raindrops in the H case (before the 50th minute) were primarily distributed above 2 km (Fig. 6d), indicating that they were mainly formed by the autoconversion of cloud droplets. After about 50 minutes, the raindrop in the H case can be formed from the ice particle melting process. However, the raindrops in the I and C cases below 2 km were produced from melting ice particles (Figs. 6e and 6f). These can be explained that the immersion freezing in the I and C cases led to an earlier production of ice crystals, yielding a host of large sizes of ice particles, and thus, the melting of ice particles in the region where the temperature was above 0℃ was produced. It is evident from Figs. 6b and 6c that the mixing ratio of cloud droplets in the I and C cases was less than that in the H case. This effect from small cloud droplets led to the restriction of the autoconversion of cloud droplets. Therefore, the raindrops were not produced above 2 km. Moreover, the raindrop production in the I case was slightly larger than that in the C case (Figs. 6e and 6f). This can be explained by the fact that more small ice particles (ice crystals) in the C case led to stronger cloud water competition, resulting in the production of graupel particles in the C case being less than that in the I case, and thus the melting of the large size of ice particles was limited in C case.

Figs. 6e–6i show the temporal evolution of the maximum graupel mixing ratio and number concentration for the three cases. Graupel was firstly produced by auto-conversion of ice-graupel. In the H case, the reduction in the size of the ice crystal arose from the enhancement of ice crystal number concentration, and therefore, the ice crystal was harder to convert to graupel particles. The graupel particles in the H case emerged at about the 30th minute, which was later than those in the I and C cases (about 10 min). Additionally, the number concentration of graupel particles in the H case was 1.6×10⁴ kg^–1, much less than those in the cases of I and C (3.0×10⁵ and 2.9×10⁵ kg^–1, respectively). However, the maximum mixing ratio of graupel in the three cases was 9.1 g kg^–1, 9.4 g kg^–1, and 8.1 g kg^–1, respectively. Hence, it can be inferred that the size of the graupel in the H case was larger than those in the cases of I and C. Furthermore, the size of the graupel in the C case was smaller than that in the I case, probably attributing to the water vapor competition in the C case. It seems reasonable to conclude that the H case which only considered the homogeneous freezing process produced large graupel particles, but the number concentration was small. The immersion freezing caused a large size of ice crystal production, contributing to the higher number concentrationof graupel. The water vapor competition caused by considering two freezing processes was responsible for a slight reduction of graupel growth. A similar result can be found in the study by Diehl and Grützun ^[48].

Figure 7 shows the time evolution of non-inductive charging rates by ice crystal and inductive charging rates by graupel. It can be seen from Figs. 7a–7c that the positive non-inductive charging rates were the highest at altitudes of 4–8 km, while the negative non-inductive charge mainly resided at 2 to 4 km. Therefore, it is clear that ice crystals charged positively at lower temperatures (< –15℃), and ice crystals gained a negative charge in the regions where the temperature was higher (> –15℃). As smaller ice crystal production in the H case caused the delay times for graupel formation, the appearance time of the non-inductive charging process in the H case was much later. In contrast, the immersion freezing process in the I and C cases supported the early ice particle production, leading to the early non-inductive charge separation between ice crystal and graupel (see Figs. 7a–7c). Table 1 gives the maximum non-inductive charging rate in three cases. The H case predicted almost six times stronger positive non-inductive charging rates than those in the cases of I and C due to the enhancement of ice crystal number concentration and the large size of the graupel. Because heterogeneously nucleated ice crystals in the cases of I and C resided in the low levels of the storm where the temperature was above –15℃, ice crystals charged negatively between the height of 2.5–4.5 km (Figs. 7b and 7c). The negative non-inductive charging rate in the I case was much larger than that in the C case (see Table 1). This was attributed to stronger ice particle production. Furthermore, the water vapor competition between two nucleation processes led to speeding the cloud water consumption, and after that, the ice crystals were charged positively by non-inductive charging under the condition of low cloud water content. Therefore, it can be found from Figs. 7b and 7c that the time interval of the negative non-inductive charging process in the I case was between about 12–43 minutes, while the time evolution of the negative non-inductive charging process was roughly between 12–35 minutes in the C case.

Figure 7.  Time-height plots showing the non-inductive and inductive charging rate evolution in the three cases. (a–c) Non-inductive charging between ice crystal and graupel in three cases with contour intervals of ±1, ±5, ±10, ±50, ±100, ±600 pC m^–3 s^–1. (d–f) Inductive charging between graupel and cloud droplets in three cases with contour intervals of ±0.1, ±1, ±10 pC m^–3 s^–1. Isotherm (horizontal black lines at 0℃, –15℃, and –40℃) is in the panel. See Table 1 for maximum and minimum values and their space-time distribution.

The inductive charging rates in three cases is listed by the time-height plots in Figs. 7d–7f. The inductive charging rates roughly resided between 2–8 km. Similar to the non-inductive charging process, the later development of the inductive charging process (starting at about the 35th minute) can be found in the H case. It can be found from Figs. 7d–7f that the positive and negative inductive charging rates were greater than those in the cases of I and C. Two factors that affected the inductive charging process were electric field and hydrometer particle content. An increase in vertical electric field arising from non-inductive charging enhancement in the H case can promote the charge separation between graupels and cloud droplets, and greater production of cloud droplets and graupels below the altitude of 5 km led to stronger inductive charge separation in the H case. Therefore, the inductive charging rate in the H case was significantly greater than those in cases I and C (see Table 1).

The three sensitivity tests exhibited similar charge structures. Fig. 8 shows the charge structure in the developing stage at the 25th minute, the mature stage at the 40th minute, and the dissipate stage at the 55th minute. It can be found in Figs. 8a–8c that the dynamic structures in the thunderclouds were roughly similar, and the cloud top reaches about 5.5 km. In the developing stage, the strong updraft resided in a height range of 2–5.5 km, while the downdraft emerged at 4 km. Since the ice particles in the H case were not produced at the 25th minute, the electrification process was still not found in Fig. 8a. Figs. 8b and 8c show that both the I and C cases had an inverted dipole charge structure consisting of a strong positive charge region at about 2 to 4 km (higher in the updraft region) and strong negative charge above 4 to 6 km. The cause of the inverted charge structure from the I and C cases was clearly the positive non-inductive charging of graupel and negative non-inductive charging of ice crystal (see Figs. 9b and 9c). It also can be seen from Figs. 9b and 9c that the charge obtained by ice particles in the I case was stronger than that in the C case, because of the enhancement of ice particle production. Moreover, because the inductive charging processes at the 25th minute in the in three cases were not profound, the charging of cloud droplets was less than that of graupel and ice crystal (See Fig. 9a–9c).

Figure 8.  Charge structures of the three cases at the 25th, 40th and 55th minutes. (a, d, g) the H case, (b, e, h) the I case, and (c, f, i) the C case. Storm-relative wind vectors are shown only where the speed is greater than 1 m s^–1, and the cloud water content is over 0.01 g kg^–1. Thick black lines show the contour structure characteristics of thunderclouds.

Figure 9.  Vertical distribution of charged hydrate particles in the three cases (nC m^–3). Blue, red, and green lines represent the charge of ice, graupel, and cloud droplet, respectively.

The thunderstorms reached the mature stage at about the 40th minute. The convection strength and charge structure in the three cases showed a significant difference. The H case developed a strong vertical draft. The storm in the mature stage depicted a triple charge structure with an upper positive charge region (Fig. 8d). As shown in Fig. 8d, the main positive region was located near the cloud top with the strong wind shear, which corresponded to the positive charging by ice crystal (Fig. 9d). Furthermore, a strong negative charge region resided between 5–7 km where the wind structure was the dominant mode of updraft pulse, composed of graupel charged negatively by non-inductive charging (Fig. 9d). It also can be seen from Fig. 8d that a positive region (the so-called lower positive charge or LPC) also appeared below 5 km. The downdraft below 5 km in the H case was greater than that in the I and C cases, which was attributed to the cloud water contents. Under this condition, the inductive graupel-cloud droplet charge separation had the effect of enhancing the main negative and the lower positive charge regions. Moreover, due to the profound graupel production in the I and C cases, the I and C cases both showed a larger region of strong electrification than that in the H case (Figs. 8e and 8f). The homogeneous freezing process in the C case produced higher amounts of ice crystal at the height of 6–9 km. The greater positively charged ice particles in the C case contributed to the main positively charged region than in the I case (Figs. 9e and 9f). In addition, the lower positive charge region in the C case was "missing, " and thus, the charge structure in the C case was a "normal" polarity charge structure consisting of a strong negative region at 2 to 5 km altitude and a strong positive charge above 5–8 km (Fig. 8f). This can be explained that the water competition in the C case resulted from both homogeneous freezing process and heterogeneous freezing process led to the rapid cloud droplet consume. The inductive charging process in the C case at 2 to 4 km was too weak to contribute to the lower positive charge region.

At the 55th minute of simulation, the cases of I and C with larger ice particle production predicted the great downdraft below 4 km, while the strong upper-level updraft can be found above 4 km in the H case. It can be seen from Fig. 8g that the H case at the 55th minute had a triple charge structure consisting of a weak lower positive charge region. The charge structure was mainly developed due to the negatively charged by graupel and the positively charged by ice crystal (Fig. 9g). Moreover, the I case developed a lower dipole structure (see Fig. 8h). The main positive charge region and negative charge region were mainly located about 2–4 km where large size of ice particles led to profound collision between graupel and ice particles (Fig. 9h). When the homogeneous freezing process was considered in the C case, ice particles contributed by homogeneous freezing near cloud tops was charged positively between 6–8 km (Figs. 8i and 9i). Furthermore, the enhanced downdrafts led to stronger evaporation or melting of precipitating particles, and thus, the graupel particles charged negatively in the lower layer became weaker (Fig. 9i).

Figure 10 illustrates the spatial and temporal distribution of the charge structure of thunderstorms in the three cases. Compared to the homogeneous freezing process, the heterogeneous freezing process led to earlier electrification. In the H case, the charge structure developed a normal triple with a strong positive charge region from the 40th to 50th minute. Additionally, at the simulation time of 60th to 70th minute, the storm in the H case had a complex charge structure (four layers) (Fig. 10a), probably attributing to the stronger inductive charging between the graupel and the cloud droplet. After that, the storm in the H case developed a normal dipole structure. The electrification in the I and C cases began at about the 17th minute, with an inverted dipole structure first produced. By the 30th minute, the storms in the I and C cases exhibited a triple structure. It can be seen from Figs. 10b and 10c that the immersion freezing process in the I case produced a stronger lower positive charge region than that in the C case, while the homogeneous freezing process in the C case led to a stronger main upper positive region (resided between 6–8 km) production than that in the I case. In the dissipating stage, both the I and C cases produced a normal dipole structure. As the larger size of ice particles produced from the immersion freezing process in the I case tended to be decent, the dipole structure in the I case only resided below 5 km, which was significantly lower than 8 km in the C case.

Figure 10.  The evolution of charge structure in the three cases. (a) H case, (b) I case, and (c) C case. The maximum and minimum values are provided and height of their achievement are in the following brackets (units: min, km).

In this paper, the influence of different ice nucleating processes on dynamics, microphysics, and electrification in thunderstorms has been studied using a two-dimensional cloud model. Three model runs have been carried out driven by three ice nucleating processes. The analyses of results demonstrate that the different ice nucleating processes have a significant influence on the thunderstorm cloud microphysical processes and electrification.

It is found that the heterogeneous freezing process led to quick convective development because of more latent heat release. When both nucleating freezing processes were considered, more latent heat was released in the thunderstorm; thus, the maximum updraft was the strongest among the three cases. However, the storm, including homogeneous and heterogeneous freezing processes, developed a weaker downdraft than the storm.

Homogeneous freezing was dominant at the level above –35℃, while immersion freezing produced frozen hydrometeors in advance between the isotherms of 0℃ and –20℃. Homogeneous freezing was the major source of ice crystals in the cloud anvil, so the cloud anvil was composed of a large number of small size of ice crystals, while the heterogeneously formed was a small number of large size of ice crystals. Additionally, water vapor competition can be found between two nucleation freezing processes, which might inhibit both homogeneous and heterogeneous ice nucleation. Water vapor was depleted due to the enhancement of condensation on the heterogeneously formed ice crystals. Therefore, the formation of cloud droplets was stronger for the homogeneous freezing case than for the two other cases. This led to more autoconversion of cloud droplets in the homogeneous freezing case, while the heterogeneous ice nucleation, in general, can lead to less cloud droplet production, which resulted in only the melting of ice particles (rain formation) for the storm simulation with the heterogeneous freezing process. Moreover, as the small ice crystal was harder to convert to graupel particles, the graupel particle production in the homogeneous ice freezing case was less than that in the other two cases.The heterogeneous freezing process caused a large size of ice crystal production, contributing to the higher number concentration of graupel. The water vapor competition caused by considering two freezing processes was responsible for a slight reduction of graupel growth.

The electrification in thunderstorms was dramatically affected by freezing processes. It has been found that homogeneous freezing can maintain the ice particles near cloud tops, supporting the important role in the main positive charge regions. Since smaller ice crystal production from homogeneous freezing caused the delay times for graupel formation, the appearance time of the electrification process was much later. The heterogeneous freezing led to fast ice particle development and sufficient supplies of early non-inductive and inductive charging processes. Because heterogeneously nucleated ice crystals resided in the low levels of the storm where the temperature was above –15℃, ice crystals charged negatively, and the charge structure in the developing stage tended to be an inverted dipole. Additionally, the heterogeneous freezing process played a role in the development of lower charge regions, and the lower dipole charge structure can be found in the dissipating stage when only the heterogeneous freezing process was considered. The competition between the two freezing processes led to less cloud water content, leaving weaker electrification. Under this condition, the graupel particles were charged negatively by non-inductive charging under the condition of low cloud water content (lower than 0.1 g kg^–1) in the low levels of the storm, and thus, two freezing processes were also found to be capable of reverting the charge structure toward dipole structure in the early simulation time.

With a weak initial humidity and updraft velocity in a simulated thunderstorm, the competition between freezing processes was obvious in the electrical development of clouds. Future research need to consider modifying other parameters, such as humidity, aerosols, and the ascent of the air, to provide further insights. Additionally, heterogeneous freezing consists of more than one type of nucleation scheme and more heterogeneous freezing nucleation schemes might be presented in the forthcoming study.

where u and w are the velocity components in the x and z directions, respectively; C_p is the specific heat of air at constant pressure; θ is the potential temperature; θ_v0 is the virtual potential temperature; P is the nondimensional perturbation pressure from the initial state; g is the gravitational acceleration, and the six hydrometeors mixing ratio (water vapor, cloud droplet, rain, graupel, hail and ice crystal) are represented by Q_v, Q_c, Q_r, Q_g, Q_h, and Q_i. Furthermore, the term D represents turbulent diffusion and is evaluated with a prognostic equation for turbulent kinetic energy (Klaassen et al. ^[49]).

where β is the quasi-elastic coefficient; R_d is the gas constant for dry air; ρ₀ is the air density; C_v is the specific heat of air at constant volume; ρ₀ is the nondimensional pressure.

where T is the atmospheric temperature.

The prognostic equations denote the mixing ratio (Q_x ) of water vapor, cloud droplet, rain, graupel, hail, ice crystal and the number concentration (N_x) of cloud droplet, rain, graupel, hail, ice crystal. S_Mx denotes the source and sink terms for each hydrometeor and V_x is the average fall speeds.

The cloud droplet spectral width (F_c ) is a parameter that is only used to calculate the conversion of cloud droplets to rain, and it has no relationship with the mixing ratio and number concentration of cloud droplets.

where Q_ex is the charge density carried by each hydrometeor.

where D is diameter; N₀, and λ are two parameters for different hydrometeors: (1) cloud droplet, α=2; (2) rain and graupel, α=0; (3) ice crystal, α=1; (4) hail, α=0.

where A_mr=0.524gcm^–3, A_mg=0.065g cm^–3, A_vr=2100cm^0.2 s^–1, A_vg=500 cm^0.2 s^–1, and α₁ =0.286. V_r(g) is the average velocity of raindrop and graupel.

where A_mi=0.001 g cm^–2, A_vi=70 cm^2/3 s^–1, and α₂ =0.3.

The hail is assumed to have an interceptive gamma function distribution of diameter:

where D_* is assumed to 0.5 cm denotes the minimum diameter of hail.

where A_mh=0.471 g cm^–3, and A_vh=810 cm^0.2 s^–1.

where N_ccn is the number of activated cloud condensation nuclei (CCN) and S is the supersaturation of the cloud, k is a constant and depends on the chemical composition and physical properties of the aerosol. This paper, like Wang ^[50], refers to k as 0.7. For simplicity, C_o is the concentration of CCN activated numbers under 1% supersaturation and is used to indicate the initial aerosol concentration (Li et al. ^[51]) in each numerical experiment. On this basis, this paper adds a diagnostic process to ensure that the model conforms to common sense.

where N_c^new is the calculated activation cloud droplet number concentration within a new time step (∆t), N_c^old is the cloud droplet number concentration at the former time step, a new cloud droplet forms when N_c^new > N_c^old, and the activation rate is N_c.

In this paper, the CCN type is assumed to be sulfate.

where Q₀ and N₀ are the mixing ratio and number concentration of cloud droplets, respectively. D represents its diameter, L_v is the latent heat in the process of phase transition, Q_v−Q_sw is supersaturation of water vapor. Besides, k_d, and k_T are functions of temperature. They represent the water vapor diffusion coefficient, and thermal conductivity in the air, respectively.

where A_mi=0.001 g cm^–3, Q_i and N_i are the mixing ratio and number concentration of ice crystals, respectively. Q_v−Q_si represents water vapor supersaturation on the ice surface. a₁ and a₂ are functions of temperature.

where A_mr=0.524 g cm^–3, and A_vr=2100 cm^0.2 s^–1. L is the phase transformation latent heat, L_f, L_v, and L_s are the the latent heat by freezing, condensing and sublimating, respectively. Besides, k_d, k_T, and μ are the functions about temperature. They represent the water vapor diffusion coefficient, thermal conductivity and kinetic viscosity coefficient in the air, respectively.

where Amg=0.065 g cm^–3, and Avg=500 cm 0.2 s^–1. Q_v − Q_s0 is water vapor supersaturation corresponding to freezing point. S_vgw and S_vgd are the wet and dry growth of graupels, respectively. Q_s0 is the saturated specific humidity at freezing point. C_cg, and C_rg represents collision between cloud droplet and graupel, and collision between rain and graupel, respectively.

where D_*=0_vh=810 cm^0.2s^–1. S_vbw and S_vbd are the wet and dry growth of hail, respectively. C_ch, and C_rh represents collision between cloud droplet and hail, and collision between rain and hail, respectively.

where D₁^*=15 μm, $A_{\mathrm{vi}}=70 \mathrm{~cm}^{\frac{2}{3}} \mathrm{~s}^{-1}$, and A_mi=0.001 g cm^–2. E_ci is the collision coefficient. E_ci can be calculated by using the mass median diameter $D_{\mathrm{i} 0}=\frac{3.67}{\lambda_{\mathrm{i}}}$ of the ice crystal.

where $A_{\mathrm{mr}}=\frac{\pi}{6} \rho_{\mathrm{w}}=0.524 \mathrm{~g} \mathrm{~cm}^{-3}$, A_vr=2100 cm^0.2 s^–1, and α₁ =0.286.

where A_mg=0.065 g cm^–3, A_vg=500 cm^0.2 s^–1, and α₁ =0.286.

where A_mh=0.471 g cm^–3, A_vh=810 cm^0.2 s^–1, D_*=0.5 cm.

where E_ri = 0.8.

where the coefficient of coalescence between hail and ice crystals is below:

where kk is an indicator of hail growth status. When the cloud droplets and raindrops totally collide and freeze into hail, kk is set to 0. When the cloud droplets and raindrops partly collide and freeze into hail, kk is set to 1.

where the coefficient of coalescence between graupel and ice crystals is related to the temperature and the surface state of graupel. For simplicity, E_ig = 0.1.

where the coefficient of coalescence between hail and ice crystals is below:

where the first part on the right of the upper formula is the secondary raindrop which is produced by raindrops colliding with each other and breaking up. S_n is the average secondary drop produced by each collision and breakage and S_n=3.

where