# Some explicit microphysical and electrical aspects of a cloud resolving model: description and thunderstorm case study

Gilles Molinié, Jean-Pierre Pinty, Frank Roux

<figure><img src="image_1.png" /></figure>

## To cite this version:

Gilles Molinié, Jean-Pierre Pinty, Frank Roux. Some explicit microphysical and electrical aspects of a cloud resolving model: description and thunderstorm case study. Comptes Rendus. Physique, 2003, 3 (10), pp.1305-1324. 10.1016/S1631-0705(02)01409-3. hal-00137537

HAL Id: hal-00137537

https://hal.science/hal-00137537v1

Submitted on 8 Jul 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Some microphysical and electrical aspects of a Cloud Resolving Model: description and thunderstorm case study

Gilles Molinié, Jean-Pierre Pinty, Frank Roux

Laboratoire d'aérologie, Observatoire Midi-Pyrénées, 14, avenue E. Belin, 31400 Toulouse, France

**Abstract** An electrification scheme, consistent with the mixed-phase microphysical parameterization, has been developed for the French cloud resolving model MésoNH. There are four successive steps: (i) charge separation is assumed to result only from non-inductive processes; (ii) electrical charges carried by the different hydrometeor species are transported along the air flow and redistributed according to the microphysical processes; (iii) the electric field is deduced from the integration of a modified Poisson equation; (iv) a lightning parameterization simulates triggering, propagation and pseudo-fractal branching of the flashes and associated charge neutralization. Two numerical experiments are conducted firstly to evaluate the performances of the lightning scheme, secondly to test the simulated evolution of the electrical characteristics of a idealized supercellular storm.

numerical modeling / convective clouds / storm electrification / lightning parameterization

## Aspects microphysiques et électriques d'un modèle explicite de nuage : description et étude d'un orage académique

**Résumé** Un schéma d'électrisation, cohérent avec la paramétrisation microphysique en phase mixte, a été développé pour le modèle numérique de nuage MésoNH. Il y a quatre étapes successives : (i) la séparation de charges est supposée résulter uniquement des processus non-inductifs ; (ii) les charges électriques emportées par les différents types d'hydromé-téores sont transportées par le flux atmosphérique et redistribuées par les processus microphysiques ; (iii) le champ électrique est déduit de l'intégration d'une équation de Poisson modifiée ; (iv) une paramétrisation des éclairs simule le déclenchement, la propagation et les branchements pseudo-fractals des décharges et les neutralisations de charge associées. Deux expériences numériques sont conduites, d'abord pour évaluer les performances du schéma d'éclairs, ensuite pour tester l'évolution simulés des caractéristiques électriques d'un orage supercellulaire idéalisé.

modélisation numérique / nuages convectifs / électrisation des nuages / paramétrisation des éclairs
# 1. Introduction

The use of a Cloud Resolving Model to study electrification processes is a powerful tool to investigate the electrical activity inside thunderstorms and the characteristics of lightning flashes (e.g., [1–5], among others). The reason is that many intriguing aspects of electrical phenomena are now thought to be closely related to the dynamics and to the microphysical nature of the storms. So probably, only a physically consistent, three dimensional (3D), time evolving picture of numerical storms will give credit to the representation of the lead mechanisms that control the complex space charge structure and the rate, location, length, and polarity of the flashes. In that sense, electrification of simulated storms has a growing importance for cloud modellers [6] because modern observing lightning devices, either spaceborne detectors (Lightning Imaging Sensor (LIS), see [7], radio frequency and optical sensors aboard FORTÉ [8]) or ground based VHF receivers (Lightning Mapping Array (LMA), [9] or ONERA’s VHF-interferometer (ITF), [10]) can now give many details on the electrical characteristics of the thunderstorms.

This preliminary work aims at presenting an explicit electrification scheme which is closely linked to a bulk microphysical mixed-phase scheme routinely used in the French non-hydrostatic mesoscale model [11]. The additional lightning scheme and charge neutralization treatment borrow many ideas recently expressed by [5,12]. Because the present development has been made in the framework of a mesoscale model which is supplied with the grid-nesting technique and with a meteorological initialization scheme (see [13] on this topic), one issue of this work is to prepare the way for simulating the electrical aspects of actual cloud situations at fine (kilometer) scale. In short, the project of parameterizing cloud electrification and lightning flash occurrence in MésoNH tries to favor the full coupling of the dynamics, microphysics, and electrical characteristics within clouds.

As previously done in many schemes, an exclusive non-inductive charge separation mechanism based upon ice-particle interactions in the presence of supercooled water is assumed. Several existing parameterizations quantatively describe the sign and magnitude of this charge transfer and most of them can be selected in the scheme for the purpose of sensitivity tests. The electrical charges are carried by the five different condensate species (cloud droplets, rain drops, primary ice crystals, aggregates and graupeln) and so they are transported by the dynamics of the storm as well. The charges are also exchanged between each category of water substance at the rate dictated by the microphysical interactions that continuously occur. For this reason, dynamics and microphysics play an important role in cloud electrification process through their contribution in the continuity equation associated to each type of charge (see [14] for a comprehensive description).

In complement to the dynamical and microphysical processes responsible for the early electrification, a lightning parameterization is required to constrain the growth of the intense in-cloud electric field within an observed range of values. The electric field is the result of the multipolar structure of a storm due to charge accumulation as organized by the flow and by the sorting of the hydrometeors due to their different sedimentation rates [15,16]. At present, the true mechanism for lightning initiation, the electric breakdown, is still uncertain [17], but enough observational evidence indicates that sufficiently low electric fields (less than 100 to 150 kV·m⁻¹) are able to trigger lightning flashes [18]. There are several reasons that explain why an accurate simulation of the onset, location, and structure of individual lightning flashes is a very difficult task. Among them, the coarse model resolution as compared to the filamentary aspect of the ionized channels and the still poorly understood physics of the discharges and their propagation in cloudy air ([19, 5, 12] are major obstacles. Nevertheless, recent criteria based upon isolines of maximum electric field with possible connections to nearby regions of high space charge density seem to draw a realistic picture of the lightning flash morphology [12]. Note that the charge neutralization along the whole lightning path with respect to total charge conservation is the prime result of a lightning scheme. Therefore, the efficiency of neutralization is linked to the lightning frequency because the charging mechanisms are comparatively slow.
Although there are several shortcomings in the representation of the electrical charge cycle with steady separation phases (electrification) and sporadic neutralization phases (lightning), this work attempts to test whether a current 'state-of-the-art' treatment of cloud electrification and lightning is able to catch the salient features of the electrical activity in 3D thunderstorms. In the following, Section 2 summarizes the MésoNH model characteristics and the mixed-phase microphysical scheme is described in Section 3. Section 4 is devoted to the electrification scheme and to the lightning scheme. Results of the numerical simulations are presented in Section 5. It is shown that the model reproduces some of the electrical characteristics of 3D thunderstorms.

## 2. The MésoNH model

The MésoNH code is a non-hydrostatic mesoscale model which results from a joint development of Laboratoire d'Aérologie and Météo-France [11]. The model integrates an anelastic system of equations that is able to simulate academic and real atmospheric flows ranging from the large eddy turbulent motion to the synoptic scale. A recent description of some model characteristics can be found in [13]. It is worth mentioning that the code is fully vectorized and parallelized so the 3D evolution of clouds with electric effects can be simulated on large grid domains.

## 3. The mixed-phase microphysical scheme

### 3.1. Generalities

In its essence, the scheme follows the approach of [20] in that it is a three-class ice parameterization coupled to a Kessler's scheme [21] for the warm processes. The scheme predicts the evolution of the mixing ratios of six water species: $r_v$ (vapor), $r_c$ and $r_r$ (cloud droplets and rain drops) and $r_i$, $r_s$ and $r_g$ (pristine ice, snow/aggregates and frozen drops/graupeln defined by an increasing degree of riming). The concentration of the pristine ice crystals, here assumed to be plates, is diagnosed. The concentration of the precipitating water drops and ice crystals is parameterized according to [22], with the total number concentration $N$ given by:

$$N = C\lambda^x, \qquad (1)$$

where $\lambda$ is the slope parameter of the size distribution, $C$ and $x$ are empirical constants drawn from observations. The size distribution of the hydrometeors is assumed to follow a generalized $\gamma$-law:

$$n(D) dD = Ng(D) dD = N \frac{\alpha}{\Gamma(\nu)} \lambda^{\alpha\nu} D^{\alpha\nu-1} \exp(-(\lambda D)^{\alpha}) dD, \qquad (2)$$

where $g(D)$ is the normalized form which reduces to the Marshall-Palmer law when $\alpha = \nu = 1$. Finally, simple power laws are taken for the mass-size ($m = aD^b$) and velocity-size ($v = cD^d$) relationships to perform useful analytical integrations using the moment formula:

$$M(p) = \int_{0}^{\infty} D^{p} g(D) dD = \frac{\Gamma(\nu + p/\alpha)}{\Gamma(\nu)} \frac{1}{\lambda^{p}}, \qquad (3)$$

where $M(p)$ is the $p$th moment of $g(D)$. A first application of (3) is to compute the mixing ratio $r_x$ as:

$$\rho r_x = a N M_x(b). \qquad (4)$$

Table 1 provides the complete characterization of each ice category and of the cloud droplets/raindrops.
<table><caption>**Table 1.** Characteristics of each hydrometeor category. $\alpha$ and $\nu$ are used in Eq. (2). The other coefficients are related to power law relationships for the mass ($m = aD^b$), the fall speed ($\nu = cD^d$), the charge ($q = eD^f$) where D is the particle size, and for the concentration in Eq. (1).</caption><thead><tr><th>Parameters</th><th>r<sub>i</sub></th><th>r<sub>s</sub></th><th>r<sub>g</sub></th><th>r<sub>c</sub></th><th>r<sub>r</sub></th></tr></thead><tbody><tr><td>α</td><td>3</td><td>1</td><td>1</td><td>3</td><td>1</td></tr><tr><td>ν</td><td>3</td><td>1</td><td>1</td><td>3</td><td>1</td></tr><tr><td>a</td><td>0.82</td><td>0.02</td><td>196</td><td>524</td><td>524</td></tr><tr><td>b</td><td>2.5</td><td>1.9</td><td>2.8</td><td>3</td><td>3</td></tr><tr><td>c</td><td>800</td><td>5.1</td><td>124</td><td>3.2 × 10<sup>7</sup></td><td>842</td></tr><tr><td>d</td><td>1.00</td><td>0.27</td><td>0.66</td><td>2</td><td>0.8</td></tr><tr><td>f</td><td>0.5</td><td>1.3</td><td>2.0</td><td>0.5</td><td>1.3</td></tr><tr><td>C</td><td></td><td>5</td><td>5 × 10<sup>5</sup></td><td></td><td>10<sup>7</sup></td></tr><tr><td>x</td><td></td><td>1</td><td>-0.5</td><td></td><td>-1</td></tr></tbody></table>

## 3.2. A short overview of the scheme

The pristine ice category is initiated by homogeneous nucleation (HON) when $T \le -35$ °C, or more frequently by heterogeneous nucleation (HEN), so the small ice crystal concentration is a simple function of the local supersaturation over ice. These crystals grow by water vapor deposition (DEP, see below) and by the Bergeron–Findeisen effect (BER). The snow phase comes by autoconversion (AUT) of the primary ice crystals; it grows by deposition (DEP) of water vapor, by aggregation (AGG) through small crystal collection and by the light riming produced by impaction of cloud droplets (RIM) and of raindrops (ACC). The graupels are formed as a consequence of the heavy riming of snow (RIM and ACC) or by rain freezing (CFR) when supercooled raindrops come in contact with pristine ice crystals. Distinction between light and heavy riming is made on the basis of a critical size of the snowflake (droplets) or by estimation of the mean density of the resulting particle (raindrops). According to the heat balance equation, graupel can grow more efficiently in the (WET) mode than in the (DRY) mode when riming is very intense (as for hailstone embryos). In the latter case, the excess of non-freezable liquid water at the surface of the graupels is shed (SHD) to form raindrops. When $T \ge 0$ °C, pristine crystals immediately melt into cloud droplets (MLT) while snowflakes are progressively converted (CVM) into graupels which melt (MLT) as they fall. The other processes are those described by the Kessler scheme: autoconversion of cloud droplets (AUT), accretion (ACC) and rain evaporation (EVA). Cloud droplets excepted, each condensed water species has a non-zero fallspeed. See Fig. 1.

## 3.3. Some details of the scheme

One difficulty inherent to mixed-phase clouds is the possible coexistence of cloud droplets and small ice crystals, which necessitates a special treatment of the fast vapor exchanges (DEP and CND). It is important to keep attention to this question because the amount of supercooled cloud water has a strong influence in the parametrization of the charge separation mechanisms (see next section). As is usually done, the ‘floating’ water vapor saturation pressure, $r_{v,c,i}^{sat}$, is defined by a barycentric formula using the vapor saturation curves over water and ice and the mass amounts $r_c$ and $r_i$, respectively. In the parameterization,
<figure><img src="image_3.png" /><figcaption>Figure 1. Microphysical processes included in the mixed-phase scheme (see text for the acronyms and explanations).</figcaption></figure>

the DEP and CND terms result from an implicit adjustment relative to $r_{v_{c,i}}^{sat}$, but with an original closure where any deficit/excess of $r_v$ due to the adjustment, is compensated/absorbed by each phase in proportion to its actual amount. This is in contrast to other schemes where the closure is based upon an artificial linear function of temperature for $0 \le T \le -40$ °C. The adjustment algorithm is non-iterative and 2nd order accurate. Note that the BER effect is treated independently in an explicit way.

The other family of process that require careful treatment are the collection processes. When non- ($r_c$) or very slowly ($r_i$) precipitating categories are involved, the collection rates are computed analytically using the geometric sweep-out concept of the collection kernels. When both particles are precipitating, an analytical integration over the spectra is no longer possible and the use of pre-tabulated kernels have been preferred to approximate integrations. However in any case of ice-ice interaction, a major point of concern resides in the tuning of the sticking efficiencies which are still poorly understood functions of temperature in most cases. After a series of experiments, the last set of coefficients retained by [23] has been adopted in the scheme. Note that these efficiencies intervene also in the computation of the rebounding collision rates (inefficient microphysical mass transfers) which are involved in the primary charge separation mechanisms.

# 4. The electrification scheme

## 4.1. Electrical variables

The storm electrification scheme aims to describe complex electrical phenomena that originate from the microphysical state of the storm. Once separated by ice-ice elastic collisions, the electrical charges are transported with the hydrometeors and are redistributed between the hydrometeors according to the various microphysical conversions. A system of prognostic equations is numerically solved for the temporal and spatial changes of the charge density carried by each of the five microphysical species. Only a single charge type (positive or negative) is considered at each grid point for a given type of microphysical particle. Electrical charges attached to water vapor are deemed to simulate the formation of ions when particles evaporate or sublimate. This additional reservoir of charges must be considered to keep the total charge within the storm. The charge density $q$ of individual particle is assumed to obey a power law relationship with the characteristic size $D$ of the particle:

$$q(D) = eD^f. \qquad (5)$$
The above coefficients *f* are chosen as recommended by [24], their values are recalled in Table 1. The bulk charge density of each microphysical specy *x* is a prognostic variable. It results from integration of Eq. (5) over the number concentration of the hydrometeors *n*<sub>*x*</sub> defined by Eq. (2) using the moment definition of Eq. (3):

$$q_x = \int_0^\infty q(D)n_x(D) dD = eNM_x(f). \qquad (6)$$

The parameters *e* are deduced from equation (6) at the beginning of each time step.

## 4.2. Charge separation mechanisms

As a result of many laboratory studies [25–28], the non-inductive (NI) charge transfer after rebounding collisions between ice crystals and graupels is likely the dominant mechanism which must be taken into account in a first step. The different parameterizations of NI processes that are available in the scheme are described now. Note that the sign and magnitude of charge transfer $\delta q_{NI}^P$ for parameterization *P*, are those of the biggest (graupel) particle in the following schemes.

Historically, laboratory work by [25] leads to the first formula for graupel that may charge positively or negatively depending on the cloud water content (LWC) and temperature,

$$\delta q_{NI}^T = f_T(T, \text{LWC}) \text{Min} \left[ 10, 5 \left( \frac{D_{i,s}}{D_0} \right)^2 \frac{|v_g - v_{i,s}|}{v_0} \right], \qquad (7)$$

where $f_T(T, LWC)$ results from a linear ($T$) and logarithmic (LWC) interpolation of tabulated values from Fig. 8 of [25] as in [4,29]. $D_{i,s}$ is the ice crystal or snowflake diameter and $v_g$ and $v_{i,s}$ are the terminal velocities of graupels and crystals. The scaling values of $D_0$ and $v_0$ are $10^{-4}$ m and $8 \text{ m} \cdot \text{s}^{-1}$ respectively.

In order to retain the essential features of many charge separation mechanisms in a very simple formula, [2] suggest that one considers solely a polarity reversal at a fixed temperature, $T_C = -10 \text{ }^{\circ}\text{C}$ with $T_C = T - 273 \text{ K}$, and an associated fixed quantity of separated charge per collision leading to

$$\delta qNI^{HF} = \begin{cases} 2 \times 10^{-15} & \text{if } T_C > -10 \text{ }^{\circ}\text{C}, \\ -2 \times 10^{-15} & \text{if } T_C \le -10 \text{ }^{\circ}\text{C} \end{cases} \qquad (8)$$

for graupel/ice crystals interactions and,

$$\delta qNI^{HF} = \begin{cases} 2 \times 10^{-13} & \text{if } T_c > -10 \text{ }^{\circ}\text{C}, \\ -2 \times 10^{-13} & \text{if } T_c \le -10 \text{ }^{\circ}\text{C} \end{cases} \qquad (9)$$

for graupel/snowflakes interactions.

Later and using laboratory results of [26] and in situ measurements, Gardiner et al. [30] indicate that the charge transferred to graupel particle could be parameterized by

$$\delta qNI^G = 73D_{i,s}^4|v_g - v_{i,s}|^3 (\text{LWC} - \text{LWC}_{\text{crit}}) f_G(T_C), \qquad (10)$$

where $f_G(T_C) = -1.7 \times 10^{-5}T_C^3 - 0.003T_C^2 - 0.05T_C + 0.13$ with $T_C < 0 \text{ }^{\circ}\text{C}$ and $\text{LWC}_{\text{crit}} = 10^{-4} \text{ kg} \cdot \text{m}^{-3}$ is the critical value of the liquid water content.

Following the same approach, [31] proposes a similar expression but with several sets of parameters (prefactor and exponents) dependent on the range of $D_{i,s}$:

$$\delta q_{NI}^{S91} = F D_{i,s}^m |v_g - v_{i,s}|^n f_{S1}(T_C, \text{ELWC}). \qquad (11)$$

Now ELWC is the effective liquid water content, the efficient portion of LWC for riming on graupels. However, the multiple discontinuous charging regimes and ‘anomalies’ seen by [31] render the use of (11)
very delicate (see [29] for a useful synthesis about the constant values and the definitions of $f_{S1}(T_C, ELWC)$ in (11), but also for some appropriate modifications brought to the original scheme). The many charge transfer measurements operated in the UMIST cloud chamber lead [32] to recast $f_{S1}(T_C, ELWC)$ in a much compact form

$$f_{S2}(T_C, ELWC) = \begin{cases} 6.74 \times \text{RAR} + 1.36 \times T_C + 10.05 & \text{if } \text{RAR} > \text{RAR}_{\text{crit}}, \\ 2.95 \times \text{RAR}^2 - 10.59 \times \text{RAR} + 3.02 & \text{if } \text{RAR} \le \text{RAR}_{\text{crit}}, \end{cases} \quad (12)$$

where RAR is assimilated to the rime accretion rate with a critical value RARcrit, a polynomial function of T_C, that determines the charge sign reversal.

The full integration of the different charge separation parametrisations over the particle distribution law (2) is performed with the same technique as for the microphysical scheme (see, for instance, [1]). Three types of ice-ice interactions are considered in the scheme: cloud ice-snow (AGG), cloud ice-graupel and snow-graupel (both referred to DRY). The first two charging mechanisms involving cloud ice are integrated analytically but the last one needs look-up tables because the impact velocity $|v_g - v_{i,s}|$ term precludes an exact integration.

### 4.3. Charge transfers

Charge transfer rates are basically computed by integrating mass transfer rates multiplied by the charge-dimension relationship [24]. Some microphysical processes involving sizeable hydrometeors are not directly integrable because of the absolute value of the hydrometeor fall speed difference. These cases are treated by means of look-up tables as in [33]. The specific case of water vapor adjustment over the cloud water and cloud ice (DEP and CND assumes that the charge carried by the flux of water vapor is proportional to the two third power of the exchanged mass).

### 4.4. Integration of the electric field

The electric potential V is solution of a Poisson's equation forced by the local net charge density, $q_{total}$. This charge density is the algebraic sum of the electric charges carried by the hydrometeors and water vapor (pseudo ions) inside the cloud and the fair weather charges, derived from the standard electric field profile, outside the cloud. The electric field $\vec{E}$ is the gradient of V. However, in order to draw profit from the existing elliptic pressure solver of MésoNH (see Section 2.1), we choose to define a pseudo electrical potential by

$$\vec{\nabla} V = -\frac{\vec{E}}{\tilde{\rho}}, \quad (13)$$

where $\tilde{\rho}$ is the reference density of the anelastic system of equations (see [11]). The MésoNH solver can then be used to integrate directly a modified Poisson equation which reads

$$\nabla \cdot (\vec{E}) = \frac{q_{\text{total}}}{\epsilon_0} = -\nabla \cdot (\tilde{\rho} \vec{\nabla} V), \quad (14)$$

where $\epsilon_0$ is the dielectric constant of the air. In MésoNH, solution of Eq. (14) is obtained with appropriate boundary conditions in any curvilinear non-orthogonal coordinate system when following-terrain coordinates and/or map factors are considered in case of orography and/or geographical projections, respectively. Top and bottom boundary conditions are constituted by two perfect conductors (ionosphere and earth surface) so $\vec{E}$ is orthogonal to these surfaces (mirror effect). Note that the hypothesis of ionospheric conditions would need a last level model as high as 80 km, so these conditions are approximate when the last model level lies in the lower stratosphere. The lateral boundary conditions are of Neuman type

$$\vec{n} \cdot \vec{\nabla} V = \vec{n} \cdot E_{fw}, \quad (15)$$
where $\vec{n}$ is the vector normal to the surface and $\vec{E}_{fw}$ is the standard fair weather electric field given by

$$ \vec{E}_{fw} = \vec{E}_0 \exp \left( -K_E \left( H \frac{z - z_s}{H - z_s} \right) \right). \qquad (16) $$

$\vec{E}_0$ is vertical, directed downward and its module is of 100 V·m⁻¹. $K_E = 292 \times 10^{-6}$ m⁻¹, H is the height of the last model level, $z_s$ the maximal relief height. Although the relief is known to enhance the ambient electric field due to the electrode effect [34], no modification of (16) is available for complex terrain. It is believed that the electric fields are sufficiently high in clouds to overcome this effect. Conversely the assumption of far enough lateral boundary conditions with respect to the magnitude of the electric field found in the thunderstorms is probably less satisfied.

## 4.5. The lightning parameterization

### 4.5.1. Introduction

The lightning flashes result in the propagation of electrical discharges in a hot (~30 000 K) narrow channel. The flashes are triggered by the breakdown of intense electric field of order 100 to 300 kV m⁻¹ due to corona emission or to electron avalanche. The emerging theories for the discharge propagation in thundercloud is the bi-leader theory [35] and the dielectric breakdown model [36] with a recent application to the simulation of 2D lightning flashes [5].

At the origin of the bi-leader theory, two segments of opposite polarity propagate in two opposite directions from an initial point where $|E|$ is the highest. The positive (negative) flash elements propagate toward negative (positive) values of $\vec{E}$. The propagation of a discharge occurs when ionisation processes at the tip of the lightning channel induce a sufficient charge current toward the channel to enhance the electric field at its tip. The propagation of the negative segment occurs by steps of 20–90 m at 20–100 µs intervals while positive segments propagate in a continuous way [14].

In a dielectric breakdown model, the channel steps are chosen at random according to a probability of bond generation. Once a new grid point is added to the discharge channel, the Poisson's equation (14) is solved to update the electric potential V because the growing lightning channel behaves as a perfect conductor. Despite its heavy computational cost due to a repetitive resolution of (14), the strength of this approach is that complex branching structures on the flashes can be simulated in a very realistic way [5].

Note that at the end and in any case, the quantity of charge neutralized in the multiple channels is exactly balanced by an equal quantity of charge of opposite polarity unless the lightning does not hit the ground.

### 4.5.2. The lightning channels

For the sake of efficiency, the present lightning proposal follows the parameterization of [12] which mimics a first 'bi-leader' phase completed by a 'multiple streamer' phase. As this model of discharge development is purely diagnostic, no time sequence nor time scale are assumed.

The procedure starts when a lightning flash is triggered at a grid point where the electric field magnitude reaches a threshold value given in [18]:

$$ E_{\text{init}}(z) = \pm 167 \rho(z), \qquad (17) $$

where $E_{\text{init}}$ has units of kV·m⁻¹. Apart from this grid point, a flash traces $\vec{E}$ in both parallel and antiparallel directions until the magnitude of the ambient electric field at each end falls below $E_{\text{stop}} = 15$ kV·m⁻¹. The channel is drawn by segments linking grid points that are in the line of sight to which the local $\vec{E}$ is pointing in a 3D geometry. When the lower end of the channel falls below 1000 m AGL, the lightning is directly connected to the ground. The next expansion phase is initiated to account for the tortuous aspect of the flash propagation from some selected forking points where the net charge density $|q_{\text{total}}(x, y, z)|$ is larger than
0.5 nC·m⁻³ and where the local electric potential |V(x, y, z)| is larger than the mean potential of the leader channel. In order to grossly introduce a true fractal dimension to the simulated lightning [36], an iterative algorithm is set to generate a stochastic stepwise branching with a number of branch obeying the scaling law [37]

$$ \mathcal{N}(d) \sim d^{\chi^{-1}}, \qquad (18) $$

where N(d) is the number of discharge lines that crosses a sphere of radius d centered at the origin of the lightning flash (the initial grid point where |E| is maximum). The exponent χ is the Hausdorff or fractal dimension with 1 < χ < 2 (2 < χ < 3) in a 2D (3D) geometry. The algorithm works as follows. First it seeks for the contiguous grid points of the leader channel that fulfills the above charge and potential requirements. Then a connection probability, taken as inversely proportional to the number of possible bonds that can be formed from these grid points, is computed. The adjunction of new bonds at a distance d is explored by progressing along each of the bi-leader path. The N new branches estimated from (18) are selected at random from the subset of grid points having the highest probability. This growth mode attempts to reproduce qualitatively the screening effect of the dielectric breakdown model. Once the bi-leader paths are explored, the procedure of constructing new sub channels is repeated for the leaders and the first streamers and so on until no more branching possibility exists.

### 4.5.3. Charge neutralization

Because of the very high temperatures inside the lightning channels, air molecules are ionized like a plasma, and ions migrate at the surface to maintain the electrostatic equilibrium. These mobile ions further diffuse away from the channel and so partially neutralize the charged hydrometeors of opposite polarity. The net result is that lightning flashes eliminate an equal amount of positive and negative charges and hence damp the electric field in the cloud.

The charge neutralization along the complex lightning channel essentially follows [12]. The net charge densities in excess of qneut are neutralized with a simple scheme that preserves an exact polarity balance. The charge density induced by the lightning is then redistributed over each hydrometeor category in proportion to their total surface area.

# 5. Numerical simulations

## 5.1. Static tests of the lightning scheme

Setting-up the lightning scheme is a complex operation in the context of 3D parallelized code so beginning with a simple 2D test case is very helpful to test the algorithm and to illustrate the merit of the original pseudo fractal description of the streamers.

A realistic charge density distribution must be prescribed to evaluate the lightning scheme so the simple bipolar structure reproduced in Fig. 2 is adopted. The domain covers 25 km on the horizontal and 12 km on the vertical with a resolution of 1 km and 500 m, respectively. Two oblate ellipsoidal zones with semi-axes of 8.5 km and 2 km and filled with a uniform electrical charge of ±0.7 nC·m⁻³ (so with opposite polarity), are prescribed in the 2D domain. In the lower 6 km, the negative charge zone is filled by a second pocket of charge of −1 nC·m⁻³ while in the upper 6 km, the positive charge zone is filled by two adjacent pockets of charge, each containing +1 nC·m⁻³. The maximum electric field obtained in this configuration is 240 kV·m⁻¹. Three experiments are performed with the scheme. They differ by the choice of the fractal dimension χ.

First we concentrate on the results of Fig. 2, for χ = 1.3. The thick line traces the path of the two leaders. The bi-leader is initiated at random among the grid points where |E| > 150 kV·m⁻¹. The ignition point is marked by the gap between the two branches. The descending leader propagates vertically toward the negatively charged zone. It stops when |E| < 15 kV·m⁻¹. The ascending leader follows a quasi symmetrical path which passes between the two positive pockets of positive charges. The slight asymmetry at the end of
<figure><img src="image_4.png" /><figcaption>**Figure 2.** 2D lightning branching test with varying fractal dimension $\chi$ (Eq. (18)) with: (a) $\chi = 0.3$; (b) $\chi = 0.5$; and (c) $\chi = 0.7$ (see the text for the charge configuration). The two leaders are plotted with a bold solid line.</figcaption></figure>

the branch (tilting to the left) is due to the left shift of the main charge centers which makes them slightly closer of the left boundary of the domain. The result is that, contrarily to the positive leader case, the negative leader cannot reach the main positive charge centers because $\vec{E}$ is too weak there. The consequence is that the ‘bi-leader phase’ does not neutralize a large quantity of charge so $\vec{E}$ remains high and several lightning flashes are necessary to reduce the magnitude of $\vec{E}$.

The usefulness of fractal streamers is illustrated by the branches with dashed lines in Fig. 2. These tortuous segments considerably lengthen the bi-leader path. The positive leader develops symmetric branches which progress through the pocket of negative charge center. The same feature occurs for the negative leader with a complex structure emerging from the leader tip to reach each of the two positive charge pockets. Note that the stochastic nature of the growing streamer segments leads to non-symmetrical branching trees and to a remarkable space occupancy in the pockets. Hence, the pseudo fractal representation of the streamers clearly allows for a better lightning flash development because the main charge centers can be fully connected, so that a much efficient charge neutralization is possible.

The effect of varying $\chi$ is illustrated by the three plots in Fig. 2. Increasing the $\chi$ value in (18) favors an early branching along the leader path. This feature is well illustrated for the positive leader by comparing cases with $\chi = 1.3$ and with $\chi = 1.5$. With $\chi = 1.5$, the development of branches at long distance from the lightning initiation point simulates a screening effect because no branch can grow from the tip of the positive
<table><caption>**Table 2.** Total length of the lightning and quantity of charge along the lightning paths corresponding to Fig. 2. The total charge is the sum of the positive and negative charges.</caption><thead><tr><td>Lightning element</td><td>χ</td><td>Length (km)</td><td>Charge (C)</td></tr></thead><tbody><tr><td>bi-leader (only)</td><td>(0)</td><td>7.7</td><td>2.3</td></tr><tr><td>bi-leader and streamers</td><td>1.3</td><td>49.6</td><td>16</td></tr><tr><td>bi-leader and streamers</td><td>1.5</td><td>78.0</td><td>21</td></tr><tr><td>bi-leader and streamers</td><td>1.7</td><td>83.7</td><td>21</td></tr></tbody></table>

leader as observed when $\chi = 1.3$. Looking now to the (upward) negative leader, one can notice also a rapid quasi-symmetrical branching toward the negative charge centers in the case of $\chi = 1.5$. This structure develops because the local charge is larger than the threshold value $(0.7 \text{ nC} \cdot \text{m}^{-3})$ against $0.5 \text{ nC} \cdot \text{m}^{-3}$). However the presence of this arborescence hampers the development of the upper branch which grows from the tip of the leader (screening effect). The case with $\chi = 1.7$ is close to the case with $\chi = 1.5$ except for the increased densification of the branching around the negative leader.

Table 2 recalls the lightning length and the quantity of charges that are attained by the lightning segments as $\chi$ varies. Roughly speaking, Table 2 reveals that the sole bi-leader represents 10–15% of the total length and charge. This means that the stochastic model of streamer is very efficient to extend the lightning, particularly in the horizontal direction. Increasing $\chi$ results in longer lightning paths and more charge neutralized. Large values of $\chi$ do not modify substantially the lightning characteristics because the charge criterium and the screening effect limit the lightning expansion. However, this lightning parameterization suffers of the poor knowledge of many adjustable parameters in the ‘bi-leader phase’ so more sensitivity experiments are needed to get a robust tuning of the scheme.

## 5.2. 3D case study of a severe storm

A first application of this work is to investigate the charging processes and subsequent lightning characteristics in a variety of deep convective systems from 2D squall lines to 3D isolated storms. Here we choose to concentrate on the simulation of an prototype of 3D supercellar storm [38] which presents numerous ingredients for a sustained electrical activity. The case study is characterized by a well organized dynamics with an intense updraft ($> 20 \text{ m} \cdot \text{s}^{-1}$) surrounded by several evaporative downdrafts. The microphysical processes that take place produce a copious amount of ice particles with different degree of riming but also maintain a sufficient amount of supercooled cloud water to produce rimed particles. Thus, these environmental conditions appear to be very favourable for the formation and accumulation of electrical charges.

The initial sounding is given in Fig. 3. It is characterized by a strong instability in a sheared environment with veering and by a very dry atmosphere aloft. Convection is initiated by placing a warm 1.5 K bubble of radius 10 km in the planetary boundary layer. The domain of simulation is $40 \times 40 \times 30$ points with a resolution of 1 km in the horizontal and 0.5 km in the vertical. The timestep is 2 s.

### 5.2.1. General characteristics

In Fig. 4, horizontal cross sections of vertical velocity and horizontal wind at 5 km altitude, and precipitation rate at ground level are shown at three moments of the simulation: lightning initiation stage (40 min), mature stage (60 min), dissipation stage (80 min). Convection organizes from the buoyant bubble placed in the southwestern quadrant of the domain because the thunderstorm develops toward the northeast during its lifetime of more than 90 min. An important consequence of the directional wind shear is the three-dimensional organization of the updraft and downdraft, so that falling precipitation that partly evaporates
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Synthetic sounding of the storm environment adapted from Klemp and Wilhelmson [38].</figcaption></figure>

in low levels at the rear of the storm, supplies additional moisture to the air inflow, leading to stronger convective motions during the simulation. Results are now analysed along the vertical cross-sections marked in Fig. 4 which approximately pass through the updraft core in the direction of the main flow.

At 40 min of simulated time (Fig. 5), upward velocities reach $18 \text{ m·s}^{-1}$ at 7 km altitude and compensating downward motions of $-7 \text{ m·s}^{-1}$ start to form at 8.5 km altitude on the northeastern side of the updraft. The condensation of water vapor in the updraft produces non-precipitating cloud droplets and ice crystals with maximum mixing ratios of 2.2 and $2.6 \text{ g·kg}^{-1}$. At this stage, precipitating snow, graupel and rain are already formed with colocalized small ice crystals and large ice hydrometeors (snow and graupel). Furthermore and as deduced from the freezing level (273 K isoline in Fig. 5(a)), some supercooled cloud and rain exist above 4 km altitude, a favourable condition for storm electrification to occur according to most of the charge separation parameterizations. At 60 min (Fig. 6), the updraft vertical velocity still exceeds $18 \text{ m·s}^{-1}$ at 6 km altitude. Several downdrafts with $|w| < 5 \text{ m·s}^{-1}$ are present. The evaporative downdraft of $\sim -3 \text{ m·s}^{-1}$ at 2 km altitude separates two updraft branches. The hydrometeor contents keep their peak value of $2.5 \text{ g·kg}^{-1}$ for cloud water, $2 \text{ g·kg}^{-1}$ for cloud ice, $3 \text{ g·kg}^{-1}$ for rain, $0.8 \text{ g·kg}^{-1}$ for snow and $5 \text{ g·kg}^{-1}$ for the graupel. Above 6 km altitude, hydrometeors are carried in the anvil which extends in the direction of the upper flow At 80 min (Fig. 7), the updraft has splitted and weakened with velocities lower than $12 \text{ m·s}^{-1}$, and the low-level downdraft is increased approaching $-5 \text{ m·s}^{-1}$. Non precipitating cloud water and ice contents are lower ($<1.5 \text{ g·kg}^{-1}$) due to the decrease of the water vapor supply by the 7 updraft. This leads to a reduction of the graupel content ($<4.5 \text{ g·kg}^{-1}$) which however remain relatively high aloft due to the time necessary for these hydrometeors to fall. As a consequence, the anvil still spreads northeastwards while sinking between altitudes of 5 and 7 km.

### 5.2.2. *Storm electrification*

Among the parameterizations of the non-inductive charging processes available in the model (see Section 3.2), the simplest one of [2] has been chosen to illustrate the model capabilities. The temperature of charge reversal (TCR) is set to $-10 \text{ }^{\circ}\text{C}$. This parameterization is independent of the supercooled liquid water content, so charge exchanges occur when both graupel and snowflakes or ice crystals coexist. The influence of charge neutralization by lightning on the storm electrical structure will be discussed later.
<figure><img src="image_6.png" /><figcaption>(a)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

**Figure 4.** Horizontal cross section of the wind vector (horizontal component in arrows) and of the wind vertical velocity (hatched area) at 5000 m height and the instantaneous precipitation rate at ground level (isolevels 2, 10, 20, 30, 40, 50 mm·h⁻¹) at: (a) *t* = 40 min; (b) *t* = 60 min; and (c) *t* = 80 min. The straight lines indicate the location of the vertical cross sections in Figs. 5, 6 and 7. Only a portion of the computational domain is shown.

A significant amount of charges is observed after 30 min of simulated time, with maximum charge separation rates due to graupel–snow and graupel–crystal bouncing collisions which reach 20 and 2 nC·m⁻³·s⁻¹, respectively. The main charging zone is located in the upper part of the updraft, between 4 and 8 km altitudes where the temperature varies from 0 to –20 °C (across the TCR). A few minutes are then necessary for the electric charges acquired by the different hydrometeors to separate spatially. At 40 min but not shown here, graupel particles predominantly have negative (resp. positive) charges, snow and small ice gain positive (resp. negative) charges, where temperature is colder (resp. warmer) than –10 °C. As charge distribution is dominated by the charge carried by the graupels (–1.5 nC·m⁻³ at high altitude and
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Vertical cross section along the line of Fig. 2(a). (a) Horizontal wind (vector field), isotherms at 273 K, 263 K and 253 K, and vertical velocity (hatched areas); (b) mixing ratios: $r_r$ ([0.1, 1, 2, 3, 4, 5 g·kg⁻¹] solid lines), $r_s$ ([0.1, 0.2, 0.3, 0.4, 0.5, 0.6 g·kg⁻¹] grey scale) and $r_g$ ([0.1, 1, 2, 3, 4, 5 g·kg⁻¹] dashed lines); (c) mixing ratios: $r_c$ ([0.1, 0.5, 1, 1.5, 2 g·kg⁻¹] solid lines) and $r_i$ ([0.1, 0.5, 1, 1.5, 2, 2.5 g·kg⁻¹] dashed lines); (d) total charge $q_{\text{tot}}$ [-2, -1, -0.5, -0.1, 0.1, 0.5, 1, 2 nC].</figcaption></figure>

+2 nC·m⁻³ at mid levels), the early dipolar structure becomes rapidly tripolar (Fig. 5(d)) as more graupels are found at an altitude higher than TCR's altitude. This structure which persists until the storm dissipates, is qualitatively consistent with the analysis of [39] and with the numerical results of [3]. The region of negative charge below this dipole is due to snow. At 60 min (Fig. 6(d)), the charges above 5 km altitude result mainly from those associated with graupel, with weaker contributions from snow and ice crystals. The
<figure><img src="image_10.png" /><figcaption>Figure 6. As in Fig. 5, but the vertical cross section is along the line of Fig. 2(b).</figcaption></figure>

negative zone below is explained by the charging of snow. A major difference with the previous situation (Fig. 5(d)) is the weaker positive pole in the mid-levels and the stronger one below 3 km altitude. This results from the melting of graupels which transform into raindrops. These graupels keep their dominant positive charge, while falling in region of downward motions (see Fig. 5(a)). The charge stratification is clearly visible in the anvil with apparently a large horizontal extension on the $q_{tot}$ plot in Fig. 5(d) (compared to the microphysical fields in Figs. 5(b) and (c)) which is a consequence of the chosen isocontour levels. At 80 min (Fig. 7(d)), the region of positive charges between 5 and 6 km altitude has considerably reduced and is
<figure><img src="image_11.png" /><figcaption>**Figure 7.** As in Fig. 5, but the vertical cross section is along the line of Fig. 2(c).</figcaption></figure>

now restricted to the very upper portions of the remaining weak updrafts (see Fig. 7(a)). In the low levels, relatively important positive charges are still observed and result from the melting of positively charged graupeln aloft. In contrast, the negative charges between 1.5 and 4 km height and carried by rain are now produced by the melting of negatively charged snow. This means that later in the storm history, the polarity of the precipitation at the ground will reverse. In the upper levels, the negative charges are further advected with ice particles in the overhanging anvil.
### 5.2.3. Lightning activity

The rate of simulated lightning discharges is shown in Fig. 8. The displayed lightning frequencies (Fig. 8(a)) and triggering altitudes (Fig. 8(b)) are 1-min averages. Lightning activity begins around 40 min of simulated time with the first intra-cloud (IC) discharges. Cloud-to-ground (CG) flashes occur about 10 min later. The lightning flash frequency increases rapidly after 60 min and remains high until the dissipation of the cloud. The maximum frequency is $0.87 \text{ s}^{-1}$, and the mean frequency during the period of maximum activity (between 60 and 90 min) is $0.6 \text{ s}^{-1}$ approximately. CG flashes are found to be about ten times less numerous than IC ones, with the most intense rate between 65 and 80 min. The simulated flash rates are in a similar range to those currently observed in intense thunderstorms [14]. In Fig. 8(b), the triggering altitude of both IC and CG lightnings tends to decrease with time with a large scatter between 2 and 6 km altitude. The triggering altitude lies between 5 and 8 km before 60 min and spans in the 1 and 6 km range between 60 and 90 min. At the earliest times, lightnings are initiated at the edge of region of positive charges aloft. Later, the presence of a low-level zone of positive charges leads to a more frequent triggering in other parts of the storm. Charge neutralization due to the lightning flashes limits the intensification of the local charges. During the dissipating phase after 80 min, less CG flashes are produced while they are triggered at slightly lower altitudes. Meanwhile, the horizontal extension of IC lightning increases in the anvil.

An illustration of the induced charges in the lightning channels at 60 min of time simulation is given in Fig. 9. In the 3D plot, the trace of IC and CG lightnings can be easily depicted and located with respect to the storm. A remarkable feature of the lightning channels in Fig. 9 is their relatively large horizontal extension which is a consequence of the fractal branching algorithm of the lightning scheme.

### 5.2.4. Discussion

The electrification scheme is applied to simulate the electrical behaviour of a severe storm in MésoNH. With the simple parameterization of [2], the model efficiently generates a multipolar electrical structure of this storm [16,15]. The formation of zones with increasingly positive and negative charges produces an intense enough electric field to trigger lightning flashes.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** (a) Total and CG lightning flash frequency; (b) lightning triggering altitude.</figcaption></figure>
<figure><img src="image_13.png" /></figure>

**Figure 9.** 3D contour of the storm (very light shaded surface) after 1 hour. The light gray shaddied surface near 8 km height is the area where the lightning carries a negative charge density less than $-0.5 \ \text{nC} \cdot \text{m}^{-3}$. The dark gray surfaces show regions where the flashes carry a positive charge higher than $0.5 \ \text{nC} \cdot \text{m}^{-3}$ (note that the vertical scale is exaggerated).

Although the present storm simulation is no more than a coarse simulation of an idealized supercell, many indicators show that many model results are sustainable. The lightning flash activity begins 40 min after the initiation of the convection and holds up almost 40 min. This is typical of a high lightning activity with a maximum flash rate of 0.8 per second (with about 10% of CG) as measured by [40] during the EULINOX field campaign. Also the simulated delay of 10 to 15 min. between total and CG activity has been often reported [39]. The altitude of lightning initiation plotted in Fig. 8(b), occurs between 1 and 8 km and shows that CGs could be triggered at altitude as high as the triggering altitude of ICs. This was reported by [41] and confirmed by others including the recent observations made by [42] who mapped lightning flashes crossing the cloud from 8 to 10 km high before striking the ground. The lightning trajectories (with an example given in Fig. 9) fill the regions of high charge density displaying the multipolar charge structure in the cloud as reported by [43] earlier. Moreover, the new time-of-arrival technique [9] used to map lightning discharges inside storms reveals the case of horizontally expanded discharge ended by a stroke to the ground. Such CGs (such as the one of Fig. 9) initiate at the cloud base. The horizontal part of the discharge joins up with two precipitation shafts of opposite polarity before reaching the ground.

There are, however, some weaker points that certainly need further improvements. The upper part of the cloud remains negatively charged which is not usually observed. This may be a consequence of the abrupt polarity change at TCR when using the separation charge scheme of [2]. Indeed, charge neutralization by lightning, a function of the equivalent surface of the hydrometeors, seems to be more efficient for snow and ice crystal than it is for graupel. As a consequence, charges carried by graupels largely control the storm electrical structure during its lifetime. This point will be clarified next in coming studies with the model.

**6. Conclusions**

The objective of this work is to describe a complete cloud electrification scheme which is part of a 3D mesoscale model with non-hydrostatic dynamics and with explicit microphysics involving five liquid and solid water species. Each cloud particle and hydrometeor follows a size distribution governed by a generalized $\gamma$-law and carries its own charge density with a specific size-charge relationship. The electrification scheme contains several parameterizations of the non-inductive charge separation mechanism due to ice-ice collisions in presence of supercooled cloud water. Charge transfers that result from the numerous microphysical processes occuring in mixed-phase clouds are carefully taken into account by integrating the transfer rates of individual charges as done in the microphysical scheme. The electric field is computed after solving the Poisson's equation for the potential with appropriate boundary conditions. Finally an original lightning scheme is devised to relax the electrical stress when the electric field becomes disruptive locally. In particular, the scheme incorporates recent views on the stochastic behavior of the lightning with a fractal geometry.
In conclusion, the development of a detailed electrification scheme has been accomplished in a non-hydrostatic mesoscale model with other high level capabilities. The evolution of the charge densities attached to each particle type remains by this way closely coupled with the dynamical and microphysical state of the host model. This is a strong argument when simulating with some uncertainties the charging mechanisms and lightning activity in response to the complex internal organization of the storms. Furthermore and because thunderclouds or real mesoscale flows have 3D dynamics by nature, this 3D electrification scheme (with a considerable effort to parallelize the code) suggests many perspectives for future work, e.g., test of more elaborated charge separation parametrizations, deeper analyse of the lightning scheme, and above all simulations at high spatial resolution to study of correlation between precipitation, flash rate, updraft, and graupel echo volume.

**Acknowledgements.** We would like to thank Professor S. Chauzy (Laboratoire d'Aérologie) for helpful discussions, D. Gazen and J. Escobar (Laboratoire d'Aérologie) for their precious advices when coding the parallelized version of the electrical module in the MésoNH code. J. Duron (also from Laboratoire d'Aérologie) is gratefully acknowledge for the graphic support. Gilles Molinié warmly acknowledge CNES (Centre National d'Etudes Spatiales) for a three year financial support. Computer resources on Massively Parallel machines were provided by CINES (Centre Informatique National de l'Enseignement Supérieur, Montpellier, France, under project: c20000122004-CP1) and by CALMIP (CALcul MIdi-Pyrénées, Toulouse, France, under project: P0121).

References

[1] C.L. Ziegler, P.S. Ray, D.R. MacGorman, Relations of kinematics, microphysics and electrification in an isolated mountain thunderstorm, J. Atmos. Sci. 43 (1986) 2098-2115.

[2] J.H. Helsdon Jr., R.D. Farley, A numerical modeling study of a montana thunderstorm. 2. Model results versus observations involving electrical, J. Geoph. Res. 92 (1987) 5661-5675.

[3] C.L. Ziegler, D.R. MacGorman, Observed lightning morphology relative to modeled space charge and electric field distributions in a tornadic storm, J. Atmos. Sci. 51 (1994) 833-851.

[4] S.A. Randell, S.C. Rutledge, R.D. Farley, J.H. Helsdon Jr., A modeling study on the early electrical development of tropical convection: Continental and oceanic (monsoon) storms, Mon. Wea. Rev. 122 (1994) 1852-1877.

[5] E.R. Mansell, Electrification and lightning in simulated supercell and non-supercell thunderstorms, Ph.D. thesis, Univ. Oklahoma, 2000.

[6] M.B. Baker, A.M. Blyth, H.J. Christian, J. Latham, K.L. Miller, A.M. Gadian, Relationships between lightning activity and various thundercloud parameters: satellite and modelling studies, J. Atmos. Res. 51 (1999) 221-236.

[7] H.J. Christian, R.J. Blakeslee, S.J. Goodman, D.A. Mach, M.F. Stewart, D.E. Buechler, W.J. Koshak, J.M. Hall, W.L. Boeck, K.T. Driscoll, D.J. Bocippio, The lightning imaging sensor, in: Int. Conf. on Atmospheric Electricity, Guntersville, AL, Amer. Meteor. Society, 1999, pp. 746-749.

[8] A.R. Jacobson, S.O. Knox, R. Franz, D.C. Enemark, Forte observations of lightning radio-frequency signatures: Capabilities and basic results, Radio Sci. (1999).

[9] W. Rison, P.R. Thomas, R.J. Krehbiel, T. Hamlin, J. Harlin, A gps-based three-dimensional lightning mapping system: Initial observations, Geoph. Res. Lett. 26 (1999) 3573-3576.

[10] E. Defer, P. Blanchet, C. Théry, P. Laroche, J.E. Dye, M. Venticinque, K.L. Cummins, Lightning activity for the July 10, 1996, storm during the stratosphere-troposhere experiment: Radiation, aerosol, and ozone-a (STERAO-A) experiment, J. Geoph. Res. 106 (2001) 10151-10172.

[11] J.P. Lafore, J. Stein, N. Asencio, P. Bougeault, V. Ducrocq, J. Duron, C. Fischer, P. Hereil, P. Mascart, J.P. Pinty, J.L. Redelsperger, E. Richard, J. Vila-Guerau de Arellano, The meso-NH atmospheric simulation system. Part I: Adiabatic formulation and control simulations, Ann. Geoph. 16 (1998) 90-109.

[12] D.R. Mac Gorman, J.M. Straka, C.L. Ziegler, A lightning parameterization for numerical cloud models, J. Appl. Meteor. 40 (2001) 459-478.

[13] J. Stein, E. Richard, J.-P. Lafore, J.-P. Pinty, N. Asencio, S. Cosma, High-resolution non-hydrostatic simulations of flash-flood episodes with grid-nesting and ice-phase parameterization, Meteorol. Atmos. Phys. 72 (2000) 203-221.

[14] D.R. Mc Gorman, W.D. Rust, The Electrical Nature of Storms, Oxford University Press, 1998.

[15] M. Stolzenburg, T.C. Marshall, W.D. Rust, B.F. Smull, Horizontal distribution of electrical and meteorological conditions across the stratiform region of a mesoscale convective system, Mon. Wea. Rev. 122 (1994) 1777-1797.

[16] E.R. Williams, The tripole structure of thunderstorm, J. Geoph. Res. 94 (1989) 13151-13167.

[17] R. Solomon, V. Schroeder, M.B. Baker, Lightning initiation-conventional and runaway-breakdown hypotheses, Quart. J. Roy. Meteor. Soc. 127 (2001) 2683-2704.
[18] T.C. Marshall, M. Mc Carthy, W.D. Rust, Electric field magnitudes and lightning initiation in thunderstorms, J. Geoph. Res. 100 (1995) 7097-7103.

[19] V. Mazur, L.H. Ruhnke, Model of electric charges in thunderstorms and associated lightning, J. Geoph. Res. 103 (1998) 23299-23308.

[20] Y.-L. Lin, R.D. Farley, H.D. Orville, Bulk parameterization of the snow field in a cloud model, J. Clim. Appl. Meteor. 22 (1983) 1065-1092.

[21] E. Kessler, On the distribution and continuity of water substance in atmospheric circulation, in: Meteor. Monogr., No 32, Amer. Meteor. Society, 1969.

[22] G. Caniaux, J.L. Redelsperger, J.P. Lafore, A numerical study of the stratiform region of a fast-moving squall line. Part I: General description and water and heat budgets, J. Atmos. Sci. 51 (1994) 2046-2074.

[23] B.S. Ferrier, W.-K. Tao, J. Simpson, A double-moment multiple phase four-class bulk ice scheme. Part II: Simulations of convective storms in different large-scale environments and comparisons with other bulk parameterizations, J. Atmos. Sci. 52 (1995) 1001-1033.

[24] K.V. Beard, H.R. Ochs, Charging mechanisms in clouds and thunderstorms, in: The Earth's Electrical Environment, Studies in Geophysics, Nat. Academic Press, Washington, DC, 1986, pp. 114-130.

[25] T. Takahashi, Riming electrification as a charge generation mechanism in thunderstorms, J. Atmos. Sci. 35 (1978) 1536-1548.

[26] E.R. Jayaratne, C.P.R. Saunders, J. Hallet, Laboratory studies of the charging of soft-hail during ice crystal interactions, Quart. J. Roy. Meteor. Soc. 103 (1983) 609-630.

[27] E.E. Avila, G.G. Aguirre Varela, G.M. Caranti, Temperature dependence of static charging in ice growing by riming, J. Atmos. Sci. 52 (1995) 4515-4520.

[28] E.E. Avila, G.M. Caranti, N. Castellano, C.P.R. Saunders, Laboratory studies of the influence of cloud droplet size on charge transfer during crystal-graupel collisions, J. Geoph. Res. 103 (1998) 8985-8996.

[29] J.H. Helsdon Jr., W.A. Wojcik, R.D. Farley, An examination of thunderstorm-charging mechanisms using a two-dimensional storm electrification model, J. Geoph. Res. 106 (2001) 1165-1192.

[30] B. Gardiner, D. Lamb, R.L. Pitter, J. Hallet, Measurements of initial potential gradient and particle charges in a Montana summer thunderstorm, J. Geoph. Res. 90 (1985) 6079-6086.

[31] C.P.R. Saunders, W.D. Keith, R.P. Mitzeva, The effect of liquid water on thunderstorm charging, J. Geoph. Res. (1991).

[32] C.P.R. Saunders, S.L. Peck, Laboratory studies of the influence of the rime accretion rate on charge transfer during crystal/grauplen collisions, J. Geoph. Res. 103 (1998) 13949-13956.

[33] B.S. Ferrier, A double-moment multiple phase four-class bulk ice scheme. Part I: description, J. Atmos. Sci. 51 (1994) 249-280.

[34] Reiter, Precipitation and cloud electricity, Quart. J. Roy. Meteor. Soc. 91 (1965) 60-72.

[35] H.W. Kasemir, A contribution to the electrostatic theory of a lightning discharge, J. Geoph. Res. 65 (1960) 1873-1878.

[36] L. Niemeyer, L. Pietronero, H.J. Weismann, Fractal dimension of dielectric breakdown, Phys. Rev. Lett. (1984).

[37] N.I. Petrov, G.N. Petrova, Physical mechanisms for intracloud lightning discharges, Technical Phys. 44 (1993) 472-475.

[38] J.B. Klemp, R.B. Wilhelmson, The simulation of three-dimensional convective storm dynamics, J. Atmos. Sci. 35 (1978) 1070-1096.

[39] E.R. Williams, M.E. Weber, R.E. Orville, The relationship between lightning type and convective state of thunderclouds, J. Geoph. Res. 94 (1989) 13213-13220.

[40] N. Dotzek, H. Höller, C. Thèry, T. Ferh, Lightning evolution related to radar-derived microphysics in the 21 July 1998 eulinox supercell storm, J. Atmos. Res. (2001).

[41] D.E. Proctor, Regions where lightning flashes began, J. Geoph. Res. 96 (1991) 5099-5112.

[42] R. Mardiana, Z.I. Kawasaki, T. Morimoto, Three-dimensional lightning observations of cloud-to-ground flashes using broadband interferometers, J. Atmosph. Solar-Terr. Phys. (2002).

[43] X.M. Shao, P.R. Krehbiel, The spatial and temporal development of intracloud lightning, J. Geoph. Res. 101 (1996) 26641-26668.