# Corona emission from raindrops in strong electric fields as a possible discharge initiation: comparison between horizontal and vertical field configurations

Sylvain Coquillat, Bruno Combal, Serge Chauzy

► To cite this version:

Sylvain Coquillat, Bruno Combal, Serge Chauzy. Corona emission from raindrops in strong electric fields as a possible discharge initiation: comparison between horizontal and vertical field configurations. Journal of Geophysical Research Space Physics, 2003, D7 (108), pp.4205. 10.1029/2002JD002714. hal-00137520

HAL Id: hal-00137520

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

Submitted on 23 Aug 2021

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Corona emission from raindrops in strong electric fields as a possible discharge initiation: Comparison between horizontal and vertical field configurations

Sylvain Coquillat, Bruno Combal, and Serge Chauzy

Laboratoire d'Aérologie, UMR CNRS/UPS 5560, Observatoire Midi-Pyrénées, Toulouse, France

[1] The purpose of this work is to determine which electric field configuration (vertical or horizontal) produces corona emission from raindrops for the lowest ambient electric field. For that, a numerical modeling of the distortion of uncharged raindrops falling at terminal velocity in quiescent air in a horizontal electric field is performed. The results are compared with previous numerical study involving a vertical electric field. It is shown that the fall velocity is quite unaffected by ambient field intensities lower than 200 kV/m. The disruption and the corona onset fields are lower than those corresponding to the vertical field configuration and the larger the drop the larger the difference. For a given altitude the difference between the corona onset fields in both configurations can reach about 100 kV/m; meanwhile for a given ambient field intensity, the difference in altitude of corona emission can rise to about 1 km. These onset field intensities are, though, too high for allowing drop breakup and corona emission from an uncharged drop in permanent field conditions below 12 km height. An estimation of the critical onset field (disruption and/or corona emission) is carried out for charged drops in the horizontal field configuration. For a drop 2 mm in spherical equivalent radius carrying the quarter of its Rayleigh maximum net charge, the critical onset field is approximately equal to 290 kV/m below 7.4 km height and decreases down to 110 kV/m at 10 km. Given that present modeling does not take into account the effect of turbulence that could induce drop oscillations, the critical field intensities calculated here may be considered as upper limits.

## 1. Introduction

[2] In order to account for the electrostatic processes in meteorological models of cloud convection, the various interactions between cloud microphysics and thunderstorm electrification should be investigated. The final aim is to find relevant parameterizations for phenomena such as lightning initiation, collision and collection efficiencies, riming effect, raindrop distortion and disruption, and falling velocity changes. Each of them must therefore be analyzed and quantified. As regards to the lightning initiation, many parameters still remain unknown: cloud volume of intense electric field, onset ambient electric field, and microphysics involved. Several mechanisms likely to produce corona, which is the first step to a propagative discharge, have been studied experimentally [Dawson, 1969; Richards and Dawson, 1971; Griffiths and

Latham, 1972, 1974; Griffiths et al., 1976]. They involve various types of hydrometeors like water drops, ice particles, or melting hailstones. The corresponding electric fields required to trigger corona, whatever the altitude level is, range above about 400 kV/m which is the strongest electric field that could have been measured in a thundercloud [see MacGorman and Rust, 1998]. However, Crabb and Latham [1974], Blyth et al. [1992], and Blyth et al. [1998] showed that glancing collisions involving at least one raindrop can produce long liquid filaments at the tips of which corona is easily emitted in electric fields down to 150 kV/m [Blyth et al., 1998]. Even if the probability of occurrence of such interactions needs to be evaluated, this type of mechanism provides the lowest thresholds of corona emission from hydrometeors ever measured experimentally.

[3] Models by Coquillat and Chauzy [1994] and Georgis et al. [1995] indicate that corona can also be triggered at realistic altitudes from stable raindrops in vertical electric fields of the same intensities, provided either their electric
charge is high, or pairs of drops interact. On the other hand, if most of the cloud-to-ground lightning flashes are supposed to be initiated in a vertical ambient field, the horizontal structure of intracloud discharges is often predominant [MacGorman et al., 1981], especially in storms with large stratified regions [Krehbiel et al., 2000]. Several experiments [Nolan, 1926; Macky, 1931; Ausman and Brook, 1967; Kamra and Ahire, 1989; Kamra et al., 1993; Georgis et al., 1997] were performed to study the influence of a horizontal electric field on the behavior of raindrops. They show that the disruption field at sea level pressure is clearly lowered in such a field configuration for uncharged drops since both aerodynamical and electrical effects act in the same direction when disturbing the drop stability. In this way, we wonder whether pure corona at higher altitudes would also be easier to trigger in the horizontal field configuration, especially for charged drops. It can be pointed out from the experimental work by Dawson [1969] that corona emission from water drop surfaces is favored by an intense surface electric field and by a large radius of curvature. Previous modeling of raindrop behavior in a vertical ambient field downward directed [Coquillat and Chauzy, 1993, 1994] shows that corona is emitted from the bottom of positively charged drops and from the top of negatively charged drops where the surface field is rather low and the radius of curvature large. In a horizontal ambient field, the combination of electric and aerodynamic distortions leads to a corona emission in the ambient field direction from the equator of the drop where the surface field is rather high but the radius of curvature small (see pictures in Kamra et al., 1993). However, we still do not know which field configuration is the most favorable to corona triggering.

[4] The goal of present study is to compare the influence of the field direction (vertical or horizontal) on the triggering of corona from raindrops. We carried out this work in a numerical way in order to overcome the practical constraints of a laboratory experiment almost impossible to conduct at variable pressure and with charged raindrops. In order to achieve this comparison, we performed a modeling of the behavior of uncharged raindrops falling at terminal velocity in quiescent air in a horizontal electric field - not modeled, the case of charged raindrops will however be discussed. The results are directly comparable to those of Coquillat and Chauzy [1993, 1994] since their model is based on the same assumptions. The equilibrium shape is computed in a step-by-step method by distorting an initially spherical drop in order to make the internal overpressure uniform. As noted by Kamra et al. [1993], not only the terminal velocity and the direction of the disruption field but also the drop oscillations are key parameters of the breakup characteristics of the drops. Therefore, a raindrop modeling that neglects the oscillations provides the average shape of raindrops. The altitude of corona occurrence computed subsequently must be regarded as an average altitude as it is in Coquillat and Chauzy [1994] and in the present work.

[5] The main problems arising from this modeling are first, the computation of the surface electric field and second, the determination of the aerodynamic pressure necessary to calculate the terminal velocity. The former is

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Geometry of the pseudo-ellipsoidal shape of a raindrop elongated by the horizontal ambient field $E_a$ in the x direction. $v_f$ features the airflow direction.</figcaption></figure>

solved by an integral method and the latter is realized by adapting the semiempirical method by Beard and Chuang [1987] to the shape of horizontally stretched drops.

## 2. Distortion Modeling

[6] As observed and analyzed by Kamra et al. [1993] in a wind tunnel experiment, raindrops undergo oscillations and even rotation around their vertical minor axis when falling at terminal velocity. However, the turbulence to which the drops are submitted inside thunderclouds may not exhibit the same characteristics than that produced in a wind tunnel. Therefore, the shape of the drops could oscillate around the equilibrium shape in a way that seems difficult to accurately model. This is one reason for which we do not take into account the oscillations in the present modeling. The other reason arises from the comparison between the present work, the numerical study by Coquillat and Chauzy [1993, 1994] as indicated in section 1, and the experimental work by Georgis et al. [1997]. As a matter of fact, the disruption electric field deduced from the present modeling will be compared and scaled with the results of this laboratory experiment in which uncharged raindrops fell at terminal velocity in quiescent air in a horizontal electric field. Georgis et al. [1997] have shown that the oscillations undergone by the drops are largely damped by their fall through the still air in a 17 meters height column so that their shape can be compared to that in an equilibrium state. Therefore, we focus this modeling on the equilibrium shape of the drops in order to allow a correct comparison with this experiment and the numerical study in the vertical field configuration by Coquillat and Chauzy [1993, 1994].

[7] In a horizontal electric field, the drops are stretched along the field direction due to induction. Consequently, their equilibrium shape can be described as two half-ellipsoids having two common horizontal axes, the longer being aligned with the field direction, as shown in Figure 1. $E_a$ features the ambient electric field, a and b the perpen-
dicular horizontal semiaxes, c<sub>t</sub> the vertical semiaxis of the top half of the drop, c<sub>b</sub> the vertical semiaxis of the bottom half, and v<sub>f</sub> the direction of the airflow velocity. Also displayed in Figure 1 are the spherical coordinates (r, θ, φ) of a surface point M used for the calculation of the surface electric field and of the aerodynamic pressure introduced in section 3.

[8] The basic procedure of the present modeling consists in adjusting, in a quasi-static way, the length of each semiaxis in order to make the internal overpressure ΔP uniform at several test points labeled from 1 to 4 in Figure 1. This procedure is performed for increasing ambient electric fields as far as the convergence of the ΔP values is possible. The maximum electric field compatible with an equilibrium shape is considered as the disruption field. On the basis of the local processing of Laplace's equation by Taylor [1964], developed by Abbas and Latham [1969] and Coquillat and Chauzy [1993], the computation of the internal overpressure arises from the local pressure balance. This balance involves the hydrostatic pressure P<sub>h</sub>, the aerodynamic pressure P<sub>a</sub>, the electrostatic pressure P<sub>e</sub>, and the pressure due to surface tension P<sub>s</sub>:

$$ \Delta P + P_h + P_a + P_e = P_s, \quad (1) $$

where:

$$ P_h = \rho_w g h, \quad (2) $$

ρ<sub>w</sub> being the water density, g the gravitational acceleration, and h the height of water above the considered point;

$$ P_a = \frac{\rho_a v^2}{2} P_d, \quad (3) $$

ρ<sub>a</sub> being the air density, v the terminal velocity of the drop (see section 4), and P<sub>d</sub> the dimensionless aerodynamic pressure (see section 3) at the considered point;

$$ P_e = \frac{\epsilon_0 E^2}{2}, \quad (4) $$

ε<sub>0</sub> being the dielectric constant of the air, and E the electric field near the surface at the considered point. This equation is valid as long as the electric field is perpendicular to the drop surface, that is, as long as it can be assumed that the drop behaves as a conductor, which is stated in the present paper. At last, the pressure due to surface tension is given by:

$$ P_s = \frac{2T_s}{R_c}, \quad (5) $$

T<sub>s</sub> being the coefficient of surface tension, and R<sub>c</sub> the equivalent curvature radius at the considered point. R<sub>c</sub> is analytically calculated from both principal curvature radii determined at each test point from the corresponding ellipses. For the test points located at the equator (i.e., labeled 1 and 2), the principal curvature radius in the vertical plane is averaged since the top and bottom ellipses have different semiaxes in the vertical direction.

[9] As far as the surface electric field E is considered, an integral method based on a sampling of the surface charge density was performed. It is described in the following. The uncharged raindrop is considered as a conductor having a uniform electric potential equal to zero and the presence of the ambient electric field E<sub>a</sub> induces a surface charge separation by electrostatic influence. The drop is divided into n surface samples. An unknown value σ<sub>i</sub> of the surface charge density is associated with each sample i. The electric potential is then written at the center of the ith sample as follows:

$$ V(i) = \sum_{j=1}^{n} \frac{\sigma_j}{4\pi\epsilon_0} \int \int \frac{dS_j}{d_{ij}} - E_a x_i = 0, \quad (6) $$

where d<sub>ij</sub> is the distance from a given point of the jth sample to the center of the ith sample, dS<sub>j</sub> is the elementary surface of the jth sample, E<sub>a</sub> is the ambient electric field, and x<sub>i</sub> is the abscissa of the center of the ith sample. If the electric potential is expressed at each center of the n surface samples, we obtain a system of n linear equations, which is set in the following matricial form:

$$ \begin{pmatrix} a_{11} & a_{12} & \cdots & a_{1n} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{pmatrix} \begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \vdots \\ \sigma_n \end{pmatrix} = E_a \begin{pmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{pmatrix}, \quad (7) $$

where the coefficients a<sub>ij</sub> are the multiplicands of σ<sub>j</sub> in equation (6). However, since the electric charge distribution displays a plane of symmetry defined by y = 0 and a plane of antisymmetry defined by x = 0, the values of σ<sub>j</sub> are the same (except for the sign) on each quarter of the drop. It is therefore possible to factorize them so that the number of equations required to solve the problem is reduced by a factor 4. In this way, equation (6) is written only on one quarter of the drop and the new coefficients a<sub>ij</sub> take into account the contribution of four samples. A specific process is applied in the case where i = j corresponding to the electric potential calculated on the considered sample itself. The coefficient a<sub>ii</sub> is numerically computed by a double integral routine like all coefficients a<sub>ij</sub> except on a small element, located at the center of the sample. The potential on this element, of sides 2A and 2B, is approximated by the potential due to a plane quadrilateral of same sides 2A and 2B according to Durand [1964]:

$$ V_{ii} = \frac{\sigma_i}{\pi \epsilon_0} \left[ A \ln \left( \frac{B + \sqrt{A^2 + B^2}}{A} \right) + B \ln \left( \frac{A + \sqrt{A^2 + B^2}}{B} \right) \right]. \quad (8) $$

At last, equation (7) corresponding to the quarter of the drop is solved by the Gauss method and the surface electric field E(i) near the surface of any point i is given by:

$$ E(i) = \frac{\sigma_i}{\epsilon_0}. \quad (9) $$
<table><caption>**Table 1.** Reynolds Number of Precipitating Raindrops of Spherical Equivalent Radius R for Three Atmospheric Conditions<sup>a</sup></caption><thead><tr><th>Radius, mm</th><th>0.5</th><th>1.0</th><th>1.5</th><th>2.0</th><th>2.5</th></tr></thead><tbody><tr><td>T = -10°C; P = 500 hPa</td><td>202</td><td>668</td><td>1277</td><td>1877</td><td>2446</td></tr><tr><td>T = 0°C; P = 700 hPa</td><td>234</td><td>802</td><td>1515</td><td>2249</td><td>2915</td></tr><tr><td>T = 20°C; P = 1013 hPa</td><td>267</td><td>862</td><td>1605</td><td>2345</td><td>3018</td></tr></tbody></table>

It can be pointed out that the present method can easily be made suitable to the case of a charged raindrop. The electric potential of the drop being unknown, an additional equation is required which is given by the equation of the charge conservation:

$$ \sum_{i=1}^{n} \iint \sigma_i \, dS_i = Q. \qquad (10) $$

### 3. Computation of the Aerodynamic Pressure Distribution

[10] The main problem encountered in the modeling of raindrop distortion concerns the determination of the aerodynamic pressure distribution at the drop surface, and therefore the calculation of the flow with detached wake about nonspherical bodies. For raindrops of spherical equivalent radii ranging from 0.5 mm to 2.5 mm, the corresponding Reynolds numbers $N_{Re}$ lie between about 200 and 3000 as indicated in Table 1 for three different atmospheric conditions. We suppose that the flow regimes are similar to that of a rigid sphere in a laminar steady free-stream velocity as described by Beard [1976]. The flow should therefore be characterized by the appearance of steady vortex trails ($N_{Re} < 300$) and the formation of vortex loops by roll-up and partial detachment of wake vortex ring ($300 < N_{Re} < 450$). At large Reynolds numbers ($N_{Re} > 450$) a wake detachment takes place and becomes turbulent. In the following, we consider that the general feature of the flow exhibits a detached and turbulent wake whatever the size is. The induced error for small drops is negligible since for radius lower than about 0.7 mm ($N_{Re} = 481$), the contribution of the aerodynamic pressure never exceeds 7% of the pressure due to surface tension in the pressure balance (equation (1)).

[11] Given that the pressure distribution is to be introduced in an iterative model of raindrop distortion, we deliberately avoided carrying out any 3D computation of the airflow with turbulent detached wake since it would have been much more expensive in computation time. Thus, we determine the pressure distribution in the semiempirical way proposed by Beard and Chuang [1987] by adjusting the pressure distribution measured around a sphere to the shape of a nonaxisymmetrical body.

#### 3.1. Description of the Method

[12] The choice of the pressure distribution around a sphere was very restricted. No distribution corresponding to Reynolds numbers lower than $N_{Re} = 60 \times 10^3$ like that

measured by Maxworthy [1969] has been found in the literature because of experimental constraints that arise from this type of measurement [see Beard and Chuang, 1987]. Therefore we base our calculations on the nondimensional distribution of Maxworthy [1969] whose characteristics are described hereafter. The angle at which the local pressure is the same as that upstream is close to $\theta = 135^\circ$ (actually $45^\circ$ from the upstream direction), $\theta$ being defined in Figure 2. The minimum pressure before separation appears at $\theta = 109^\circ$ ($71^\circ$ from the upstream direction) with a magnitude of about $-0.56$. The detachment occurs close to the equator with a maximum pressure after separation of about $-0.38$. At last, the pressure keeps close to $-0.4$ for the whole downstream half of the sphere. Maxworthy measured all these parameters in function of Reynolds numbers ranging between $60 \times 10^3$ and $200 \times 10^3$. The asymptotic behavior of his measurements indicates that only the location of the minimum pressure should change, though slightly from $109^\circ$ to $110^\circ$ ($71^\circ$ to $70^\circ$ from the upstream direction, respectively) for Reynolds numbers lower than $10^3$. In this way, the pressure distribution by Maxworthy [1969] is expected to be a good reference for the flow around a sphere in the range of Reynolds numbers corresponding to that of raindrops.

[13] Since the minimum pressure location is not very sensitive to the Reynolds number value in the case of a sphere, we make the same hypothesis in the case of a distorted raindrop. In further calculations about nonaxisymmetrical bodies, we will thus consider that the minimum pressure keeps close to the tangent angle $\chi = 71^\circ$, which implicitly corresponds yet to different zenith angles $\theta$ in each $\phi$-direction. The tangent angle $\chi$, which is defined in Figure 2, allows normalizing the pressure distribution for spheroids [Beard and Chuang, 1987]. The actual location of the detachment is governed by the flow in the boundary layer about which we make no hypothesis. Hence, the detachment is supposed to appear at the equator even though a turbulent airflow (like updrafts in thunderclouds) could induce turbulence in the boundary layer and consequently a backward motion of the detachment. Therefore, no correction is made on the location of these two particular points. According to Beard and Chuang [1987], the cor-

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Definition of the tangent angle $\chi$ in a vertical plane $\phi$.</figcaption></figure>

<footnote><sup>a</sup>Calculations were performed using the terminal velocity of raindrops measured at see level pressure by Gunn and Kinzer [1949], and the method by Beard [1976] to take into account the influence of pressure and temperature.</footnote>
rected pressure for the drop P<sub>d</sub>(χ) in a plane φ is computed
from:

$$
P_d(\chi) = 1 - k(\chi)[1 - P_s(\chi)] \quad (11)
$$

by using the stagnation point as a reference, P<sub>s</sub>(χ) being the
pressure distribution for a sphere measured by Maxworthy
[1969], and k(χ) a shape dependent factor. Assuming that
the fractional deviation of P<sub>d</sub>(χ) from P<sub>s</sub>(χ) would be the
same as the fractional deviation in potential flow of a
pseudo-ellipsoid from a sphere, k(χ) is given in the unse-
pared flow region (0° < χ < 71°) by:

$$
k(\chi) = \frac{1 - P_d(\chi)}{1 - P_s(\chi)} = \frac{1 - P_{\text{potd}}(\chi)}{1 - P_{\text{pots}}(\chi)}, \quad (12)
$$

where P<sub>pots</sub>(χ) = 1 - 9/4 sin<sup>2</sup>(χ) is the analytical solution for the potential flow around a sphere, and P<sub>potd</sub>(χ) is the pressure distribution for the potential flow around a pseudo-ellipsoidal body presented in the following section 3.2. In the downstream half of the drop (90° < χ < 180°), the comparison between the pressure distributions for a sphere by Fage [1937], Maxworthy [1969], and Achenbach [1972] at different Reynolds numbers shows that the nondimensional pressure undergoes minor variations and remains close to about -0.4 in the whole turbulent wake. Thus, the pressure does not seem to be sensitive neither to the location at drop surface nor to the Reynolds number. Therefore k(χ) is held at the uniform value k<sub>wake</sub> which is expected to be close to unity since the pressure in the turbulent wake should not be very dependent on the drop shape unlike in the unseparated flow region. At last, the value of k(χ) between the minimum pressure and the equator is a linear transition from k(χ) at the minimum pressure to k<sub>wake</sub>.

**3.2. Calculation of the Potential Flow**

[14] Milne-Thomson [1938] found a solution for the
potential flow around ellipsoidal bodies by solving the Lap-
lace equation. He gave a specific expression for the potential
velocity using an ellipsoidal coordinate system. This kind of
solution has been used in particular topics like atmospheric
turbulence measurements [Oost, 1991] or aeronautical prob-
lems [Band and Payne, 1980]. For axisymmetrical shapes
like those of freely falling raindrops, Grover and Beard
[1974] determined a series solution for the potential flow.
Beard and Chuang [1987] used an analytical solution, based
on the stream function given by Happel and Brenner
[1965]. They calculated the potential velocity v<sub>p</sub> at the
surface of an oblate spheroidal shape (i.e., the revolution
shape around the small vertical axis of an ellipse) featuring
the lower half of the drop. This velocity is given in oblate
spheroid coordinates (ξ, η) by:

$$
v_p = \frac{v \sin \eta}{[(\lambda^2 + 1) - \sin^2 \eta]^{1/2} [(\lambda^2 + 1) \cot^{-1} \lambda - \lambda]}, \quad (13)
$$

where η is defined by tan η = α tan θ, and λ is the ratio of
the axis ratio α to the eccentricity.

[15] Nevertheless, none of these solutions is suitable to
the pseudo-ellipsoidal shape of the horizontally stretched
drops studied here. Therefore we performed a numerical
computation of the potential velocity v<sub>p</sub> (see Appendix A)

from which the dimensionless pressure is deduced accord-
ing to the Bernoulli theorem:

$$
P = 1 - \frac{v_P^2}{v^2}, \quad (14)
$$

where v is the terminal velocity of the drop. This
computation led us to point out very convenient results.
As far as the potential flow is considered, the pressure
distribution versus the tangent angle χ is independent of
the azimuth angle φ. Furthermore, this pressure distribu-
tion around an ellipsoidal body of vertical semiaxis c and
of horizontal semiaxes a and b is exactly the same than
that around a spheroidal body of same vertical semiaxis c
and of horizontal semiaxis equal to (a + b)/2. In this way
it is possible to apply the analytical solution of the
potential velocity v<sub>p</sub> for spheroids (equation (13)) used by
Beard and Chuang [1987] to the top and bottom
ellipsoids of the drops studied here. The pressure
distribution of the potential flow around a pseudo-
ellipsoidal body is thus a combination of both analytical
solutions P<sub>t</sub>(χ) and P<sub>b</sub>(χ) deduced from equations (13)
and (14) for the top and bottom spheroids of horizontal
and vertical semiaxes [(a + b)/2, c<sub>t</sub>] and [(a + b)/2, c<sub>b</sub>],
respectively. Actually it is built up as a weight average in
which each analytical distribution is associated to a
weight function depending on the tangent angle. The
weight function w<sub>b</sub>(χ) of the bottom spheroid distribution
decreases linearly from 10 to 1 between χ = 0° and χ =
90° and from 1 to 0 until χ = 180°. The weight function
w<sub>t</sub>(χ) of the top spheroid distribution increases linearly
from 0 to 1 between χ = 0° and χ = 90° and from 1 to
40 until χ = 180°. At each tangent angle χ the pressure
distribution of the potential flow is given by:

$$
P_{\text{potd}}(\chi) = \frac{w_b(\chi)P_b(\chi) + w_t(\chi)P_t(\chi)}{w_b(\chi) + w_t(\chi)}. \quad (15)
$$

Figure 3 displays the comparison between P<sub>potd</sub>(χ) and the numerical solution calculated by spherical harmonics (see Appendix A) in three different φ-planes for one given shape. Also displayed are portions of the analytical pressure distributions P<sub>t</sub>(χ) and P<sub>b</sub>(χ) corresponding to the associated top and bottom spheroids from which P<sub>potd</sub>(χ) is deduced. The maximum weights 10 and 40 were chosen to ensure a good enough fitting of the numerical solution whatever the shape is. This figure confirms that the pressure distribution is the same regardless of the φ-direction and one can see the fairly good agreement between the fast calculation of P<sub>potd</sub>(χ) presented above and the numerical computation.

**3.3. Determination of k<sub>wake</sub>**

[16] The determination of k<sub>wake</sub> arises from the balance
between the pressure drag D<sub>p</sub> of the drop, computed by
integration of the pressure distribution P<sub>d</sub> to the whole
surface:

$$
D_p = -\frac{1}{2} \rho_a v^2 \int_{\text{drop}} P_d n_z \, dS, \quad (16)
$$

where ρ<sub>a</sub> is the air density and n<sub>z</sub> is the vertical component
of the local normal to the surface, and the theoretical
<figure><img src="image_4.png" /><figcaption>Figure 3. Comparison between the dimensionless pressure computed with the spherical harmonics method and the analytical solutions for spheroids of horizontal semiaxis equal to the average of both semiaxes of the drop: a = 1.00, b = 0.90, c_t = 0.85, c_b = 0.65.</figcaption></figure>

pressure drag which is given by the following equation if the buoyancy force is neglected:

$$D_p = D \frac{C_{D_p}}{C_D} = M_d g \frac{C_{D_p}}{C_D}, \quad (17)$$

where D is the total drag on the drop (i.e., the aerodynamic drag force), C_Dp and C_D are the pressure drag and the total drag coefficients respectively, M_d is the drop mass, and g is the gravitational acceleration (9.8 m/s²). According to the assumption of *Beard and Chuang* [1987] and *Chuang and Beard* [1990], the ratio C_Dp/C_D is considered to be the same for a raindrop as for a sphere. It is calculated from the following empirical formula based on the measurements of *Achenbach* [1972]:

$$\frac{C_{D_p}}{C_D} = 1 - n N_{Re}^{-m}, \quad (18)$$

where the equivalent radius is used to define the Reynolds number N_Re. The relation between C_Dp/C_D and N_Re is displayed in Figure 4 with n = 13.4 and m = 0.58 corresponding to the interpolation performed by *Beard and Chuang* [1987] of *Achenbach's* [1972] data with the numerical result of *LeClair et al.* [1970] for N_Re = 400.

[17] In order to test the determination of k_wake we performed a calculation for a raindrop of radius R = 2 mm falling at its terminal velocity v = 8.83 m/s in still air [Gunn and Kinzer, 1949]. Its shape is defined by the geometrical parameters a = 1, b = 1, c_t = 0.96, and c_b = 0.62 [Coquillat and Chauzy, 1993] corresponding to the actual semiaxes of

the drop normalized by the semiaxis in the x direction. We found k_wake with a relative accuracy of 0.3% corresponding to an accuracy of 1% on the pressure drag balance. It can be pointed out that this balance is rather sensitive to the value of k_wake but as expected, k_wake is close to unity (0.977). Consequently the corrected pressure is close to the measured distribution by *Maxworthy* [1969] in the turbulent wake as it is shown in Figure 5. The pressure distribution exhibits an artificial break due to the change in the coefficient k(χ) at the equator but this region plays a negligible part in the pressure drag since the local normal to the surface is horizontally directed.

## 4. Computation of the Terminal Fall Speed

[18] Knowing the aerodynamic pressure distribution for a given shape, it is possible to evaluate the change in terminal fall speed due to the change in raindrop shape. In the following, the subscript 0 corresponds to the raindrop in a zero field; meanwhile no subscript corresponds to the distorted drop in the horizontal electric field. If we neglect the buoyancy effects, we can say that the drop is in a vertical equilibrium state under the influence of its weight and of the aerodynamic drag force when it reaches its terminal velocity (see *Gay et al.* [1974] in the case of a vertical electric field). The presence of an electrical field does not add a new force in this balance since its acts in the horizontal direction and above all, the drop is uncharged. Therefore, the sole influence of the electric field is to change the shape, the velocity, and the drag coefficient of the drop in such a way that the aerodynamic drag force
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Ratio of the pressure drag coefficient to the total drag coefficient in function of the Reynolds number.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Corrected dimensionless pressure for the real flow over a raindrop (a = 1.00, b = 1.00, c<sub>t</sub> = 0.96, c<sub>b</sub> = 0.62) falling at terminal velocity in a zero electric field.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 6. Axis ratio in the ambient field direction versus spherical equivalent radius for various ambient field intensities. For measurements of Best, see Lane and Green [1956].</figcaption></figure>

remains equal to the weight of the drop. The total drag being unchanged, one can write:

$$C_{D0} \frac{1}{2} \rho_a v_0^2 \pi a_0 b_0 = C_D \frac{1}{2} \rho_a v^2 \pi a b, \quad (19)$$

where the total drag coefficients can be expressed in function of the pressure drag coefficients and the Reynolds numbers according to equations:

$$\left\{ \begin{aligned} C_{D0} &= \frac{C_{Dp0}}{1 - n N_{Re0}^{-m}} \\ C_D &= \frac{C_{Dp}}{1 - n N_{Re}^{-m}} \end{aligned} \right. . \quad (20)$$

At last, the pressure drag coefficients are calculated from the dimensionless pressures:

$$\left\{ \begin{aligned} C_{Dp0} &= \frac{1}{\pi a_0 b_0} \int_{\text{drop}_0} P_{d0} n_{z0} dS_0 \\ C_{Dp} &= \frac{1}{\pi a b} \int_{\text{drop}} P_d n_z dS \end{aligned} \right. . \quad (21)$$

The combination of equations (19), (20), and (21) leads to the following expression for the terminal fall speed of the distorted raindrop:

$$v = v_0 \left( \frac{1 - n N_{Re}^{-m}}{1 - n N_{Re0}^{-m}} \frac{\int_{\text{drop}_0} P_{d0} n_{z0} dS_0}{\int_{\text{drop}} P_d n_z dS} \right)^{\frac{1}{2}} . \quad (22)$$

[19] Since the Reynolds number $N_{Re}$ is calculated in function of the actual fall speed and horizontal section, $v$ is computed in an iterative way. In the first iteration the unknown $N_{Re}$ is replaced by $N_{Re0}$ to evaluate the fall speed. In the following iterations, $N_{Re}$ is calculated in function of $v$ previously evaluated. This procedure converges very quickly to the solution since the change in $N_{Re}$ is very weak and $v$ is computed with an accuracy of about 0.01% in 3 iterations. It can be pointed out that the decrease in fall speed from $v_0$ to $v$ is of the same order as the increase in horizontal section. However, several computations with various arbitrary distorted shapes showed us that there is no linear relation between both variations.

## 5. Results

[20] In order to compare the present modeling with the previous one in vertical field configuration [Coquillat and Chauzy, 1993, 1994], we performed the calculations with the same pressure $P = 700$ hPa and temperature $T = 0°C$. These atmospheric parameters can be considered as average pressure and temperature corresponding to the range of altitudes where large liquid water drops reside inside thunderclouds. This low cost procedure does not introduce a significant error in the results since the distortion is only slightly affected by the couple of parameters $P$ and $T$. As a matter of fact, the decrease in $P$ at higher altitudes tends to distort the drops via the increase in fall velocity, while the combined decrease in $T$ counteracts this process via the increase in the surface tension.

### 5.1. Raindrop Distortion

[21] As expected, the shape of raindrops distorted by a horizontal electric field is characterized by a horizontal stretching in the field direction. Figure 6 displays the axis ratio in the field direction $\alpha = (c_t + c_b)/2a$ as a function of
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Influence of the ambient field intensity on the terminal velocity: comparison between vertical and horizontal field configurations in relative units.</figcaption></figure>

drop size. As expected, the higher the ambient field, the higher the horizontal stretching, and the lower the axis ratio. If we compare the relative decrease of $\alpha$, in a horizontal electric field, to the relative increase of the axis ratio, in a vertical field [see *Coquillat and Chauzy*, 1993], they are of the same order of magnitude for sizes equal to 0.5 and 1 mm. But the distortion becomes markedly more pronounced for larger drops for which the aerodynamic pressure becomes high. For example, a drop of radius $R = 1.5$ mm falling in a horizontal field of 800 kV/m is characterized by a decrease in its axis ratio of about -22% meanwhile in a vertical field of the same intensity, its axis ratio increases only of about 11%. This illustrates the combination of both aerodynamic and electric distortions, which act together in a horizontal field whereas they counteract in a vertical field. Nevertheless, the distortion is rather weak for ambient fields lower than 200 kV/m since the relative decreasing of $\alpha$ does not exceed about -2.5% for the largest drop studied here ($R = 2$ mm) in such an ambient field. Furthermore, the raindrops keep the same silhouette in the vertical plane perpendicular to the field direction since the axis ratio in the plane perpendicular to the field direction defined by $\alpha_{\perp} = (c_t + c_b)/2b$ is almost unaffected whatever the field intensity and the drop size are.

[22] The present modeling is however not fully satisfactory for large drops. As a matter of fact, no equilibrium shape was found in the case where $R = 2.5$ mm, even in a zero field. This problem could originate from the simplified geometry that would not be suitable for the larger drops or, more probably, from the neglecting of the internal circulation. As analyzed by *Brazier-Smith* [1992], this phenomenon tends to limit the distortion, and therefore to reinforce the stability of the drop. In order to understand the behavior of the present modeling, we compared it with that in a vertical field configuration performed by *Coquillat and Chauzy* [1993]. Both of them are based on a local pressure

balance of the internal overpressure, use the same geometry, and neglect the internal circulation. Therefore they should give comparable results in a zero field when *Coquillat and Chauzy* [1993] found an equilibrium shape for $R = 2.5$ mm. In fact, in this previous study the terminal velocity was computed according to the method by *Gay et al.* [1974] adjusted to the distortion under the influence of a vertical electric field. This procedure led to underestimate the velocity at the ground in a zero field as compared with *Gunn and Kinzer* [1949] experimental data. In the present modeling, the terminal velocity at the ground in a zero field fits the data by *Gunn and Kinzer* so that the neglecting of the internal circulation seems not to be balanced by a reduction of the velocity.

[23] Given the above analysis, we reduced the terminal velocity for the larger drops. In doing so, we assume that this change in the velocity (from 7% in a zero field to 20% in the disruption field for $R = 2.5$ mm) represents, via the aerodynamic pressure, the contribution of the internal circulation. As a matter of fact, this correction solves the problem since the internal overpressure $\Delta P$ at the bottom of large drops was systematically higher than all other overpressures (equator and top of the drop), and the contribution of the aerodynamic pressure to the internal overpressure is positive at this sole test point.

## 5.2. Velocity Change

[24] To evaluate the change in terminal fall speed due to the change in raindrop shape, we used the method presented in section 4. The results are displayed in Figure 7 in relative units, dividing the terminal velocity by the velocity in a zero field. Both horizontal and vertical field configurations are displayed. Whatever the configuration is, the velocity variation versus the increasing electric field exhibits an increase for small drops and a decrease for large drops. This variation is shifted toward faster relative
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the disruption field computed by the present model with Taylor's theoretical criterion and three previous experimental results. The equivalent spherical radius ranges between 0.5 and 3 mm.</figcaption></figure>

velocities in the vertical electric field configuration since, in this case, the vertical stretching due to electric effects counteracts the horizontal flattening due to the aerodynamic distortion, leading to better shaped drops. The unexpected velocity behavior for large drops is tightly connected to the evolution of their horizontal section. In increasing horizontal fields, this one markedly increases leading to a velocity reduction. This behavior is also undergone by large drops in a vertical field. In this case, the higher the field, the higher the vertical stretching. This leads to a higher velocity and a higher aerodynamic distortion, which induces, in turn, a velocity reduction because of a strong flattening of the drop base. Nevertheless, the relative velocity just barely diverges for typical electric field values observed in thunderclouds. As far as charged raindrops are concerned, they undoubtedly should be accelerated in the field direction so that their shape should be asymmetrically distorted by the horizontal induced airflow. But the change in the vertical component of their velocity should not be as large as it is in a vertical field for which the Coulomb force acts in the same direction as the weight [see *Coquillat and Chauzy*, 1993].

**5.3. Disruption Field and Critical Field Profile**

[25] The computed disruption field is displayed in Figure 8 and compared with previous experimental studies for raindrops whose equivalent spherical radius ranges between 0.5 and 3 mm. For R = 0.5 mm, the disruption field roughly agrees with the theoretical criterion of Taylor [1964] since neither the hydrostatic pressure nor the aerodynamic pressure are key parameters for the distortion of small raindrops. For increasing sizes until 1.5 mm, the disruption field tends to the lower bound of all data that stem from experiments in quiescent air. It is logically lower

than the experimental results by *Macky* [1931] and *Ausman and Brook* [1967] since the drops were not allowed to fall at terminal velocity in these experiments. Furthermore, the curve seems to be in such a good agreement with experimental results by *Georgis et al.* [1997] that we based the modeling correction of the terminal velocity (see section 5.1) in such a way that the disruption field matches these experimental data for the largest drops. It can be pointed out that the disruption field measured by *Kamra et al.* [1993] are markedly lower than the whole set of data since on one hand, the drops fall at terminal velocity, and on the other hand they undergo oscillations and vibrations due to the continuous turbulence induced in the wind tunnel. The comparison between the data by *Kamra et al.* [1993] and *Georgis et al.* [1997] clearly shows that the turbulence is a key parameter in determining the disruption field. Therefore, the actual disruption field in thundercloud conditions should be closer to the results by *Kamra et al.* [1993]. However, one must keep in mind that the turbulence of intracloud conditions does probably not display the same features as those of a wind tunnel experiment.

[26] The critical horizontal field, that is the horizontal ambient field that causes either disruption at lower altitudes or pure corona emission at higher altitudes, has subsequently been calculated. The procedure is the same as that used by *Coquillat and Chauzy* [1994]. For a given ambient field, the local electric field and the corresponding radii of curvature of the drop surface are determined by the present modeling. Thus the corona occurrence altitude is calculated from the corona experimental results obtained by *Dawson* [1969]. Let us remind that this altitude should be lower for an oscillating and vibrating drop [*Kamra et al.*, 1993] than that presently calculated for a stationary drop, whatever the field direction is. Here, the comparison between horizontal
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Critical field profiles established for 4 values of the equivalent spherical radius (0.5 to 2.0 mm) of the drop. Comparison of the vertical and the horizontal field configuration.</figcaption></figure>

and vertical field configuration remains valid since both results derive from equilibrium shapes of raindrops. The vertical profile of the critical horizontal field is displayed in solid lines in Figure 9; the dashed lines correspond to the vertical field configuration. As usual, each field profile displays a vertical straight line up to the transition altitude. This line corresponds to the drop disruption, which is roughly altitude independent. Above this altitude, the critical field decreases since it corresponds to pressure dependent pure corona. The comparison between the vertical and the horizontal field configurations shows that, for a given drop size, the latter provides lower critical field intensities and allows corona triggering at lower altitudes than the former. Furthermore, the larger the drop, the more important the lowering of the corona occurrence altitude. For example, the difference between the altitude of corona occurrence from a small drop (0.5 mm) at a given field tends to zero. On the contrary, it may become larger than 1 km for a drop of 2 mm in radius: it reaches 2.2 km for an ambient field of 400 kV/m. In this way, we can state that corona emission is more easily triggered from raindrops falling in a horizontal rather than in a vertical ambient field. Obviously, the critical field profiles in Figure 9 should be disregarded for altitudes above about 11–12 km that correspond to temperatures lower than –40 °C in which water drops cannot remain in a liquid state. This limitation of the height of the –40 °C isotherm is in agreement with most of the thundercloud observations [see MacGorman and Rust, 1998], in which the –40°C isotherm is most often located higher than 10 km in the convective core of thunderclouds. For example, *Weinheimer et al.* [1991] reported a temperature of –26°C at 9.8 km.

## 6. Case of Charged Raindrops

[27] The behavior of charged raindrops in a horizontal electric field appears to be very difficult to observe in a

laboratory experiment since a charged particle submitted to an electric field is instantaneously ejected from the observation volume because of the Coulomb force action. Even a numerical study is difficult to perform because there is no symmetry that could simplify the treatment, which is already rather complex for uncharged raindrops as it can be seen in the present paper. As a matter of fact, the shape of the drop is: (1) flattened at the bottom because of the airflow induced by its fall, (2) elongated in the field direction because of the induced charge separation, and (3) flattened at one tip in the field direction since the velocity gets a horizontal component that disorganizes the whole shape. However, as it has been shown by Coquillat and Chauzy [1993, 1994] in the vertical field configuration, the presence of electric charge on a drop leads to a pronounced reduction of the disruption field and/or the corona occurrence altitude. Given that the critical field is lower, first, for charged raindrops than for uncharged in the vertical field configuration, and second in the horizontal than in the vertical field configuration for uncharged raindrops, it should logically reach the lowest values for charged raindrops in a horizontal electric field. As far as the triggering of corona is concerned, it has been shown by Dawson [1969] that both the surface electric field and the curvature radius are the key parameters. In a horizontal ambient field, the corona must be triggered at the drop surface where the local field is higher, that is, at the horizontal tips of a raindrop. If the raindrop is charged, this local field is enhanced by the excess of charge at one of these tips. Moreover, the horizontal component of the fall velocity creates a flattening at the same tip that increases the corresponding curvature radius, which is known to favor the corona emission at equal local field [Dawson, 1969]. Even if this flattening reduces in turn the surface field, the combination of the large radius of curvature and the somewhat lower surface field leads to even better conditions for corona triggering because the
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Critical field profiles for a 2 mm radius drop in vertical and horizontal field configurations. The profile for Q/Qr = 1/4 in the horizontal configuration is an estimation.</figcaption></figure>

volume of high field is large enough to ensure that a minimum number of free electrons can trigger the avalanche process. Therefore it seems sound to extrapolate the results presented in the previous sections to charged drops in order to better approach the real mechanisms and therefore enhance the quality of the lightning triggering parameterizations introduced into thundercloud models at this time.

[28] According to Coquillat and Chauzy [1993], the decrease of the vertical disruption field versus net charge normalized to the Rayleigh's limit QR is quite linear. Furthermore, the slope of this decrease is a weak function of the radius except for the larger size studied (2.5 mm). If we omit this size because of the uncertainties arising from the present modeling, the vertical ambient field Ea (V/m) that causes disruption can be expressed in function of the radius R (m) and the net charge Q by:

$$E = E_0 - (k_1 R + k_2) \frac{Q}{Q_R},$$

where Ea0 is the disruption field of the uncharged raindrop, k1 = 289 × 10⁶ V m⁻², and k2 = 1.43 × 10⁶ V m⁻¹. Assuming that this decrease is of the same order in the horizontal field configuration, one can estimate the horizontal disruption field for charged raindrops. As far as pure corona is concerned, the lowering of the corona onset field from Q/QR = 0 to Q/QR = 1/4 in the horizontal field configuration - at each altitude above the transition altitude - is assumed to be the same than that in the vertical field configuration. Following these assumptions, we calculated the critical field profile for a 2 mm in radius raindrop carrying a net charged Q = 450 pC equal to the quarter of the Rayleigh's limit (QR = 1800 pC at 20°C), which allows the comparison with previous results by Coquillat and

Chauzy [1993]. This electric charge appears high, indeed unlikely to exist in thunderclouds. Actually, we cannot definitely state about the charge carried by raindrops in thunderclouds because there are very few simultaneous measurements of charge, size, and nature of hydrometeors in the literature. At the present time, only the kind of observations performed by Weinheimer et al. [1991] and Takahashi et al. [1999] can provide this information. The former reported net charges of about -460 pC but they were unable to associate the corresponding particle shape, and the data of the later concern only winter storms. On the other hand, several measurements of charge and size measurements have been performed by soundings [Marshall and Winn, 1982; Marsh and Marshall, 1993; Bateman et al., 1995; Bateman et al., 1999] but nature of the observed particles remains unknown. Among them, Marshall and Winn [1982] reported net charges up to 400 pC with diameters ranging between 1 and 3 mm. The same kind of measurements have been carried out from aircraft and benefit from the information of an observer [MacCready and Proudfit, 1965; Latham and Stow, 1969; Gaskell et al., 1978; Christian et al., 1980]. MacCready and Proudfit [1965] reported several charges between 500 and 1000 pC at altitudes ranging from 5000 to 5500 m in supercooled clouds. They found charges up to -300 pC in melting drops, between -150 and -300 pC in rain, of about +300 pC in a mixture of graupel and rain, more than -1000 pC at an altitude of about 3500 m, and a maximum negative charge of about -1500 pC on a liquid water drop the radius of which was supposed to be around 3 mm. These measured charges remain by far the highest available in the literature.

[29] The critical field profile is plotted in Figure 10. For a drop 2 mm in spherical equivalent radius carrying the quarter of its Rayleigh maximum net charge, this estimated critical field reaches low values as compared to the dis-
ruption fields for uncharged raindrops. It is approximately equal to 290 kV/m below the altitude of 7.5 km and decreases down to 110 kV/m at 10 km. These field intensities remain within the range of intense fields already measured inside thunderclouds [see MacGorman and Rust, 1998]. Furthermore, it must be reminded that these results relate to nonoscillating raindrops so that all the values calculated here may be considered as upper limits.

## 7. Conclusion and Perspectives

[30] The present modeling of raindrop behavior in a horizontal electric field is far from perfect because it does not allow us to describe the hydrodynamic oscillations and vibrations undergone in thundercloud conditions. Nevertheless, it is suitable to reach the objectives: determining if the raindrops located in a horizontal electric field can be more efficient for discharge triggering than the raindrops located in a vertical electric field. For doing so, we made a comparison with a previous numerical study by Coquillat and Chauzy [1994] based on the same assumptions and in which the conditions of corona emission from drop surface where calculated for a vertical field configuration. Since both studies determine the average shape of raindrops, the corona occurrence altitude calculated is therefore an upper limit of the actual corona occurrence altitude. In the case of uncharged raindrops, the results clearly show that the horizontal configuration is much more efficient for corona triggering than the vertical one. For example, the altitude of corona emission can be more than 2 km lower in a horizontal field than in a vertical field, that is, for an uncharged raindrop of radius R = 2.0 mm in a 400 kV/m field intensity.

[31] As far as charged raindrops are concerned, an estimation of the corona and disruption onset fields shows that the critical horizontal field can be reduced to intensities that remain within the range of intense fields actually measured. Thus, the present results sustain the conclusions by Kamra et al. [1993] according to whom the corona triggering from large raindrops in horizontal electric field may be responsible for lightning discharge initiation. The transition altitude from disruption to pure corona emission being rather high (7.4 km for Q/QR = 1/4 in a horizontal field) this process would be more efficient in deeply convective cells. Furthermore, given that the oscillations and the vibrations of the drops induced by the intracloud turbulence [Kamra et al., 1993], or the interaction between drops with or without collision [Blyth et al., 1992; Georgis et al., 1997], can reduce these onset fields, the triggering of corona discharges by raindrops located in an horizontal field region of a thundercloud is a process that should be taken into account in the parameterization of the electrical processes in the atmospheric models.

[32] Further comparisons with in situ observations are required in order to evaluate the probability of occurrence of such a mechanism. These observations should gather together nature, size, and charge of precipitation particles, electric field soundings, and detection of lightning flash origins, especially in tropical regions. Few studies have been performed to localize the lightning sources and to simultaneously observe the microphysics in the corresponding cloud volumes. Proctor [1991] showed a bimodal distribution of the median heights of the origins of 773

lightning flashes located between 4.4 and 5.7 km above mean sea level for a lower group, and between 6.7 and 9.7 km for an upper group. The corresponding temperatures ranged from +1°C to -8.5°C and from -21.1°C to -33°C, respectively. However, the microphysical data were reported in terms of radar reflectivity. The use of data available in literature would give hypothetical conclusions since the multiplicity of parameters required to state about this kind of mechanism (net charge, size, nature of hydrometeors, electric field, altitude, temperature, lightning flash origins) are related to different storms, even to different phases of storm evolution.

[33] The subsequent perspective that arises from such a numerical study of raindrop behavior is the coupling with a model of discharge propagation. As a matter of fact, the computed conditions of corona emission can only state if corona is triggered or not from a given raindrop in a given ambient field. The further development of the corona discharge still remains unknown. Does it propagate and lead to a cloud or ground discharge, or does it participate in a local dissipation of the ambient field? This question remains open at this time. For ensuring a propagation at large scale, the discharge channel should include a hot thermalized part of high conductivity associated with a current flow. Griffiths [1977] suggested that several streamer systems could propagate in rapid succession. The whole mechanism should probably involve more than one triggering particle. On the other hand, Baker [2001] proposed an alternative mechanism based on the acceleration of high energy electrons produced by cosmic rays that could lead to propagative discharges. The coupling could help us either to determine the most favorable conditions for natural triggering of lightning or to evaluate the amount of charge emitted during an aborted discharge, which is of importance for a good parameterization of the electrical processes in a thundercloud model. This latter evaluation could also be performed through a laboratory experiment provided we are able to control the ambient pressure.

## Appendix A: Computation of the Potential Flow by Use of Spherical Harmonics

[34] As far as potential flows are considered, the absence of viscosity implies that the velocity field remains irrotational. Consequently, the condition:

$$ \vec{\nabla} \times \vec{v}_p = \vec{0} \qquad (A1) $$

corresponds to a velocity field which derives from a potential function $\psi_p$:

$$ \vec{v}_p = \vec{\nabla} \psi_p. \qquad (A2) $$

Given the continuity equation:

$$ \frac{d\rho}{dt} + \rho \vec{\nabla} \cdot \vec{v}_p = 0, \qquad (A3) $$

the incompressible nature of the potential flow imposes that the velocity field is nondivergent, therefore it yields:

$$ \vec{\nabla} \cdot \vec{v}_p = \vec{\nabla} \cdot \vec{\nabla} \psi_p = \Delta \psi_p = 0 \qquad (A4) $$
<figure><img src="image_12.png" /><figcaption>**Figure A1.** Dimensionless pressure distribution for the potential flow over oblate spheroids. Comparison between the spherical harmonics method and the analytical solutions of Beard and Chuang [1987].</figcaption></figure>

where $\Delta\psi_p = 0$ is the Laplace's equation. As for the limit conditions, they bear on the normal component of the velocity, with respect to the surface of the pseudo-ellipsoid, which is zero since no fluid flux exists throughout that surface. Therefore, the fluid velocity of a potential flow is purely tangential to the surface:

$$ \vec{\nabla} \psi_p \cdot \vec{n} = 0, \qquad (A5) $$

and there is no condition on the tangential component of the velocity because of the absence of viscosity.

[35] The present methodology consists in solving the Laplace's equation (A4) with condition (A5) on the whole surface of the pseudo-ellipsoidal body corresponding to the raindrop shape, the velocity field being deduced from equation (A2) and the pressure distribution from equation (14). The mathematical form of the velocity potential $\psi_p$ is the sum of the velocity potential $\psi$ of the uniform flow at an infinite distance from the raindrop, with a product of three functions of separated variables. These ones are used to take into account the decrease in distance $r$ of the flow disturbance and its dependence on both angles $\theta$ and $\phi$. The latter is expressed as a linear combination of powers of $1/r$ multiplied by spherical harmonics $Y_{\ell}^m$ depending on $\theta$ and $\phi$. Consequently, the solution that we are looking for is:

$$ \psi_p = \psi + \sum_{\ell=0}^{\infty} \sum_{m=0}^{\ell} C_{\ell}^{m} \frac{Y_{\ell}^{m}}{r^{\ell+1}}, \qquad (A6) $$

where the coefficients $C_{\ell}^{m}$ are scalar real values to determine. Assuming that $t = \cos\theta$, the spherical harmonics $Y_{\ell}^m$ are given by:

$$ Y_{\ell}^{m} = \frac{1}{2^{\ell}\ell!} \sqrt{\frac{2\ell + 1}{4\pi} \frac{(\ell - m)!}{(\ell + m)!}} (1 - t^2)^{\frac{m}{2}} \frac{d^{m+\ell}}{dt^{m+\ell}} (t^2 - 1)^{\ell} e^{im\phi}. \qquad (A7) $$

[36] Furthermore, the above general expression of the velocity potential can be simplified according to the three following items:

[37] (1) The velocity being a real value vector, the potential $\psi_p$ should take only real values. Thus, only the real part of $e^{im\phi}$ is considered in equation (A7).

[38] (2) The flow is symmetrical with respect to the planes defined by $x = 0$ and $y = 0$, this implies $m$ is even.

[39] (3) The fluid flux throughout a surrounding sphere of given radius $R$ large compared to the drop size is zero since there is no fluid source [Guyon et al., 1991]. Therefore the first term of the development ($\ell = 0$) is removed from equation (A6).

[40] In order to finely account for the pseudo-ellipsoidal shape, the mesh on the surface is built with about $N_{\phi} = 30$ longitudinal angles $\phi$, and with at least $N_0 = 30$ zenithal angles $\theta$ for weaker distortion areas to about $N_0 = 80$ for higher distortion areas. The $N_0$ zenithal angles are determined so that the surface points are equidistant in a given plane $\phi$ since in the case of the most distorted bodies studied here, the accuracy is sensitive to angle $\theta$ discretization.

[41] Introducing the potential definition (A6) in the equation (A5) and taking into account the simplifications pre-
viously mentioned, the problem can be expressed for the
whole grid of points as a set of linear equations:

$$
AX = B, \tag{A8}
$$

where A is a rectangular matrix with N₀ × Nₚ lines
corresponding to the number of points on the mesh, and Nℓm columns corresponding to the number of terms in ψₚ. The
A-matrix elements are:

$$
\vec{\nabla} \left( \frac{Y_{\ell}^{m}}{r^{\ell+1}} \right) \cdot \vec{n}, \quad (A9)
$$

while the column matrix B contains the gradient of the
potential ψ multiplied by the local normal. The unknown in
equation (A8) is the column matrix X whose Nℓm elements
are the Cℓm values to determine. Given that the minimum
value for Nθ × Nφ is about 900 and the maximum value for
Nℓm is 271 (see below), the set of linear equations (A8) is
strongly overdetermined. Its resolution needs the least
squares method which consists in transforming the set of
equations in a square Nℓm × Nℓm set, multiplying (A8) by
the transpose of the rectangular matrix A:

$$
A^t A X = A^t B, \quad (A10)
$$

where A† denotes the transpose of matrix A. Equation (A10)
is solved with a classical Gauss' pivot method.

[42] The higher the distortion, the higher the number N<sub>ℓm</sub> of terms necessary to reach a given accuracy on ψ<sub>p</sub> in equation (A6). For example, in the case of a sphere, only one term (C<sub>1</sub><sup>0</sup> = 0.5) is required so that the potential function is:

$$
\psi_p = \psi + \frac{1}{2r^2} \cos \theta, \quad (A11)
$$

when, if more complex spheroidal shapes are considered,
one needs $\ell$ varying from 1 to 5 for $\alpha = 0.9$, from 1 to 9 for
$\alpha = 0.8$, and from 1 to 27 for $\alpha = 0.7$ if an accuracy at least
equal to 1% is required. For higher distortions, the
summation on $\ell$ in equation (A6) is arbitrarily limited to a
maximum value $\ell_{max} = 31$ because of computer capacity
constraints, so that the maximum value of $N_{\ell m}$ is equal to
271 since m takes successive even values in the range $[0, \ell]$.
Consequently, this method is here no more valid for
distorted ellipsoids characterized by an axis ratio lower than
0.6 as it can be seen in Figure A1. This figure displays the
dimensionless pressure distribution versus the tangent angle
$\chi$ for the potential flow around various oblate spheroids
computed with the spherical harmonics method, and the
corresponding analytical solutions (see equations (13) and
(14)) by *Beard and Chuang* [1987] for comparison. One can
observe in Figure A1 that the agreement between the two
kinds of solution is fairly good for axis ratio ranging from
1.0 to 0.7. The error is about 2% for $\alpha = 0.6$ at the equator
but it rises about 20% for $\alpha = 0.5$ whose pressure
distribution curve also exhibits oscillations for tangent
angles close to $\chi = 0^\circ$ and $\chi = 180^\circ$.

[43] In the case of ellipsoidal bodies, the results found
here are in very good agreement with *Band and Payne*
[1980] who found that the dimensionless pressure is con-
stant around the equator of an ellipsoid (of axes a ≠ b ≠ c).

For pseudo-ellipsoidal bodies we found a very fast solution
(see section 3.2) based on the analytical solutions by *Beard
and Chuang* [1987] by using the present computation as a
calibration available for axis ratio higher than 0.6.

References

Abbas, M. A., and J. Latham, The disintegration and electrification of charged water drops falling in an electric field, *Q. J. R. Meteorol. Soc.*, **95**, 63–76, 1969.

Achenbach, E., Experiments on the flow past spheres at very high Reynolds numbers, *J. Fluid Mech.*, **54**, 565–575, 1972.

Ausman, E. L., and M. Brook, Distortion and disintegration of water drops in strong electric fields, *J. Geophys. Res.*, **72**, 6131–6135, 1967.

Baker, M., Physical processes of thunderstorm electrification, *The Physics of Thundercloud and Lightning Discharges, Onera Scientific Day*, November 6th, Chatillon, France, 2001.

Band, E. G. U., and P. R. Payne, The pressure distribution on the surface of an ellipsoid in inviscid flow, *Aeronaut. Q.*, **31**(1), 70–84, 1980.

Bateman, M. G., W. D. Rust, B. F. Smull, and T. C. Marshall, Precipitation charge and size measurements in the stratiform region of two mesoscale convective systems, *J. Geophys. Res.*, **100**, 16,341–16,356, 1995.

Bateman, M. G., T. C. Marshall, M. Stolzenburg, and W. D. Rust, Precipitation charge and size measurements inside a New Mexico mountain thunderstorm, *J. Geophys. Res.*, **104**, 9643–9653, 1999.

Beard, K. V., Terminal velocity and shape of cloud and precipitation drops aloft, *J. Atmos. Sci.*, **33**, 851–864, 1976.

Beard, K. V., and C. Chuang, A new model for the equilibrium shape of raindrops, *J. Atmos. Sci.*, **44**, 1509–1524, 1987.

Blyth, A. M., H. J. Christian, and J. Latham, Corona initiation in thunderclouds, *Proc. 9th Int. Conf. Atmos. Electr.*, **1**, 104–107, 1992.

Blyth, A. M., H. J. Christian, and J. Latham, Corona emission thresholds for three types of hydrometeor interaction in thunderclouds, *J. Geophys. Res.*, **103**, 13,975–13,977, 1998.

Brazier-Smith, P. R., On the shape and fall velocities of raindrops, *Q. J. R. Meteorol. Soc.*, **118**, 749–766, 1992.

Christian, H., C. R. Holmes, J. W. Bullock, W. Gaskell, A. J. Illingworth, and J. Latham, Airborne and ground-based studies of thunderstorms in the vicinity of Langmuir Laboratory, *Q. J. R. Meteorol. Soc.*, **106**, 159–174, 1980.

Chuang, C., and K. V. Beard, A numerical model for the equilibrium shape of electrified raindrops, *J. Atmos. Sci.*, **47**, 1374–1389, 1990.

Clift, R., J. R. Grace, and M. E. Weber, *Bubbles, Drops and Particles*, 380 pp., Academic, San Diego, Calif., 1978.

Coquillat, S., and S. Chauzy, Behavior of precipitating water drops under the influence of electrical and aerodynamical forces, *J. Geophys. Res.*, **98**, 10,319–10,329, 1993.

Coquillat, S., and S. Chauzy, Computed conditions of corona emission from raindrops, *J. Geophys. Res.*, **99**, 16,897–16,905, 1994.

Crabb, J. A., and J. Latham, Corona from colliding drops as a possible mechanism for the triggering of lightning, *Q. J. R. Meteorol. Soc.*, **100**, 191–202, 1974.

Dawson, G. A., Pressure dependence of water-drop corona onset and its atmospheric importance, *J. Geophys. Res.*, **74**, 6859–6868, 1969.

Durand, E., *Electrostatique, Tome 1: Les distributions*, pp. 245–246, Masson et Cie, Paris, 1964.

Fage, A., Experiments on a sphere at critical Reynolds numbers, *Aero. Res. Comm., England, Rep. and Memo. No. 1766*, 20 pp., 1937.

Gaskell, W., A. J. Illingworth, J. Latham, and C. B. Moore, Airborne studies of electric fields and the charge and size of precipitation elements in thunderstorms, *Q. J. R. Meteorol. Soc.*, **104**, 447–460, 1978.

Gay, M. J., R. F. Griffiths, J. Latham, and C. P. R. Saunders, The terminal velocities of charged raindrops and cloud droplets falling in strong electric fields, *Q. J. R. Meteorol. Soc.*, **100**, 682–687, 1974.

Georgis, J. F., S. Chauzy, and S. Coquillat, Computed conditions of corona emission from two interacting raindrops, *Q. J. R. Meteorol. Soc.*, **121**, 1853–1866, 1995.

Georgis, J. F., S. Coquillat, and S. Chauzy, Onset of instability in precipitating drops submitted to horizontal electric fields, *J. Geophys. Res.*, **102**, 16,793–16,798, 1997.

Griffiths, R. F., A hypothetical cloud-seeding method for facilitating the occurrence of lightning, *J. Appl. Meteorol.*, **16**, 658–661, 1977.

Griffiths, R. F., and J. Latham, The emission of corona from falling drops, *J. Meteorol. Soc. Jpn*, **50**, 416–422, 1972.

Griffiths, R. F., and J. Latham, Electrical corona from ice hydrometeors, *Q. J. R. Meteorol. Soc.*, **100**, 163–180, 1974.

Griffiths, R. F., J. Latham, and R. L. Reed, Electrical corona from melting hailstones, *J. Meteorol. Soc. Japan*, **54**, 123–125, 1976.
Grover, S. N., and K. V. Beard, A numerical computation of the collision efficiency of raindrops colliding with aerosol particles and cloud drops, *Preprints Conf. on Cloud Physics*, pp. 120–121, Amer. Meteorol. Soc., Tucson, 1974.

Gunn, R., and G. D. Kinzer, The terminal velocity of fall for water drops in stagnant air, *J. Meteorol.*, **6**, 243–248, 1949.

Guyon, E., J. P. Hulin, and L. Petit, *Hydrodynamique Physique*, 506 pp., InterEditions Editions du CNRS, Paris/Meudon, 1991.

Happel, J., and H. Brenner, *Low Reynolds Number Hydrodynamics*, 553 pp., Prentice-Hall, Englewood Cliffs, N. J., 1965.

Kamra, A. K., and D. V. Ahire, Wind-tunnel studies of the shape of charged and uncharged water drops in the absence or presence of an electric field, *Atmos. Research*, **23**, 117–134, 1989.

Kamra, A. K., R. V. Bhalwankar, and A. B. Sathe, The onset of disintegration and corona in water drops falling at terminal velocity in horizontal electric fields, *J. Geophys. Res.*, **98**, 12,901–12,912, 1993.

Krehbiel, P. R., R. J. Thomas, W. Rison, T. Hamlin, J. Harlin, and M. Davis, GPS-based mapping system reveals lightning inside storms, *Eos Trans.*, **81**, 21–25, AGU, Washington, D. C., 2000.

Lane, W. R., and H. L. Green, The mechanics of drops and bubbles, in *Surveys in Mechanics*, edited by G. K. Batchelor and R. M. Davies, pp. 183–184, Cambridge Univ. Press, New York, 1956.

Latham, J., and C. D. Stow, Airborne studies of the electrical properties of large convective clouds, *Q. J. R. Meteorol. Soc.*, **95**, 486–500, 1969.

LeClair, B. P., A. E. Hamielec, and H. R. Pruppacher, A numerical study of the drag on a sphere at low and intermediate Reynolds numbers, *J. Atmos. Sci.*, **27**, 308–315, 1970.

MacCready, P. B., Jr., and A. Proudfit, Observations of hydrometeor charge evolution in thunderstorms, *Q. J. R. Meteorol. Soc.*, **91**, 44–53, 1965.

Magono, C., On the shape of water drops falling in stagnant air, *J. Meteorol.*, **11**, 77–79, 1954.

Marsh, S. J., and T. C. Marshall, Charged precipitation measurements before the first lightning flash in a thunderstorm, *J. Geophys. Res.*, **98**, 16,605–16,611, 1993.

Marshall, T. C., and W. P. Winn, Measurements of charged precipitation in a New Mexico thunderstorm: Lower positive charge centers, *J. Geophys. Res.*, **87**, 7141–7157, 1982.

Maxworthy, T., Experiments on the flow around a sphere at high Reynolds numbers, *J. Appl. Mech., Trans. A.S.M.E.*, **36**, 598–607, 1969.

MacGorman, D. R., and W. D. Rust, *The Electrical Nature of Storms*, 422 pp., Oxford Univ. Press, New York, 1998.

MacGorman, D. R., A. A. Few, and T. L. Teer, Layered lightning activity, *J. Geophys. Res.*, **86**, 9900–9910, 1981.

Macky, W. A., Some investigations on the deformation and breaking of water drops in strong electric fields, *Proc. R. Soc. London, Ser. A*, **133**, 565–587, 1931.

Milne-Thomson, L. M., *Theoretical Hydrodynamics*, 552 pp., MacMillan, Old Tappan, N. J., 1938.

Nolan, J. J., The breaking of water drops by electric fields, *Proc. R. Ir. Acad., Sect. A*, **37**, 28–39, 1926.

Oost, W. A., Flow distortion by an ellipsoid and its application to the analysis of atmospheric measurements, *J. Atmos. Oceanic Technol.*, **8**, 331–340, 1991.

Proctor, D. E., Regions where lightning flashes began, *J. Geophys. Res.*, **96**, 5099–5112, 1991.

Pruppacher, H. R., and K. V. Beard, A wind tunnel investigation of the internal circulation and shape of water drops falling at terminal velocity in air, *Quart. J. Roy. Met. Soc.*, **96**, 247–256, 1970.

Pruppacher, H. R., and R. L. Pitter, A semi-empirical determination of the shape of cloud and rain drops, *J. Atmos. Sci.*, **28**, 86–94, 1971.

Richards, C. N., and G. A. Dawson, The hydrodynamic instability of water drops falling at terminal velocity in vertical electric fields, *J. Geophys. Res.*, **76**, 3445–3455, 1971.

Takahashi, T., T. Tajiri, and Y. Sonoi, Charges on graupel and snow crystals and the electrical structure of winter thunderstorms, *J. Atmos. Sci.*, **56**, 1561–1577, 1999.

Taylor, G., The disintegration of water drops in an electric field, *Proc. Roy. Soc. London*, **A280**, 383–397, 1964.

Weinheimer, A. J., J. E. Dye, D. W. Breed, M. P. Spowart, J. L. Parrish, T. L. Hoglin, and T. C. Marshall, Simultaneous measurements of the charge, size, and shape of hydrometeors in an electrified cloud, *J. Geophys. Res.*, **96**, 20,809–20,829, 1991.

S. Chauzy, B. Combal, and S. Coquillat, Laboratoire d'Aérologie, UMR CNRS/UPS 5560, Observatoire Midi-Pyrénées, 14, Avenue Edouard Belin, 31400 Toulouse, France. (serge.chauzy@aeroobs-mip.fr; bcom@scot.cnes.fr; sylvain.coquillat@aeroobs-mip.fr)