Charge Structure of a Summer Thunderstorm in North China: Simulation Using a Regional Atmospheric Model System

Investigations into the charge structure and discharge processes inside thunderstorms have significantly improved in recent years, due to improvements in observational technologies. Observations by jets and balloons in thunderstorms can only reveal the properties along their routes; the direct detection of the charge structure remains difficult to achieve. Hence, numerical simulations play an important role in investigating the electrification process of thunderstorms. In the past few decades, numerical models with a parameterized electrification scheme have provided increasingly accurate simulations of thunderstorm charge structure (Takahashi, 1984; Ziegler et al., 1991; Helsdon et al., 2001; Mansell et al., 2005; Fierro et al., 2008; Barthe et al., 2012; Li et al., 2012).

The models used to simulate charge structure are mostly based on a cloud-scale model with electrification parameterization (Norville et al., 1991; Ziegler et al., 1991; Zhang et al., 1999; Helsdon et al., 2001; Barthe et al., 2005; Barthe and Pinty, 2007), which is suitable for small-scale thunderstorms. However, for mesoscale convective systems (MCSs) and squall lines, which frequently occur in summer, simulations of large-scale thunderstorms are also necessary. Using a 3D model, the electrification and lightning activity of thunderstorms (Mansell et al., 2002, 2005), an idealized hurricane (Fierro et al., 2007), and a squall line (Fierro et al., 2008) have been analyzed. Based on a 3D mesoscale model, RAMS (Regional Atmosphere Model System), (Altaratz et al., 2005) simulated the charge structure of a winter thunderstorm in Israel, considering a non-inductive electrification scheme without a discharge parameterization, and the simulation was ended before the first occurrence of lightning. The results indicated a tripole charge structure with Takahashi (Takahashi, 1978) non-inductive charging parameterization and a dipole charge structure with the Saunders scheme (Saunders et al., 1991). (Barthe et al., 2012) coupled electrification and lightning parameterization into a non-hydrostatic Meso-NH model to simulate the charge structure of a storm over large computational domains. (Xu et al., 2012) utilized the Weather Research and Forecasting (WRF) model electrification and lightning parameterization and simulated an idealized supercell, which presented a tripole charge structure.

In this study, the RAMS model coupled with a charging mechanism is used to simulate the charge structure of a summer thunderstorm. The simulated electrical activity is then compared with observed lightning characteristics. (Altaratz et al., 2005) examined the two non-inductive charging mechanisms and obtained the charge structure of a thunderstorm. However, the simulation was ended before the first occurrence of lightning. This study employs a RAMS-coupled electrification and discharge model and discusses the relationship between the simulated charge structure and the measured lightning features. The relationships between the dynamics, microphysics and electrification process are explored, and the simulation results are compared with the observed lightning activity of the thunderstorm. In addition, the study combines the simulated results and observed data to provide a meaningful electrification and discharge simulation. Although the case simulated in this study is an isolated storm, the overall objective is to develop an electrical model that can be used for large-scale convective systems in large grids.

The model used for the simulation is based on a mesoscale model of RAMS (Tripoli and Cotton, 1982; Tremback, 1990), which was developed by Colorado State University. It is capable of simulating thunderstorms efficiently and effectively. More details of the dynamics and microphysics of the model are described in the works of (Pielke et al., 1992) and (Cotton et al., 2003).

The microphysical parameterization in the model is based on a bulk two-moment microphysical scheme, whereby the hydrometeor particles are categorized as cloud, rain, graupel, ice, snow, hail and aggregate. The microphysical processes parameterized in the model include condensation, deposition, evaporation, sublimation, nucleation of ice crystals, melting, collection, freezing, secondary ice production, auto-conversion, collision and coalescence of hydrometeor particles, breakup, and sedimentation.

A warm moist bubble is introduced as an initiation mechanism into the model, which is a warm, moist disturbed field (relative to the local environment), located within the low-level range of the simulation region. The initial convection is governed by the vertical equations of motion, as follows:

Here, (x _c,y _c,z _c) are the x, y, and z values of the center coordinates, respectively; (x _r,y _r,z _r) are the three-dimensional radii; θ stands for potential temperature; ∆θ is the maximum disturbed potential temperature of the center; Q _v stands for the mixing ratio of water vapor; and Q _vs is the saturation specific humidity; Q _v0 is the initial environmental field; β is the arbitrary factor. Considering the characteristics of the influence of ice on radiation, the Harrington radiation scheme is adopted. Further details can be found in (Altaratz et al., 2005).

2.2.1 Inductive electrification parameterization

According to (Chiu, 1978), the charge transfer that occurs between small particles with charge q _se and radius r _s, and larger particles with charge q _Le and radius r _L are described as follows:

Here, ∆ q=4πεΓ₁|E|cos(E,r _Ls)r _s²+Aq _Le-Bq _se, where E stands for the electrical field, ε is the permittivity of air, Γ₁ is the complete gamma function, A and B are dimensionless coefficients:

Γ₂ is dependent on the radius ratio between the collision particles in the formula. When the value of r _s/r _L is between 0.01 and 1.0, the ranges of Γ₁ and Γ₂ are π²/6≤Γ₁≤π²/2 and 1.0≤Γ₁≤π²/6, respectively. When r _s« r _L, then Γ₁=π²/2 and Γ₂=π²/6.

Large particles with radius R _L and charge q _Le collide with smaller particles. The formula describing the charge separation rate of large particles is given below:

Here, E _Ls is the collision efficiency, V _Ls is the fall velocity, S(θ) is the separation probability function, θ is the linear angle between the centers of the small and large particles and the z-axis, N _s stands for the number concentration for small particles and A is the collision cross-section area. The integration of Eq. (7) can be written as:

where cos(E,V _Ls)=EV _Ls/(|E||V _Ls|) and cosθ are the average collision angles.

2.2.2 Non-inductive electrification mechanism

The non-inductive electrification mechanism is considered as the most important mechanism involved in the charging of a thunderstorm. When choosing different non-inductive charging parameterizations mentioned in the literature (Takahashi, 1978; Altaratz et al., 2005), the adopted charge transfer function is accordingly different.

The charge transfer of graupel particles on liquid water content and temperature is adopted from (Takahashi, 1978). The temperature ranges from 0^°C to -30^°C, and cloud water content (CWC) values range from 0.01 to 30 g m^-3. For temperatures lower than -30^°C, the charge separation values for -30^°C are used. (Takahashi, 1984) introduced a revised factor of per collision charge transfer based on laboratory experiments, as follows:

where D _i is the diameter of the ice crystal, and V _g and V _i stand for the terminal falling velocities of graupel particles and ice crystals (m s^-1), respectively. The magnitude of D₀ is 100 μm and V₀ is defined as 8 m s^-1. The value of α is limited to 10, and stands for the apparent saturation value of large snow crystals in electrical charge (Marshall et al., 1978). The non-inductive electrification mechanisms of graupel/hail-ice and graupel/hail-snow are considered in the model.

To simulate the entire electrification process of the storm, the lightning parameterization is taken into account. However, the discharge process of a thunderstorm is highly complicated and includes various electrical processes. It is difficult to reflect the real lightning channel in simulations with lightning parameterization. Usually, the discharge parameterization can be separated into two categories. One category is the whole discharge parameterization, and the other category is the lightning scheme, which includes the propagated branch channel. (Takahashi, 1984) defined the breakdown threshold as 340 kV m^-1 and assumed that the same positive and negative charge is neutralized in each hypothetical discharge. In the second category, the development of lightning channels and the bidirectional extension of discharge following a step-by-step process are assumed (Mansell et al., 2005; Tan et al., 2006).

The simple lightning parameterization used in the present study's model is in accordance with (MacGorman et al., 2001), wherein lightning occurs when the magnitude of the electric field exceeds the breakdown value:

where E _init stands for the initial electric field, z is the height in km, and ρ _a is the air density in kg m^-3. Lightning parameterization is initiated when the electric fields (E) at the grid point satisfy the condition of E>0.9E _init. The range of discharge is influenced by a small initial value and random effects.

A sketch of the simple discharge scheme is shown in Fig. 1. If the grid points satisfy the condition that the electric field exceeds the breakdown value, the discharge scheme will work. The initial point of lightning is arbitrarily chosen from those meeting the condition that more than three grid points satisfy a breakdown of the electric field around the point at greater than 0.9E _init(E>0.9E _init). The simulated lightning frequency in the RAMS model obeys the rule that a single lightning flash should be identified in at least seven or more adjacent grid points to satisfy the condition that the corresponding electric field reaches the breakdown value. The discharge mechanism neutralizes the positive and negative charges previously stored in the discharge channels, resulting in a charge redistribution of the thunderstorm. In order to reflect this discharge function in the model, it is assumed that the concentration of hydrometeor particles, and the charge density of positive and negative charge in the associated grid points decrease by a certain percentage of 40% after the discharge ends, resulting in a reduction in the net charge concentration. Consequently, the net charge density decreases and is redistributed by a gradual change process. When the discharge is finished, the SOR-type (successive over relaxation) method is used for recalculating the electric potential with one or multiple iterations.

Figure 1.  Sketch of a simple grid showing the discharge scheme. Open circles represent grids not taking part in discharge. Solid black circles stand for the initial point satisfying (E>0.9E _int), and the solid red circles stand for the adjacent seven consecutive grid points. Solid lines indicate the grids participating in the discharge scheme as single lightning bolts without a discharge path.

A 3D simulation region with dimensions of 120×120×14 km³ is represented in the model, which has a horizontal resolution of 500 m, a vertical resolution of 100 m, and a vertical stretching ratio of 1.11. At a low-level center of the simulation region, the warm moist bubble is initialized with a temperature of 4^°C, a humidity ratio of 1.3 g kg^-1, a horizontal dimension of 4 km, and a vertical extent of 2 km.

The sounding data at 1200 UTC 15 July 2007 after the thunderstorm dissipated were obtained from the Beijing Meteorological Station. Thus, the sounding data are modified in the code to reflect the real-time situation to make the simulated spurious convection closer to the observed reality. The revised sounding applied to the initial profiles is shown in Fig. 2, indicating a ground temperature of nearly 30^°C, a dew point of 22^°C, a convective available potential energy (CAPE) of 4800 J kg^-1, and a K index of 34.2^°C.

Figure 2.  Modified sounding data according to the initial one at the 1200 UTC 15 July 2007 at Beijing meteorological station during the storm. Solid line stands for the modified temperature; dashed line for the modified dew point.

The storm forms approximately 10 minutes after the simulation begins, when the bubble is elevated to the lifting condensation layer. Figure 3 displays the observed and simulated radar CAPPI (Constant Altitude Plan Position Indicator) of the thunderstorm with a height of 3.5 km, which moves slowly from the northwest to southeast, approaching the quasi-stationary state with high radar reflectivity in front of the storm. The position of the convective leading line is slightly different between the simulation and the observation, whereas the profile of the thunderstorm is roughly consistent between both datasets. Overall, the simulated radar echo is smoother than the observed radar echo. Furthermore, the trailing rearward region of the simulated storm is larger than that of the observed situation.

Figure 3.  Simulated radar reflectivity and observed radar reflectivity at the 3.5-km height level at 45 min: (a) simulation; (b) observation.

The evolution of the storm observed by the Doppler radar in different stages is shown in Fig. 4. The thunderstorm develops vertically at the beginning, with the radar reflectivity greater than the simulated reflectivity. In the developing stage, the top of the thunderstorm extends to 10 km MSL, and the horizontal range gradually enlarges.

Figure 4.  Cross section of the radar echo as the stom develops: (a)10min;(b)45min;(c)60min;(d)80 min

With the thunderstorm entering the mature stage, the strong radar echo reaches a maximum value of 60 dBZ. An anvil region forms in front of the storm, and a trailing rearward region forms at the rear, similar to in the simulation. In the dissipation stage, both the intensity of radar reflectivity and the range of the thunderstorm decrease.

According to the diameter and concentration of hydrometeor particles, the radar reflectivity is calculated. Figure 5 illustrates a cross section of the simulated radar reflectivity of the thunderstorm. A comparison between the simulation and the observation (Fig. 3) of radar echo suggests that the distribution and intensity of the radar reflectivity are largely consistent. Ten minutes after convection initiation (Fig. 5a), the maximum radar reflectivity reaches 55 dBZ, the horizontal extent of the thunderstorm (with radar echo of 5 dBZ) is approximately 20 km, and the top of the thunderstorm extends to 8.5 km MSL. Thirty minutes after the beginning of the simulation, the scale of the radar reflectivity greater than 40 dBZ gradually enlarges. When the updraft reaches the top of the thunderstorm, the airflow transported from the convective region moves downward to the trailing stratiform region. As the simulated thunderstorm enters the mature stage at 45 minutes, a mushroom-shaped storm forms, and the thunderstorm is characterized by strong central updraft. At the same time, the strong radar echo reaches its largest range, and extends to the top of the thunderstorm from the bottom layer. With the developing convection activity, the horizontal range of the thunderstorm gradually extends to between 60 km MSL and 10 km MSL in the vertical direction at the simulation time of 45 min. At this moment, the maximum radar reflectivity reaches 60 dBZ, and is mainly associated with the updraft zone. At 60 min, the convective activity enters into the most vigorous stage, gradually weakens, and then the decaying convective cell turns into the trailing rearward region. When the thunderstorm enters the dissipating stage, evaporative cooling is generated, driving a strong downdraft, and both the intensity of radar reflectivity and the spatial extent are reduced.

Figure 5.  Cross section (x = 30 km) of the radar echo in the thunderstorm, taken at four simulation time slices: (a) 10 min; (b) 45 min; (c) 60 min; (d) 80 min.

Once the convection processes initiates, both the concentration of hydrometeor particles and the updraft strengthens rapidly. With the maximum updraft increased, the center of the thunderstorm lifts to about 5-km MSL at 30 min. Cloud droplets (not shown here) are mainly located in the lower level of the storm where the temperature is -10^°C. Rain drops are automatically converted by cloud droplets and are mostly located where temperatures are above 0^°C. Snow and ice crystals have already formed and are mainly located in the upper level, where the temperature is -25^°C. Graupel particles are mainly generated by glaciation processes between the ice and cloud droplets, which are roughly centered at 5 km MSL, where the temperatures are -15^°C. At this stage, hail particles are located at the level of 2-6 km MSL formed by the freezing and riming mechanism in the thundercloud. The high radar reflectivity in the lower level is mainly formed by raindrops and graupel particles. Different hydrometeor categories are interconverted via melting, nucleation, evaporation, and other microphysical processes.

With the aggregation process commencing, both the riming process and the growth of ice crystals accelerate. Graupel and hail particles are mainly generated in the mixed-phase region of the thunderstorm. Due to the wind shear effect, the bulk of the mass is located in the upwind zone of the thunderstorm, whereas the riming mechanism is more efficient in the region with the strong updraft. The microphysical process of coalescence and the collision between graupel and hail particles with raindrops creates the core of hail particles. It should be noted that the concentration of aggregates is relatively less than other hydrometeors, and the influence of electrification can be ignored in the present case.

Figure 6.  Cross section (x = 30 km) of the charge density (units: nC m^-3) at 60 min simulation time: (a) cloud; (b) ice + rain; (c) graupel; (d) snow; (e) hail. (f) Legend for (a-e), in which CD stands for charge density of different hydrometeor particles. Solid black lines show the peripheral profile of the storm.

The simulated charge density carried by the different categories of hydrometeor particles is shown in Fig. 6. The charge separation process starts at 10 min, after the convection becomes active. As the cloud droplets are mainly distributed in the lower layer of the thunderstorm, the maximum positive charge region is located at a height of approximately 4 km and a negative charge is centered at 6 km MSL mainly carried by cloud drops. Ice crystals and rain gain positive charge with a maximum density of 2.2 nC m^-3 in the upper level and negative charge with a minimum density of -1.6 nC m^-3 in the lower level. The charge density of raindrops is mainly located in the bottom of the storm, with some of the hydrometeor particles melting and the charge transferring to rain. Graupel particles carrying the negative charge are located at the 8-km level, while the positive charge region is located at 4 km MSL. Snow crystals with positive charge dominate at approximately 8 km MSL (Fig. 6d), whereas negatively charged crystals predominate in the lower level at approximately 5 km MSL, with the light mass of snow particles being elevated to the upper level via updraft. Figure 6e shows the distribution of hail charge density with a positive charge in the lower levels of 2-4 km and a negative charge region in the 4-6-km levels. Some of the ice crystals are transferred via downdraft from the upper layer to the lower layer of the thunderstorm.

Via the non-inductive charging mechanism function, which plays an important role in the entire electrification simulation, the magnitude of the electric field increases gradually, and then the effect of the inductive electrification mechanism progressively strengthens. When the simulation enters the medium-term period, the function of the inductive mechanism strengthens, and the electrification activity becomes intense, partly due to the continuous increase of the inductive charge transfer rate.

Figure 7 displays the charge density distribution of different time slices under the electrification and discharge mechanisms. From the initial simulation to 10 min (Fig. 7a), with the charge separation having started, the graupel particles are mainly positively charged, and the ice crystals are negatively charged. Due to the different sedimentation rates between the graupel particles and the ice crystals, the thunderstorm is characterized by an inverted dipole charge structure, a low-layer positive charge region with a temperature of -10^°C and an upper negative charge region with a temperature of 0^°C. The maximum positive charge density reaches 0.18 nC m^-3, and the minimum negative charge density reaches -0.24 nC m^-3. However, no lightning occurs during this stage because of the small charge density. The non-inductive electrification mechanism is primarily influenced by the charge separation between the ice crystals and the graupel particles.

Figure 7.  Charge structure distribution of the thunderstorm at the grid of x=30 km. Red shading indicates a positive charge density (nC m^-3); blue shading indicates a negative charge density (nC m^-3). The peripheral profile is the mixing ratio of hydrometeors with 0.01 g kg^-1: (a) 10 min; (b) 45min; (c) 60 min; (d) 80 min. The legend is the same as in Fig. 6.

As the simulation continues and the convection develops, the electric field gradually intensifies. The effect of the inductive electrification parameterization then becomes considerable. With the charge separation process represented more clearly, the inverted dipole charge structure has already disappeared. The net charge density distribution of the thunderstorm displays a normal tripole charge structure, with an upper positively charged region between the temperatures of -20^°C and -30^°C, a significant negatively charged region in the middle level of the storm between -10^°C and -0^°C, and a smaller positively charged region at the lower level with a temperature higher than 0^°C. Most of the charge polarities of the graupel particles shift from negative to positive polarity. In the region where the graupel particles and supercooled water coexist, the inductive mechanism is actively manifested, and the net charge carried by the graupel particles gradually increases.

After 45 min, the electric field reaches a large value, and the inductive electrification scheme becomes increasingly significant. The inductive-charging parameterization is important to form the lower-level positive charge layer and for strengthening the intensity of the main negatively charged region, which is conducive to the generation of CG lightning. The storm exhibits a tripole charge structure, with a negative maximum charge density of -2.8 nC m^-3 in the middle of the negatively charged region and 2.6 nC m^-3 in the upper layer of the positively charged region at 60 min (Fig. 7c). When the simulation enters the dissipating stage at 80 min, both the intensity and the range of the charge regions reduce. As the downdraft gradually enhances, the graupel particles drop to the lower level and melt, forming precipitation particles with negative polarity.

The charge structure of the upper-negative charge and the lower-positive charge forms in the region of a thunderstorm cloud-anvil at 45 min, and the charge structure from the main charge region extends to the back of the storm. At 60 min, the charge structure of the cloud-anvil extends further and joins into the main charge region. As the simulation reaches 90 min, the charge structure in the cloud-anvil region of the storm disappears. The charge separation of the cloud-anvil is mostly formed by the forward region itself, as caused by the collision of hydrometeor particles that mainly originated from the core of the thunderstorm by the flow. Furthermore, with the convection and the updraft weakened, the collision opportunities and the charge separation of hydrometeor particles are reduced, resulting in the dipole charge structure of the cloud-anvil eventually disappearing.

It is assumed in the model that the convective region corresponds to the region where precipitation is greater than 5 mm h^-1, whereas the region with precipitation (not shown herein) lower than 5 mm h^-1 defines the trailing rearward region with weak precipitation. According to the position inferred in Fig. 3, Fig. 8a1 plots the distribution of the mixing ratio of the various hydrometeors with the height at the grid point of (x=50 km, y=20 km) that is located in the convective region. Cloud droplets are primarily distributed in the lower levels at 2-7 km MSL, with a maximum value of 1.2 g kg^-1 at 4 km MSL. Raindrops are concentrated in the lower levels, formed mainly by the melting of ice-phase particles. Graupel particles are located mainly in the middle of 6-9 km MSL, and the mixing ratio has a peak value of 1.1 g kg^-1 at a height of 6 km MSL. There is a region of graupel particles at the height of 2-4 km MSL, and some of them melt to form rain drops. Ice and snow crystals are mainly located in the higher levels, and the mixing ratio reaches a maximum at 7 km. The mixing ratio of hailstones is mainly located in the 4-6-km MSL level, with a maximum of 0.3 g kg^-1. Most of the negative charge density in the middle level is carried by rimed cloud droplets and graupel particles, while the positive charge in the lower levels depend mainly on larger graupel particles and hailstones. Lower cloud droplets contribute to an increasing mixing ratio of rain drops. The value of the charge carried by the hailstones is approximately 20- to 100-fold smaller than that carried by graupel particles (Fierro et al., 2007).

Figure 8.  Evolution of the mixing ratio of different hydrometeor particles (a1, b1), total charge density (a2, b2),the distribution of various hydrometeor particles in charge density (a3, b3) with height of the convective region (x = 50 km, y = 20 km) and of the trailing rearward region (x = 30 km, y = 40 km) at the simulation time of 60 min. Panels (a1–a3) represent the convective region; (b1–b3) represent the trailing rearward region.

The distribution of hydrometeors changing with altitude in the grid point of (x=30 km, y=40 km) in the trailing rearward region is plotted in Fig. 8b1. In general, the mixing ratio of different categories of hydrometeors in the trailing rearward region is apparently smaller than that in the convective region, which is the reason why the scaling used in Fig. 8a1 is different than that of Fig. 8b1. Cloud droplets are located mainly in the lower level at 3.5 km MSL and display a peak value of 0.22 g kg^-1, which is significantly smaller than that in the convective region. Raindrops are mainly distributed in the bottom of the storm with a peak value of 0.2 g kg^-1. The content of ice crystals is greater than other categories of hydrometeors, with a maximum value of 0.25 g kg^-1, peaking at 6 km MSL. Owing to the light mass of ice and snow crystals, both of them can be more easily transported backward by airflows from the top of the convective region into the transition area of the stratiform region (Carey et al., 2005), whereas graupel particles are less active in the transfer process due to mass detrainment. Hence, graupel particles are mainly concentrated in the middle level of the thunderstorm because of the rapid descent. The graupel-height in the rearward region is consistent with that in the convective region. The mixing ratio of graupel reaches a peak value of 0.19 g kg^-1 at the 8-km MSL level, which is caused by the riming processes of ice, snow, and supercooled water in the upper level of the thunderstorm. A small section of graupel is associated with updraft toward the back of the thunderstorm (Braun and Houze, 1994), indicating that the updraft in the stratiform region is strong enough to produce graupel particles. These results suggest that the in situ non-inductive charging mechanism is a crucial reason for most of the charge manifested in the stratiform region. It is more likely that the remnant charge in the stratiform region originated from the strong convective zone of the thunderstorm. The greater charge separation occurs in the convective mature zone in severe convective systems, which can be transported rearward to promote the occurrence of lightning. In the trailing rearward region, the negative charge is mainly carried by graupel particles, whereas the positive charge is associated with ice and snow crystals. (Schuur and Rutledge, 2000) demonstrated that, in the non-inductive charging mechanism, non-graupel particles play a major role in the electrification process of the stratiform region in thunderstorms.

Figures 8a2 and b2 show the distribution of the charge density changing with height in the convective area and the trailing rearward region. The comparison of the charge density in two regions indicates that the net charge in the convective area is larger than in the trailing rearward region. At the lower level of 3 km MSL in the convective region, there is a positive charge region with a maximum value of 0.9 nC m^-3. At higher levels, the charge density switches from positive to negative polarity, dropping to -1.3 nC m^-3 at the height of 6 km MSL. The positively charge region is located at the upper level of the thunderstorm with a peak charge density of 0.5 nC m^-3. Due to the weak updraft in the rearward region of the storm, the upper-positive charge density is far lower than in the convective region. In summary, the rearward region displays a small positive charge with a charge density peak value of 0.55 nC m^-3 at the lower level of 2-4 km MSL, a negative charge at approximately 6 km, which reaches a minimum value of -0.3 nC m^-3, and a positive charge at the upper level, with a center of charge density equal to 0.25nC m^-3 at the height of 9 km MSL.

Generally, the quantity of different categories of hydrometeors in the convective region is greater than that in the trailing rearward region. Some particles are transferred into the rearward region from the convective region by the airflow. Although the rearward region of the thunderstorm is weak and small (in the present case), it involves the charging mechanism via two separate electrification mechanisms: a dominant process, generated by the rearward region itself, and another process involving the charging of particles transferred from the convective region to the back of the thunderstorm, where the charge structure forms. With the stable charge structure, when the electric field reaches the threshold value, the discharge process is activated in the grid that satisfies the condition.

Figures 8a3 and b3 display the charge density of different hydrometeors for the convective region (x=30 km, y=20 km) and for the trailing rearward region (x=50 km, y=40 km) at the simulation time of 60 min. Without considering ions in the model, the charge release is barely tracked by the evaporation of rain and ice sublimation. Therefore, the total charge conservation cannot be established. In the convective region, the graupel particles are positively charged in the lower layer and negatively charged in the middle level, and hailstones are carried with positive charge in the lower level of 2-4 km MSL and negative charge in the higher level of 4-6 km MSL. Charge separation is generated between ice-phase hydrometeor particles, such as graupel particles and ice crystals, colliding with each other in the middle of the thunderstorm, enhancing the impact opportunity created by updraft. Ice crystals carrying negative charge are confined to the middle level of the thunderstorm, whereas positively charged particles are present in the upper level. Under the inductive parameterization, the charging mechanism of the cloud droplets is slightly weakened. Consequently, the hydrometeor particles are positively charged at the lower level of 3 km MSL and negatively charged in the core of the thunderstorm at approximately 7 km MSL. Raindrops are weakly positively charged between the ground and 4 km. Ice crystals have a charge density of negative polarity between 4.5 km and 6.5 km, and a positive one at the height of 6.5-11 km. Snow crystals are negatively charged in the middle layer of 4-5 km and positively charged at the level of 5-8 km. The charge density of the various hydrometeors fluctuates from the altitudinal range of 4 km MSL to 8 km MSL, where the maximum variation rate occurs between ice crystals and other hydrometeor particles. In the convective region, the charge is mainly obtained by the electrification of graupel and reaches a value of -1.15 nC m^-3 at 6.5 km MSL.

In the rearward region, the particles are mainly located in the lower level due to weakening updraft and the drag effect from the convective region. The cloud drops are positively charged in the lower level and negatively charged in the middle level. Graupel particles mainly carry a negative charge in the middle level, reaching a minimum charge density of -0.35 nC m^-3, and are positively charged in the upper level, with a charge density of 0.3 nC m^-3.

To examine whether the simulated charge structure distribution agrees with the actual situation, the lightning radiation sources distribution of the thunderstorm is investigated. As Fig. 9 shows, the lightning radiation sources distribution occurs at 30 min along the tangent line within a 2 km MSL interval from the convective region to the rearward region, as observed by the SAFIR 3000 lightning detection system located in Beijing area. More details of the SAFIR 3000 lightning network and data processing method are discussed by (Liu et al., 2011). The lightning radiation sources mainly display a three-layer distribution with an upper level centered at 9 km MSL, a middle level centered at 5 km MSL and a sparse density of lightning radiation source located at the level of 2 km that almost merges into the middle level. According to the results of (Carey et al., 2005) and (Ely et al., 2008), thunderstorms demonstrate a tripole charge structure with a negative charge region in the middle level of 6-8 km and a positive charge region in the lower level of 2-5 km, plus an additional positive charge region centered at the level of 9 km MSL, located at the top of the thunderstorm (Stolzenburg et al., 1998). To some extent, the simulated charge structure is successful in simulating the conditions as compared with the observed lightning activity of the present thunderstorm.

Figure 9.  Density of the lightning radiation source distributed with height of the thunderstorm at 1030 UTC within half an hour, as observed by the SAFIR 3000 lightning detection system. The range of lightning radiation sources occupy 2-km intervals.

Figure 10 depicts the simulated and observed evolution of the lightning frequency. The observed lightning data were obtained by the SAFIR 3000 lightning detection system, and the storm occurred and disappeared in the center of the lightning network. For determining the lightning number, the lightning radiation sources that occurred in one second within a 5 km distance are considered as one single lightning flash. As Fig. 10 shows, the simulated lightning increases gradually from the discharge starting, with a maximum value of 84 flashes per minute, which corresponds to the period from 0930 to 1110 UTC in the real thunderstorm, and the lightning frequency increases progressively to reach a peak value of 96 flashes minute. The simulated and observed maximum value of lightning frequency appear at the same moment, and both of them decrease sharply. The aforementioned results demonstrate that the simulated lightning frequency correlates well with observations. The tendencies of the simulated and observed lightning frequency curve basically coincide with each other. With the inclusion of the inductive mechanism in the model, the electric field increases rapidly and the discharge process begins at 20 min with an electric field of 148 kV m^-1 (not shown). From 30 min to 80 min, the electric field varies considerably, corresponding to a significant increase in the lightning frequency. After 80 min, both the amount of particles and the updraft decrease, implying a reduction in particle collisions and a gradual decline in charge density. It is anticipated that simulated storms produce slightly less total lightning activity compared with observed data. Without considering the lightning channel propagation, the lightning flash is identified by seven or more contiguous grids exceeding the breakdown value at the same time under the selected simple discharge parameterization. This threshold is prone to not being achieved in the simulation. As the thunderstorm enters the dissipating stage, the magnitude of the electric field and the lightning frequency decrease gradually, but the total simulated lightning rate is still higher than the observed rate.

Figure 10.  Temporal evolution of lightning frequency during the lifetime of the thunderstorm: black solid line stands for the simulation; gray solid line stands for the observation.

Using a 3D mesoscale model (RAMS) coupled with electrification and simple discharge parameterization schemes, the simulated charge structure of a summer thunderstorm presented a tripole feature, which is similar to previous results achieved using different models (Takahashi, 1984; Mansell et al., 2005). The results indicate that the main negatively charged region was located between the levels of 4 km and 7 km MSL, whereas the upper positively charged region extended to 10 km MSL, and a further positively charged region was located between the 2 km and 4 km MSL levels. In this study, the observational results of lightning radiation sources detected by the SAFIR 3000 lightning network indicated that the thunderstorm had a tripolar charge structure, with the negative charge situated in the middle level at 6-8 km MSL and a positive charge located in the upper level at 9 km MSL and lower levels at 2-5 km MSL. As a result, the simulated charge structure of the thunderstorm displays a tripolar charge structure, indicating good consistency with the observed situation.

This numerical simulation improves our understanding of the dynamics, microphysics and electrical processes within a thunderstorm, even in terms of the intensity of cosmic ray muons when they pass through thunderstorms to the ground (Wang et al., 2012). In this work, the electrical characteristics of a thunderstorm that occurred in Beijing were investigated. By the non-inductive Takahashi electrification parameterization, the inductive charging mechanism and discharge parameterization, the simulated charge structure of the storm was found to be in agreement with observations, and is likely applicable to most thunderstorms in the Beijing area. (Krehbiel et al., 2000) observed a tripolar charge structure with the main negative charge region at the 5-6-km MSL height and the positive charge region in the upper levels (Shao and Krehbiel, 1996). In addition to the tripole charge structure, other observations suggest a different type of charge structure in thunderstorms, converting from a simple dipole (Chauzy et al., 1985) to an inverted dipole (Rust et al., 2005) or a tripole structure (Williams et al., 1989; Qie et al., 2005, 2009; Zhao et al., 2009), and an even more complex multi-polarity (Stolzenburg et al., 1998). An examination of only one thunderstorm is insufficient to draw conclusions on the electrification and charge structure of thunderstorms. More cases should be investigated to confirm the charge structure of thunderstorms. Further research should focus on the simulation of the charge structure in combination with electric field sounding.