A Numerical Study of Aerosol E ﬀ ects on Electriﬁcation with Di ﬀ erent Intensity Thunderclouds

: Numerical simulations are performed to investigate the e ﬀ ect of varying CCN (cloud condensation nuclei) concentration on dynamic, microphysics, electriﬁcation, and charge structure in weak, moderate, and severe thunderstorms. The results show that the response of electriﬁcation to the increase of CCN concentration is a nonlinear relationship in di ﬀ erent types of thunderclouds. The increase in CCN concentration leads to a signiﬁcant enhancement of updraft in the weak thunderclouds, while the high CCN concentration in moderate and severe thunderclouds leads to a slight reduction in maximum updraft speed. The increase of the convection promotes the lift of more small cloud droplets, which leads to a faster and stronger production of ice crystals. The production of graupel is insensitive to the CCN concentration. The content of graupel increases from low CCN concentration to moderate CCN concentration, and slightly decreases at high CCN concentration, which arises from the profound enhancement of small ice crystals production. When the intensity of thundercloud increases, the reduction of graupel production will arise in advance as the CCN concentration increases. Charge production tends to increase as the aerosol concentration rises from low to high in weak and moderate thundercloud cases. However, the magnitude of charging rates in the severe thundercloud cases keeps roughly stable under the high CCN concentration condition, which can be attributed to the profound reduction of graupel content. The charge structure in the weak thundercloud at low CCN concentrations keeps as a dipole, while the weak thunderclouds in the other cases (the CCN concentration above 100 cm − 3 ) change from a dipole charge structure to a tripole charge structure, and ﬁnally disappear with a dipole. In cases of moderate and severe intensity thunderclouds, the charge structure depicts a relatively complex structure that includes a multilayer charge region.


Introduction
In recent years, considerable progress has been made in understanding the interaction between aerosol and lighting. Forest fires inject aerosol particles into the troposphere, allowing an investigation of the response of lightning activities on aerosol. It has been found that smoke aerosol released from fire likely led to the noticeable increasing percentage of positive cloud to ground (CG) lightning [1][2][3][4][5][6][7]. The thermodynamic and dynamic contrast between continental and maritime can be largely traced to the difference of lightning production between maritime and continental thunderclouds, and the CCN is also hypothesized to be an important factor [7,8]. The CCN hypothesized for the land-ocean lightning contrast has been developed by a well-established contrast in boundary layer aerosol concentration [9]. However, Williams and Standfill [8] suggested that the surface Bowen ratio is responsible for the land-ocean lightning contrast. Recent observations and theoretical studies have where N 0 is the initial CCN concentrations on the ground, N ccn (z) is the CCN concentration at a different height, z s (CCN elevation) is set to 2 km in this study [32]. In addition, we assume that the CCN concentration at some height layer is uniformly distributed when CCN concentrations are initialized in each grid point. Therefore, the CCN concentrations' distribution for the whole space can be calculated by given CCN concentrations on the ground. Soluble CCN particles (sulfates, etc.) can grow into cloud droplets, which affects the microphysical development of thundercloud. This paper couples a classical CCN activation parameterization in the improved 2D cloud model: where N a represents activated CCN number concentration and S represents the water vapor supersaturation within the cloud. k of Equation (2) is related to the chemical composition and physical properties of aerosols, where k is 0.7 [25]. C 0 is used to represent the initial CCN concentrations [26]. Thus, it can be found that the number of cloud droplets activated by CCN is not only related to the initial concentration of CCN but also depends on the water vapor supersaturation within the cloud. On this basis, we add a diagnostic procedure given by Tan et al. [23] to determine whether there are new cloud droplet generations.
where N c represents the activation rate of cloud droplets, N c old is cloud droplets concentration at the former time step, N c new represents activated cloud droplets calculated within a new time step (∆t).
Therefore, new cloud droplets are generated when N c new > N c old in this study.
Nucleation of ice crystals occurs both through homogeneous-freezing of cloud droplets and by heterogeneous nucleation caused by IN (ice nuclei). However, it must be clarified that an empirical formula of ice nucleation based on Fletcher [33] is considered in the original model, while the freezing process of cloud droplets is not be included. Therefore, this work makes an improvement for the existing model, and the process of cloud droplets freezing into ice crystals is added into the model. The freezing scheme of cloud droplets developed by Sun et al. [34] is expressed as: Atmosphere 2019, 10, 508 4 of 20 where Q cif represents the nucleation rate of cloud droplets, that is the mixing ratio of ice crystals formed by cloud droplets freezing, Q c is the cloud-water mixture ratio, the parameter of pc is the percentage of frozen droplets, and pc varies by temperature, pc = 0.008 × (1.274) T 0 −T . T is the temperature within the cloud, and T 0 (−20 • C) is the temperature at which the water droplets begin to freeze. The water droplets have partially frozen from −20 • C and completely frozen at −40 • C. The heterogeneous nucleation scheme is based on Cooper [35]: where N in is the heterogeneous nucleation rate, and the number of ice crystals formed by heterogeneous nucleation is only related to temperature. This study also considered the influence of water vapor supersaturation referring to Hu et al. [31]. Therefore, Equation (5) is rewritten as: When dT dt < 0 and Q v > Q si 0 When dT dt ≥ 0 and Q v ≤ Q si , (6) where k = 5, dT dt ≈ w ∂T ∂t , w represents wind speed, and S i is for the water vapor supersaturation on the ice surface. Therefore, the heterogeneous nucleation rate is determined by both temperature and water vapor supersaturation.
Raindrops can be produced from auto-conversion of cloud droplets or the melting of ice particles. The microphysical conversions include condensation and evaporation of five types of particles. The secondary ice multiplication depends on the number of cloud droplets with a diameter greater than 24 µm, which is based on Mossop [36]. Since a large amount of supercooled water resides in cumulus cloud, ice crystals with the diameter over the critical value (0.03 cm) will grow up quickly during colliding with supercooled water, and the growth speed is larger than that in the process of sublimation. Under this condition, as soon as the diameter of ice crystal exceeds 0.03, ice crystal should become graupel. Similarly, after a significant growth of graupel (the diameter exceeds 0.5 cm), the graupel should be assigned to hail. In addition, the microphysical processes include various collisions between two types of particles (the collision between cloud droplet and ice crystal, rain, graupel, and hail; the collision between rain and ice crystal; the collision between rain and graupel, and hail; the collision between ice crystal and graupel, and hail).
The initial perturbation on the convective cloud simulation that is with a horizontal uniform horizontal flow is used to start the convection. In this article, the way to onset convection is setting temperature and specific humidity perturbations. This setup adds a limited range of warm-wet bubble perturbation in the lower layer of the model domain, as a result, the humidity and temperature of this warm-wet bubble are higher than that of the surround. Then the buoyancy in the atmospheric equation of vertical motion is used to start the initial convection. Perturbation function is as follows: here, (x c , y c ) is the central coordinate of perturbation, x r , y r is the radius of two coordinate directions on the warm-wet bubble. ∆θ is the central maximum disturbance potential temperature. Q vn is saturation specific humidity under the condition of temperature perturbation. Therefore, we modulate the intensity of the convective system by changing center potential temperature and saturation specific humidity on disturbance.

of 20
Electrification via induction in the model occurs when graupel collides with cloud droplets/ice crystals. Inductive collision charging parameterization is based on Ziegler et al. [37]. The inductive graupel-cloud droplet charging rate can be given by: where Q eg is the individual charge from graupel, D c/i and D g are the diameter of cloud droplets and graupel, respectively. V g is the falling speed of graupel, N c/i and N g are the cloud droplet/ice crystal and graupel concentrations, respectively. N 0g is the number concentration intercept for graupel. C is the complete gamma function, and E z is the vertical electric field. The symbols E gc/i and E rc/i denote graupel-cloud droplet/ice crystal collision efficiency and rebound probability, respectively. It is the polar collision angle. According to Mansell et al. [38], coefficients for inductive graupel cloud droplet charging (E rc = 0.01 and cos θ = 0.4) in this study are between the moderate to strong values which range from E rc = 0.007 to 0.015 and cos θ = 0.2 to 0.5. In addition, coefficients for inductive graupel-ice crystal charging E ri = 0.7 and cos θ = 0.2. The non-inductive graupel charging rate takes the form: where D g and D i are diameters of the colliding particles (graupel and ice crystal). E r abounds probability, which is assumed to be 0.2, while E gi is graupel-ice crystal collision efficiency (set to be 0.8).
N is the number concentration. V i and V g are the mass-weighted mean terminal speeds for ice crystal and graupel. An arbitrary factor β is used to limit charging at low temperature and given is similar to Mansell et al. [38] by: The scheme of individual graupel charging rate (δq) we adopted is modified from the scheme of Gardiner et al. [39] based on the experiment results of Pereyra et al. [40]. The reversal temperature in the modified parameterization scheme, replacing the original fixed value, considers the functional dependence on cloud water content (CWC), so not only positive charges but also negative charges may be acquired by graupel when CWC is less than 1 g m −3 . The modified parameterization scheme is denoted as follows.
The expression of δq following Gardiner et al. [39] is approximated as: where D i is the diameter for ice, δL is a parameter related to CWC, modified as: where τ = (−21/T r ) (T − 273.16) is the scaled temperature used by Ziegler et al. [37] to allow the reversal temperature T r to be varied. It must be stated that T r (−15 • C) is in • C and T is in K. Electrification parameterizations follow Tan et al. [23,41] including non-inductive and inductive charge separation processes. The non-inductive charging parameterization is modified from Gardiner et al. [39]. It mainly takes into account the collision between graupel and ice crystals, which is related to the size of ice crystals and graupel, drops velocity, relative collision rate, and the coefficient of the inversion temperature. Lightning discharges are parameterized based on Tan et al. [42]. Lightning initiation uses the runaway electron threshold for break-even field and thereafter bidirectional channels are propagated in a stochastic step-by-step fashion. The leaders of IC (intracloud) lightning do not reach the ground, and a height threshold (1.5 km or six grid points above ground) is used to define a flash to be a CG (cloud to ground) lightning (includes +CG and −CG).

Model Initialization and Method in Numerical Experiments
A mountain thundercloud SEET (Studies of Electrical Evolution in thunderclouds) case occurred on 31 July 1999 and was used for numerical simulation in this study. The environmental temperature, humidity stratification, and vertical wind profile at the corresponding time are shown in Figure 1. A humid and warm bubble of horizontal radius 5 km and vertical radius 1 km was given at the grid point with a height of 1 km. Three kinds of thunderclouds were set by changing the environmental humidity and temperature stratification of the SEET case at an initial time. The temperature disturbance and humidity disturbance of weak thundercloud were set as 2.5 k and 60%, that of moderate thundercloud were set as 4.5 k and 80%, and 6 k and 100% in Equations (7) and (8). The different CCN concentrations were characterized by changing N 0 in Equation (4). Four numerical simulation experiments (N 0 was 100, 500, 1000, and 3000 cm −3 respectively) were conducted. Simulations were performed in a 76 km × 20 km domain with a constant grid spacing of 250 m in the horizontal and 250 m in the vertical from the surface to 20 km. Moreover, all simulations were carried out for 80 min, and the time step was 2 s. Electrification parameterizations follow Tan et al. [23,41] including non-inductive and inductive charge separation processes. The non-inductive charging parameterization is modified from Gardiner et al. [39]. It mainly takes into account the collision between graupel and ice crystals, which is related to the size of ice crystals and graupel, drops velocity, relative collision rate, and the coefficient of the inversion temperature. Lightning discharges are parameterized based on Tan et al. [42]. Lightning initiation uses the runaway electron threshold for break-even field and thereafter bidirectional channels are propagated in a stochastic step-by-step fashion. The leaders of IC (intracloud) lightning do not reach the ground, and a height threshold (1.5 km or six grid points above ground) is used to define a flash to be a CG (cloud to ground) lightning (includes +CG and −CG).

Model Initialization and Method in Numerical Experiments
A mountain thundercloud SEET (Studies of Electrical Evolution in thunderclouds) case occurred on 31 July 1999 and was used for numerical simulation in this study. The environmental temperature, humidity stratification, and vertical wind profile at the corresponding time are shown in Figure 1. A humid and warm bubble of horizontal radius 5 km and vertical radius 1 km was given at the grid point with a height of 1 km. Three kinds of thunderclouds were set by changing the environmental humidity and temperature stratification of the SEET case at an initial time. The temperature disturbance and humidity disturbance of weak thundercloud were set as 2.5 k and 60%, that of moderate thundercloud were set as 4.5 k and 80%, and 6 k and 100% in Equations (7) and (8). The different CCN concentrations were characterized by changing N0 in Equation (4). Four numerical simulation experiments (N0 was 100, 500, 1000, and 3000 cm −3 respectively) were conducted. Simulations were performed in a 76 km × 20 km domain with a constant grid spacing of 250 m in the horizontal and 250 m in the vertical from the surface to 20 km. Moreover, all simulations were carried out for 80

Dynamic and Microphysical Processes
Three types of thunderclouds with 12 cases were carried out until the thundercloud dissipated (usually 80 min). The time evolution of maximum and minimum vertical velocities (see Figure 2) in 12 cases is roughly similar, and the variations in updraft and downdraft experience from quick

Dynamic and Microphysical Processes
Three types of thunderclouds with 12 cases were carried out until the thundercloud dissipated (usually 80 min). The time evolution of maximum and minimum vertical velocities (see Figure 2) in 12 cases is roughly similar, and the variations in updraft and downdraft experience from quick enhancement to slow reduction. The maximum updraft velocity increases greatly when the intensity of thunderclouds increases from weak to severe (see Figure 2), primarily because of the dramatic convection activities. Similarly, the downdraft speed increases with thundercloud intensity increasing Atmosphere 2019, 10, 508 7 of 20 because of a larger mass loading of hydrometeors to decrease buoyancy in the stronger intensity thundercloud case [26,43,44]. In addition, in the weak thundercloud, the peak of updraft velocity arises in advance when CCN concentration is increasing (Figure 2a). With more latent heat released by the process of cloud droplets condensation, the maximum updraft speed increases from low CCN concentration to high CCN concentration, and similar characterization also can be found in the maximum downdraft speed. However, in moderate and severe cases, the updraft velocities in the case of N 0 = 1000 cm −3 is stronger than that in the other cases. As the CCN concentration further increases to 3000 cm −3 , the maximum updraft velocity keeps stable. This unusual instance only occurs in the condition of high CCN concentration and strong convolution, which is primarily because of the extreme conditions that are closely related to less latent heat released by different microphysical processes [14].
Atmosphere 2019, 10, x FOR PEER REVIEW 7 of 21 enhancement to slow reduction. The maximum updraft velocity increases greatly when the intensity of thunderclouds increases from weak to severe (see Figure 2), primarily because of the dramatic convection activities. Similarly, the downdraft speed increases with thundercloud intensity increasing because of a larger mass loading of hydrometeors to decrease buoyancy in the stronger intensity thundercloud case [26,43,44]. In addition, in the weak thundercloud, the peak of updraft velocity arises in advance when CCN concentration is increasing (Figure 2a). With more latent heat released by the process of cloud droplets condensation, the maximum updraft speed increases from low CCN concentration to high CCN concentration, and similar characterization also can be found in the maximum downdraft speed. However, in moderate and severe cases, the updraft velocities in the case of N0 = 1000 cm −3 is stronger than that in the other cases. As the CCN concentration further increases to 3000 cm −3 , the maximum updraft velocity keeps stable. This unusual instance only occurs in the condition of high CCN concentration and strong convolution, which is primarily because of the extreme conditions that are closely related to less latent heat released by different microphysical processes [14].  Figure 3 presents the variation of cloud droplet number concentration and mixing ratio with time in weak, moderate, and severe thunderclouds, respectively. With the increase of CCN concentration, the cloud droplet mixing ratio and number concentration both increase. Especially in the cases of 3000 cm −3 (see Figure 4a-c), the cloud droplets are produced in a large amount at the lower layer. The increase in cloud droplets number concentration is much stronger than the enhancement of the mixing ratio. For example, in the moderate thundercloud, when N0 = 100 cm −3 , the average mixing ratio of cloud droplet is 0.56 g kg −1 , the mean number concentration only reaches 4.75 × 10 6 kg −1 , but when N0 = 3000 cm −3 , the average mixing ratio of cloud droplet is 0.97 g kg −1 , and the mean number concentration is 9.83 × 10 7 kg −1 . In other words, a higher concentration of CCN will lead to a smaller scale of cloud droplets. This is consistent with previous research [21,23,[45][46][47][48]. When the thundercloud development becomes more intense, the cloud droplet begins to decrease and disappears significantly in advance. In a weak thundercloud (Figure 3a,d,g,j), cloud top height is low and the cloud dissipates after 70 min, while in a severe thundercloud (Figure 3c,f,i,l and Figure 4c,f), the cloud top height is about 8 km and the cloud droplets dissipate after 45 min. The primary cause is that the stronger convective can lift the small cloud droplets to the colder region, promoting the conversion of cloud droplets to ice phase particles. Therefore, cloud droplets will be consumed during this process.  The increase in cloud droplets number concentration is much stronger than the enhancement of the mixing ratio. For example, in the moderate thundercloud, when N 0 = 100 cm −3 , the average mixing ratio of cloud droplet is 0.56 g kg −1 , the mean number concentration only reaches 4.75 × 10 6 kg −1 , but when N 0 = 3000 cm −3 , the average mixing ratio of cloud droplet is 0.97 g kg −1 , and the mean number concentration is 9.83 × 10 7 kg −1 . In other words, a higher concentration of CCN will lead to a smaller scale of cloud droplets. This is consistent with previous research [21,23,[45][46][47][48]. When the thundercloud development becomes more intense, the cloud droplet begins to decrease and disappears significantly in advance. In a weak thundercloud (Figure 3a,d,g,j), cloud top height is low and the cloud dissipates after 70 min, while in a severe thundercloud (Figure 3c,f,i,l and Figure 4c,f), the cloud top height is about 8 km and the cloud droplets dissipate after 45 min. The primary cause is that the stronger convective can lift the small cloud droplets to the colder region, promoting the conversion of cloud droplets to ice phase particles. Therefore, cloud droplets will be consumed during this process.
It can be found from Figure 5 that the thundercloud with strong violent activity promotes the formation of raindrops significantly, and the precipitation formation also arises in advance. In the case of 100 cm −3 , the formation time of raindrops advances from 24 min in the weak thundercloud to 18 min in the severe thundercloud. This is explained by a stronger convection intensity that leads to a quicker warm rain production by auto-conversion of cloud droplet-rain and the collection of cloud droplets by raindrops [47,[49][50][51][52]. Additionally, a higher concentration of CCN reduces raindrop number concentration and mixing ratio, and stronger thundercloud leads to a more obvious reduction (see Figure 5g,j-l). In high CCN cases, raindrops are mainly distributed below the melting layer. Most of the raindrops are produced in the mature phase of the thundercloud, which is caused by the melting of ice particles. The high concentration of CCN reduces the collision efficiency of cloud droplets, the automatic conversion rate of cloud to rain decreases, and reduces warm rain coalescence. Therefore, most of the precipitation comes from the later melting of ice particles. It also can be seen from Figure 6 that the production of the raindrop is reduced with increasing CCN, and especially the raindrops residing in the high region are hindered by the inefficiency conversion of cloud droplets to raindrops. Many studies have shown that the decreased precipitation of warm clouds with increasing CCN support this idea [22,47,48].    . Spatial and temporal distribution of the mixing ratio and the number concentration of cloud droplet for initial aerosol concentrations: (a-c) 100 cm −3 , (d-f) 500 cm −3 , (g-i) 1000 cm −3 , (j-l) 3000 cm −3 ; the number represents the average mixing ratio (Qc Ave) and the average number concentration (Hc Ave) , and below is the same; the red solid line represents the isotherm (−40 and 0 °C); the black contour respectively represents: 10 6 , 10 7 , 10 8 , and 10 9 kg −1 .  It can be found from Figure 5 that the thundercloud with strong violent activity promotes the formation of raindrops significantly, and the precipitation formation also arises in advance. In the case of 100 cm −3 , the formation time of raindrops advances from 24 min in the weak thundercloud to 18 min in the severe thundercloud. This is explained by a stronger convection intensity that leads to a quicker warm rain production by auto-conversion of cloud droplet-rain and the collection of cloud droplets by raindrops [47,[49][50][51][52]. Additionally, a higher concentration of CCN reduces raindrop number concentration and mixing ratio, and stronger thundercloud leads to a more obvious reduction (see Figure 5g,j-l). In high CCN cases, raindrops are mainly distributed below the melting layer. Most of the raindrops are produced in the mature phase of the thundercloud, which is caused by the melting of ice particles. The high concentration of CCN reduces the collision efficiency of cloud droplets, the automatic conversion rate of cloud to rain decreases, and reduces warm rain coalescence. Therefore, most of the precipitation comes from the later melting of ice particles. It also can be seen from Figure 6 that the production of the raindrop is reduced with increasing CCN, and especially the raindrops residing in the high region are hindered by the inefficiency conversion of cloud droplets to raindrops. Many studies have shown that the decreased precipitation of warm clouds with increasing CCN support this idea [22,47,48].  As illustrated in Figure 7, the thundercloud in three types with higher CCN concentration produced a larger ice crystals content, which increases from the low-intensity case to the strongintensity case. From Figure 8c,f it can be seen that the convective intensity increases the number concentration of ice crystals, and the high concentration center rises from 6 km in the weak thundercloud to 8 km in the severe thundercloud. Some ice crystals are initiated from vapor deposition onto ice nuclei, and some ice crystals are produced during graupel riming collection of cloud droplets with a diameter greater than 24 μm. The homogenous droplet freezing is also a major contributor to ice crystal production. The high CCN concentration case is associated with more cloud water content lifted by stronger updrafts, and thus it is favorable for ice nucleation, droplet freezing, and the collision between ice crystals and cloud droplets. Therefore, the content of ice crystals increases significantly with increasing CCN concentration. In addition, the reduction of ice crystals' size in high CCN cases is primarily associated with the process that a host of small cloud droplets can be frozen to ice crystal, thus significant cloud water competition hiders the growth of ice crystals. It should be borne in mind that the response of ice crystals' microphysical characteristics to CCN concentrations is similar in the different intensity thunderclouds, as found Khain et al. [53] and Yang et al. [54]. As illustrated in Figure 7, the thundercloud in three types with higher CCN concentration produced a larger ice crystals content, which increases from the low-intensity case to the strong-intensity case. From Figure 8c,f it can be seen that the convective intensity increases the number concentration of ice crystals, and the high concentration center rises from 6 km in the weak thundercloud to 8 km in the severe thundercloud. Some ice crystals are initiated from vapor deposition onto ice nuclei, and some ice crystals are produced during graupel riming collection of cloud droplets with a diameter greater than 24 µm. The homogenous droplet freezing is also a major contributor to ice crystal production. The high CCN concentration case is associated with more cloud water content lifted by stronger updrafts, and thus it is favorable for ice nucleation, droplet freezing, and the collision between ice crystals and cloud droplets. Therefore, the content of ice crystals increases significantly with increasing CCN concentration. In addition, the reduction of ice crystals' size in high CCN cases is primarily associated with the process that a host of small cloud droplets can be frozen to ice crystal, thus significant cloud water competition hiders the growth of ice crystals. It should be borne in mind that the response of ice crystals' microphysical characteristics to CCN concentrations is similar in the different intensity thunderclouds, as found Khain et al. [53] and Yang et al. [54].
Since the liquid water lifted by strong updraft led to the rapid growth of ice crystals, graupel is firstly produced by the auto-conversion of ice-graupel. Thereafter, the collisions between the liquid droplets and ice particles play a signification role in the development of graupel growth [21,22]. It can be seen from Figure 9 that the graupel particles are produced in advance with increasing thundercloud intensity. For example, the occurrence time of graupel in the cases (the CCN concentration is 100 cm −3 ) is from 48 min in weak thundercloud to 18 min in severe thundercloud (Figure 9a-c). When the CCN concentration increases, the production of graupel is strong at first and then tends to become weak. For example, the peak content of graupel arises at about 42 min in the case of 1000 cm −3 (weak thunderclouds) (Figure 9g,j), and the vertical distribution of maximum mixing ratio and number concentration is similar to those in the case of 1000 and 3000 cm −3 (Figure 10a,d). While in moderate thunderclouds, the average concentration (1.89 × 10 4 kg −1 ) and the average mixing ratio (1.93 g kg −1 ) of graupel reach a peak in the case of 500 cm −3 (Figure 9e) which is similar to the case with CCN concentration of 1000 cm −3 (Figure 9h). From Figure 10b,e, it also can be found that the production of graupel in the case of 500 and 1000 cm −3 is significantly larger than that in the other cases. In severe thunderclouds, the content of graupel reaches a peak in the case of 500 cm −3 (Figure 9f) and it is obvious that the residence area of graupel in the severe case of 500 cm −3 is much larger than that in other cases. Overall, deeper convection and less CCN concentration can accelerate the development of graupel, but an extremely high concentration of CCN would result in weakening the formation of graupel. In the initial state, the CCN particles are fully activated under the sufficient water vapor condition, and a large number of ice crystals are formed through vapor deposition and condensation freezing. This promotes the auto-conversion of ice crystal to graupel and the collision coalescence between the liquid droplets and ice particles, leading to more graupel produced. However, when the concentration of CCN reaches a certain extent, due to the consumption of water vapor by precipitation, the water vapor supply is reduced, which drives the competition between ice crystal and graupel to available water vapor. Therefore, a high concentration of CCN reduces the size of ice crystals, hindering the conversion of ice crystals to graupel. In addition, the variation of graupel is very similar to that of the maximum updraft. The development is the most vigorous among three cases, where the CCN concentration of weak thundercloud is 3000 cm −3 , the CCN concentration of moderate thundercloud and severe thundercloud is 500-1000 cm −3 . Therefore, the process that leads to the decrease of maximum updraft is likely to be the microphysical development of graupel. The formation of graupel is reduced by high concentration CCN and deep convection so that the latent heat released decreases and the maximum updraft becomes weaker.  Since the liquid water lifted by strong updraft led to the rapid growth of ice crystals, graupel is firstly produced by the auto-conversion of ice-graupel. Thereafter, the collisions between the liquid droplets and ice particles play a signification role in the development of graupel growth [21,22Error! Reference source not found.]. It can be seen from Figure 9 that the graupel particles are produced in advance with increasing thundercloud intensity. For example, the occurrence time of graupel in the cases (the CCN concentration is 100 cm −3 ) is from 48 min in weak thundercloud to 18 min in severe thundercloud (Figure 9a-c). When the CCN concentration increases, the production of graupel is strong at first and then tends to become weak. For example, the peak content of graupel arises at about 42 min in the case of 1000 cm −3 (weak thunderclouds) (Figure 9g,j), and the vertical distribution of maximum mixing ratio and number concentration is similar to those in the case of 1000 and 3000 cm −3 (Figure 10a,d). While in moderate thunderclouds, the average concentration (1.89 × 10 4 kg −1 ) and the average mixing ratio (1.93 g kg −1 ) of graupel reach a peak in the case of 500 cm −3 (Figure 9e) which is similar to the case with CCN concentration of 1000 cm −3 (Figure 9h). From Figure 10b,e, it also can be found that the production of graupel in the case of 500 and 1000 cm −3 is significantly larger than that in the other cases. In severe thunderclouds, the content of graupel reaches a peak in the case of 500 cm −3 (Figure 9f) and it is obvious that the residence area of graupel in the severe case of 500 cm −3 is much larger than that in other cases. Overall, deeper convection and less CCN concentration can accelerate the development of graupel, but an extremely high concentration of CCN would result in weakening the formation of graupel. In the initial state, the CCN particles are fully activated under the sufficient water vapor condition, and a large number of ice crystals are formed through vapor deposition and condensation freezing. This promotes the auto-conversion of ice crystal to graupel and the collision coalescence between the liquid droplets and ice particles, leading to more graupel produced. However, when the concentration of CCN reaches a certain extent, due to the consumption of water vapor by precipitation, the water vapor supply is reduced, which drives the competition between ice crystal and graupel to available water vapor. Therefore, a high concentration of CCN reduces the size of ice crystals, hindering the conversion of ice crystals to graupel. In addition, the variation of graupel is very similar to that of the maximum updraft. The development is the most vigorous among three cases, where the CCN concentration of weak thundercloud is 3000 cm −3 , the CCN concentration of moderate thundercloud and severe thundercloud is 500-1000 cm −3 . Therefore, the process that leads to the decrease of maximum updraft is likely to be the microphysical development of graupel. The formation of graupel is reduced by high concentration CCN and deep convection so that the latent heat released decreases and the maximum updraft becomes weaker.   Figure 9. Spatial and temporal distribution of the mixing ratio and the number concentration of graupel for initial aerosol concentrations: (a-c) 100 cm −3 , (d-f) 500 cm −3 , (g-i) 1000 cm −3 , (j-l) 3000 cm −3 ; the red solid line represents the isotherm (−40 and 0 °C); the black contour respectively represents: 10 3 , 10 4 , 10 5 , and 10 6 kg −1 .

Charging Rate
The electrification process of the thundercloud is based on the non-inductive electrification mechanism of ice crystals and graupel particles and the inductive charging between graupel and cloud droplets [42]. Therefore, the concentration of graupel is one of the key factors affecting the electrification process [55]. From the analysis in the previous section, we know that the cloud droplets, ice crystals, and graupel are sensitive to CCN concentration, so the change of those hydrometer particles should lead to a different electrification process with different CCN concentrations. The time evolution of the non-inductive charging rate is given in Figure 11. The negative non-inductive charging rate of all cases roughly resides between 3 and 5 km (0 to −15 • C), while the stronger production of ice particles makes a greater contribution to actively positive non-inductive charging. Non-inductive charge separation in the weak-case starts at 40 min, while in the severe-case charge separation occurs at an earlier time of 18 min, implying the great electrification is closely connected with heavy convection. In addition, the non-inductive charging rate demonstrates a complex trend with CCN concentration increasing. In four cases of the weak thundercloud, higher CCN concentration makes the vertical distribution of non-induced electricity wider, while in the moderate thundercloud and severe thundercloud, the time evolution of the non-inductive charging rate in the case of 3000 cm −3 shows the complex trend, which is roughly similar to that in the 1000 cm −3 case. The complicated evolution of the charging process is because of the distribution of mixed-phase hydrometer particles produced from strong convection.
To further show the characteristics of induction and non-induction charging rates, the induction and non-induction charging rates of the 12 cases are listed in Table 1. The increase non-inductive charging rate with CCN concentration increases is well correlated to the stronger convection. The charge separation parameterization depends on graupel-ice collisions in mixed-phase updrafts [22]. Therefore, the change in charging rate means that stronger production of small ice crystal and graupel related to CCN concentrations ranging from 100 to 1000 cm −3 and deep convection can enhance non-induction electrification and lightning processes. Note that under extreme conditions (3000 cm −3 cases of severe thunderclouds), the positive non-inductive charging keeps steady, but the negative non-inductive charging rate is reduced from −1724.95 to −1289.69 pC m −2 s −1 . From the analysis in the previous section, it can be found that the production of graupel particles at higher CCN concentration (1000−3000 cm −3 ) and deep convection is reduced by vapor competition, but more ice crystals are still able to improve the collision efficiency, thus the non-induction charging process is gradually enhanced. When the CCN concentration is over 3000 cm −3 , the decrease in ice crystals' size and the small supply of graupel is likely to hold the answer to the decrease of the non-induction charging rate. This also shows that the amount of charge separated is associated with ice crystals' concentration and size spectrum, which means stronger particle formation led to a higher non-induction charging rate [22,33,56].
while the stronger production of ice particles makes a greater contribution to actively positive noninductive charging. Non-inductive charge separation in the weak-case starts at 40 min, while in the severe-case charge separation occurs at an earlier time of 18 min, implying the great electrification is closely connected with heavy convection. In addition, the non-inductive charging rate demonstrates a complex trend with CCN concentration increasing. In four cases of the weak thundercloud, higher CCN concentration makes the vertical distribution of non-induced electricity wider, while in the moderate thundercloud and severe thundercloud, the time evolution of the non-inductive charging rate in the case of 3000 cm −3 shows the complex trend, which is roughly similar to that in the 1000 cm −3 case. The complicated evolution of the charging process is because of the distribution of mixedphase hydrometer particles produced from strong convection. To further show the characteristics of induction and non-induction charging rates, the induction and non-induction charging rates of the 12 cases are listed in Table 1. The increase non-inductive charging rate with CCN concentration increases is well correlated to the stronger convection. The charge separation parameterization depends on graupel-ice collisions in mixed-phase updrafts [22].  The collision of the non-induction electrification process between ice particles enhances the environmental electric field, and induction electrification begins to occur with strong electric fields. Figure 12 shows the time evolution of inductive charging rates by graupel. The inductive charging rate is smaller than the non-induction charging rate, and greater production of cloud droplet and graupel below the attitude of 6 km leads to stronger inductive charge separation in higher CCN cases. The maximum positive and negative inductive charging rates in Table 1 depends on the amount of charge acquired from graupel. From the table and Figure 12, one can conclude that the charge separation reduces sharply from 100 to 1000 cm −3 and then increases in the moderate and severe thundercloud. As we all know, the charge transfer processes caused by the collision of graupel and cloud droplets can be significantly influenced by the non-induction electric field. With CCN concentration increasing from 100 to 500 cm −3 (four cases in moderate and severe thundercloud), the production of cloud droplet and graupel is increased, and the enhanced electric field is explained by the stronger non-induction electrification process, which causes the steady increases of inductive charging rates. However, the supply of graupel is restricted at higher CCN concentration (1000 cm −3 ), resulting in the weak inductive charge separation. Figure 12g,k shows a smaller distribution of positive inductive charging rate. The studies on the electrification mechanism point out that the inductive charging process is weak, and the non-induction charging process plays a major role in the whole electrification process [22,37]. Therefore, when the concentration of CCN increases to 3000 cm −3 , the influence of the vertical electric field from non-inductive charge enhancement is dominant, implying that the inductive charge separation process exhibits an obvious enhancement.

Electrification Charge Structure
The evolution of charge structure that depends only on charge separation at each minute is shown in Figure 13. In the very early electrical development stage, all storms depict a dipole charge structure (positive charge above negative), which maintains a positive charge at higher altitudes (5-8 km), with negative charge raised by enhancement of convection at 4-5 km. However, with the increase of CCN concentration, the enhancement of convective intensity and the process of lightning, the charge structure becomes more and more complex. The weak thundercloud with no lightning and low CCN (100 cm −3 ) shows a stable normal dipole structure ( Figure 13a). As shown in Figure 2, the maximum updraft speed in the weak thundercloud is less than 10 m s −1 , so the cloud height is limited, and the positive charge region is corresponding to ice crystals which are charged positively during the non-inductive electrification process. The negative charge region of a weak storm is composed of graupel that is negatively charged above the anti-rotation temperature layer and ice crystal which is negatively charged under the anti-rotation temperature layer. However, other cases depict a tripole structure with a relatively weak lower positive charge before first lightning. The bottom secondary positive charge region is coincident with the non-inductive graupel carrying positive charge below the inversion temperature layer. In these cases, the secondary positive charge

Electrification Charge Structure
The evolution of charge structure that depends only on charge separation at each minute is shown in Figure 13. In the very early electrical development stage, all storms depict a dipole charge structure (positive charge above negative), which maintains a positive charge at higher altitudes (5-8 km), with negative charge raised by enhancement of convection at 4-5 km. However, with the increase of CCN concentration, the enhancement of convective intensity and the process of lightning, the charge structure becomes more and more complex. The weak thundercloud with no lightning and low CCN (100 cm −3 ) shows a stable normal dipole structure ( Figure 13a). As shown in Figure 2, the maximum updraft speed in the weak thundercloud is less than 10 m s −1 , so the cloud height is limited, and the positive charge region is corresponding to ice crystals which are charged positively during the non-inductive electrification process. The negative charge region of a weak storm is composed of graupel that is negatively charged above the anti-rotation temperature layer and ice crystal which is negatively charged under the anti-rotation temperature layer. However, other cases depict a tripole structure with a relatively weak lower positive charge before first lightning. The bottom secondary positive charge region is coincident with the non-inductive graupel carrying positive charge below the inversion temperature layer. In these cases, the secondary positive charge region at the bottom of the thundercloud with low CCN concentration and shallow convective (Figure 13b-d,g) gradually disappears, and then the charge structure tends to be dipole, in which the main positive charge region is larger than the main negative charge region. It is such a difference that the positive charge can get freedom from the bondage of the negative charge region and reach the ground, thus causing triggering the positive CG lightning [52,57]. For higher CCN concentration and deeper convective (weak thundercloud of 3000 cm −3 and 500-3000 cm −3 cases in moderate thundercloud and severe thundercloud), with the development of thunderclouds turning to maturity, the development of multipolarity charge structure shows a characteristic as follows: (1) the main negative charge resides in the lower region where the temperature is higher. (2) In addition to the three middle charge regions formed by the non-induction mechanism, the secondary positive charge region contains some small negative charge regions, which can be found in weak thunderclouds at N0 = 3000 cm −3 and moderate or severe thunderclouds at N0 ≥ 500 cm −3 . (3) The increasing CCN concentration leads to more negative charge production on a small scale. For example, in moderate storms of 3000 cm −3 cases, charge structure with six layers appears in 40−50 min. This feature is highly similar to non-inductive charge separation and inductive charge separation. Our results are consistent with the studies that are conducted to simulate severe thunderclouds, as they found that the charge generation is dominated by non-inductive charge separation between graupel-ice crystals on the updraft region and inductive charging processes on the periphery of updraft region [57,58]. The non-inductive electrification was mainly produced in the For higher CCN concentration and deeper convective (weak thundercloud of 3000 cm −3 and 500-3000 cm −3 cases in moderate thundercloud and severe thundercloud), with the development of thunderclouds turning to maturity, the development of multipolarity charge structure shows a characteristic as follows: (1) the main negative charge resides in the lower region where the temperature is higher. (2) In addition to the three middle charge regions formed by the non-induction mechanism, the secondary positive charge region contains some small negative charge regions, which can be found in weak thunderclouds at N 0 = 3000 cm −3 and moderate or severe thunderclouds at N 0 ≥ 500 cm −3 . (3) The increasing CCN concentration leads to more negative charge production on a small scale. For example, in moderate storms of 3000 cm −3 cases, charge structure with six layers appears in 40−50 min. This feature is highly similar to non-inductive charge separation and inductive charge separation. Our results are consistent with the studies that are conducted to simulate severe thunderclouds, as they found that the charge generation is dominated by non-inductive charge separation between graupel-ice crystals on the updraft region and inductive charging processes on the periphery of updraft region [57,58]. The non-inductive electrification was mainly produced in the altitude of 3−8 km, and the negative non-inductive charging rate decreases with the development of storm. At the highest CCN concentration (3000 cm −3 ), the non-inductive process remains stable, while the inductive charging rate tends to increase, so it is an obvious contribution to the local area showing a small range of negative charge region. In addition, compared with the variation of graupel in Figure 9, it can be seen that the content of graupel above the inversion temperature layer in weak thundercloud rises slightly with the increase of CCN concentration. However, in moderate and severe thunderclouds, high CCN concentrations result in a decreasing range of graupel center. For example, in severe thundercloud of the 500 cm −3 case (Figure 9f), the graupel center is located from 6 to 8km, but the center of graupel in the 3000 cm −3 case (Figure 9h) was reduced to 4−6 km, and the distribution was greatly reduced, which is quite consistent with the evolution of the main negative charge region.

Conclusions
The influences of aerosol concentration on dynamics, microphysics, and electrification of different intensity thunderclouds were evaluated by using a two-dimensional thundercloud model. Due to the significant difference in the production and development of graupel caused by different CCN concentrations in different intensity thunderclouds, the response of dynamics, electrification, and charge structures is a nonlinear relationship. Our results demonstrate that the pronounced increase in CCN concentration leads to a significant enhancement of updraft in the weak thunderclouds. At high CCN concentration in the moderate and severe thundercloud cases, however, it can be found that a slight reduction in maximum updraft is closely connected with the microphysics of graupel. In addition, the most noticeable effect of increasing aerosol in three types of the thundercloud is the increase of cloud droplet number concentration and mixing ratio, but a decrease in the size of cloud droplets. The increase of the convection promotes the lift of small cloud droplets, which leads to a faster and stronger production of ice crystals. Additionally, small cloud droplets reduce the auto-conversion rate of cloud droplets and raindrops so that the formation of raindrops is suppressed. Since the raindrops are produced from the cloud-rain auto-conversion to the melting of the ice-phase particles, the formation of the raindrops is delayed, and raindrop production decreases very rapidly in the case of high-intensity thunderclouds. The production of graupel is insensitive to the CCN concentration. The content of graupel increases from low CCN concentration to moderate CCN concentration and slightly decreases at high CCN concentration. The production graupel in high CCN concentration cases is hindered by the formation of large amounts of small ice crystals, which enhance the competition of available water vapor. When the intensity of thundercloud increases, the reduction of graupel production will arise in advance as the CCN concentration increasing.
Charge separation tends to increase as the aerosol concentration rises from low to high in the weak and moderate thundercloud cases. Although a slight reduction of graupel production in the high CCN concentration case was observed, the enhancement of ice crystal and cloud droplet formation still arises in the increase of non-inductive and inductive charge separation rates. However, the magnitude of charging rates in the severe thundercloud cases remain roughly stable under high CCN concentration condition, which can be attributed to the profound reduction of graupel content. The charge structure in the weak thundercloud at low CCN concentrations (100 cm −3 ) remains as a dipole, while the weak thunderclouds in the other cases (the CCN concentration above 100 cm −3 ) experience change from a dipole charge structure to a tripole charge structure, and finally disappear with a dipole. The time evolution of charge structure under the low CCN concentration condition in the moderate and severe thunderclouds is roughly similar to that in the weak thunderclouds with high CCN condition. In cases of strong intensity thunderclouds, the charge structure depicts a relatively complex structure that includes a multilayer charge region with alternating polarity. Those small-range and short-lived weak negative charge regions under the lower positive charge region are formed by graupel charging during the inductive charge separation process.
It is believed that aerosol can act as ice nuclei (IN). As cloud-ice nuclei (IN) interaction is increasingly recognized as one of the factors influencing the microphysical structure of clouds, and ice crystal (nucleation of IN) is of crucial importance in thundercloud electrification, this process may have an important impact on electrification properties in thunderclouds [59][60][61]. The coalescence between water drops is enhanced due to the inclusion of giant CCN, resulting in the early development of raindrops. Therefore, the giant CCN perhaps has a profound influence on microphysics and electrification in thunderstorms. Further aspects of this problem will be addressed in forthcoming studies. Due to the limitations of the 2D model, the frequency of lightning obtained from the 2D cloud model is relatively small as compared to some observations. The influence of aerosol on the lightning discharges process in thunderclouds will be presented by a 3D model in the future, but we still believe that the inherent physical mechanisms and lightning discharge characteristics in the two model categories are similar.