CELLS v 1 . 0 : updated and parallelized version of an electrical scheme to simulate multiple electrified clouds and flashes over large domains

The paper describes the fully parallelized electrical scheme CELLS which is suitable to simulate explicitly electrified storm systems on parallel computers. Our motivation here is to show that a cloud electricity scheme can be developed for use on large grids with complex terrain. Large computational domains are needed to perform real case meteorological simulations with many independent convective cells. The scheme computes the bulk electric charge attached to each cloud particle and hydrometeor. Positive and negative ions are also taken into account. Several parametrizations of the dominant non-inductive charging process are included and an inductive charging process as well. The electric field is obtained by inverting the Gauss equation with an extension to terrain-following coordinates. The new feature concerns the lightning flash scheme which is a simplified version of an older detailed sequential scheme. Flashes are composed of a bidirectional leader phase (vertical extension from the triggering point) and a phase obeying a fractal law (with horizontal extension on electrically charged zones). The originality of the scheme lies in the way the branching phase is treated to get a parallel code. The complete electrification scheme is tested for the 10 July 1996 STERAO case and for the 21 July 1998 EULINOX case. Flash characteristics are analysed in detail and additional sensitivity experiments are performed for the STERAO case. Although the simulations were run for flat terrain conditions, they show that the model behaves well on multiprocessor computers. This opens a wide area of application for this electrical scheme with the next objective of running real meterological case on large domains.


Introduction
The ground detection of the electrical activity inside convective systems revealed the strong links with the dynamics (Goodman et al., 1988;Wiens et al., 2005), the cloud microphysics and even the atmospheric chemistry through the formation of nitrogen monoxide, an ozone precursor (Schumann and Huntrieser, 2007).As a result, cloud discharges were related to the presence of precipitating ice in deep clouds (Blyth et al., 2001;Petersen et al., 2005;Prigent et al., 2005;Deierling et al., 2008;Barthe et al., 2010), to the intensification of tropical cyclones (Cecil and Zipser, 1999;Squires and Businger, 2008;Price et al., 2009) and to useful nowcasting index of severe hail-bearing storms (Darden et al., 2010;Emersic et al., 2011).
In contrast, modeling the electrical activity of a storm is still a very difficult task owing to the large number of physical mechanisms to represent and to the poor knowledge and parameterization of basic processes.To reproduce the electric charge cycle in a thunderstorm, the following issues must be considered: a micro-scale charge separation mechanism, the transfer and the transport of the electric charges according to the evolution of the hydrometeors at cloud-scale, the computation of the electric field, the propagation of the lightning flashes and a partial neutralization of the charges.
Most of the modeling studies in cloud electricity have concentrated on the microphysical aspects of the charge cycle or attempted to reproduce the bulk effect of lightning flashes.Altaratz et al. (2005) introduced an electrification scheme in the Regional Atmospheric Modeling System (RAMS) model, but in absence of lightning scheme, they restricted Published by Copernicus Publications on behalf of the European Geosciences Union.
their study to the charging processes until the first lightning flash.In the same way, Hou et al. (2009) limited the study of the charge structure during the pre-lightning stage of five thunderstorms since their model solely integrated an electrification scheme.Rawlins (1982) was the first to introduce a lightning parameterization in a numerical cloud model.This simple scheme reduced arbitrarily the net charge density at each grid point when the electric field exceeded a threshold.Takahashi (1987) and Ziegler and MacGorman (1994) also treated the bulk effect of lightning flashes in a 2-D axisymmetric model, and in a 3-D model, respectively.Solomon and Baker (1996) developed a complete electrical scheme with individual lightning flash treatment, but in a simplified 1.5dimension kinematic cloud model which prevented real case studies.On the other side, Mazur and Ruhnke (1998) and Riousset et al. (2007) discussed on the propagation of cloud discharges, but with an idealized charge distribution, and did not integrate their lightning scheme in a cloud model.
Today, only a few models attempt to simulate the structure and the evolution of the electric charges in a thunderstorm.The Storm Electrification Model with an electrification and lightning flash scheme was pioneered by Helsdon and Farley (1987) and Helsdon et al. (1992).It was used to investigate the charge structure and the Maxwell currents in an idealized storm (Helsdon et al., 2001) and then to study the lightning-produced NO x (Zhang et al., 2003a,b).Sun et al. (2002) adopted the electrical scheme of Helsdon et al. (1992) to simulate the feedbacks of cloud electricity on convection.This was done by adding the three components of the electric force acting on all the charges to the momentum equation of their host model.Nevertheless Sun et al. (2002) found a strenghening of the convection in their thunderstorm case study but the validity of this result is also questionable since the flash scheme they used was too much simplified.
The first complete electrification scheme coupled to a realistic but expensive lightning flash scheme, with leader and branches, was developed by Mansell et al. (2002).It was widely used to study the sensitivity of the lightning activity to the non-inductive parameterization (Mansell et al., 2005), and to analyze the lightning activity in a idealized tropical cyclone (Fierro et al., 2007) and in the 29 June 2000 Severe Thunderstorm Electrification and Precipitation Study (STEPS) storm case (Kuhlman et al., 2006).Independently but in the same vein, Barthe et al. (2005) developed another electrical scheme.They introduced an original fractal approach in their lightning scheme which was developped in the framework of the french mesoscale model Meso-NH.The fractal law was introduced to estimate the degree of branching of the discharge when expanding from the bidirectional leader.This leads to a complex code with a probabilistic search of new lightning segments to add to the growing lightning structure.The model was used to investigate the sensitivity of the charge structure to the noninductive parameterizations and to the sensitivity of lightning flash parameters non amenable to direct observations (Barthe and Pinty, 2007a).The robustness of the full electrical scheme was then demonstrated by Pinty and Barthe (2008).Besides, a first direct modeling of the production of nitrogen oxides by lightning flashes was realized for the 10 July 1996 STERAO storm case (Barthe et al., 2007;Barth et al., 2007).
Until now and despite their success to simulate isolated electrified storms, a number of difficulties prevented the last two models (Mansell et al., 2005;Barthe et al., 2005) from being used over large computational domains or for real meteorological applications.There are several reasons for that.First, the commonly shared view of sequential and stepwise propagation of the flashes makes the lightning path algorithm not well adapted to massively parallel computing.It is a difficult task to parallelize and to check a lightning flash algorithm in the context of domain decomposition but even so, an acceptable multiprocessor computing efficiency cannot be achieved as long as the spatial growth of a branched structure is based on an iterative process.Second, several isolated cells can trigger flashes in convective systems during a single time step.Consequently, the lightning flash scheme must apply to all cells at once or needs to be repeated in a determined order to explore carefully each of the electrified convective cells present in the domain of simulation.Finally, one can expect numerical difficulties linked to the distortion of the curvilinear vertical coordinate due to orography in real case studies.This problem arises when computing the electric field but solutions exist to invert the key elliptical equation with extra metric terms (see below).In addition, one can expect also serious complications due to terrain-following coordinates in the description of the filamentary structure of the flashes in case of uneven locations of the grid points.
However and in the context of the next Hydrological cycle in Mediterranean Experiment (HyMeX) (http://www.hymex.org/)during which several lightning sensors will be deployed, it is intended to perform Meso-NH simulations of three-dimensional (3-D) electrified cloud systems on a very large computational domain at kilometer scale resolution with the grid-nesting technique to downscale the meteorological analyses.To this aim, the Meso-NH lightning scheme must be revised while keeping as realistic as possible the electrical behavior of the flashes, mostly the horizontal and vertical extensions of the intra-cloud (IC) discharges, and the quantity of neutralized charge per flash.
This paper describes the new lightning flash scheme developed in the cloud-resolving model Meso-NH with the perspective of running real case electrified storms over complex terrain at high resolution on multiprocessor computers.The electrical scheme labelled CELLS for "Cloud ELectrification and Lightning Scheme" is detailed in Sect. 2 with a focus on the lightning flash scheme.Section 3 is dedicated to a first evaluation of the new lightning scheme on the 10 July 1996 STERAO (Stratospheric-Tropospheric Experiment: Radiation, Aerosols and Ozone) storm with a sensitivity study of the lightning flash parameters.Finally, the 2 Description of the electrical scheme 2.1 The cloud electrification scheme in Meso-NH model

Generalities
The Meso-NH model (Lafore et al., 1998) is able to simulate idealized precipitating systems at high resolution and real meteorological events on large domains with complex terrain.In the later case, Meso-NH needs meteorological analyses for the initialization and the open boundary conditions while high resolution, typically the kilometer scale, is achieved automatically via the grid nesting facility.Since the code is fully vectorized and efficiently parallelized (Jabouille et al., 1999), the 3-D evolution of any cloud system is currently simulated on large grids with hundreds of points in each horizontal direction.
The cloud electrification scheme of Meso-NH has been already described in Barthe et al. (2005) and Barthe and Pinty (2007b).However, due to the sequential algorithm of the flash scheme and to the numerical cost induced by the frequent communications between processors, simulations of electrified storms in Meso-NH were mostly performed on a single processor.
In the scheme, the mass charge densities (q x in C kg −1 of dry air) are the bulk prognostic electrical state variables to fit with the conservation law of the scalar fields in Meso-NH.They are closely related to the mixing ratio (r x in kg kg −1 ) of the microphysical species x (cloud droplets, rain, pristine ice crystals, snow/aggregates, graupel and hail).For instance and similarly to the mass of individual particles, a chargeparticle size power law relationship is assumed as explained in Barthe et al. (2005).The bulk charges q x are evolving according to: where U is the 3-D air velocity and ρ dref a fixed, dry air density reference state (Meso-NH integrates an anelastic system of equation).The source terms S q x include the turbulent diffusion, the charging mechanism rates, the charge sedimentation by gravity and the charge neutralization by the lightning flashes.The transfer rates due to the microphysical evolution of the particles are collected in T q x .Each microphysical process T r x is associated to an electrical tranfer rate in proportion of the mixing ratio and electric charge, i.e.T q x = (q x /r x ) × T r x where T r x is provided by the microphysical scheme.

Charge separation mechanisms
Even if the physical explanations are still unclear, laboratory studies (Takahashi, 1978;Jayaratne et al., 1983;Saunders et al., 1991;Avila et al., 1995;Saunders and Peck, 1998, among others) show indeed that the non-inductive (NI) charging mechanism after rebounding collisions between small unrimed and big rimed ice particles is likely to be the dominant process for charge separation which must be considered at first.Four different parameterizations of the noninductive mechanism are available in Meso-NH.They result from the published work of Takahashi (1978), Gardiner et al. (1985), Saunders et al. (1991) and Saunders and Peck (1998).For each colliding event, the polarity and the quantity of separated charge is given as a function of the temperature and the liquid water content or riming rate.This concerns only three types of collision: pristine ice-snow, pristine icegraupel and snow-graupel.Hail is not efficient to generate electric charges in Meso-NH because these particles are supposedly wrapped by a film of water (Saunders and Brooks, 1992).The analytical expressions of the charging rates relies heavily on the microphysical scheme: with D x and D y the diameter for species x and y, respectively.|V x − V y | is the relative fall speed, n x and n y are the number concentrations of species x and y, respectively, and E xy is the collection efficiency.The collection efficiency depends on the temperature and follows Kajikawa and Heymsfield (1989) for ice-snow and snow-graupel collisions, and Mansell et al. (2005) for ice-graupel collisions.
As in Mansell et al. (2005), the charge exchanged per rebounding collision δq is limited to prevent unreasonable charging rate.Based on Keith and Saunders (1990), it is assumed that the charging rate of the pristine ice crystal with D max ∼ 100 µm is the most limiting one, that is 30(10) fC per collision with graupel (aggregate) particles.We take a larger value (100 fC) for the graupel-snow collisions because it corresponds roughly to an average of the saturation levels when the particle sizes reach ∼1 mm (see Keith andSaunders, 1990 or Fig. 3.13 in MacGorman andRust, 1998).This limitation is introduced in the computation of the bulk charging rates which result from an integration over the size spectrum of the ice particles.
The inductive charging that results from graupel-droplets collisions in a preexisting electric field is also taken into account following Ziegler et al. (1991).

Small ions
In order to close carefully the electric charge budget when the cloud particles and hydrometeors (the main electric charge carriers) evaporate and to simulate the screen charges, it is necessary to integrate two conservation equations for the positive (n + ) and for the negative (n − ) ion concentrations (Helsdon and Farley, 1987).Assuming that all the ions have an elementary charge, the condensed form of the ion governing equation writes: where S att , S evap , S light , and S pd are source/sink terms corresponding to ion attachment to charged hydrometeors (sink), release of ion when hydrometeors evaporate (source), production by lightning flash and by point discharge current from the surface, respectively.The term G is the ion generation rate by cosmic rays, and α is the ion-ion recombination coefficient.The first three terms are given in this order: ion transport by the mean flow, electrostatic drift motion with parameterized mobilities µ ± and turbulent mixing with K the eddy diffusivity.The ion attachment is a complex term with a combination of free ion diffusion to the particle surface (an electrical attraction due to the presence of a net particle charge) and ion conduction (due to ion motion in the presence of an electric field).The analytical case-dependent expressions of the ion attachement were first given by Chiu (1978) and then were adopted by Helsdon et al. (2002) and Mansell et al. (2005).Fair weather conditions for the mean current density and for the vertical decrease of the electric field profile, are used to initialize the positive and negative free ion concentrations as proposed by Helsdon and Farley (1987).Then assuming steady state conditions, the intensity of the constant cosmic ray source, G, can be estimated from a balance involving the ion drift and ion recombinaison (see also MacGorman and Rust, 1998).The fair weather ion concentrations are used to treat inflow conditions on the lateral boundaries during the model integration.Furthermore, because the downward drift motion enables the ions to cross the top of the domain, it is necessary to relax the ion concentrations to their fair weather value in upper levels to avoid their accumulation.
The ion generation source G is height dependent as in the previous studies, but it should reflect the changes in ionization intensity along the solar cycle.In the following two case studies, the same profile is used since the events occurred at the late afternoon (in local time) over a short period.However, it is probably more consistent to use time-variable height profiles if the convective systems have a longer lifetime.

Electric field computation
The electric field (E) is diagnosed each time step and after the charge rearrangement following a flash.E is solution of the Gauss equation forced by the total charge volume density with a = 8.85 pF m −1 , the permittivity of air and |e| = 1.602 × 10 −19 C, the elementary electric charge.In order to compute E, it is useful to introduce a pseudo electric potential V such as ρ∇V = −E so that a diagnostic "pressure equation" analog in Meso-NH (Lafore et al., 1998) is recovered: GDIV is the generalized divergence operator in the nonorthogonal curvilinear coordinates system, and ρ = ρ dref × J is the mass of dry air (J , the Jacobian of the coordinate transform, corresponds to the local gridbox volume).As a result, the optimized elliptic standard pressure solver of Meso-NH can be employed to get V with Neumann boundary conditions in Eq. ( 5).Finally the electric field E is derived by applying a gradient operator on the V field.

Lightning flash scheme
The objective of the new lightning flash scheme is to reproduce some morphological characteristics of the lightning flashes as in Barthe and Pinty (2007b), but for electrified storms growing over large grids and complex terrain.This is achieved by simplifying the original algorithm to get a parallel code as explained below.Details about the parallelization of Meso-NH are available from http://mesonh.aero.obs-mip.fr/mesonh/.Good efficiency of the parallelization is provided by a library of high level functions which greatly helps the coding for scientific end users.Because the nature of most of the calculations involves only a local knowledge of the global 3-D fields (with storage on distributed memory), each processor can easily work independently on its side.In our case however, building the filamentary structure of a lightning flash path is leading ineluctably to frequent communications between processors which must be optimized.
In the following, the variables suffixed by ll refer to global variables with a single updated value available to all processors.It is hypothesized that the domain is divided into N proc subdomains, with N proc being the number of working processors.
The different steps of the lightning flash scheme are sketched in Fig. 1 and described in this section.

Electrified cells identification
Previous simulations of storm lightning activity were performed in an idealized framework, over a limited domain The cross and the + show the grid points where the flash is initiated and where the leader propagates, respectively.(d) Regions of possible branching: the blue (red) contours limit the regions with negative (positive) total charge density where the negative (positive) branches can propagate.(e) Choice of the branches: the grey points represent the grid points with possible branching between two successive circles (or spheres in 3-D).N ll(r) is the maximum number of branches between two circles of radius r and r + dr (see Eq. 7).Here, the maximum number of branches is given for a 2-D framework and a mean grid size of 1450 m.(f) Lightning flash: the resulting lightning flash path made up with the triggering point (black cross), the bidirectional leader (black +) and the branches (grey +).
area and for a rather short duration time.Therefore, only a few electrified cells were present in the simulation domain at the same time, but the goal now is to treat the individual flashes of several cells simultaneously.
An iterative algorithm is first developed to identify all the electrified cells in the domain of simulation.In the following, a cell is termed "electrified" if conditions to trigger and to propagate a flash inside it are fulfilled.The electric field module is first multiplied by a factor to get free of the height effect.It is noted E 0 and corresponds to the electric field module reduced to the ground level.The peak value of E 0 , E max , is sought in each subdomain.Then, the global value of the maximum electric field E max ll = MAX(E max ) is determined and the processor number (IP cell ) where E max ll = E max is identified.If E max ll is higher than 200 kV m −1 , the electric field threshold for flash triggering at the ground level (see Sect. 2.2.2), a first electrified cell is detected.The maximum electric field is a natural marker of lightning-triggering cells since a flash is triggered only if E max > 200 kV m −1 .The point where E max ll > 200 kV m −1 is hereafter called the cell center.The local coordinates of the cell center and IP cell are then broadcast to all processors.
The next step explores the vertical and horizontal extensions of the selected cell (Fig. 1a).The domain volume is scanned from the bottom to the top.The cell center is www.geosci-model-dev.net/5/167/2012/Geosci.Model Dev., 5, 167-184, 2012 projected onto the horizontal plane of the running level and contiguous grid points are tagged if they meet the following conditions: r tot > 1 × 10 −5 kg kg −1 to restrict a flash propagation to a single cloud, -at least one hydrometeor category has | qx | > qcell where qcell is a given threshold to isolate individual storm cells.qx is the volume charge density (C m −3 ) for species x.The process is repeated along the horizontal until no more grid points can be added to the cell volume.Updates in the halo zones (in a parallel architecture, a "halo" zone contains the overlapping grid-points which are exchanged with the neighbor processors) are necessary because electrified cells may span over several neighboring subdomains.Then the algorithm loops to analyse the electric field out of the electrified cell to find out if another disjuncted electrified cell exists in the whole domain.

Flash triggering
The local electric field condition which initiates a flash, follows MacGorman et al. (2001) and Barthe and Pinty (2007b).The triggering electric field, E trig decreases with altitude as observed by Marshall et al. (1995): where z is the altitude (km) and E trig is given in kV m −1 .To account crudely for grid scale uncertainty, a flash is triggered where the electric field exceeds a slightly smaller value than E trig (such as k × E trig , with k = 0.9).If more than one grid point per convective cell meets the condition E > 0.9×E trig , then the triggering point is chosen at random (Fig. 1b).
The processor IP trig containing this point is identified.The value of the triggering electric field, the coordinates of the flash origin and the sign of the vertical component of the electric field at this point are broadcast from IP trig to all processors.This procedure is repeated for each cell.Then, if several cells exist in the domain several flashes can be treated simultaneously since there is a mask that discriminate the different cells in the domain.
Once the characteristics (center and extension) of all electrified cells are available, the lightning flash stage follows.The treatment of the flashes is broken down into two parts with a "leader" phase that precedes a phase that generates the branches.

Bidirectional leader
The approach follows Helsdon et al. (1992) that relies on the bidirectional leader theory of Kasemir (1960).Kasemir assumes that the flash leader propagates bi-directionally from the triggering point, in the parallel and anti-parallel directions of the ambient electric field.The propagation is stopped once the electric field drops below a threshold value.As previously done by Helsdon et al. (1992), Barthe and Pinty (2007b) simplified this concept since they used the ambient electric field to control the leader propagation instead of the total electric field.They acknowledged that it was a shortcoming, but argued that computing the local electric field at the tip of each segment added to the leader was computationally expensive.In the present scheme, a new simplification is considered, still for a sake of reducing the computational cost.The bidirectional leader is allowed to propagate along the vertical axis only and not slantwise along E as in the previous scheme, to avoid communication between the processors each time a new segment is added at the tip of the leader.The two branches of the leader propagate until the ambient vertical electric field (E z ) at the tips of the last segments falls below ∼15 kV m −1 or when the sign of the vertical component of the electric field reverses (Fig. 1c).
As in other studies (MacGorman et al., 2001;Mansell et al., 2002;Barthe et al., 2005;Mansell et al., 2005Mansell et al., , 2010)), a flash is categorized as "cloud-to-ground" (CG) when the lower end of the leader reaches the bottom of the cell which altitude is below 2 km above ground level (AGL).CG flashes are artificially prolonged to the ground.
Only processor IP cell is in charge of the bidirectional leader.The coordinates of the leader channel and the flash type are broadcast to all processors.

Horizontal extension of the flash
VHF mapping systems have highlighted the extensive horizontal structure of lightning flashes in two distinct layers (Shao and Krehbiel, 1996;Rison et al., 1999;Thomas et al., 2001;Wiens et al., 2005;Bruning et al., 2007), with a single vertical channel connecting the two layers.Therefore, in this context, the new lightning flash scheme must reproduce this feature but in an economical way, since a physically consistent representation of the discharges is too expensive and would be technically impraticable on powerful massively parallel computers.
According to VHF observations, a positively and a negatively charged region must be delineated (propagation is not allowed in a third region in case of a tripole structure of charges).The positively and negatively charged regions where the flash can propagate are explored separately.From the positive part of the leader, the region with negative total charge density where the positive branches can propagate is explored.First, a 3-D mask M(:,:,:) is initialized with value equal to 1 where grid points are reached by the positive leader.Then the neighboring grid points matching the following conditions are selected as part of the negative charge pocket: -the grid point belongs to the electrified cell of the leader -the total charge density must be negative The fields in the halo zones are updated to allow a continuous extension of the pocket of negative charge to nearby subdomains.This step is resumed until no more point can be added.The positive charge pocket is build the same way around the negative leader, leading to two regions of opposite charge embedded in the electrified cell contour (blue and red contours in Fig. 1d).Williams et al. (1985) initiated discharges through plastic slabs with regions of stronger and weaker negative charge density.They observed that the discharges tend to propagate toward regions of high charge density, which underlined the importance of the charge density in discharge propagation.Niemeyer et al. (1984) showed that a stochastic dielectric breakdown model naturally leads to a fractal structure of the discharge.Tsonis and Elsner (1987) analyzed a set of lightning pictures and deduced an average fractal dimension of the lightning projection.The dielectric breakdown model has been widely used in the past to simulate different types of lightning discharges (Wiesmann and Zeller, 1986;Wiesmann, 1988;Petrov and Petrova, 1993;Sañudo et al., 1995;Kawasaki and Matsuura, 2000, among others).However, none of these studies simulated lightning discharges in a real storm context.Mansell et al. (2002) first introduced the dielectric breakdown concept to simulate the lightning flashes in a cloud model.Then Barthe and Pinty (2007b) developed a probabilistic branching algorithm adapted from the dielectric breakdown concept to mimic the horizontal extension of the flash toward regions of high charge density.The present scheme keeps the idea of charge density criterion to build a 3-D branched discharge (Williams et al., 1985) and to monitor the fractal nature of the flash (Niemeyer et al., 1984) as previously highlighted.

Distribution of the branches
The electrified cell domain is divided into concentric spheres with radius r centered on the triggering point (Fig. 1e).The global number of branches N ll at a distance r from the triggering point is assumed to follow a fractal law (Niemeyer et al., 1984): with L mean , the mean mesh size (m), L χ , a characteristic length scale (m), and χ, the fractal dimension (2 < χ < 3 according to Petrov and Petrova (1993)).The running integer i, computed as i = NINT(r/L mean ), varies from i min = 0 to i max where i max corresponds to the maximum distance where branching is possible, i.e. for gridpoints belonging to mask M. NINT is a function returning the nearest integer of a real number.So a local array A that contains all the nearest integer distances i between the triggering point and each grid point passing mask M(:,:,:) = 1 is filled in each subdomain.The minimum (i min ) and maximum (i max ) distances are checked so that the next steps are iterated for i min ≤ i ≤ i max .
On each processor, the number of grid points belonging to mask M and located at distance i (N poss (i)) is computed and summed over all processors to get the global number of possible locations N poss ll(i).Then N poss ll(i) is compared to N ll(i) of Eq. ( 7): -if N poss ll(i) ≤ N ll(i): all the possible grid points at distance i from the triggering point are selected, -if N poss ll(i) > N ll(i): too much grid points are found so a subset must be selected at random.
In order to randomize properly the selection of the grid points which are dispersed on several subdomains, two 1-D working arrays V(:) and V (:) of size N poss ll(i) are allocated to each subdomain.Each processor packs the 3-D array A into a 1-D array V under a running mask control defined by A(:,:,:) = i.V is initially set to 0.
The number of grid points, N poss (i), at a given distance i is gathered by each processor and the result is broadcast to all processors.Consequently, each of the IP flash processor knows the proportion of grid points which is granted in its physical subdomain since the total lightning flash path should contain at most N ll(i) grid points.A random integer n is then taken in the interval 0 and N poss ll(i).The grid point number V(n) is extracted and element V (V(n)) is turned to 1.This step is repeated until N ll(i) points are chosen.Finally, the elements of V are scattered (equivalent to an unpack operation) in a 3-D array S(:,:,:) under the same mask condition A(:,:,:) = i.As a result, the sparse array S(:,:,:) obtained at the last iteration i = i max , marks the full path of the lightning flash (Fig. 1f).The branches coordinates are then broadcast to all processors.
Compared to the original scheme of Barthe et al. (2005), the new lightning scheme ignores deliberately any rules of connectivity that should apply to the grid points tracing the flash path.The higher randomization and consequently the lack of structuration of the branching phase is the price to pay to get a parallel source code even if it is done essentially for technical reasons.We believe however that running simulations of electrified storms on large grid systems, and so accessing to parallel computing, is very promising for real meteorological applications of cloud electricity.

Neutralization
The total charge in excess of qneut along the lightning channel is neutralized in the lightning flash.It is distributed to the ions of opposite sign at each grid point of the flash.For intracloud flashes, a charge correction is applied to all grid points of the flash to ensure the total charge neutrality (MacGorman et al., 2001) before the redistribution of net charge at grid points.For cloud-to-ground discharges, the charge neutrality constraint does not apply.Once the charge is neutralized, the electric field is updated.If at least one new triggering point is found in the domain, the procedure (from Sect.2.2.3 to Sect.2.2.6) is repeated.Thus, in a single time step, each cell can generate several flashes.

Technical aspects
The whole scheme is embedded in the cloud-resolving mesocale model Meso-NH in which the resolved-scale and turbulent transport terms (Eqs. 1 and 3) are computed.The code is written in Fortran 90 and parallelization is obtained by using MPI functions.The bulk electrical charge scheme is developped in the framework of the single moment microphysical scheme of Meso-NH to take advantage of calculating the mass microphysical rates (to get the transfer rates T q x in Eq. 1).The inversion of Eq. ( 5) to extract V is done by direct and inverse FFT on the horizontal and by resolving a tridiagonal system along the vertical.The computation of E is exact in absence of orography.In the other case (orography is present), a conjugate gradient method or equivalent which implies an iterative loop, is necessary to include cross derivative terms.The lightning flash scheme excepted, the computation of the charges and of the electric field are not easily adaptable to other models.CELLS v1.0 is then available in the release of Meso-NH v4.9.The use of the Meso-NH model by groups other than the developers is subject to a licence agreement.Additional details on the Meso-NH model appear in http://mesonh.aero.obs-mip.fr/mesonh/with technical and scientific documentation.

Simulation of the lightning activity during the 10 July 1996 STERAO storm
A first evaluation of the new flash algorithm is undertaken for the well-known idealized test case of the 10 July 1996 STERAO-A storm that occurred along the Wyoming-Nebraska-Colorado border (Dye et al., 2000).The storm started to develop as a multicell and then moves to a supercell ∼90 min later.In addition to the dynamical aspects with a multi-to supercell transition to capture, this storm case is interesting for two reasons at least.The presence of several convective cells in the domain of simulation allows testing of the cell identification algorithm and a unique set of data (Defer et al., 2001(Defer et al., , 2003)), including IC and CG flash characteristics, is an opportunity to tune and to evaluate the sensitivity of some critical flash parameters.

Numerical set-up
The numerical simulation was performed with the nonhydrostatic mesoscale model Meso-NH version 4.8.4.The environment was assumed to be homogeneous, then a single profile was used for initialization (Skamarock et al., 2000).
As in previous numerical studies of this storm (Skamarock et al., 2000(Skamarock et al., , 2003;;Barth et al., 2001Barth et al., , 2007;;Barthe et al., 2007;Barthe and Barth, 2008), three warm bubbles (+3 • C) were oriented in a north-west to south-east line (Fig. 2a).A 160 × 160 km 2 horizontal domain was used with a 1-km resolution.The vertical grid had 51 levels up to 23 km with a level spacing of 50 m close to the surface stretching to 700 m at the top of the domain.The time step was 2.5 s, and the simulation lasted 3 h.
The physics of the model included a mixed-phase microphysics scheme (Pinty and Jabouille, 1998) and a 3-D turbulence scheme (Cuxart et al., 2000).The parameterization of Takahashi (1978) is used to describe the NI processes as in the previous simulation of this storm by Barthe et al. (2007).To prevent unreasonably large charging and flash rates, the magnitude of charge separated per rebounding collision is limited to 100 fC, 30 fC and 10 fC for snow-graupel, ice crystal-graupel and ice crystal-snow collisions, respectively.The branching parameters χ and L χ are set to 2.3 and 1500 m, respectively.The impact of these parameters on the flash rate and lightning characteristics will be investigated in Sect.3.4.

Dynamics, microphysics and charge structure
As observed (Dye et al., 2000;Skamarock et al., 2000), the 10 July 1996 STERAO storm simulated with Meso-NH evolved from a multicell to a supercell.Figure 2 shows the composite reflectivity (maximum radar reflectivity in a grid cell column, Z max , in dBZ) at different stages of the storm.At 20 min (Fig. 2a) and 70 min (Fig. 2b), the storm was in a growing stage with remanent individual cells oriented along a NW-SE axis.At 120 min (Fig. 2c), the simulated storm started to move from a multicell to a supercell.The supercell stage is represented in Fig. 2d, with maximum radar reflectivity reaching 50 dBZ in the convective core.From radar and lightning data observations, the periods during which the storm was in a multicell stage, in transition to a supercell, and expanded as a supercell were evaluated to be 23:00-00:30 UTC, 00:30-01:15 UTC and 01:15-02:30 UTC, respectively (Skamarock et al., 2000;Barthe et al., 2010).
www.geosci-model-dev.net/5/167/2012/Geosci.Model Dev., 5, 167-184, 2012 The transition occurred after 150 min of simulation while Skamarock et al. (2000) estimated the real system transitionned after 3 h of active convection so the model is fairly successful in reproducing the evolution of the storm.Vertical cross sections of the total charge density along the wind axis at 5 km altitude were performed at different stages of the storm to analyze the charge structure and its evolution.At 20 min (Fig. 3a), the south-eastern cell was the first one to become electrified.Charges were separated by the non-inductive mechanism at ∼7.5 km altitude, generating a negative dipole (negative charge layer above positive charge layer) centered at 7.5 km altitude above sea level (ASL).As a result, the electric field increased above 10 kV m −1 at the interface of the two layers of opposite charge.50 min later (Fig. 3b), the total charge density exhibited a more complex structure, alternating between tripoles in the convective cores and positive dipoles in the cloud anvils.A negative screening layer due to negative ions was visible at the top of the cloud.At 120 min (Fig. 3c), during the transition stage, the charge structure looked like a tripole.Positively charged ice crystals were responsible for the upper positive charge layer, while snow and graupel particles were involved in the main negative layer.The positive charge layer above the ground was made up with positive graupel which fell below the 0 • C isotherm and melted into positive raindrops.During the supercell stage (170 min; Fig. 3d), the charge structure evolved between a dipole and a tripole, depending on the location in the storm.Negative ice crystals were advected in the anvil (not shown).

Lightning characteristics
The flash rate is displayed in Fig. 4a for the 3 h of simulation.During the multicell stage, the flash rate does not exceed 20 fl.min −1 which is significantly lower than the observations.Defer et al. (2001) reported indeed a maximum of 50 fl.min −1 at ∼23:45 UTC from the ONERA (Office National d'Etudes et Recherches Aérospatiales) VHF interferometric mapper (ITF).However, during the multicell stage, a large portion of the flashes detected by the ITF were identified as "short-duration" flashes.At 23:30 UTC, the "shortduration" flash ratio per 5-min period raised up to 0.46 (as estimated from Fig. 10 of Defer et al., 2001).The simulated flash rate decreases at 110 min, and reaches a plateau (∼5-10 fl.min −1 ) until 150 min.The flash rate detected by the ITF (∼10-35 fl.min −1 ) is still larger than the simulated one.During this stage of the storm, the "short-duration" flash ratio reaches 0.2.Then, the flash rate increases up to 44 fl.min −1 during the supercell stage, while in the observations, it peaks at 45 fl.min −1 .
Figure 4b-d shows the 1-min averaged number of segments, triggering altitude, triggering electric field and positive and negative charge neutralized.One can note that the average number of segments is less than 400 segments per flash except a peak of 700 segments per flash during the The number of segments, the triggering electric field, the triggering altitude and the charge neutralized are averaged over a 1-min interval.
multicell stage (Fig. 4b).The histogram showing the number of segments per flash confirms that 91 % of the flashes have less than 400 segments (Fig. 5a).
In the first part of the multicell stage (30-60 min), the average triggering altitude is 10 km a.s.l., then it decreases to 6-8 km altitude (Fig. 4c). Figure 5b confirms that two different layers exist for flash triggering.Most of the flashes are triggered between 5.5 and 7 km a.s.l., between the main negative charge and the lower positive charge (Fig. 3).Note that the triggering altitude and the triggering electric field have an opposite evolution (Fig. 4c) as expected from Eq. ( 6).
The total electric charge neutralized per flash lays between 0.4 and 20.23 C. The temporal evolution of the 1-min averaged neutralized charge shows that a peak occurs between 45 and 60 min which is associated with the larger number of segments.More than 60 % of the flashes neutralize between 1 and 4 C.These results are in agreement with data reported in the literature.Indeed, from modeling studies, the charge transfer during intra-cloud discharges was estimated to be between ∼2.7 C and ∼52.4 C (MacGorman et al., 2001;Mansell et al., 2002;Riousset et al., 2007).Similar values were found from observational studies (Krehbiel, 1981;Shao and Krehbiel, 1996;Rakov and Uman, 2003).
This simulation was performed on the University de la Réunion cluster (hereafter referred to as CCUR; Bull Novascale R422 with 20 compute nodes each with 2 × 2 Intel Xeon processors quadcore 2.26 Ghz and 24 GB of memory).About 30 % of the computing time is attributed to the electrical scheme for this simulation on 16 processors.

Sensitivity analysis
Since the electrical scheme results rely on several thresholds (for charge separation, cell detection, branching and charge neutralization), it is important to investigate the impact of varying the value of these thresholds.

Non-inductive parameterization
Several studies (Helsdon et al., 2001;Mansell et al., 2005;Altaratz et al., 2005;Barthe and Pinty, 2007a) have investigated the sensitivity of the flash rate to the NI parameterization which is recognized as the process mainly responsible for charge separation (Reynolds et al., 1957;Williams and Lhermitte, 1983;Dye et al., 1989;Latham et al., 2007).To prevent unreasonably large charging and flash rates, the magnitude of the charge exchanged per rebounding collision (δq) is generally limited.Mansell (2000) limited δq to 200 fC for graupel-snow collisions and 2 fC for graupel-crystal interactions.Mansell et al. (2005) revised these δq values to 50 fC and 20 fC for rebounding graupel-snow and graupel-crystal collisions, respectively.
Increasing δq leads to an increase of the number of flashes as shown in Table 1.With the recommended setting of Mansell et al. (2005) (NI1 case), 752 flashes are triggered, while if the charge exchanged per collision is unbounded (NI3 case), 10 times more flashes are produced (7534 flashes).The NI2 case is intermediate between NI1 and NI3.From NI1 to NI3, the number of segments is increased by a factor 3 leading to a higher quantity of charge neutralized per flash (2.98 C vs. 4.36 C on average).Consequently, the total charge neutralized during the storm is ∼14 times higher in the NI3 test (32813.6C) than in the NI1 test (2242.5C).Note that the first flash in NI3 is triggered 13 min before the first flash in NI1 as a result of allowing a larger charging rate.The flash rate evolution keeps the same trend whatever the peak value δq is set to.
The structure of the total charge density does not depart from the tripole (Fig. 3) when the limitation of the amount of neutralized charge per rebounding collision is changed (not shown).The total charge density never exceeds ±3 nC m −3 and ||E|| remains under 110 kV m −1 for the NI1, NI2 and NI3 cases.If the charge exchanged per collision is unlimited, the maximum charging rate is 288 pC m −3 s −1 for snowgraupel interactions and never exceeds 130 pC m −3 s −1 for ice crystal-snow and ice crystal-graupel elastic collisions.
In conclusion, even if the flash rate is sensitive to the setting of δq, the charge structure which depends on the temperature and on the microphysical composition (NI charging processes) is not impacted.Moreover, the new lightning scheme is efficient to limit the total charge density and the electric field module to acceptable values.In the following, the NI2 setting of the NI charging parameters which leads to the best flash data comparison with Defer et al. (2001), is preferred.
When the threshold for cell detection qcell is decreased in Q1, the size of the detected cells is larger.At the begining of the simulation when the three cells are well separated, this results in few differences.However after 30 min of simulation when the cells start to merge, the cell detection algorithm recognizes a single "big cell".Consequently, the lightning flashes are allowed to spread horizontally over the artificially "big cell", while the three original cells are still distinguishable.Decreasing qcell leads to slightly less flashes (1727 vs. 1849), but with a significantly wider extension (438 vs. 212 segments) and charge neutralization efficiency (4.32 C vs. 3.72 C).
If the charge neutralization threshold qneut is increased to 0.2 nC m −3 in Q2, it is expected that the charge neutralized per flash decreases as well (2.06 C vs. 3.72 C).As a consequence, the number of flashes is increased (3361 vs. 1849) but a similar total charge amount is neutralized (6933.74 C vs. 6880.18 C).Note that for all the experiments reported in Table 2, the mean triggering altitude lies between 7.21 km and 7.39 km, except for simulation Q2 where it is higher (8.68 km).This increase of the mean triggering altitude is due to a larger proportion of flash triggered between the upper positive and the middle negative charge layers in Q2 compared to REF.In Q2, more than 50 % of the flashes are triggered between 10 and 12 km altitude, while this proportion falls down to 20 % in REF (not shown).A possible explanation for this difference can be found in Figs 3b and 3c.The total charge density has a higher value in the lower positive charge layer (0.1-1 nC m −3 ) than in the upper positive charge layer (0.1-0.5 nC m −3 ).Then, if the charge neutralization threshold is increased, the flashes triggered in the upper part of the cloud will neutralize less charges than the ones initiated in the lower part of the cloud.So because they are less efficient, more flashes are needed to neutralize the same amount of charge in the 10-12 km layer.
In agreement with Barthe and Pinty (2007b), if χ is increased, the number of flashes decreases while the number of segments per flash (given by Eq. 7) and consequently, the amount of neutralized charge per flash increases.Keeping identical the NI parameterization and δq in the C1-C4 simulations, the total charge neutralized differs by only 13.5 % between C1 (χ = 2.1) and C4 (χ = 2.9).When χ is held constant (χ = 2.3) and L χ is increased from 500 m to 2500 m in L1-L2, the average number of segments per flash increases by a factor 3.4.As a consequence, the mean neutralized charge is higher for L χ = 500 m than for L χ = 2500 m by a factor 2.5.The L1 simulation is the only one that produces CG flashes for the STERAO storm.All are negative CGs.The first CG is produced at 94 min, and the three others are triggered between 157 and 159 min.The main difference with the other simulations lies in the relatively low number of segments per flash in L1 (91.1 ± 60.5), leading to a lower average neutralized charge per flash (2.01 ± 1.11 C).Some indicators are not significantly influenced by a variation of χ, L χ , qcell or qneut .In all the sensitivity tests performed here, the first flash still occurs at 1367.5 s because the triggering time of the first flash is only determined by the charge structure of the storm (and therefore by the electrification processes).Also, the total neutralized charge amount does not vary too much between all the sensitivity tests related to the lightning flash parameters.It is simulation L1 where globally the lowest amount of charge is neutralized.In this case, the highest number of flashes (3156) with the lowest number of branches (91.1) neutralizes on average the lowest amount of charge per flash (2.01 C).On the contrary, the C4 simulation produces the lowest number of flashes ( 1325) with the highest number of segments ( 463) and the highest amount of neutralized charge per flash (5.75 C).However, among the tests performed here but strictly for a same storm dynamics and microphysics (and so for the same electrical charging rate conditions), the quantity of neutralized charge by the flashes varies as much as ∼18 %.This discrepancy is mainly attributable to the wide range of explored parameters which brings to light the thresholds and non-linearities of the electrification scheme.(Huntrieser et al., 2002) is simulated to evaluate the tuning of the new lightning flash parameters for a severe event that developped on the evening, West of Munich, Germany.After a first period of intensification, the storm split into two distinct cells.The northernmost cell became multicellular in structure and was observed to decay soon after the cell-splitting event, while the southern cell strengthened and developed supercell characteristics.Cloud-to-ground discharges were recorded by an LPATS (Lightning Position and Tracking System) system, while the total lightning activity (IC + CG) was mapped by the ONERA interferometer.This storm has been previously simulated with the Penn State/NCAR Mesoscale Model 5 (MM5) (Fehr et al., 2004) and the 3-D Goddard Cumulus Ensemble (GCE) Model (Ott et al., 2007) to investigate lightning NO x production and transport.

Numerical set-up
The simulation performed is similar to Fehr et al. (2004).
A single thermodynamic composite sounding is used to initialize the storm (see Fig. 4 in Fehr et al., 2004).The convection is initiated with a warm (3 • C perturbation) bubble.
The model domain is chosen to be 180 × 180 × 23 km 3 in x, y, and z directions with 1 km horizontal grid spacing and 51 grid points in the vertical direction, with a variable resolution beginning at 250 m at the surface and stretching to 500 m at the top of the domain.The simulation is integrated with a 2.5 s time step for a 3 h period.As in the 10 July 1998 STERAO simulation, the physics of the model included a mixed-phase microphysics scheme (Pinty and Jabouille, 1998) and a 3-D turbulence scheme (Cuxart et al., 2000).The parameterization of Takahashi (1978) is used to describe the non-inductive processes.The lightning scheme parameters are the same as for the STERAO storm simulation.The maximum magnitude of charge separated per rebounding collision is limited to 100 fC, 30 fC and 10 fC for snow-graupel, ice crystal-graupel and ice crystal-snow collisions, respectively.The branching parameters χ and L χ are set to 2.3 and 1500 m, respectively.

Results: electrical activity
At the begining and as expected from the sounding, the storm is composed of a single vigorous cell which starts to split after the first precipitation at 60 min.This leads to a southern cell that first developed rapidly into a supercell and then decayed after 150 min (Fig. 6a-c  The electrical structure is shown in Fig. 6b and 6d.At 60 min, the storm shows a forward tilted tripole, except for the anvil with a slightly positive dipole.The highest charges are found in the updraft core at 5-6 km altitude.At 120 min, the updraft charge structure of the southern cell is much complex and an inverted dipole tends to emerge in the anvil (Fig. 6b and d).
The flash rate simulated by Meso-NH for the 21 July 1998 EULINOX storm is shown in Fig. 7.The first flash is triggered after 19 min of simulation.Then, during 75 min, the total flash rate ranges between 5 and 15 fl.min −1 .This corresponds to the first stage of weak lightning activity observed and reported by Fehr et al. (2004) (see their Fig.6b) with a slow increase of the total flash rate close to 10 fl.min −1 .After 75 min, the flash rate increases and reaches a peak of 58 fl.min −1 at 128 min.Then, the flash rate remains in the range 35-55 fl.min −1 until the end of the simulation.This sequence is slightly different in the observations.First the total lightning activity achieved an abrupt increase up, to 52 fl.min −1 around 17:45 UTC, i.e. ∼95 min after the first flash was detected by the ITF.Then, the observed total flash rate decreased gradually until 19:00 UTC (∼20 fl.min −1 ).At this time however, Fehr et al. (2004) indicated that the lightning flash detection efficiency was possibly degraded when the supercell moved over the interferometer antenna.Meso-NH produced 4521 flashes, 33 of them being cloudto-ground flashes, representing 0.7 % of the total lightning activity.Fehr et al. (2004) reported that 3321 flashes were detected, and 7 % were cloud-to-ground flashes.So one can estimate that Meso-NH is able to reproduce the characterics (amplitude and time evolution of the total flash rate) of the EULINOX storm approximately but the cloud-to-ground proportion is clearly underpredicted by a factor 10.This is not surprising since the criterion to decide the formation of cloud-to ground flashes must be refined.
The histograms of the number of segments per flash (Fig. 7b) shows that almost 85 % of the lightning flashes are made up of less than 600 segments, the average number of segments per flash being 348 segments per flash.As in the STERAO storm, lightning flashes are preferentially triggered in two layers (Fig. 7c): 5.5-6.5 km and 10-12 km, which is consistent with the charge structure displayed in Fig. 6b and  6d.Each flash neutralizes 5.02 C on average.
In conclusion, Meso-NH succeeded in reproducing the EULINOX storm and its lightning activity in case of highly simplified environmental conditions.This is a first step since only a full real case simulation followed by a direct comparison with detailed lightning measurements (Rison et al., 1999) allows an unambiguous validation of the electrical scheme.
For the EULINOX simulation performed on 32 processors of the CCUR, the electrical scheme computation represents 40 % of the total computing time.The computing time due to CELLS is expected to vary with more or less lightning activity in the domain.

Conclusions
In this paper, we report the recent improvements brought to the electrification and lightning flash scheme in Meso-NH.The governing equations of the electric charge carried by the condensed species (including hail) are completed by two equations treating the physics of positive and negative ions.Ions ensure a continuity of the electric charge out of the clouds and the precipitation in particular when particules evaporate.The ions are also responsible of the screen charge layers which form on the edge of electrified clouds.
The most important change however concerns the lightning flash scheme which has been heavily revised in order to be run in a multi-processor environment.The treatment of the flashes is the bottleneck of an electrification scheme because the filamentary aspect of the channel which needs to be resolved at high resolution, imposes to develop an algorithm which is not suitable to parallelization.Previously, the growth of the lightning discharge was based on a recursive description of the flash propagation into positive and negative pockets of charges.Here recursion comes from the rule that at a given stage of the flash, new segments are added only after a random selection of all the possible grid points that can be connected to the current structure.The new scheme simplifies this view by relaxing the connectivity criterion in order to select at once all the grid points of the flash structure.
The new flash algorithm identifies first all the independent electrified cells in order to prepare a parallel treatment of all flashes that propagate inside the cells.In most of the cases, an electrified cell is split over several subdomains on a distributed-memory computer.A bidirectional leader phase of the flash is defined for each cell with an upward and downward vertical tracing from the flash triggering grid point.The density of the branches is assumed to follow a fractal law so that the number of grid points reached by the flash increases with the distance from the triggering point.These distances, filtered by an electric charge density conditions, are computed by the processors attached to the subdomains where the flash propagates and the corresponding end grid points are stored.Then a subset of grid points at an equal distance from the triggering point are selected at random in order to fit the number of points deduced by the fractal law for this distance.All these operations can be parallelized owing an exchange of messages between all processors.
The electrical scheme was tested for two idealized storm cases drawn from the STERAO and EULINOX field experiments in order to compare first of all, the observed and the simulated flash rates.Three-hour simulations were performed over relatively large computational domains on clusters with up to 32 processors.Using the same setting for the electrification scheme, the model succeeded in reproducing www.geosci-model-dev.net/5/167/2012/Geosci.Model Dev., 5, 167-184, 2012 favorably the thousands of flashes in both cases but with a few number of CG flashes, obtained for the EULINOX case only.A sensitivity study carried out for the STERAO case helped to limit some excessive NI charging rates and to estimate non-measurable flash parameters.The next step is to run the electrical scheme for storms developing over complex terrain in order to verify that the new lightning scheme supports well a high coordinate distortion as anticipated.This objective is part of the HyMeX experiment planned in 2012 with the purpose of studying the heavy rainfalls produced by electrified orographically-forced storms in the South of France.

Fig. 1 .
Fig. 1.Different steps of the lightning flash scheme illustrated on a two-dimensional cross section of the total charge density (colored area; nC m −3 ).The black line represents the cloud contour.(a) Cell definition: the star locates the grid point where ||E|| is maximum (E max ).The dashed area delineates the electrified cell.(b) Flash triggering: the grey contours represent ||E|| with 10, 30, 50 and 70 kV m −1 contours.The crosses show the grid points where ||E|| > E trig .The black cross locates the origin of the flash.(c) Bidirectional leader: the gray contour shows locations where E z > E prop .The cross and the + show the grid points where the flash is initiated and where the leader propagates, respectively.(d) Regions of possible branching: the blue (red) contours limit the regions with negative (positive) total charge density where the negative (positive) branches can propagate.(e) Choice of the branches: the grey points represent the grid points with possible branching between two successive circles (or spheres in 3-D).N ll(r) is the maximum number of branches between two circles of radius r and r + dr (see Eq. 7).Here, the maximum number of branches is given for a 2-D framework and a mean grid size of 1450 m.(f) Lightning flash: the resulting lightning flash path made up with the triggering point (black cross), the bidirectional leader (black +) and the branches (grey +).

Fig. 2 .
Fig. 2. STERAO storm: Simulated composite radar reflectivity (Z max , in dBZ) at (a) 20 min, (b) 70 min, (c) 120 min and (d) 170 min.The + symbols indicate the origin of the lightning flashes in a 10-min interval from the time of the cross section.The black line segment [AB] corresponds to the location of the vertical cross sections of Fig. 3.

Fig. 3 .
Fig. 3. STERAO storm: Vertical cross sections of the total charge density (colors; in nC m −3 ) along the [AB] segment defined in Fig. 2 at (a) 20 min, (b) 70 min, (c) 120 min and (d) 170 min.The black solid line corresponds to the cloud contour.Dashed gray contours show the electric field module (10 and 50 kV m −1 contours).

Fig. 4 .
Fig. 4. STERAO storm: Time evolution of (a) the flash rate (fl.min −1 ), (b) the average number of segments per flash, (c) the average triggering electric field (black curve; kV m −1 ) and triggering altitude (gray curve ; km), and (d) the average positive (black curve) and negative (gray curve) charge neutralized per flash (C).The number of segments, the triggering electric field, the triggering altitude and the charge neutralized are averaged over a 1-min interval.

Fig. 5 .
Fig. 5. STERAO storm: Histograms of (a) the number of segments per flash, (b) the triggering altitude (kV m −1 ) per flash, and (c) the charge neutralized (C) per flash.Since all flashes were IC flashes, the same amount of negative charges was neutralized.

Fig. 6 .
Fig. 6.EULINOX storm: (a/c) Simulated composite radar reflectivity (Z max , in dBZ) and (b/d) vertical cross section of the total charge density (colors ; in nC m −3 ) along the [AB] segment defined in (a/c) at 60/120 min.The + symbols indicate the origin of the lightning flashes in a 10-min interval from the time of the cross section.

Fig. 7 .
Fig. 7. EULINOX storm: Lightning flash characteristics as simulated by Meso-NH with (a) temporal evolution of the total (black curve) and CG (blue curve) flash rates and histograms of (b) number of segments per flash, (c) triggering altitude (kV m −1 ) per flash, and (d) charge neutralized (C) per flash.Due to the presence of negative CG flashes, slighlty more negative charges were neutralised.
EULINOX (European Lightning Nitrogen Oxides Project) golden case of the 21 July 1998 in Germany is investigated and simulated flash statistics are provided in Sect. 4. The paper concludes on the improvements brought to the electrical scheme of Meso-NH and gives the perspective of accurate calibration when used in real case simulations.

Table 1 .
Summary of the sensitivity tests related to the limitation of the quantity of charge separated per rebounding collision.

Table 2 .
Summary of the sensitivity tests related to the flash propagation.In all these simulations, the first flash occurrence is 1367.5 s.