Simulated Impacts of Thundercloud Charge Distributions on Sprite Halos Using a 3D Quasi-Electrostatic Field Model

Abstract

Sprite halos are transient luminous phenomena in the lower ionosphere triggered by tropospheric lightning. The effect of removed charge distributions on sprite halos has not been sufficiently discussed. A three-dimensional (3D) quasi-electrostatic (QES) field model was developed in this paper, including the ionospheric nonlinear effect and optical emissions. Simulation results show that, for a total charge of 150 C removed within 1 ms with different spatial distributions, higher altitudes of charge removal lead to stronger electric fields and increase sprite halos’ emission intensities. The non-axisymmetric horizontal distribution of charge affects mesospheric electric fields, and the corresponding scales and intensities of emissions vary with observation orientations. Considering the tilted dipole charge structure due to wind shear, the generated electric field and the corresponding position of sprite halos shift accordingly with the tropospheric removed charge, providing an explanation for the horizontal displacement between sprite halos and the parent lightning.

1. Introduction

Transient luminous events (TLEs) are upper atmospheric electrical discharges triggered by lightning discharges in tropospheric thunderstorms, including elves, sprite halos, sprites, starters, blue jets, and gigantic jets [1,2]. Understanding TLEs contributes to our knowledge of their potential influences on global circuits and climate, offering insights into their significant impact on the chemistry of the upper atmosphere and their capacity to disrupt subionospheric radio signals [3,4,5]. Since the first image of sprites captured by Franz et al. [6], TLEs have attracted extensive interest from researchers worldwide and have received increased attention for investigation over the past three decades.

Sprite halos and sprites are two types of TLEs that occur in the mesosphere/lower ionosphere; they are closely interlinked phenomena driven by the quasi-electrostatic (QES) fields produced by intense cloud-to-ground (CG) lightning discharges [7,8,9,10,11,12]. Sprites are one of the most morphologically complex types of TLEs, occurring at altitudes ranging from 40 to 90 km. They are composed of vertical filamentary streamers, referred to as sprite streamers [12], and are often accompanied by brief descending diffuse glows at higher altitudes, termed sprite halos. Sprite halos, or simply halos, are short-lived (~1–2 ms) and pancake-shaped descending glows at ~70–85 km altitudes, with lateral scales typically less than 100 km. The optical emissions of sprite halos are diffuse and less structured than sprites, and they often accompany or precede sprites, but they can also occur as an isolated event [13,14]. Sprites are predominantly produced by the quasi-electrostatic field generated by positive cloud-to-ground (+CG) lightning, while sprite halos are caused by both positive and negative CG lightning [15,16]. The typical brightness of sprites ranges from several to tens of mega-Rayleigh (MR) [17,18], and the brightness at the sprite streamer head can reach up to 1 giga-Rayleigh (GR) [19]. However, the average brightness of sprite halos is clearly lower than that of sprites. Kuo et al. [20] analysed TLEs based on the data set of the Imager of Sprites and Upper Atmospheric Lightning (ISUAL) experiment and found that the spatially averaged brightness of sprites and halos is 1.5 and 0.3 MR, respectively. Kuo et al. [21] reported five sprite halos without visible sprite streamers with a maximum brightness of 230–720 kR detected by the ISUAL satellite with an exposure time of ~29 ms. The results of optical observations may be affected by many uncertain factors such as band percentage uncertainties, altitude errors, and lightning contamination [18]. It is possible that in some sprite events, only sprite streamers were observed; sprite halos may indeed exist but they are not observed due to their diffuse emissions being too weak to be detected [22]. Due to their lower optical intensities and shorter durations, observations of the evolution of sprite halos need optical measurements with a higher time resolution than that of sprites.

Previous observational and simulation studies have concluded that the parent lightning parameters and the initial ambient conditions of the lower ionosphere (such as density inhomogeneities, gravity waves, and meteoritic dust) play key roles in the formation of TLEs. The vertical charge moment change (CMC), defined as the amount of charge removed during lightning discharge multiplied by the altitude from which it was removed, is a reliable metric to evaluate the capacity of sprite halo and sprite production. Another similar metric is the impulsive charge moment change (iCMC) [8,9,23,24,25]. For example, Hu et al. [23] analysed many sprite-producing strokes of summer thunderstorms and indicated that +CG lightning associated with CMC ≥ 1000 Ckm has >90% probability of producing sprites, while CMC ≤ 600 Ckm only has <10% probability of generating sprites. Yair et al. [26] found that the causative +CG lightning of winter thunderstorms has a mean CMC of 1400 ± 600 Ckm, with some extreme events exceeding 3500 Ckm. Based on satellite observation data, Chen et al. [27] found that the minimum critical CMC of +CG lightning capable of initiating sprites is as low as 63 Ckm, while the mean CMC value is ~1480 Ckm. Based on numerical simulations, Liu et al. [28] found that downward positive sprite streamers (produced by positive CGs) require a smaller electric field to initiate than negative sprite streamers. Qin et al. [22] estimated the threshold CMC of positive and negative sprites under typical night-time conditions are ~320 and ~500 Ckm, respectively. Qin et al. [29] further found that the minimum CMC for positive and negative sprite production could drop to ~200 and ~320 Ckm, respectively, under favourable ionospheric conditions, such as relatively high altitude of ionosphere and the pre-existing plasma inhomogeneities in the D-region ionosphere, which have been evidently observed by Qin et al. [30]. The presence of plasma inhomogeneities may be a necessary condition for the initiation of sprite streamers [30,31]. The lightning-induced electric fields driven by causative lightning discharges also influence sprite morphology. Qin et al. [32] demonstrated that carrot sprites with upward negative streamers are produced for lightning with CMC > 500 Ckm, while column sprites dominated by downward positive streamers occur for lightning with CMC < 500 Ckm under typical night-time conditions. The formation process of sprite halos is similar to that of sprites, except that the duration of the lightning-induced electric field perturbation exceeding the critical electric field is relatively short [22,33].

The existing research has discussed the roles of the lightning parameters such as polarity, peak current and multiplicity, the magnitude and duration of the continuing current, and the resultant vertical CMC. However, the influence of thundercloud charge distributions is still poorly understood. Although the data provided by radar, rocket, and balloon measurements of electric fields and radio remote sensing of lightning electromagnetic fields can be combined to estimate charge distributions and evolutions inside a thundercloud, it is still difficult to obtain a detailed spatial charge distribution of causative lightning [34,35,36,37]. Numerical simulation remains an effective tool in the study of TLEs, and many numerical models have been established so far. For example, Pasko et al. [8,9] developed a 2D QES model for previous sprite studies, which was later demonstrated to be suitable for simulating sprite halos. Barrington-Leigh et al. [13] presented a 2D electromagnetic pulse (EMP) model to study elves and halos. The simulation of sprite initiation, in particular the formation and evolution of sprite streamers, requires the plasma fluid models with higher resolution, such as the literature introduced by Pasko et al. [38], Luque and Ebert [39], Liu et al. [28,31,40], and Qin et al. [22,29,30,32,41]. The majority of existing models have been developed in a 2D cylindrical coordinate system to reduce the computational cost. Consequently, only symmetric charge distributions located on the axis of symmetry can be simulated, such as the dipole charge structure with spherical or disc-shaped charge distributions. In reality, the distribution of neutralised charges by lightning is complex, and the existence of vertical wind shear may ‘tilt’ the typically upright positive dipole; the resulting tilted dipole charge structure is considered favourable for CG lightning discharges [42,43,44,45]. To model non-symmetric charge distributions and tilted dipole charge structure, a fully 3D model is needed. For example, Asano et al. [46] presented a 3D electromagnetic pulse (EMP) model and simulated a horizontal channel in the model. Their results indicate an obvious correlation between the horizontal channel length and the horizontal displacement of potential sprite inception regions. Kabirzadeh et al. [47] established a 3D QES model for the first time and simulated the effects of the geomagnetic field dip angle on mesospheric electric fields. Zhang et al. [48] established a 3D QES model and simulated the effects of atmospheric gravity waves on the initiation and optical emission of sprite halos. Recently, Haspel et al. [49] established a simple 3D electrostatic model based on the method-of-images solution of Poisson’s equation to explore the influence of complex charge configurations on sprite initiation, but they ignored atmospheric conductivity and its time evolution, which is important for determining the possible regions of sprite initiation. Haspel et al. [50,51] improved their model by taking into account the time evolution of the free charge and the ambient conductivity, which varies in space but is still independent of time, and simulated the region of possible sprite inception in the mesosphere above winter thunderstorms under wind shear. Zhang et al. [52] employed a 3D QES model to investigate the influence of the ionospheric electron density profile and its inhomogeneities on sprite halos.

In summary, the effect of complex thundercloud charge distributions on sprite halos has not been sufficiently discussed in the literature and models. To evaluate the effect of thundercloud charge distributions on sprite halos, a 3D QES model is developed in this paper, based on the 2D QES heating model developed by Pasko [8]. It is noteworthy that this 3D QES model is a typical large-scale electrostatic coupling model; it cannot simulate the dynamics of sprite streamers, although its results can be used to evaluate the possible regions of sprite initiation in some literature [49,50]. The simulation results will facilitate a more profound comprehension of the correlation between charge distribution and sprite halos, potentially offering a plausible explanation for the horizontal displacement between sprites and their parent lightning. The rest of this paper is organised as follows: Section 2 introduces the 3D QES heating model. Section 3 presents the simulation results and discussion on the effect of the thundercloud charge distributions. Finally, the conclusions are given in Section 4.

2. Methods

2.1. The 3D QES Heating Model

A 3D QES heating model is developed to simulate thundercloud electrostatic interaction with the lower ionosphere and the resulting optical emissions, extending the 2D QES heating model proposed by Pasko [8] to 3D Cartesian coordinates, as shown in Figure 1. The simulation domain is 180 × 180 × 90 km³, which is divided into 1 × 1 × 2/3 km³ grid cells in the x, y, and z directions, respectively. The boundary condition is assumed to be perfectly conductive, and the topography is flat ground [8]. In the simulation, the charge density ρ and electrostatic potential φ satisfy the charge conservation equation and Poisson’s equation [48]:

$$\left\{\begin{array}{c}{\displaystyle \frac{𝜕\rho }{𝜕t}}-\nabla \sigma \cdot \nabla \phi +{\rho }_s\sigma /{\epsilon }_0=0\\ \hspace{0.33em}\nabla \cdot E=-{\nabla }^2\phi =(\rho +{\rho }_s)/{\epsilon }_0\end{array}\right.,$$

(1)

where E is the electrostatic field, and ρ_s is the thundercloud source charge density. ε₀ and σ are the permittivity and the atmospheric conductivity, respectively. The details of the updated equations for resolving Equation (1) are presented in Appendix A. The thundercloud charge model used in this paper is described in detail in Section 2.2. The self-consistent evaluation of the atmospheric conductivity σ is presented in Section 2.3.

2.2. Thundercloud Charge Model

The thundercloud charging and discharging process mainly includes two stages, namely the charge accumulation phase of the storm (τ_f) and the charge removal during a lightning discharge (τ_s). The total source charge is Q(t) = ∫ρ_s(x,y,z,t)dV, and the evolution of the charges is modelled using the formulation proposed by Pasko [8]:

$$Q(t)=\left\{\begin{array}{c}\hspace{1em}{Q}_0{\displaystyle \frac{tan\mathrm{h}(t/{\tau }_f)}{tan\mathrm{h}(1)}}\hspace{1em}\hspace{1em},0<t<{\tau }_f\\ Q_0\left[1-{\displaystyle \frac{tan\mathrm{h}((t-{\tau }_f)/{\tau }_s)}{tan\mathrm{h}(1)}}\right],{\tau }_f<t<{\tau }_f+{\tau }_s\end{array}\right.,$$

(2)

where Q₀, τ_f = 0.5 s, and τ_s = 1 ms are the total removed charge and the charging and discharging times, respectively. For a +CG discharge, the thundercloud positive charge is removed with a magnitude of Q₀ and the lightning duration is set to be 1 ms. The thundercloud charge density has a Gaussian spatial distribution:

$${\rho }_s(x,y,z,t)={\rho }_0(t)\mathit{exp}\left[-\left({\displaystyle \frac{(x-x_{\pm }{)}^2}{a^2}}+{\displaystyle \frac{(y-y_{\pm }{)}^2}{b^2}}+{\displaystyle \frac{(z-z_{\pm }{)}^2}{c^2}}\right)\right],$$

(3)

where ρ₀(t) = Q(t)/V = Q(t)/(4πabc/3), and (x_±, y_±, z_±) is the location of the positive and negative charge centres. By default, a dipole charge structure is considered with the positive and negative charge regions located at z₊ = 10 km and z₋ = 5 km altitude, respectively. The sizes of the charge distribution are parameterised by a, b, and c in the x, y, and z directions, respectively, as shown in Figure 2. In our 3D model, both symmetric and asymmetric charge distributions can be simulated.

2.3. Computation of Atmospheric Conductivity

Similarly to the atmospheric parameterization scheme proposed by Pasko [8], the atmospheric conductivity σ in Equation (1) consists of ion conductivity σ_i and electron conductivity σ_e (i.e., σ = σ_i + σ_e). Below 60 km, ion conductivity is dominant, expressed by σ_i = 5 × 10⁻¹⁴e^z/6km S/m [53]. Above 60 km, electron conductivity dominates the total conductivity and is calculated by σ_e = q_eN_eμ_e [54], where q_e is the elementary charge, μ_e is the electron mobility, and N_e is the electron density.

The electron mobility μ_e, as shown in Figure 3b, is a nonlinear function of the electric field E and the neutral density N, given by the following expression:

$$\left\{\begin{array}{c}{\mathit{log}}_{10}({\mu }_{e}N)={\sum }_{i=0}^2{a}_i{x}^i,{\displaystyle \frac{E{N}_0}{N}}\ge 1.62\times {10}^3 \mathrm{V}/\mathrm{m}\\ {\mu }_{e}N=1.36N_0\hspace{1em},{\displaystyle \frac{E{N}_0}{N}}<1.62\times {10}^3 \mathrm{V}/\mathrm{m}\end{array}\right.,$$

(4)

where x = log₁₀(E/N), a₀ = 50.97, a₁ = 3.026, and a₂ = 8.4733 × 10⁻². N₀ = 2.688 × 10²⁵ m⁻³ is the neutral density on the ground. The neutral density N at the actual altitude is taken from the NRLMSISE-00 atmospheric model [55], as shown in Figure 3a.

The background electron density profile N_e (m⁻³) at altitude z is calculated by N_e(z) = 1.43 × 10¹³exp(−0.15z′)exp[(β − 0.15)(z − z′)], where β = 0.5 km⁻¹ and z′ = 85 km are used as the typical electron density profile of the night-time ionosphere [56], as shown in Figure 3a. The electron density N_e follows the continuity equation:

$${\displaystyle \frac{d{N}_e}{dt}}=({\nu }_i-{\nu }_a)N_e,$$

(5)

where ν_i and ν_a shown in Figure 3c are the ionization and attachment rates, respectively.

The analytical expression of ionization rate ν_i is given by Papadopoulos et al. [57]:

$$v_i=7.6\times {10}^{-13}N{x}^{2}f\left(x\right)e^{-4.7(1/x-1)},$$

(6)

where f(x) = (1 + 6.3e^−2.6/x)/1.5, x = E/E_k, and E_k = 3.2 × 10⁶N/N₀ is the conventional electrical breakdown field.

The analytical form of the attachment rate ν_a is expressed by Equation (7):

$$v_a=\left\{\begin{array}{c}{\displaystyle \frac{N}{N_0}}{\sum }_{i=0}^2{a}_i{x}^i,{\displaystyle \frac{E{N}_0}{N}}\ge 1.628\times {10}^6 \mathrm{V}/\mathrm{m}\\ {\displaystyle \frac{N}{N_0}}{10}^{{\sum }_0^3{b}_i{y}^i},{\displaystyle \frac{E{N}_0}{N}}<1.628\times {10}^6 \mathrm{V}/\mathrm{m}\end{array}\right.,$$

(7)

where x = EN₀/N and y = log₁₀(EN₀/N), and the coefficients are a₀ = −2.41 × 10⁸, a₁ = 211.92, a₂ = −3.545 × 10⁻⁸, b₀ = −1073.8, b₁ = 465.99, b₂ = −66.867, and b₃ = 3.1970.

Given the 1 ms timescale of lightning discharges employed in our simulations, the primary chemical reactions are considered, which include the electron impact ionization of N₂ and O₂, as well as the electron dissociative attachment to O₂. Numerous slower chemical reactions are ignored, such as three-body detachment [40]. In addition, the influence of the geomagnetic field is neglected due to electron collisions dominating in the D-region ionosphere [8].

2.4. Computation of Optical Emissions

In a +CG discharge, the positive charge region within the thundercloud is neutralized, and the resulting electric field triggers nonlinear effects in the electric field of the lower ionosphere. Concurrently, the generated quasi-electrostatic fields will facilitate the acceleration of electrons. Once sufficient energy has been absorbed, the electrons will transition to an excited state. Upon returning to the ground state, the excited electrons release energy in the form of photons, thereby generating an optical emission. In this paper, the optical emission intensities I_k (in Rayleighs) generated by the nitrogen molecule (${\mathrm{N}}_{2}1\mathrm{P}$, ${\mathrm{N}}_{2}2\mathrm{P}$), nitrogen ion (${\mathrm{N}}_2^{+}\mathrm{M}$, ${\mathrm{N}}_2^{+}1\mathrm{N}$), and oxygen ion (${\mathrm{O}}_2^{+}1\mathrm{N}$) are calculated using the following expression:

$$I_k={10}^{-6}{\int }_L{A}_k{n}_{k}dl,$$

(8)

where the integral is taken along L (m), the line of sight through the emissions to an observer. In this paper, it is assumed that sprite halos are observed from the horizontal direction, with an integration path of 180 km corresponding to the horizontal scale of the entire simulation space. A_k is the radiation transition rate, and the number density n_k (m⁻³) in state k is calculated with its continuity equation as follows:

$${\displaystyle \frac{𝜕n_k}{𝜕t}}=-{\displaystyle \frac{n_k}{{\tau }_k}}+{\sum }_m{n}_m{A}_m+v_k{N}_e,$$

(9)

where ${\tau }_k={\left[A_k+{\alpha }_1{N}_{N_2}+{\alpha }_2{N}_{O_2}\right]}^{-1}$ is the lifetime of state k, α₁ and α₂ are the quenching rates, and $N_{N_2}$ and $N_{O_2}$ are the number densities of N₂ and O₂, respectively. The optical coefficients are given in

Table 1. Optical coefficients.

Band	${\mathit{\alpha }}_1\left({\mathbf{m}}^3{\mathbf{s}}^{-1}\right)$	${\mathit{\alpha }}_2\left({\mathbf{m}}^3{\mathbf{s}}^{-1}\right)$	${\mathit{A}}_{\mathit{k}}\left({\mathbf{s}}^{-1}\right)$
${\mathrm{N}}_{2}1\mathrm{P}$	${10}^{-17}$	0	$1.7\times {10}^5$
${\mathrm{N}}_{2}2\mathrm{P}$	0	$3\times {10}^{-10}$	$2\times {10}^7$
${\mathrm{N}}_2^{+}\mathrm{M}$	$5\times {10}^{-10}$	0	$7\times {10}^4$
${\mathrm{N}}_2^{+}1\mathrm{N}$	$2\times {10}^{-10}$	0	$1.4\times {10}^7$
${\mathrm{O}}_2^{+}1\mathrm{N}$	$4\times {10}^{-10}$	0	$8.5\times {10}^5$

[58] (p. 119). The term $\sum _m{n}_m{A}_m$ represents the contribution of the cascade from higher states to n_k, and the cascading process from N₂ 2P to N₂ 1P is considered here. v_kN_e denotes the excitation progress of electrons.

The optical excitation rate ν_k, which is shown in Figure 3d, as a nonlinear function of the electric field [59]:

$$\mathit{log}({\displaystyle \frac{v_k{N}_0}{N}})={\sum }_{i=0}^3{a}_i{x}^i,$$

(10)

where x = log₁₀(EN₀/N), and the approximation coefficients a_i for different optical bands are given in

Table 2. Approximation coefficients.

Band	a0	a1	a2	a3
${\mathrm{N}}_{2}1\mathrm{P}$	−1301.0	563.03	−80.715	3.8647
${\mathrm{N}}_{2}2\mathrm{P}$	−1877.1	814.70	−117.46	5.6554
${\mathrm{N}}_2^{+}\mathrm{M}$	−1760.0	724.70	−99.549	4.5862
${\mathrm{N}}_2^{+}1\mathrm{N}$	−1802.4	750.91	−104.28	4.8508
${\mathrm{O}}_2^{+}1\mathrm{N}$	−2061.1	870.25	−122.43	5.7668

[60,61].

2.5. Model Validation

Compared with the 2D QES model developed in a 2D cylindrical coordinate system, the 3D QES model is built in 3D Cartesian coordinates, thereby enabling the simulation of any realistic charge distributions. In order to validate the accuracy of the calculations performed by the 3D QES model, a symmetric dipole charge structure is considered for comparison with the 2D QES model. The positive and negative charge centres were located at 10 km and 5 km altitude, respectively. By default, each thundercloud charge density is assumed to have a Gaussian spherical distribution with a radius of 3 km, and all other parameters of the model are as described above. For the case of 200 C charge removed within 1 ms from an altitude of 10 km, the horizontal simulated space of the 3D model is set to 120 km here. Figure 4 compares the spatial distribution of the charge density ρ at the end of the lightning discharge (marked as t = 1 ms, the same as below) simulated by the 3D QES model with those calculated by the 2D QES model [8]; a good agreement can be seen between them.

To further validate the calculation accuracy of the optical emissions, Figure 5 illustrates the altitude profiles of the instantaneous optical emission intensities for varying charge removal and for different optical bands. The altitude profiles of the instantaneous optical emission intensities of the 1st positive band (red) N₂ (N₂ 1P) calculated by our 3D QES model, for different removed charges ranging from 50 C to 250 C, as shown in Figure 5b, are also in good agreement with those obtained by Pasko [8] using the 2D QES model. Figure 5c,d compares the altitude profiles of the instantaneous optical emission intensities corresponding to different optical bands for the Q = 200 C case between the results of Pasko [8] and our 3D model. As can be seen, the agreement between the two results is still excellent. There is a slight difference between the two results, which can be attributed to the different initial electron density profiles adopted in the respective models. Details on the model formulation, parameter set, and model accuracy verification of the 3D QES model introduced above can be found in [48].

3. Results and Discussion

In this section, the impacts of thundercloud charge distribution factors on sprite halos will be modelled using the 3D QES heating model. These factors include the altitude of the removed (neutralised) charge, the horizontal scale of the charge region, and the horizontal shift of the positive and negative charge regions. Unless otherwise specified, the model parameters are set to the default values as described above.

3.1. Effect of the Altitude of the Thundercloud Charge

Figure 6 shows the vertical profiles of the electric field E, the normalized electric field E/E_k, the electron density N_e, and the instantaneous optical emission intensity I of N₂ 1P above the charge centre after a +CG discharge at t = 1 ms, considering a total charge Q = 150 C removed from different altitudes z₊ = 5~20 km within 1 ms. It should be noted that under certain relatively rare conditions, a thundercloud charge can accumulate at ~20 km altitude [8]; thus, the case of z₊ = 20 km is used as an upper limit. For the majority of the calculations presented in this paper, the results are shown at the end of the discharge (t = 1 ms), because this is when the electric field, the electron density change, and the instantaneous optical emission intensity just above the charge centre reach their maximum [8,48]. Due to the negative charge remaining unchanged during a +CG discharge, which is essentially independent of the electric field appearing at mesospheric altitudes, the altitude of the negative charge region is set to 5 km by default.

Figure 6a,b show the vertical profiles of electric field E and normalized electric field E/E_k. The downward developing trend and amplitudes of mesospheric electric fields increase with the increase in CMC. Figure 6c shows the nonlinear change in electron density. The higher the electron density, the greater the atmospheric conductivity, which will enhance the local electric field, but it will also shield the internal electric field. Figure 6d shows that, similar to the electric field trend, the optical emission intensities increase with the increase in CMC. These results are similar to the calculated results of removing different charges from the same altitude, as documented in Pasko et al. [8,9]. For the case of z₊ = 5 km and 10 km, the maximum of the electric field and optical emission intensity are located at an altitude of ~78 km, which corresponds to the inception altitude of sprite halos, and is in good agreement with observations [13,14,62]. For the case of z₊ = 15 km and 20 km, the descending altitudes of the maximum electric field are obvious, while the electric field value in the upper space remains low, due to the nonlinear change in electron density. In general, for the same value Q of the removed (neutralized) charge, the higher altitude of the charge centre, i.e., the greater CMC, resulting in higher values of the electric field, electron density change, and optical emission intensity in the mesosphere.

As in Figure 6, but for a CMC of 1500 Ckm, Figure 7 further shows the simulation results for different charge altitudes. Positive charges with charges of 300, 150, 100, and 75 C centred at 5, 10, 15, and 20 km, respectively, are removed within 1 ms. At the end of the +CG lightning discharge (t = 1 ms), the remaining negative charges located at 5 km altitude are −300, −150, −100, and −75 C, respectively.

Figure 7a,b show that, even under the same CMC conditions, the electric field amplitude in the low ionosphere region increases slightly with higher charge removal altitudes. It is considered that the reduced distance between the charge centre and the ionosphere leads to more pronounced quasi-electrostatic heating in the D-region ionosphere and a consequent increase in the electric field amplitude [47]. Correspondingly, the electron density change and optical emission intensity in Figure 7c,d are positively correlated with the charge altitude. This result is consistent with the simulation results of Tong et al. [63], in which they found that the electron density change and optical emission intensity are slightly higher at an altitude of 10 km compared to 5 km under the same CMC conditions. The results shown in Figure 7 demonstrate that the lightning CMC is not the only determinant of the electric field in the D-region ionosphere; the altitude of the removed (neutralized) charge also affects the electric field amplitude and the corresponding optical intensity of sprite halos to some extent.

3.2. Effect of the Horizontal Scale of the Charge

As described in Section 2.2, the geometry of a thundercloud charge can be controlled by the values of a, b, and c. Specifically, when a = b = c = 3 km, a = b > c = 3 km, and a ≠ b are considered, the charge geometries are spherical, oblate, and ellipsoid, respectively. To further investigate the effects of the horizontal scale of the charge region, simulations were conducted in the 3D model with varying values of a and b. By default, the vertical dimension c is set to 3 km, the positive charge centre is located at (0 km, 0 km, 10 km), and the negative charge centre is placed at (0 km, 0 km, 5 km). For a +CG discharge, a Gaussian-distributed charge of 150 C is neutralized within 1 ms.

Figure 8 depicts the vertical profiles of the electric field E, normalized electric field E/E_k, electron density N_e, and instantaneous optical emission intensity I of N₂ 1P above the thundercloud at the end of the discharge (t = 1 ms), considering different horizontal scales a = b, ranging from 3 to 40 km. Figure 8a,b show that in the D-region ionosphere, the electric field amplitude increases with the horizontal scale of the charge region for the same CMC when the horizontal scale is less than 10 km. In contrast, when the horizontal scale exceeds 10 km, the electric field amplitude decreases as the horizontal scale expands. Figure 8c,d indicate that the variations in electron density change and optical emission intensity with a horizontal scale in the D-region ionosphere are similar to those of the electric field. These results align with some simulation results reported by Pasko [8], who compared the charge sizes with horizontal scales of 3 km and 30 km, under the same CMC conditions. Their simulation results indicate that the electric field generated by the discharge process in the charge region with a 30 km horizontal scale is weaker than that with a 3 km horizontal scale. Considering two spherical charge regions in the 2D model, Pasko [8] also compared the results of 3 km and 6 km horizontal scales and found no obvious influence on the results. Haspel et al. [49], in their simplified 3D model, compared the point charge and disk-shaped charge region with horizontal scales of 4 km and 8 km, and the results also did not find an obvious influence. However, their model neglected atmospheric conductivity and ionospheric nonlinear effects.

Figure 9 shows the relationship between the normalized electric field at t = 1 ms and the horizontal scale (a = b) of charge density distributions. Our results show that, under the same CMC conditions, the horizontal scale of the charge region significantly influences the amplitude of the electric field, consequently impacting the optical emission intensities of sprite halos.

It should be noted that the charge distribution considered above remains axisymmetric and can be simulated by the 2D model. Subsequently, we will conduct simulations for non-axisymmetric scenarios, which require accurate modelling with the 3D model. For the case of a ≠ b, Figure 10 shows the simulation results with a = 3 km and b = 30 km. For comparative analysis, the results of three different geometries are compared here: (i) a = b = 3 km, (ii) a = 3 km and b = 30 km, and (iii) a = b = 30 km. The simulation results indicate that the electric field, electron density change, and optical emission intensity in case (ii) are intermediate, being lower than those in case (i) but greater than those in case (iii).

Figure 11 further shows the spatial distribution of charge density, normalized electric field E/E_k, and the time-averaged optical emission intensity I of N₂ 1P over a 2 ms period, considering charge geometries with a = 3 km, b = 30 km, and c = 3 km. Figure 11a,b show that the vertical distribution of charge density is significantly affected, particularly below 60 km, due to the extension of the charge region along the y-direction. Above 60 km, the vertical and horizontal distribution of the normalized electric field E/E_k in the cross-section of the domain at z = 78 km, shown in Figure 11c–e, becomes asymmetric because of the non-axisymmetric structure of the removed charge, and the electric field obviously extends along the y-direction. The exposure time of the optical instruments results in the captured sprite halos being optical images averaged over their main duration (2 ms), reflecting the overall spatial distribution of sprite halos [48], as shown in Figure 11f,g. The horizontal observation results along the x-axis and y-axis directions are presented, respectively. Observation along the x-axis reveals a broad horizontal distribution of optical emission with a relatively low amplitude. In contrast, observation along the y-axis shows a more concentrated emission with higher intensity. This indicates that the morphology and intensity of the captured sprite halos’ optical emission may vary with the orientation of observation.

3.3. Effect of the Horizontal Shift of Positive and Negative Charges

As mentioned in Section 1, thunderclouds in reality are frequently affected by wind shear, resulting in a tilted dipole charge structure. As depicted in Figure 12, the upper positive charge is laterally offset from the negative charge [42,43]. To evaluate the impact of a tilted dipole charge structure on sprite halos, our 3D model uses the following parameter set here: a Gaussian-distributed charge of 150 C with scale a = b = c = 3 km, removed from 10 km within 1 ms. The positive charge centre is located at (0 km, d km, 10 km) and the negative charge centre is situated at (0 km, 0 km, 5 km). The distance d denotes the horizontal displacement of the positive charge relative to the negative charge along the y-axis direction.

As in Figure 11 but for the tilted dipole charge structure with d = 20 km, Figure 13 shows the spatial distribution of the charge density, the normalized electric field E/E_k, and the time-averaged optical emission intensity of N₂ 1P over a 2 ms period. The red solid line and dotted line mark the centre positions of the negative charge and positive charge, respectively. The vertical distribution of the spatial charge density in Figure 13a is asymmetric due to the tilted dipole charge structure. The vertical and horizontal distributions of the normalized electric field E/E_k in Figure 13b,c align with the removed charge. The electric field near the model boundary becomes asymmetric, attributed to the influence of the ideal conductor boundary condition used in this model. With the displacement of the mesospheric QES field position, the location of the excited optical emission shifts, depending on the observation orientation. Figure 13d,e present the time-averaged optical emission intensity of N₂ 1P over a 2 ms period, observed along the x-axis and y-axis directions, respectively. Variations in observation orientations do not affect the shape and intensity of sprite halos’ optical emission. However, the position of optical emission observed along the x-axis direction is shifted by 20 km from the centre of the tropospheric negative charge region, aligning with the offset of the positive charge region.

Figure 14 further illustrates that the corresponding positions of sprite halos shift in response to neutralized charges, considering different horizontal offset distances (d = 0, 10, 20, and 30 km). It is widely recognized that sprites and sprite halos are triggered by the strong QES field that results from the neutralized thundercloud charges by intense lightning. The process of charge removal can be viewed as the addition of an identical charge of the opposite sign at the exact position of the removed charge; thus, the position of the mesospheric electric field and optical emission depend on the position of the neutralized charge [8,9]. These results can provide a possible explanation for the observed horizontal displacement of sprites from their parent lightning [64,65,66]. Sato et al. [66] discovered that the locations of sprite emissions are distinctly displaced by 8–20 km from the peak positions of the parent lightning emissions and proposed two possible occurrence conditions to explain these observations. One possible occurrence condition is that the pre-existing conductivity inhomogeneity at sprite altitude could readily excite sprites even under a weak QES field displaced from the charge centre and lightning; such conductivity inhomogeneity may be closely related with the inhomogeneity of electron density, as confirmed by Qin et al. [30]. Another possible condition is that the centre of the neutralized charge deviates from the location of the lightning strike point due to the horizontal lightning channel. The strongest QES field would be generated above the centre of the neutralized charge, as confirmed by our 3D model results. Consequently, sprites would be initiated just above the charge centre but with a horizontal displacement from the lightning strike point.

4. Conclusions

In this study, a 3D QES heating model was developed in 3D Cartesian coordinates, enabling the simulation of various charge distributions. The accuracy of the 3D model was validated by comparing its results with those documented in the scientific literature. The impacts of thundercloud charge distributions on sprite halos were simulated using the 3D QES model. The main conclusions are summarized as follows:

(1)

For the case of removing the same charge from different altitudes, higher altitudes are associated with greater CMC, resulting in higher values of electric fields and optical emission intensities in the mesosphere. Even under identical CMC conditions, higher charge removal altitudes will lead to increased mesospheric electric fields and optical emission intensities in the mesosphere to some extent.

(2)

For a Gaussian oblate symmetric charge density distribution with the scale of a = b = 3~40 km, under the same CMC conditions, the electric field amplitude increases with the enlargement of the charge region’s horizontal scale up to 10 km. Beyond this scale, the electric field amplitude decreases as the horizontal extent of the charge region expands further. For a Gaussian ellipsoid nonaxisymmetric charge density distribution with scale of a ≠ b, the spatial distribution of the electric field becomes asymmetric. The morphology and intensity of the time-averaged optical emission of sprite halos also change with different observation orientations.

(3)

The effect of the horizontal shift of the positive charge relative to the negative charge due to wind shear was also simulated. For a tilted dipole charge structure, the mesospheric electric fields and the corresponding positions of sprite halos’ optical emissions shift accordingly with the tropospheric removed charge. These findings provide an explanation for the observed horizontal displacement of sprite emissions relative to their parent lightning.

Overall, the spatial distribution of a thundercloud charge removed by lightning has an important effect on the mesospheric electric field and the optical emission of sprite halos. A higher-resolution plasma fluid model, which includes more charged species and more chemical reactions, can also simulate sprite halos and even sprite streamers, providing more detailed information. However, the computational cost of the corresponding 3D model will be very high. Our 3D QES model constructed in this paper, based on the classical 2D QES model, has relatively lower spatial resolution and increased computational demands, which is suitable for the simulation of sprite halos, while the simulation of sprite streamers is beyond the scope of this model. The purpose of this paper is to evaluate the effect of thundercloud charge distributions on sprite halos, ignoring the influence of ionospheric inhomogeneities on sprite halos. The presence of plasma irregularities has been proven to be a necessary condition for the initiation of sprite streamers [30,31], which may also influence the existence of sprite halos [52]. The charge configurations of two adjacent lightning flashes in the neighbouring thunderstorm cells discussed by Haspel et al. [49,50] and Zhang et al. [67] could also affect sprite halos, which can also be modelled by our 3D QES model. Furthermore, the potential impact of terrain topography on mesospheric electric fields and optical emissions has not been discussed in the existing literature and models. However, the relevance of these effects is beyond the scope of this paper and requires further investigation.