%~Mouliné par MaN_auto v.0.32.0 2024-05-30 10:18:45
\documentclass[CRPHYS, Unicode, XML, screen, biblatex, published]{cedram}


\daterevised{2024-03-21}

\addbibresource{crphys20230534.bib}


%\usepackage{amsmath}
\usepackage{caption}
\usepackage{xcolor}
\usepackage{changes}
\usepackage{float}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\graphicspath{{./figures/}}

\newcommand*{\mk}{\mkern -1mu}
\newcommand*{\Mk}{\mkern -2mu}
\newcommand*{\mK}{\mkern 1mu}
\newcommand*{\MK}{\mkern 2mu}

%\hypersetup{urlcolor=purple, linkcolor=blue, citecolor=red}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Energy in the heart of EM waves: modelling, measurements and management}
%\def\cdr@specialtitle@french{L'énergie au c{\oe}ur des ondes électromagnétiques : modélisation, mesures et gestion }
%\makeatother{}

\title{3D Computation of Lightning Leader Stepped Propagation Inside a Realistic Cloud}

\alttitle{Calcul en 3D de la propagation par étape d'un éclair à l'intérieur d'un nuage réaliste}

\keywords{\kwd{lightning}
\kwd{electrical discharge}
\kwd{simulation}
\kwd{modelling}
\kwd{lightning propagation}
\kwd{electric field}
\kwd{cloud}}

\altkeywords{\kwd{foudre}
\kwd{décharge électrique}
\kwd{simulation}
\kwd{modélisation}
\kwd{propagation de la foudre}
\kwd{champ électrique}
\kwd{nuage}}

\author{\firstname{Philippe} \lastname{Dessante}}
\address{Université Paris-Saclay, CentraleSupélec, CNRS, Laboratoire de G\'enie Electrique et Electronique de Paris, 91192, Gif-sur-Yvette, France}
\address{Sorbonne Université, CNRS, Laboratoire de Génie Electrique et Electronique de Paris, 75252, Paris, France}
\editornote{This article follows the URSI-France workshop held on 21 and 22 March 2023 at Paris-Saclay.}



\begin{abstract}
The simulation of lightning propagation is a complex problem studied for years. Here we propose to use the information from the electric potential created from a real cloud structure to study the propagation. The electric potential and field are calculated using a realistic thundercloud structure: the typical three layers cloud structure is constructed using a real cloud photograph. The different altitudes and separations of each layer are calculated from the luminosity of the picture and the space charge values are taken from data in the literature. A model of stepped leader propagation is proposed. It consists in finding by steps the path which maximises the potential difference taking into account the cloud and leader space charge. After each step, the electric potential is recalculated, and a new iteration gives a new direction. This framework permits us to analyse diverse cloud configurations. Only positive leaders from the base layer can reach the ground if the three layers are complete. Only negative lightning reaches the ground when the bottom positive layer is reduced (typical of the middle of a storm). Finally, when the two bottom layers are reduced in size (typical of the storm's end), positive lightning from the upper positive layer can make its way into the cloud toward the ground. These simulated observations agree with the hypotheses made previously by Nag and Rakov.
\end{abstract}

\begin{altabstract}
La simulation de la propagation de la foudre est un problème complexe étudié depuis plusieurs années. Nous proposons ici d'utiliser les informations du potentiel électrique créé à partir d'une structure de nuage afin d'étudier la propagation. Le potentiel électrique et le champ sont calculés en utilisant une structure de nuage d'orage réaliste : la structure de nuage typique à trois couches est construite à partir d'une photographie de nuage réel. Le modèle de propagation consiste à trouver par étapes le chemin qui maximise la différence de potentiel en tenant compte de la charge spatiale du nuage et du leader. Ce cadre nous permet d'analyser diverses configurations de nuages qui sont présentées dans cet article.
\end{altabstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\dateposted{2024-07-10}
\begin{document}
\maketitle


\section{Introduction}


Lightning is a research domain that has largely been studied for years. As the reproduction in the laboratory of lightning phenomena are difficult~\cite{group_research_1972,group_research_1974,group_positive_1977}, more and more studies use computer modelling and simulation to understand the propagation of leaders in clouds or toward the ground and structures~\cite{syssoev_numerical_2020}. They use a mix of macro modelling~\cite{vargas_development_2008,xu_improved_2011,xu_numerical_2021,cooray_modeling_2017} and heuristic like fractals~\cite{ghaffarpour_2022,ioannidis_fractal-based_2022} to simulate the stepped leader formation and propagation.

In this article, we use a macro model and a heuristic based on the electric potential's maximisation to find a leader's direction. This type of model has already been used in the article by Lalande et al.~\cite{lalande_physical_2012}, but with a cylindrical cloud structure. Here we use the same type of cloud structure (3 layers cloud) but based on a real cloud space charge repartition.

This research article aims to verify the different hypotheses made by the authors in~\cite{nag_inferences_2009,nag_positive_2012}. They try to explain, based on theoretical study and measurement, the different possibilities for the leader propagation toward the ground (cloud to ground or CTG) or in the cloud (intra-cloud: cloud to cloud CTC) based on their polarity and the cloud structure.

A typical cloud charge structure is based on three different layers: on top of the cloud, a positive charge; in the middle, a negative charge; and at the bottom a limited positive charge~\cite{cooray_lightning_2003}. Several simulations have been made based on heuristic rules for leader propagation~\cite{xu_numerical_2021,dulzon_modeling_1999,xu_3-d_2013}, but to our knowledge, none take into account a realistic cloud.

The different cloud layers may impact the propagation of the leader as proposed by Nag and Rakov~\cite{nag_inferences_2009,nag_positive_2012,group_research_1974}: the positive base layer can, in one part, accelerate the development of cloud-to-ground lightning, but it can also shield the bottom part of the cloud from upper lightning, favouring the development of cloud-to-cloud lightning, depending on the origin point of the leader's inception.

This research uses three kinds of cloud structures: a full three layers cloud structure, a reduced base layer, and then the two bottom layers reduced. These structures can mimic the accurate distribution of electrical charges at different storm moments.

This space charge repartition is used to calculate an electrical potential distribution above the ground. A leader propagation modelling is then proposed: The leader propagates in the direction of the maximum difference in terms of electrical potential. Different altitudes for the leader inception are tested. The simulations tend to validate the hypothesis that positive cloud-to-ground lightning comes from the base layer. This positive base layer form also a shield for the negative lightning which tends to be intracloud. Negative lightning must benefit from weaknesses or holes in the positive base layer to propagate toward the ground.

A definition of the cloud model will be presented in Section~\ref{cloud_model}, and then the stepped leader macro-model will be introduced in Section~\ref{leader_propagation}. Section~\ref{electric_potential} will explain how we compute the electric potential of the cloud and the leader. The results will be shown in Section~\ref{results}, and a conclusion and perspective will close this paper.


\section{Cloud Model} \label{cloud_model}


\begin{figure}[!bth]
\centering
\includegraphics[width=13cm]{Figures/photo.jpg}
\caption{Photography of thunderclouds over Paris (Ph. Dessante)}\label{photo}
\end{figure}


We will start from a real cloud to better model the structure. Figure~\ref{photo} shows photography of thunderstorm clouds above Paris. We will use the information in this photograph to construct the geometrical structure of the cloud. The photography is cropped to take only the cloud part (Figure~\ref{modelnuage}a). The resulting image is then converted into black and white, and smoothing based on a Gaussian filter is applied (Figure~\ref{modelnuage}b).

\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{Figures/2nuage.jpg}
\caption{Photography of thunderclouds over Paris (Ph. Dessante) (a), black and white and smoothed version (b), elevation calculation based on luminosity (c).}\label{modelnuage}
\end{figure}


The image luminosity is then analysed, and its level is stored in a two dimension array. These levels are converted into altitude (more light means less altitude) to form the cloud base layer (Figure~\ref{modelnuage}c).



\begin{figure}[!htb]
\centering
\includegraphics[width=10cm]{Figures/layers.pdf}
\caption{Typical cloud structure~\cite{cooray_lightning_2003,lalande_physical_2012}, in blue the positive upper and bottom layers, in red the middle negative layer.}\label{layers}
\end{figure}

A typical cloud structure~\cite{cooray_lightning_2003} is based on three different layers: on top of the cloud a large positive charge, in the middle a negative charge, and at the bottom a limited positive charge. In previous studies, \cite{lalande_physical_2012}, clouds are in general modelled in 3D by cylindrical zones like in Figure~\ref{layers}.

The generated altitude map is repeated at each charge separation altitude: see Figure~\ref{nuage}, where the colour represents the sign of the space charge (red for positive and blue for negative charge). The different altitude are taken in the literature~\cite{lalande_physical_2012,cooray_lightning_2003}. For the longitudinal and lateral dimension of the cloud, the dimensions are arbitrarily chosen based on the picture but also to have a sufficiently large space dimension for the leader propagation.


\begin{figure}[!htb]
\includegraphics[width=13cm]{Figures/charge_complete.pdf}
\caption{\emph{Cloud structure modelling: red negative space charge, blue to orange positive space charge}.}\label{nuage}
\end{figure}

Based on the literature~\cite{lalande_physical_2012,cooray_lightning_2003}, we have taken a positive space charge density of $c_1=0.21\times10^{ -9}$~C/m$^3$ for the positive bottom layer, a negative space charge density of $c_2=-0.18\times10^{-9}$~C/m$^3$ in the middle layer and a positive space charge density of $c3=1.3\times10^{ -9}$~C/m$^3$ for the upper layer.


\begin{figure}[!htb]
\includegraphics[width=10cm]{Figures/charge_axe.pdf}
\caption{Space Charge along the line in the middle of the cloud between the ground ($z=0$) and the upper part of the cloud.}\label{vert}
\end{figure}



When calculating the electric potential using the finite element method, it is advisable to prevent discontinuities between the layers of the cloud. Utilising a Gaussian distribution is one straightforward approach to achieving this; however, an error function (erf) or hyperbolic tangent (tanh) function could have also been used for the transition. Nevertheless, owing to the difficulty in measuring experimentally this transition, it cannot be definitively determined which option is most effective. Future studies may benefit from exploring various transition types and implementing different ones in distinct parts of the cloud, potentially to examine layer-specific properties.

The cloud space charge density vertical distribution is represented in Figure~\ref{vert}, the separation into three different layers is done by smoothing the transitions by three different Gaussian distributions:
\begin{equation}
\begin{split}
\rho(x, y, z) =
\frac{c_1}{\sqrt{2\pi}}\times\exp\left[\frac{ -(z - m_1(x,y))^2}{2 \sigma_1^2}\right]\\
+\frac{c_2}{\sqrt{2\pi}}\times\exp\left[\frac{ -(z - m_2(x,y))^2}{2 \sigma_2^2}\right]\\
+\frac{c_3}{\sqrt{2\pi}}\times\exp\left[\frac{ -(z - m_3(x,y))^2}{2 \sigma_3^2}\right]
\end{split}
\end{equation}
where $\sigma_{1,2,3}$ are respectively the standard deviation of each Gaussian's ($300$m, $2000m$ and $750$m), they are chosen arbitrarily to ensure a smooth transition between the positive and negative layers. $m_{1,2,3}$ are the middle of each cloud layers (4000m, 7500m and 11000m).


\section{Stepped Leader Propagation Model} \label{leader_propagation}

A start point $X_0$ for the leader propagation is chosen in the cloud. For each iteration $n$ of the propagation, we first find with the help of a deterministic gradient optimisation algorithm the maximal potential difference in a sphere of $r_c=50m$ radius around the leader's head. For this search radius the leader step should be less than 50m, which is consistent with the literature review in~\cite{rizk_modeling_2024}. In the case of a positive leader, this leads to the problem:
\begin{equation}\label{equ:1}
\begin{aligned}
& \underset{X}{\text{maximize}} & & V(X_n) - V(X) \\
& \text{subject to} & & |X -X_n|<r_c
\end{aligned}
\end{equation}
and in the case of a negative leader to the optimisation problem:
\begin{equation}
\begin{aligned}
& \underset{X}{\text{maximize}} & & V(X)-V(X_n) \\
& \text{subject to} & & |X -X_n|<r_c
\end{aligned}
\end{equation}

These two formulations can be reformulated respectively as:
\begin{equation}\label{optim1}
\begin{aligned}
&\underset{X}{\text{minimize}} & & V(X) \\
&\text{subject to} & & |X -X_n|<r_c
\end{aligned}
\end{equation}
and:
\begin{equation}\label{optim2}
\begin{aligned}
& \underset{X}{\text{maximize}} & & V(X)\\
& \text{subject to} & & |X -X_n|<r_c
\end{aligned}
\end{equation}

To represent the stochastic behaviour of the leader propagation, a random point is chosen around this maximum (in a sphere of 10m radius) which becomes the new leader's head. The leader is propagated with a constant radius $r_l=5m$ and a space charge density equals to $50\times10^{ -6}/(\pi r_l^2)$ C/m$^3$ if the leader originates from a positive region and $ -150\times10^{-6}/(\pi r_l^2)$ if it originates from a negative region~\cite{lalande_physical_2012,LALANDE20021375}. With $r_l=5$m the radius of the leader is somewhat overestimated due to mathematical constraints for the Poisson equation~\eqref{poisson} resolution). This process is illustrated in Figure~\ref{exp} for a negative leader in a two dimensions case for better clarity. The blue spheres represent the area where the optimisation is conducted (constraint $|X -X_n|<r_c$). The centre of the red sphere is the optimum found by equations~\eqref{optim1} or~\eqref{optim2}, and the red spheres represent the area where the new point $X_{n+1}$ is chosen randomly. The grey segments are two stepped leader iterations.

\begin{figure}[!htb]
\includegraphics[width=9cm]{Figures/explication.pdf}
\caption{Two iterations of a negative leader propagation.The blue spheres represent the area where the optimisation is conducted (constraint $|X -X_n|<r_c)$. The centre of the red sphere is the optimum found by equations~\eqref{optim1} or~\eqref{optim2}, and the red sphere represents the area where the new point $X_{n+1}$ is chosen randomly.The grey segments are two stepped leader iterations.}\label{exp}
\end{figure}


The potential is then finally recalculated by the equation~\eqref{poisson} taking into account the new space charge and a new refined mesh.



\section{Electric potential evaluation} \label{electric_potential}

The formation and propagation of lightning are mainly due to electrostatic phenomena~\cite{cooray_lightning_2003}, that are caused by the electric field. To compute the electric field $E$, one needs to solve the Poisson equation:
\begin{equation}
\nabla \vec E = \frac{\rho}{\varepsilon_0}
\end{equation}

where $\rho$ is the space charge, and $\varepsilon_0$ is the dielectric constant in vacuum.

With the introduction of the potential distribution $V$:
\begin{equation}
\vec E = -\nabla V
\end{equation}

The Poisson equation can be writen as:
\begin{equation}\label{poisson}
\Delta V(x,y,z)= -
\frac{\rho(x,y,z)}{\varepsilon_0}
\end{equation}

The boundary conditions are defined by a Dirichlet condition on the earth's surface:
\begin{equation}
V(x,y,0)=0
\end{equation}
and Neumann conditions on the others borders $\partial \Omega$:
\begin{equation}
\left. \frac{\partial V}{\partial n} \right |_{\partial \Omega}=0
\end{equation}


\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{Figures/potentiel_plein.pdf}
\caption{\emph{Electric Potential (iso planes), space charge (left plotted on a surface) for a full 3 layers cloud structure}.}\label{pot:1}
\end{figure}

We use a finite elements method (FEM) with the help of the Wolfram Mathematica software to solve the Poisson equation in three dimensions. FEM are often utilised to calculate electric fields when the geometry or the excitation (here the space charge) are non rectangular. This method discretises the space into a partition called the mesh. This mesh must be fine enough to be able to represent the small space charge variation inside the cloud but also in the leader. In order to solve the Poisson equation with the cloud space charge, a uniform tetrahedral mesh is used; we impose the maximum size of each mesh vertex to be 100m. The total number of hexahedron elements is 724 200, and the total number of degrees of freedom (DOF) is 2 987 552 for an order 2 formulation. The calculation time of such a resolution is about 30 minutes on a 28 cores computer and 128Gb of memory with Mathematica (a maximum of 16 cores are used).


The leader is symbolised by a 3D cylindrical shape with spherical endpoints, and its entire volume is populated with either positive or negative space charge depending on its origin point. Neglecting the discrepancy in charge density between the leader's head and body is a deliberate choice since the leader advances in large step (50m). Incorporating the difference of space charge exclusively during the recalculation of electric potential lacks realism, as the simulation does not account for it throughout the propagation phase.

One feasible solution to enhance realism in the simulation could be diminishing the step size for each leader stride. This modification might enable a more accurate representation of space charge differences during the leader's propagation but at the expense of computation time (the mesh needed for representing the space charge should be largely refined).



When the leader propagates, it deposes a certain number of space charges along its passage. The leader radius is 5m, and it is necessary to have a minimum of two mesh elements inside the leader channel, leading to a finer mesh. It is not possible to use a uniform mesh anymore, but it is possible to refine the mesh around the leader's shape. The superposition of this super fine mesh and the cloud mesh induces a mesh of 12 452 745 elements and 16 541 562 DOF for one step of the leader's propagation, which tends to be largely challenging to solve on a desktop computer. This computation must be realised at each leader propagation step; it leads to long and unpractical computing time. When the leader has several steps, the memory necessary to solve the Poisson problem exceeds the 128Gb available.

That is why we have used the superposition principle to solve the potential calculation. We separate the cloud space charge and each leader step (each grey segment in the Figure~\ref{exp}). The electrical potential distribution $V_0$ due to the cloud space charge is firstly calculated on the 100m fine mesh. Then at each step, the electrical potential $V_i$ due to the sole leader step $i$ (only the space charge of the last grey segment in the Figure~\ref{exp}) is calculated.

This could be done on a fine mesh around the leader tip and a broad mesh outside a circle twice the size of the leader step. The superposition principle gives the total electrical potential $V_n$ for the leader propagation at iteration $n$:
\begin{equation}
V_n=\sum_{i=0}^{n} V(i)
\end{equation}

This method helps us to only use a finer mesh around the leader step (mesh size of 2m in a distance of 100m around the leader) and a broad mesh (maximum of 1000m size) elsewhere. Tests have been realised to verify the effectiveness of this technic.



\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{Figures/champ_plein.pdf}
\caption{\emph{Electric Field (volume) and space charge on the cut surface for a full 3 layers cloud structure}.}\label{champ:1}
\end{figure}

The electric potential $V_0$ due to the cloud space charge is given in the Figure~\ref{pot:1}, the resulting electric field is shown in Figure~\ref{champ:1} and the Figure~\ref{champ_pot:1} gives a representation of the electric field (arrows and line plot), the electric potential as surface levels and the cloud space charge as a truncated volume (on the figure's left).

\begin{figure}[!htb]
\centering
\includegraphics[width=12cm]{Figures/champ_pot_tot.pdf}
\caption{\emph{Electric Field (vector and line plot), Electric Potential (surface levels), space charge (volume on left) for a full 3 layers cloud structure}.}\label{champ_pot:1}
\end{figure}


\section{Results}\label{results}



\subsection{Three layers cloud}

Before presenting the results, we give a flowchart of the complete simulation procedure in Figure~\ref{flowchart:1}. Firstly the cloud structure is defined from a picture and a complete space charge model is computed. The electric potential due to this space charge is then calculated. The inception altitude is chosen, defining if it's a positive or negative leader. After that the simulation enter its main loop where the next leader step is found as the optimum location of equation~\eqref{optim1} or~\eqref{optim2}. The leader deposes a space charge, and the electric potential $V(i)$ induced by this leader step is computed. The loop closes after the summation of all the electric potential (leader steps and cloud) or terminates if the leader reaches a border.

\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{organigramme.pdf}
\caption{\emph{Flow chart of the complete simulation}.}\label{flowchart:1}
\end{figure}


\begin{figure}[!htb]
\centering
\includegraphics[width=11cm]{propagation_plein.pdf}
\caption{\emph{Propagation of lightning in a full cloud structure, see Table~\ref{table1} for initialisation data, volume on back: cloud space charge}.}
\vspace*{10pt}
\label{propagation:1}
\end{figure}

\begin{table}[!h]
\caption{Initial conditions and results for the three layers cloud structure.}
\centering\label{table1}
\begin{tabular}{ccccc}
\hline Alitude (m) &Color & Charge & Direction & Type \\
\hline 3000 & Blue & + & $\downarrow$ &CTG \\
\hline 4000 &Purple & + & $\uparrow$ &CTC \\
\hline 5000 & Green & - & $\downarrow$ & CTC \\
\hline 6000 & Magenta & - & $\downarrow$ & CTC \\
\hline 7000 &Gray & - & $\downarrow$ & CTC \\
\hline 8000 & Orange & - & $\uparrow$ & CTC \\
\hline 9000 & Black & - & $\uparrow$ & CTC \\
\hline 10000 & Cyan & + & $\downarrow$ & CTC \\
\hline 11000 & Brown & + & $\uparrow$ & CTC \\
\hline 12000 & Red & + & $\uparrow$ & CTC \\
\hline
\end{tabular}
\end{table}

The first set of results is presented in figure~\ref{propagation:1}. The whole cloud structure in three complete layers is taken into account. Ten starting points are considered, they take inception in the middle of the cloud in lateral and longitudinal coordinates and altitudes from 3000m to 12000m above the ground level with 1000m steps. The results are summarised in the Table~\ref{table1}, which indicates the inception altitude, the colour of the leaders in Figure~\ref{propagation:1}, the leader's initial charge, the general upward or downward direction, and if they are cloud-to-ground (CTG) or cloud-to-cloud (CTC) leaders.

\begin{figure}[!htbp]
\centering
\includegraphics[width=12cm]{Figures/potential_leader_plein.pdf}
\caption{Value of the Electric Potential in front of the leader for the different inception altitudes. Plain: electric potential due to the cloud, Dashed: electric potential induced by the cloud space charge and the leader charge deposition.}\label{pot_leader:1}
\end{figure}

\begin{figure}[!htbp]
\centering
\includegraphics[width=12cm]{Figures/elecricfield_leader_plein.pdf}
\caption{Value of the Electric Field norm in front of the leader for the different inception altitudes. Plain: electric field due to the cloud, Dashed: electric field induced by the cloud space charge and the leader charge deposition.}\label{champ_leader:1}
\end{figure}


\break
The only leader who connects to the ground is the positive one coming from the bottom part of the base positive layer (inception at 3000m in blue). The negative leaders from the negative middle part of the cloud (between 5000m to 9000m) are trapped between the two (bottom and top) positive layers. They reach the separation altitude (either up or down) between the layers and then cannot propagate vertically anymore. These negative leaders cannot reach the ground, which is blinded and \emph{protected} from the bottom positive layer as shown theoretically by Nag and Rakov~\cite{nag_inferences_2009}.

One must note that no positive lightning coming from the upper part of the ground (inception from 10000m to 12000m) can reach the ground, the middle negative leader forms a shield toward the ground and stops the leaders to descend further.

The electric potential variation at the leader's front is shown in Figure~\ref{pot_leader:1} as a function of the step number. Positive leaders tend to follow decreasing potential values, and negative leaders tend to follow increasing potential values. One must note that taking into account the leader space charge (dashed curves in the figure) seems to be relatively important for intra-cloud lightning (inception points of 5000m to 10000m). There is an increase in the electric potential difference between the two cases when the leader reaches the layer separation as can be seen in Figure~\ref{propagation:1}.

The electric field variation at the leader's front is shown in Figure~\ref{champ_leader:1} as a function of the step number. One must note that the value of the electric field in front of the leader should be treated with care and used in qualitative terms only. The radius of the leader is too high (5m) to have a good estimation of the electric field induced by the leader space charge, even if the charge density is correct; the leader's shape in our calculation is too smooth. Another caution must be taken as we use a finite elements scheme to solve the Poisson equation: even if we have taken elements of order two, the electric field derivates from the electric potential, and a loss in precision is always present.

The general pattern of the electric field variation during the leader propagation is an increase in value followed by a decrease. When a leader connects to an opposite part of the cloud they would normally extinct themselves, this is correlated with the decreasing value of the electric field norm largely under their starting values. For further simulations, this could be implemented as a stopping condition in addition to the leader reaching one of the domain borders.

This is not true for the leader connecting to the ground (3000m) or the upper part of the simulation domain (11000m and 12000m) because the external electric field is almost constant in these areas and because the simulation stops when a leader reaches a border.


These results are reproducible with some variations on several runs.


\subsection{Two full layers cloud, reduced bottom layer}

\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{Figures/charge_bas.pdf}
\vspace*{-15pt}
\caption{Cloud structure and space charge for a reduced bottom layer.}\label{charge:2}
\end{figure}

To verify the hypothesis of Nag and Rakov~\cite{nag_inferences_2009,nag_positive_2012} in their articles, we have reduced the positive base layer around the middle of the cloud. The new cloud space charge is shown in Figure~\ref{charge:2}. One must note that as for the original cloud model, the transition to either neutral space charge or to other layers is smoothed by Gaussian functions to ensure a proper FEM resolution.

The new 3D potential (iso-potential surfaces) of the cloud space charge is shown in Figure~\ref{champ_pot:2} along the electric field (arrows and streamlines). The same inception conditions as in the previous section are taken as shown in Table~\ref{table:2}.

\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{champ_pot_bas.pdf}
\vspace*{-10pt}
\caption{Electric Field (vector and line plot), Electric Potential (surface levels), space charge (volume on left) for a reduced bottom layer cloud structure.}\label{champ_pot:2}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[width=12cm]{propagation_bas.pdf}
\vspace*{-15pt}
\caption{Propagation of lightning in a reduced bottom layer cloud structure, see Table~\ref{table:2} for initialisation data, volume: cloud space charge (cut for clarity.}\label{propagation:2}
\end{figure}


\begin{table}[!h]
\caption{Initial conditions and results for the two layers cloud structure.}
\centering
\label{table:2}
\begin{tabular}{ccccc}
\hline Alitude (m) &Color & Charge & Direction & Type \\
\hline 3000 & Blue & + & $\uparrow$ &CTC \\
\hline 4000 &Purple & + & $\uparrow$ &CTC \\
\hline 5000 & Green & - & $\downarrow$ & CTG \\
\hline 6000 & Magenta & - & $\downarrow$ & CTG \\
\hline 7000 &Gray & - & $\uparrow$ & CTC \\
\hline 8000 & Orange & - & $\uparrow$ & CTC \\
\hline 9000 & Black & - & $\uparrow$ & CTC \\
\hline 10000 & Cyan & + & $\downarrow$ & CTC \\
\hline 11000 & Brown & + & $\uparrow$ & CTC \\
\hline 12000 & Red & + & $\uparrow$ & CTC \\
\hline
\end{tabular}
\end{table}

The leader's propagations are shown in Figure~\ref{propagation:2}. This time the negative pink leader coming from the inception point at 6000m inside the negative layer can turn around the reduced base layer and reach the ground. This simulation validates the theoretical hypothesis of Nag and Rakov~\cite{nag_inferences_2009} and shows that negative lightning must profit some weaknesses or holes in the base positive layer to propagate towards the ground.

In this configuration, the positive leaders coming from the base layer are attracted to the middle negative layer, and none of them can reach the ground. This is also the case for the upper part of the clouds, as in the previous simulation.

During normal storms, the vast majority of lightning strikes reaching the ground are negative. The common model of the cloud space charge also agrees with this type of repartition~\cite{cooray_lightning_2003}, where the positive bottom layer is largely reduced in size or thickness.



\subsection{One full layer cloud and two reduced bottom layers}

\begin{figure}[!htb]
\centering
\includegraphics[width=11cm]{charge_bas_milieu.pdf}
\vspace*{-20pt}
\caption{Cloud structure and space charge for a reduced bottom and middle layers.}\label{charge:3}
\end{figure}

\begin{table}[!hb]
\caption{Initial conditions and results for the one layer cloud structure.}
\centering\label{table:3}
\begin{tabular}{ccccc}
\hline Alitude (m) &Color & Charge & Direction & Type \\
\hline 3000 & Blue & + & $\downarrow$ &CTG \\
\hline 4000 &Purple & + & $\downarrow$ &CTG \\
\hline 5000 & Green & - & $\uparrow$ & CTC \\
\hline 6000 & Magenta & - & $\uparrow$ & CTC \\
\hline 7000 &Gray & - & $\uparrow$ & CTC \\
\hline 8000 & Orange & - & $\uparrow$ & CTC \\
\hline 9000 & Black & - & $\uparrow$ & CTC \\
\hline 10000 & Cyan & + & $\downarrow$ & CTG \\
\hline 11000 & Brown & + & $\uparrow$ & CTC \\
\hline 12000 & Red & + & $\uparrow$ & CTC \\
\hline
\end{tabular}
\end{table}

Finally, we treat the representative case of the final moments of the storm. Under the effect of the different wind speeds in altitude and the horizontal development of the cloud, the upper positive layer expands laterally and becomes much larger than the bottom layers~\cite{cooray_lightning_2003}.

\begin{figure}[!htb]
\centering
\includegraphics[width=12cm]{champ_pot_bas_milieu.pdf}
\vspace*{-15pt}
\caption{Electric Field (vector and line plot), Electric Potential (surface levels), space charge (volume on left, for a reduced bottom and middle layers cloud structure.}\label{champ_pot:3}
\end{figure}

\begin{figure}[!hb]
\centering
\includegraphics[width=13cm]{propagation_bas_milieu.pdf}
\vspace*{-15pt}
\caption{Propagation of lightning in a reduced bottom and middle layers cloud structure, see Table~\ref{table:3} for initialisation data.}\label{propagation:3}
\end{figure}

The space charge resulting from this deformation is modelled and shown in Figure~\ref{charge:3}. The positive and negative bottom layers are laterally reduced with the help of a 3D gaussian centred in the middle of the cloud, ensuring that no high gradient can perturbate the Poisson equation resolution. The electrical potential (iso-potential surface) for this calculation is shown in Figure~\ref{champ_pot:3} along the electric field (arrows and streamlines). The same inception altitudes as in the previous section are taken, as shown in the Table~\ref{table:3}.




The propagation of the various leaders is shown in Figure~\ref{propagation:3}. In this situation, the only leaders connecting to the ground are those from the cloud's positive upper and bottom parts. The upper positive charge space attracts the leaders from the negative middle part.

The two leaders coming from the bottom positive layer (blue and magenta) are attracted as in the first case (3 layers scenario) by the ground: the potential difference given by equation~\eqref{equ:1} is more important in this direction as can be seen in Figure~\ref{pot_leader:3}.

The electric potential variation at the leader's front is shown in Figure~\ref{pot_leader:3} as a function of the step number. Positive leaders tend to follow decreasing potential values, and negative leaders tend to follow increasing potential values. Again, this figure shows that taking into account the leader space charge is necessary for a proper calculation of the electric potential during the propagation of intraclouds lightning.

The leader coming from the positive top layer (cyan) is attracted firstly toward the negative middle cloud, but it can pass through it as it is reduced. The bottom positive layer bends the field lines and repulses the leader toward its border. At last, it can find a path around this positive layer toward the ground. This is corroborated by the Figure~\ref{pot_leader:3} where the leader coming from the inception point 10000m (cyan) follows a potential always decreasing but with a little inflexion point at the moment of the leader bypassing the bottom positive layer.

This scenario aim at modelling the last moments of the storm. It has been observed~\cite{cooray_lightning_2003, nag_inferences_2009, nag_positive_2012} that during the last stage of the thunderstorm, positive lightning becomes the majority of the lightning that connect to the ground.


\section{Conclusions} \label{conclusions}


In this paper, we have used a lightning propagation model aimed at testing the hypotheses made by Nag and Rakov about the different possibilities of cloud-to-ground or cloud-to-cloud leader propagation. We have constructed a new cloud model from a picture of a real thunderstorm cloud for this aim. The space charge inside the cloud is composed of three layers: a positive charge at the bottom and top, and a middle negative charge layer.

This space charge is used to compute an electric potential. Introducing a stepped leader propagation model helps us characterise for different inception altitudes the behaviour of lightning inside the cloud. We have found that when the bottom positive layer covers the whole cloud in latitude and longitude, it protects the ground from negative lightning as it repulses the negative leaders inside the cloud. Only positive leaders from the bottom positive layer reach the ground in this case.

If the bottom positive layer is reduced in size, which is more representative of a typical thunderstorm cloud~\cite{cooray_lightning_2003}, the only lightning reaching the ground is a negative one that takes inception in the cloud middle layer.

Due to wind gradients, the end of a storm is characterised by reduced bottom and middle layers. In this scenario, as proposed by Nag and Rakov~\cite{nag_inferences_2009,nag_positive_2012}, only leaders coming from the two positive layers can reach the ground.

This model associated with this cloud structure has permitted us to validate these hypotheses. The photography presented in the Picture~\ref{photoroage}, where lightning reaching the ground comes from the bottom part of the cloud, corroborate these hypotheses visually.

\begin{figure}[!htbp]
\centering
\includegraphics[width=12cm]{Figures/potential_leader_bas_milieu.pdf}
\caption{Value of the Electric Potential in front of the leader for the different inception altitudes. Plain: electric potential due to the cloud, Dashed: electric potential induced by the cloud space charge and the leader.}\label{pot_leader:3}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[width=13cm]{Figures/orage2.jpg}
\caption{Photo of a storm showing bottom CTG leaders and top CTC leaders (Photography Ph. Dessante).}\label{photoroage}
\end{figure}


We think that this model, associated with other cloud structures, can improve the comprehension of the leader's propagation in different types of clouds (holes or reduction in layer thickness, for example).

In the future, we plan to introduce branching in the propagation model to have a better lightning representation, in this case, one must calculate the electric the potential of each branch, and this will induce interaction between branches.

Another improvement is to reduce the leader's radius to have a better estimation of the electric field. This could lead to a stopping condition for leader propagation.

\section*{Declaration of interests}
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

\printbibliography
\end{document}
