Articles | Volume 13, issue 11
Model description paper
12 Nov 2020
Model description paper |  | 12 Nov 2020

Modeling lightning observations from space-based platforms (CloudScat.jl 1.0)

Alejandro Luque, Francisco José Gordillo-Vázquez, Dongshuai Li, Alejandro Malagón-Romero, Francisco Javier Pérez-Invernón, Anthony Schmalzried, Sergio Soler, Olivier Chanrion, Matthias Heumesser, Torsten Neubert, Víctor Reglero, and Nikolai Østgaard

We describe a computer code that simulates how a satellite observes optical radiation emitted by a lightning flash after it is scattered within an intervening cloud. Our code, CloudScat.jl, is flexible, fully open source and specifically tailored to modern instruments such as the Modular Multispectral Imaging Array (MMIA) component of the Atmosphere–Space Interactions Monitor (ASIM) that operates from the International Space Station. In this article, we describe the algorithms implemented in the code and discuss several applications and examples, with an emphasis on the interpretation of MMIA data.

1 Introduction

Lightning flashes emit intense optical radiation that, after passing through the cloud layer, can be observed from space. This has been applied to characterize the geographic and seasonal distributions of lightning (Christian et al.2003; Cecil et al.2014), to describe and forecast the development of particular thunderstorms (Peterson2019) and to investigate phenomena related to lightning, such as transient luminous events (TLEs) (Ebert et al.2010; Pasko et al.2012) and terrestrial gamma ray flashes (TGFs) (Dwyer et al.2012; Neubert et al.2020).

Several satellite-borne instruments targeting lightning optical emissions have been commissioned over the past few years and new ones are planned for the near future. The first such device was the Optical Transient Detector (OTD) (Christian et al.2003), which was active from 1995 to 2000 onboard the Microlab-1 satellite. The OTD was a model for the longer-lasting Lightning Imaging Sensor (LIS) that operated from the Tropical Rainfall Measurement Mission (TRMM) satellite from 1997 to 2015 (Boccippio et al.2002; Mach et al.2007; Zhang et al.2019; Blakeslee2019). A second, identical LIS device has been operating on the International Space Station (ISS) since 2017 (Blakeslee et al.2016).

After the success of space-based lightning observations from the OTD and LIS, it became desirable, for operational weather prediction, to access real-time lightning data over a large geographic area. This was achieved by the Geostationary Lightning Mapper (GLM) instruments onboard the Geostationary Operational Environmental Satellites (GOES) GOES-East and GOES-West (Goodman et al.2013) and the Lightning Mapping Imager (LMI) on the FENGYUN (FY-4) satellites. Also, the Meteosat Third Generation (MTG) constellation of satellites will be equipped with lightning-dedicated instruments named Lightning Imagers (LI).

All of these instruments focus only on the detection of lightning flashes and rely on cameras with wavelength filters centered around 777.4 nm, which is a strong emission line of lightning. Their integration times are around 2 ms, and their pixel sizes cover a few kilometers on the ground.

Other instruments target not only lightning but also additional phenomena associated with lightning. The Imager for Sprites and Upper Atmospheric Lightning (ISUAL) payload on the FORMOSAT-2 satellite contained limb-pointing cameras, spectrophotometers and a photometer array designed for the observation of TLEs between the cloud tops and the lower ionosphere (Chern et al.2003). Although not its main focus, ISUAL detected routinely optical emissions from lightning. The Global Lightning and sprIte MeasurementS (GLIMS) mission, installed on the ISS, was equipped with imagers and photometers to observe lightning and TLEs towards the nadir (Sato et al.2011; Adachi et al.2016).

Our work is mainly motivated by another instrument that also has TLEs among its science objectives: the Modular Multispectral Imaging Array (MMIA), which is part of the Atmosphere–Space Interactions Monitor (ASIM) currently onboard the ISS (Chanrion et al.2019; Neubert et al.2020). Pointing towards the nadir of the ISS, MMIA is composed of two cameras and three photometers. The cameras provide images filtered at around 777.4 and 337 nm, with the latter wavelength being characteristic of nonthermal discharges that often accompany lightning flashes or are part of them. The camera pixels have a footprint on the ground of around 400 m×400 m at nadir and an integration time of 83.3 ms (yielding 12 frames per second). Of the three photometers, two are sensitive to the same wavelengths as the cameras whereas the third detects an ultraviolet band between approximately 180 and 230 nm. The photometers operate at 105 samples per second. The high spatial resolution of the MMIA cameras and the high temporal resolution of its photometers enable detailed observations of lightning flashes that are, nevertheless, conditioned by scattering within the intervening clouds. This is why this instrument demands an improved understanding of this scattering.

In order to interpret the observations of the instruments mentioned above, one must understand how the optical radiation emitted by a lightning flash is scattered inside a cloud, sometimes with an intricate shape. For example, the strong convection generally associated with active thunderclouds causes overshooting cloud turrets, and the light emerging from them is partially reflected back upwards from the lower cloud deck, as illustrated by Chanrion et al. (2017). Because the optical properties of a lightning flash as observed from above are conditioned by the properties of the cloud where it originated, they may serve to characterize thunderstorms, as recently proposed by Peterson (2019), and also to investigate the time structure of lightning discharges (Peterson and Rudlosky2019).

Quantitative models of optical lightning radiation passing through clouds have been developed in the past by researchers such as Thomason and Krider (1982), Koshak et al. (1994) and Light et al. (2001), but these models have been limited to relatively simple cloud geometries and only considered how the radiation is delayed, disregarding image composition.

Here, we present CloudScat.jl, a code that overcomes these two limitations and serves as an open-source tool for the interpretation of space-based lightning observations. CloudScat.jl deals efficiently with complex cloud shapes defined by the composition of elemental solids such as spheres, cylinders and cones that can be arbitrarily placed and deformed. It also predicts the images recorded from observing devices located anywhere relative to the lightning source.

2 Model and numerical algorithm

Our algorithm is partly based on those described by Thomason and Krider (1982) and Light et al. (2001), but it has been extended with additional physical processes and variance-reducing methods. It simulates the propagation of a set of photon packets that interact with a background that includes both cloud hydrometeors and the molecular components of air. As the types of processes that we consider do not change the photon's wavelength, we only consider a monochromatic population of photons with a certain wavelength λ in each simulation. Our model is applicable for wavelengths from around 200 nm to 1 µm.

Every photon packet propagates asynchronously and independently of all other packets; therefore, we describe the dynamics of a single packet in the following. We will also sometimes use the term “photon” as a shortened way of referring to a computational particle, which should, in reality, be considered as an ensemble of physical photons. Each of these photon packets or computational particles is assigned a statistical weight to account for the different number of physical particles that it represents.

We consider two types of processes: scattering processes that modify the propagation direction of the photons are modeled as a series of discrete, stochastic events; and molecular absorption is treated as a continuous process that dampens the statistical weight of the packet as it propagates.

Figure 1Sketch with the notation employed in this work. We account for two types of processes affecting photon propagation. Mie and Rayleigh scattering are modeled as discrete scattering events that take place at locations r1,,ri-1,ri,ri+1,. Because Mie scattering dissipates part of the photon's energy, in the ith scattering event the particle's statistical weight changes from wi- to wi+. The absorption by background molecular components of air modifies the particle's statistical weight as it propagates between consecutive scattering events. This changes the particle's weight, e.g., from wi-1+ immediately after scattering i−1 to wi- immediately before scattering i.


Regarding the discrete scattering events, our notation is that the ith collision takes place at location ri and time ti (see Fig. 1). Immediately before the collision the photon weight is wi- and its velocity is cμi−1, where c is the speed of light, and μi−1 is a unitary vector. The collision instantaneously changes these magnitudes to wi+, cμi. Due to background absorption, the weight decays through the propagation path of a photon so wi+1-wi+, with equality holding in the absence of absorption. All photons are initialized at t0=0 with weight w0+=1, with random locations r0 uniformly distributed in one straight segment and with an isotropic distribution of μ0. All of our output results are divided by the number of initial photons N, so intensities must be interpreted as intensity per physical photon in the source (equivalently, we could initialize all weights to 1∕N, but this would complicate our description).

As a photon advances, it collides with background scatterers at a total rate per unit time

(1) ν ( r ) = c p = 1 m σ ¯ p ( r ) n p ( r ) ,

where c is the speed of light, σ¯p(r) is the locally averaged cross section and np(r) the number density of the scatterers involved in the process p at location r. The cross sections are averaged locally because a full spectrum of scattering particles with different cross sections may exist around a given point (due to, for example, the different radii of cloud droplets). To sample the interaction times, we employ the null collision method, whereby we find one rate νT such that νT>ν(r) everywhere inside the computational domain. The time of the ith scattering, including null collisions, is

(2) t i = t i - 1 + ν T - 1 | log ( ξ i - 1 ) | ,

where ξi are independent random variables uniformly distributed in the interval [0,1). The type of interaction is decided with probabilities allotted proportionally to each individual rate and with the probability 1-ν(r)/νT for a null collision.

2.1 Mie scattering

In Mie scattering a photon packet with a wavelength λ interacts with a cloud droplet, modeled as a dielectric sphere with a radius Rλ (van de Hulst1981). The cross section of this interaction is

(3) σ Mie = Q ext ( R ) π R 2 ,

where Qext is called the “extinction efficiency” and is generally close to 2 (the fact that this is not 1, as the naive application of geometrical optics would suggest, is called the “extinction paradox” and is discussed at length in studies such as Bohren and Huffman1983, p. 107). The incident radiation is partially absorbed by the droplet, letting a fraction ω0 (called the “single-scattering albedo”) of the energy be reemitted. In our code we account for this absorption by updating the weight from wi- to wi+ as follows. First we compute a temporary weight wi as

(4) w i = ω 0 w i - .

If wi is larger than a prescribed wmin<1, we set wi+=wi; otherwise, we set wi+=1 with the probability wi, and wi+=0 with the probability 1-wi (i.e., in the second case, we discard the photon). Thus, we avoid computations on photons with so little weight that they do not affect the observations, but we ensure that the expected value of wi+ is wi; thus, the simulation stays unbiased (Iwabuchi2006). The default value of wmin is 10−2.

The outgoing electromagnetic wave after Mie scattering is sampled by deflecting the direction of the incident photon. The scattered intensity per unit solid angle at a given direction forming an angle θ with the incident direction is conventionally named the “phase function”, and here, following Thomason and Krider (1982), we approximate it using the Henyey–Greenstein phase function,

(5) p HG ( θ ) = 1 - g 2 ( 1 + g 2 - 2 g cos θ ) 3 / 2 .

By convention, any phase function p(θ) is normalized such that

(6) 0 π p ( θ ) sin ( θ ) d θ = 4 π .

The Henyey–Greenstein phase function is parameterized by g, known as the “asymmetry parameter” which is the average of the cosine of the scattering angle,

(7) g = 1 4 π 0 π p ( θ ) cos ( θ ) sin ( θ ) d θ .

Given an arbitrary phase function obtained, for example, by an accurate solution of the Mie scattering problem, one can compute an asymmetry parameter through Eq. (7) and then approximate the accurate phase function using the Henyey–Greenstein function with the same g. One advantage of this is that it is simple and efficient to sample the Henyey–Greenstein phase function by inverting the cumulative probability distribution of cos θ. The result is that one can draw a random number s uniformly distributed between −1 and 1 and then set

(8) cos θ = 1 2 g 1 + g 2 - 1 - g 2 1 + g s 2 .

To compute the three parameters of Mie scattering, Qext, ω0 and g, which depend on the wavelength of interest and the droplet radius, we solve the Mie problem with the open-source MieScatter.jl code (Wilkman2013). The refractive index of the medium is a required input for the Mie solver, and the code allows the user to specify this from an input file. The uppermost layers of a thundercloud are above the glaciation line, so the water content is predominantly in ice form. Nevertheless, the optical constants of water and ice are rather similar, with the main difference being that absorption in ice is negligible in the range of wavelengths from 200 to 400 nm (Warren and Brandt2008); therefore, for the sake of simplicity, we mostly employ the water optical constants, which we interpolate from Hale and Querry (1973). In most cases reported here, the difference between ice and water is rather small, and we will only consider ice clouds when studying far-ultraviolet radiation in Sect. 3.2.

In general, the droplet size spectrum of a convective cloud is inhomogeneous and assuming a constant radius everywhere may be an oversimplification. A somewhat more accurate model approximates the full spectrum of radii with a single, effective radius R but allows it to depend on the scattering position. In that case, the average cross section as well as the phase function depend on the location of the scattering event. This can optionally be implemented in the code by using a parameterized dependence of the Mie scattering features as a function of R. After inspecting plots of accurate solutions and aiming for simplicity, we settled on these expressions to which we do not attribute any physical meaning:


where a, b, g0 and R0 are fitting parameters. These functions are simple but still capture the dependence of the scattering parameters on R. The user of our code can ask to solve the Mie scattering problem for a range of radius 1µm<R<100µm, fit the results to the above functions and use the result to build an inhomogeneous model for the cloud.

2.2 Rayleigh scattering

As optical radiation travels inside the cloud or in the path to a space platform, it is also affected by interactions with the molecular components of air. For lightning emissions, this is generally a negligible correction, but we included these interactions in our code for the sake of completeness and to eventually allow for additional applications.

To compute the Rayleigh scattering cross section, we follow Bodhaine et al. (1999), which gives

(10) σ R = 24 π 3 F ( air ) λ 4 N 0 2 n s 2 - 1 n s + 2 2 ,

where N0 is the number density of air at standard temperature and pressure (the Loschmidt constant), ns is the refractive index of air under those conditions and F(air) is a so-called “depolarization term” or “King factor”, which corrects the Rayleigh theoretical estimates to account for the asymmetry of scattering molecules. We compute ns as suggested by Peck and Reeder (1972):

(11) 10 8 ( n s - 1 ) = 8060.51 + 2 480 990 132.274 - λ - 2 + 17 455.7 39.32957 - λ - 2 ,

where λ is expressed in micrometers. We calculated F(air) from the formulas provided by Bates (1984), as also compiled by Bodhaine et al. (1999) with the air composition from (Haynes2016, p. 14–19).

The phase function for Rayleigh scattering is

(12) p R ( θ ) = 3 4 ( 1 + cos 2 θ ) ,

and, after inverting the cumulative probability distribution of cos θ, one finds that the phase function can be sampled by setting

(13) cos θ = w 1 / 3 - w - 1 / 3 , w = 2 s + ( 4 s 2 + 1 ) 1 / 2 ,

where s is a random number uniformly distributed in [-1,1).

The collision rate reads νR=cN(z)σR, where N(z) is the number density of air at a given altitude z. Because Rayleigh scattering is a small correction on top of the much more relevant Mie scattering, we implemented a simple model where the air density decreases exponentially with altitude:

(14) N ( z ) = N 0 exp ( - z / H ) ,

where H is the scale height, which is configurable by the user and defaults to 7.2 km.

2.3 Background absorption

Some wavelengths of interest for the observation of lightning are significantly absorbed in air. For example, the band between 180 and 230 nm that is probed by one of the MMIA photometers is affected by molecular absorption by ozone (Molina and Molina1986) and by the Schumann–Runge band of molecular oxygen (Minschwaner et al.1992). Although, in principle, one could incorporate this into our code by an additional discrete process that randomly discards some of the computational particles, it is more efficient, in terms of a lower variance in the output for the same computational time, to continuously adapt the statistical weight of the particles along their propagation.

For a single absorbing species, in the propagation between two collision points ri and ri+1, the weight of the photon changes according to the Beer–Lambert law:

(15) w i + 1 - = w i + exp - r i r i + 1 σ abs n abs ( r ) d r ,

where the integral is along the straight line between ri and ri+1, σabs is the absorption cross section and nabs is the number density of the absorbing species.

In our code we implemented the most common case of a stratified atmosphere where nabs depends only on the altitude z. In that case, we can avoid computing the integral in Eq. (15) for every propagation step. This is implemented by computing a cumulative density during the code initialization:

(16) M ( z ) = 0 z n abs ( z ) d z .

This integral is computed within a grid of altitudes and then linearly interpolated for arbitrary z between the grid points. The integral in Eq. (15) can be rewritten as

(17) r i r i + 1 σ abs n abs ( r ) d r = σ abs | M ( z i + 1 ) - M ( z i ) | | μ i z | ,

where μiz is the z component of the unitary vector pointing from ri to ri+1 and we assumed μiz0. In the unlikely case that, within working precision, the photon propagation is horizontal, we have

(18) r i r i + 1 σ abs n abs ( r ) d r = σ abs n ( z i ) | r i + 1 - r i | .

Figure 2Outline of the local estimation method. (a) For a given packet we add the contributions to the radiated energy reaching a detector from each scattering event. (b) These contributions are affected by the optical depth of the scattering location, which involves an integral from that point to the detector. The integral is divided into segments where the collision rate is smooth.


2.4 Local estimation method

Typically, in our code we consider 106 to 107 photon packets, and we simulate the signal received by an observer hundreds of kilometers away from the source. With these numbers, the probability that a significant number of photons packets would emerge from the area covered by one image pixel and reach the detector is negligible. For this reason, we cannot simulate an image by accounting only for the final, outgoing direction of the simulated photons. We overcome this problem by means of a local estimation method as described in studies such as Iwabuchi (2006). This method is sketched in Fig. 2a and may be understood as follows.

After the ith collision, the energy of the photon packet is radiated in all directions according to the incident direction μi−1 and the phase function p(θ) (either pR or pHG). Consider a detector with a surface area A at a large distance di from the scattering event with A4πdi2. The energy that reaches this detector is

(19) δ i = A w i - p ( θ i ) 4 π d i 2 exp ( - τ i ) ,

where the energy is normalized to that of a photon of the given wavelength with unitary weight, θi is the angle between μi−1 and the line of sight Si of the event ri, and τi is the optical depth defined as

(20) τ i = S i ν all ( s ) d s ,

where s is a coordinate indexing the line of sight Si (with length di), and νall includes all processes that may affect a photon (Mie scattering, Rayleigh scattering and background absorption). During a simulation we bin all δi according to their direction of arrival to the detector to form an image or according to their arrival time ti=ti+di/c to simulate a photometer waveform. For this purpose, to set up a run the user must specify the features of one or several observing devices, including their location and the field of view and image size of their cameras as well as the frame rate of their photometers. In all photometer plots in this paper our initial time is the time of arrival of a hypothetical unscattered source photon.

The components of τi stemming from Rayleigh scattering and background absorption are computed using an analytical expression for Rayleigh scattering and with the method described in the previous section for background absorption, which is trivially adapted for the propagation towards an observer.

However, for a general cloud geometry, the integration path in Eq. (20) sometimes traverses several cloud boundaries (see Fig. 2b) where the Mie component of the collision rate (νMie) presents discontinuities. We take this into account by dividing the path Si=(sa,sb) into n subintervals (s1=sa,s2), , (sn,sn+1=sb) and rewriting the integral as

(21) τ i , Mie = k = 1 n s k s k + 1 ν Mie ( s ) d s ,

where we can assume each integrand to be smooth. Each integral, from sk to sk+1 is computed numerically via a Gauss–Legendre quadrature with an order that defaults to 3 but can be configured by the user.

The integrals in Eq. (21) are computationally expensive and take a significant portion of the computations. In some cases, however, they can be skipped: if a scattering event happens deep within the cloud, τi is so large that the contribution of δi is negligible. To account for this and save computational time, the user can select a minimum altitude zmin below which Eq. (21) is not computed and the corresponding δi is set to zero.

Note that in the local estimation method a single photon packet leaves observation traces from more than a single time and place and, thus, creates spurious correlations in the predicted observations. These correlations are visible when the simulated signal is weak, but they do not pose a practical concern.

2.5 Cloud geometry

The CloudScat code admits arbitrarily complex cloud geometries. The user constructs the cloud shape by means of Boolean combinations and affine transformations of elementary solid figures. The code is equipped with a library of figures including a sphere, a cylinder, a cone and a half-space delimited by a plane. This can be extended with user-defined figures by implementing methods that (a) compute the intersections of the figure's boundary with a straight line and (b) determine whether the figure contains a given point.

The supported Boolean operations are:

  1. union of several figures F1FM – a point is contained in the union i=1MFi if it is contained in any Fi;

  2. intersection of several figures F1FM – a point is contained in the intersection i=1MFi if it is contained in all Fi;

  3. difference of two figures F1 and F2 – a point is contained in the difference F1F2 if it is contained in F1 but not in F2.

These operations can be nested to an arbitrary depth, constructing, for example, the difference between a union and an intersection.

The figures (and any Boolean combination thereof) can also be modified by affine transformations 𝒜 consisting of a translation vector R and an invertible matrix M. A point r belongs to the transformed figure 𝒜F if r=M-1r-R belongs to F. This means linearly transforming the shape with M and afterwards translating it along R. Again, affine transformations can be combined with other affine transformations and Boolean operations up to an arbitrary depth.

3 Examples and applications

3.1 Photon diffusion model

In part as a verification of the code and in part as a modeling tool, we first consider the heavily simplified case of a point-like lightning source inside a homogeneous, infinite cloud with a planar top that we set at z=0. In this section, we also neglect Rayleigh scattering and background absorption and set a homogeneous collision rate ν=cNdQextπR2, where Nd is a droplet number density, and Qext and R are the respective extinction coefficient and droplet radius, which are also assumed to be homogeneous.

Koshak et al. (1994) applied classical methods of transport theory to the problem of radiation propagation inside a cloud. By truncating the spherical-harmonic expansion of the one-speed Boltzmann equation, they showed that the photon density ψ approximately follows a diffusion equation with an absorption term:

(22) ψ t - ( D ψ ) = - ψ τ A ,

where the diffusion coefficient D is

(23) D = c 2 3 ν ( 1 - g ω 0 ) ,


(24) τ A = 1 ν ( 1 - ω 0 )

is the photon absorption time. Equation (22) must be supplemented by the condition ψ(z=0)=0 imposed by the perfectly absorbing boundary (photons at locations z>0 never reach the cloud again).

The solution of Eq. (22) with the appropriate boundary condition and initial condition ψ(r,t=0)=δ(x)δ(y)δ(z-L), corresponding to a point source at (0,0,-L), is found by the method of images (Krapivsky et al.2010, p. 32) and reads

(25) ψ ( x , y , z ) = e - t / τ A ( 4 π D t ) 3 / 2 exp - x 2 - y 2 - ( z + L ) 2 4 D t - exp - x 2 - y 2 - ( z - L ) 2 4 D t .

From this expression, one obtains the instantaneous flux of photons emerging from the cloud at coordinates (x,y,0) and time t as

(26) f ( x , y , t ) = D ψ z z = 0 = L e - t / τ A ( 4 π D ) 3 / 2 t 5 / 2 exp - x 2 - y 2 - L 2 4 D t .

Moreover, integrating in the xy plane, this yields a total flux per unit time

(27) F ( t ) = - d x - d y f ( x , y , t ) = e - t / τ A - τ D / t π 1 / 2 τ D t τ D - 3 / 2 ,

where we have introduced the characteristic time for the photons to diffuse up to a distance L,

(28) τ D = L 2 4 D .

On the other hand, integrating Eq. (26) in time gives an expression for the number of photons that exit the cloud from a given location at anytime,

(29) Φ ( x , y ) = 0 d t f ( x , y , t ) = 1 2 π ( L 2 + ρ 2 ) L ( L 2 + ρ 2 ) 1 / 2 + 2 ϵ exp - 2 ϵ ( L 2 + ρ 2 ) L ,

where we have introduced ρ=(x2+y2)1/2 and ϵ=(τD/τA)1/2. In cases where the absorption of photons can be neglected (i.e., the time of absorption τA is much larger than the characteristic diffusion time τD) Eq. (29) reduces to

(30) Φ ( ρ ) = L 2 π ( L 2 + ρ 2 ) 3 / 2 ,

which, remarkably, leads to a spatial distribution of optical radiation that does not depend on the Mie scattering parameters but only on the source depth L.

Soler et al. (2020) fitted expression (27) to the MMIA photometer signal in order to infer the altitudes of presumed fast breakdown events (Rison et al.2016) associated with radio-detected narrow bipolar events (Le Vine1980; Smith et al.1999; Liu et al.2018). They obtained a best-fit τD, and using Eq. (28), given a plausible bracket for D, they derived an interval for L which, combined with independent measurements of the cloud-top height provides a range of source altitudes.

Another approach to estimate the source depth is by means of the spatial distribution described by Eq. (30), which offers the advantage that it does not demand assumptions about D. This is, however, limited to cases where absorption is weak and the irregularities of the cloud top do not affect the recorded image too much. Both of these conditions generally imply that the source is not too deeply buried in the cloud.

Clearly, both of these approaches to estimate the source depth also depend on approximating the source as point-like. This is acceptable for fast breakdown events which span some hundreds of meters (Rison et al.2016; Tilles et al.2019) but becomes more questionable for other types of events.

Let us apply our Monte Carlo model to this configuration. We run the code for a point source within an infinite, homogeneous cloud. The cloud top is set at 15 km and the source at 10 km. The cloud is composed of droplets with radius R=10 µm and a density Nd=100cm-3, and we consider the propagation of radiation with a wavelength λ=337 nm. The resulting Mie scattering parameters are g=0.874, ω0=1-2.82×10-6 and Qext=2+3.99×10-2.

The photons emerging from the cloud have a slight preference for the direction perpendicular to the cloud boundary, so the radiance is not perfectly Lambertian. For the purpose of comparison with the analytical expressions derived above, we must integrate the emissions in all directions. Making use of the azimuthal symmetry of the emissions, we achieve this by setting several observers, all of them at a distance of 400 km from the point in the cloud top closest to the source, but we vary the zenith angle of observation. We chose these angles to perform a five-point Gauss–Legendre quadrature of the photon flux in all directions. Figure 3 compares the Monte Carlo results and expression (27). We consider the agreement of the two curves as a verification of the code implementation as well as of the diffusion model proposed by Koshak et al. (1994).

Figure 3Total photon flux emerging from a semi-infinite cloud. The cloud extends infinitely below 15 km, and the photon source is located at 10 km. We compare a Monte Carlo (black) and an analytically solvable diffusion model (red). The photon flux is normalized with respect to the number of source photons (i.e., it must be interpreted as the flux per source photon). The parameters of the model curve are τD=0.503 ms and τA=18.5 ms. Here and in all subsequent photometer plots, the time origin corresponds to the arrival time of an unscattered photon.


Our next step is to account for the fact that the cloud is not infinitely deep and, therefore, some photons exit through the lower boundary. We now consider a cloud limited by two infinite parallel planes. The method of images is also suitable for this problem, but the solution is composed of an infinite sum of contributions from virtual images arranged in a periodic lattice. Assume that the cloud extends from z=0 to z=-b and the photon source is again at z=-L. The resulting lattice consists of positive images at z{-L+2bn:nZ} and negative images at z{L+2bn:nZ}. The contribution of each image to the integrated flux of photons at z=0 reads

(31) F n ± ( t ) = ± 2 π D e - t / τ A ( 4 π D t ) 3 / 2 ( L + 2 b n ) exp - ( L + 2 b n ) 2 4 D t .

Noting that Fn±=F-n, we write the total flux as

(32) F ( t ) = 4 π D e - t / τ A ( 4 π D t ) 3 / 2 n Z ( L + 2 b n ) exp - ( L + 2 b n ) 2 4 D t .

Notice that if we keep only the n=0 term in this series, we recover expression (27). For long time spans, Eq. (32) requires the unwieldy addition of a large number of terms. Therefore, in order to obtain the long-term behavior of F(t), we use the following identity that results from the Poisson summation formula and is valid for arbitrary u, v:

(33) n Z ( u + n v ) e - ( u + n v ) 2 = 2 π 3 / 2 v 2 k = 1 k sin 2 π k u v exp - π 2 k 2 v 2 .

This allows us to cast Eq. (32) into

(34) F ( t ) = 2 π D e - t / τ A b 2 k = 1 k sin π k L b exp - π 2 k 2 D t b 2 ,

which converges faster than Eq. (32) for long t. Indeed, from the slowest decaying term k=1 we extract a decay time

(35) τ S = b 2 π 2 D .

To check our code against these expressions, we run a simulation with the same microscopical parameters as used above for Fig. 3; however, in this case, the cloud exists only above an altitude of 7 km. The photon source is again at 10 km and, thus, closer to the lower boundary than to the cloud top. In this situation we expect Eq. (27) to perform poorly for long time spans. This is indeed the case as we show in Fig. 4, where we plot the Monte Carlo results along with expression (27), which disregards the lower boundary. Also represented are two truncations of the series (34) that keep one and two terms. We see that the long-term behavior is perfectly captured by the k=1 term in the series. At early times, the approximation is much improved by including a second term of the sum. For small t, Eq. (27) is closer to the simulation curve; this is because that expression is a one-term truncation of the series (32), which converges faster as t→0.

Figure 4Total photon flux emerging from the top of a planar homogeneous 8 km thick cloud. We plot the Monte Carlo simulation results (black) along with the analytical expression for a cloud without the lower bound (red; same as Fig. 3) and the one- and two-term truncations of the analytical approximation (blue; see text). Note that the two-term approximation reaches unphysical negative values for times shortly after the source pulse, although this is barely discernible in the plot.


3.2 Ozone absorption

As we mentioned above, the absorption of radiation by the molecular components of air cannot be neglected for certain wavelengths. A case in point is the absorption of the far-ultraviolet radiation by ozone in the range from 180 to 230 nm, to which one of the MMIA photometers is sensitive. This photometer was conceived for light emanating from TLEs in the upper atmosphere (Neubert et al.2020), but it may also be sensitive to lightning emissions. Here, we investigate this possibility by simulating the photometer response to a far-ultraviolet source inside a thundercloud.

We selected a wavelength of 200 nm, where the absorption cross section of ozone reaches a minimum σabs=3.145×10-19cm2 (Molina and Molina1986). Therefore, this is the most well-transmitted wavelength, although we do not know if strong emission lines are present around this wavelength in a lightning flash.

The ozone density inside a thunderstorm is also uncertain and may differ significantly from a fair-weather density because it is affected by the strong convection and electrical activity of the thunderstorm (Pan et al.2014; Bozem et al.2014). Nevertheless, for illustration purposes, we selected several profiles of ozone in the atmosphere corresponding to those implemented in the MODTRAN 5 code (Berk et al.2005) as previously introduced by Kneizys et al. (1980) for the LOWTRAN code. They correspond to tropical (15 N), midlatitude summer (45 N, July) and midlatitude winter (45 N, January) conditions as well as to the U.S. Standard Atmosphere (United States Committee on Extension to the Standard Atmosphere1976). We disregard subarctic conditions as they are not relevant for lightning. The ozone profiles are displayed in Fig. 5a.

Figure 5The role of ozone absorption in the photometer signal from a lightning flash. Panel (a) shows several ozone profiles employed in Monte Carlo simulations. The shaded area marks the extent of the cloud, and the horizontal dashed line marks the location of the source. In panel (b), we show the simulated photometer response to a 200 nm source observed from above at an altitude of 400 km. Except in the indicated case, all simulations use the water refractive index. The shaded areas indicate the range of exponential decay rates corresponding to the ozone density at the location of the source and at the cloud top (bear in mind that only the slopes of these curves are relevant: they are placed to pass through the peak of the Monte Carlo curves).


We implemented these ozone profiles in simulations with the same cloud configuration as in the previous section (10 µm droplets with a density of 100 cm−3 in a cloud spanning from 7 to 12 km; source at 10 km) but with a 200 nm source. Besides the role of molecular absorption, this wavelength also underlines differences between the absorption properties of water and ice: whereas water absorbs 200 nm radiation strongly, absorption by ice is negligible. However, we observe this difference only when molecular absorption by ozone does not dominate. The tropical case, for example, has the lowest ozone density and, thus, allows us to see different absorption rates between water and ice. Therefore, we run simulations for that ozone profile using both ice and water optical constants. Note also that we are now not interested in the total photon flux emerging from the cloud top; therefore, we simulated the response of a single photometer located at an altitude of 400 km directly above the source.

The results, plotted in Fig. 5b show the signal decaying quickly within a fraction of a millisecond. If the ozone density (nO3) was homogeneous within the cloud, the decay time would be

(36) τ O 3 = c σ abs n O 3 .

However, with the selected profiles, the photons explore regions with markedly different ozone densities. In Fig. 5b, we indicate, by means of a shaded region, the range of decay rates between the source location (slowest decay) and the cloud top (the fastest decay). In most cases, the decay rate is between the two extreme values predicted by the ozone concentrations at the source and at the cloud top. The exception is the tropical ozone profile, where the strong absorption by water droplets leads to a somewhat faster decay. If we use the same ozone profile but employ the ice refractive index, the decay is slower and is within the range predicted by ozone absorption.

It is important to note that although the decay in the observed optical radiation is dominated by ozone absorption in most cases in the geometry considered here, this is not always the case – as we have seen for the tropical ozone profile. Once more, the key is the interplay between the timescales defined in the previous section and the timescales for ozone absorption in the relevant range of altitudes. For example, with a thinner cloud only extending up to 12 km altitude, we clearly see the effect of τS as defined above.

3.3 The timescales involved in photometer waveforms

In the previous sections, we introduced four different timescales, summarized in Table 1, which affect the photometer response to a lightning flash. Although these timescales have a closed-form expression only within the simplified models that we have considered so far, they nevertheless provide (also in more realistic settings) a useful framework of analysis. However, it must be noted that a lightning stroke or flash possesses intrinsic timescales left aside in our discussion: the final photometer response results from the convolution of these intrinsic timescales with those introduced by in-cloud scattering.

Our experience leads us to the following broad guidelines for the analysis of photometer pulses:

  1. Absorption by ozone molecules (timescale τO3) is mostly negligible for optical wavelengths but is strong below about 300 nm. Therefore, one expects very short pulses in the ultraviolet, as shown in Fig. 5.

  2. Under most circumstances, absorption by droplet or ice particles (measured by τA) is relatively weak and is visible only in the tail of the photometer response

  3. Usually both the time for diffusion to the cloud top (τD) and the total residence time (τS) are relevant, with τD mostly affecting the rise-time of the signal and τS affecting its decay. This is visible in Fig. 4, where the rise is captured by an expression that involves only τD, whereas the decay follows an exponential with characteristic time τS.

Table 1Timescales defining the photometer response to an impulsive optical source inside a thundercloud.

Download Print Version | Download XLSX

Figure 6A satellite observes a flash-illuminated complex cloud. The observation geometry is shown in (a), where we represent the cloud geometry and the satellite location. In (b) we show a rendering of the cloud geometry as observed from the satellite's viewpoint. Finally (c) and (d) are the outcome of our Monte Carlo model simulating an image captured by the satellite's camera filtered at 337 and 777 nm, respectively. In this image we have labeled the hole in the cloud turret (1), the emissions from the lower cloud deck (2) and the surrounding, diffuse Rayleigh emissions (3).


3.4 Imaging complex geometries

Let us now turn to the imaging capabilities of the CloudScat.jl code combined with complex cloud shapes. As an example, we construct the following cloud geometry, which roughly mimics a typical convective thundercloud and is shown in Fig. 6a and b. It starts from a cylindrical base with a radius of 20 km and spans altitudes between 7 and 10 km above the ground. A turret emerges from this base that we model as a vertical truncated cone between altitudes of 10 and 15 km and respective radii of 11.25 and 7.5 km. The turret is topped by an ellipsoid centered at (x,y,z)=(0,0,15km) and semi-axes with lengths (12km,10km,2km). Finally, we dig a hole on the top of the turret by subtracting a 3 km radius ball centered at (0km,1km,17km). To aid in the visualization of the cloud geometry, CloudScat.jl generates code that can be used by the Mathematica software (newer than version 11.2) to render a three-dimensional plot of the cloud shape from a satellite's viewpoint. Using this feature, we generated Fig. 6a and b. The cloud has a homogeneous composition identical to that described in the previous section.

The figure constructed is illuminated by a vertical flash in the z axis between an altitude of 8 and 13 km and is observed by a satellite at coordinates (-200km,-200km,400km). The satellite's camera is modeled after the MMIA cameras: pointing towards the nadir with a field of view of 40 (measured as half the diagonal of a square field) and an image size of 1024 pixels×1024 pixels. We consider the two wavelengths observed by the cameras in MMIA: 337 and 777 nm.

To predict how this cloud looks when it is illuminated from the inside by a lightning flash, we run our Monte Carlo model with 108 photon packets for each of the wavelengths. The result is displayed in Fig. 6c and d, where we plot the signal normalized to the number of source photons (i.e., it shows photons in each pixel per steradian and per source photon with the considered wavelength).

Remarkably, for both wavelengths, the image is dominated by light emanating from the hole (1) at the top (note the logarithmic color scale). This stresses the counterintuitive features of observing objects illuminated from the inside, as also underlined by Peterson (2019). Quite often the standout features of a lightning flash are depressions in the cloud geometry.

Figure 7Simulated images of a cloud-scattered source at different wavelengths. Panels (a)(d) show the integrated radiance produced by a point source at 10 or 12 km within a homogenous cylindrical 20 km radius cloud spanning altitudes between 5 and 15 km with a droplet radius of 10 or 20 µm and a droplet concentration of 100 cm−3. Panel (e) shows images of a vertically extended source from 10 to 12 km. All observations are performed in the nadir, and we plot the integrated radiances in a horizontal line crossing the sub-satellite point in each panel as well as the full images in each of insets at the left. To ease the comparison, all data and images are normalized to their peak value.


Figure 8Comparison of a picture from the 337 nm filtered MMIA camera from ASIM (a) and a simulated observation from our Monte Carlo code (b) with the cloud geometry shown in panel (c). The MMIA picture was captured by the 337 nm filtered camera on 22 November 2019 at 08:43:05 at 4.45 N, 77.50 W.


Another feature of both pictures is the illumination of the lower cloud deck (2), partly due to light diffusing from inside the cloud and partly due to the reflection of light emerging from the turret's external boundary. A close look reveals a slight difference between the two wavelengths in the spatial distribution of this brightness, with the 337 nm illumination being somewhat more spread out. This difference will be analyzed in detail in the following section. The darkest areas are those not covered by the cloud, which is also somewhat contrary to our intuition.

One last feature of the simulated image is the weak, diffuse glow away from the brightest areas of the 337 nm picture (3), which is essentially absent in the 777 nm simulation. This glow results from Rayleigh scattering in the atmosphere above and around the cloud. Here, these image areas are 3 to 4 orders of magnitude dimmer than the brightest features.

3.5 Imaging at different wavelengths

Our next application is understanding differences in how a lightning source looks as it is imaged at different wavelengths. As we mentioned in the introduction, the MMIA module of ASIM contains cameras with wavelength filters at around 337 and 777 nm; thus, we restrict ourselves to these wavelengths.

Solving the Mie scattering problem shows that the extinction parameter (always close to 2) and the asymmetry parameter only weakly depend on the radiation wavelength. It is the single-scattering albedo (the absorption by cloud droplets) that dominates the differences between different wavelengths. Table 2 shows the different absorption properties of the cloud microphysical parameters that we have considered in our simulations. Remarkably, the different absorption rates are visible in a cloud-scattered image as a softer, more blurred image in the least absorbed band.

Table 2Absorption of light of different wavelengths by clouds with droplet radii of 10 and 20 µm and a density of 100 cm−3. See the text for the definition of the single-scattering albedo ω0 and expression (24) for the absorption time.

Download Print Version | Download XLSX

To illustrate this feature, we simulate observations of sources within a 20 km radius cylindrical cloud between an altitude of 5 and 15 km (the larger span of the cloud compared with other simulations limits the leakage of photons from the lower edge). The cloud is populated with a droplet density of 100 cm−3, and we test droplet radii of 10 and 20 µm to emphasize the role of absorption, with the latter value leading to much higher absorption rates. Figure 7a–d compare the outcomes of these simulations. We notice that when absorption is weaker, the emissions of a point source are more spread out and look more diffuse. The reason for this is that photons that travel radially spend more time inside the cloud and have a higher probability of being absorbed. This reduces the spread of the most absorbed radiation, which, in this case, is in the 777 nm band.

This effect is somewhat stronger if the light source extends vertically. This is shown in Fig. 7e, which contains results of a simulation similar to those in Fig. 7b and d (i.e., 20 µm droplets) but where the source is a vertical channel with a constant luminosity per unit length extending from an altitude of 10 to 12 km. In this case, we see that the difference between the two wavelengths is even more noticeable than in either Fig. 7b or d. The reason for this is that, due to a weaker absorption of the 337 nm band, we are capable of seeing emissions in this wavelength coming from deeper within the cloud and these create a wider image (Eq. 29 provides a quantitative estimate of this).

3.6 Interpretation of an MMIA observation

As a final application of the code, we analyze an actual observation recorded by the MMIA module of ASIM. We chose an event that took place on 22 November 2019 at 08:43:05 UTC at 4.45 N, 77.50 W (about 10 km off the Pacific coast of Colombia).

Figure 8a contains the image from MMIA's 337 nm filtered camera. There is a ring structure that suggests the effect of a cloud turret extending above a lower cloud surface. The cloud turret contains two differentiated regions, with the lower being somewhat brighter. Finally, the lower part of the ring surrounding the turret has two maxima to the left and right.

These considerations suggest that the discharge may have taken place below a cloud turret with two lobes and a slight depression in the cloud top to explain the two lower maxima. Therefore, we constructed the following cloud shape: the cloud base is a 15 km radius cylinder between an altitude of 5 and 12 km. From this cylinder, we subtract an ellipsoid with semi-axes (4km,6km,2km) and centered at (2km,-4km,12km) with respect to the ground point below the base center. The two lobes are represented by ellipsoids: the first one with semi-axes (4km,6km,2km) and centered at (-1km,-1km,12km) and the second with semi-axes (4km,5km,5km) and centered at (-4km,3km,12km). Moreover, the dark upper lobe in the MMIA image may indicate a strong absorption at higher altitudes, so we imposed an effective droplet radius that grows linearly from 10 µm at the lower boundary of the cloud (5 km) to 20 µm at the upper boundary of the tallest lobe (15 km). The droplet density is set to 100 cm−3 everywhere inside the cloud.

The implemented geometry is plotted in Fig. 8c, and the simulation results in Fig. 8b can be compared with the direct observation. We made no attempt to rigorously fit the geometry parameters to the observations, and, clearly, there are too many possible configurations that would produce similar observations to claim that the hypothesized shape really captures the reality of the cloud. The point of this exercise was rather to show the usefulness of the Monte Carlo model to test hypothesis about the cloud composition and configuration as well as to gain intuition about them.

4 Conclusions

We developed the CloudScat.jl code to assist in the interpretation of space-based observations of lightning-illuminated clouds. We were particularly interested in observations by the MMIA module of the ASIM instrument onboard the International Space Station, but we invite the broader research community to experiment with our code and adapt it to other observation platforms.

At present, the cloud microphysics in the code is somewhat oversimplified and does not compete in sophistication with radiative transfer models in other areas of the climate and weather research communities. Future work could address this shortcoming by combining our radiative transfer code with cloud physics models and with additional inputs from other Earth-observing spacecraft. This path, using a different code, has very recently been undertaken by Brunner and Bitzer (2019).

Currently the atmospheric electricity community regards the scattering of lightning-produced optical radiation by clouds as a hindrance that prevents a more detailed view of the lightning discharge. However, there is another point of view whereby lightning discharges are beacons that probe the cloud microphysics by allowing a deeper view into them. To realize this possibility, one needs, in addition to good space-borne detectors, reliable numerical codes to interpret the observations. It is our intention that our code serves in this role.

Appendix A: Implementation details

The code presented here is implemented in the Julia programming language. Although this language is gaining adoption in numerically intensive research areas, it is still less popular than more established languages such as Python, Fortran and C++. However, we considered it particularly suited for our application because it offers a large amount of flexibility with little or no performance overhead. This is visible in some aspects of our code that we wish to mention.

A significant feature of Julia is that each function is compiled at the time of first call with its code being specialized for the given argument types. That means that the compiler uses information contained in these types to optimize the output code.

We make use of this feature in the specification of the cloud composition, where we define different types for homogeneous and inhomogeneous clouds. For an inhomogeneous cloud, the code has to repeatedly call a function that probes the cloud composition at each collision location. In principle, one could particularize this to a homogeneous cloud by letting this function always return the same value; however, we would then needlessly call this function many times. This is avoided in Julia because even if the cloud composition is specified dynamically (i.e., at run time), the subsequent code is specialized depending on whether the cloud is homogeneous or not.

A similar benefit is obtained regarding the cloud geometry, which is defined as a complex structure of nested affine transforms and Boolean operations that act on different types of elementary shapes. We encode this structure into a data type and make it available to the compiler. In this case, to ensure that the code is fully specialized to the geometry structure, we employed another Julia feature, which is meta-programming, particularly generated functions. In these functions the code can be constructed dynamically according to the type of the function arguments. Thus, we dynamically build code tailored to the specific cloud geometry that is then compiled and fully optimized. This ensures the highest possible performance for each configuration.

These kind of optimizations are also available in other compiled languages like Fortran or C++; however, one loses the dynamical capability in those cases, leading to a slower development cycle and a poorer user experience.

Code availability

The latest version of our code is available at (last access: 11 November 2020). For the simulations presented here, we used release 1.0 (hash f5358c023a08590c72dd0488e2454f338ba0d3c0), which was archived in a Zenodo repository (Luque2020). As the code is distributed as a Julia package, the simplest way to install it and its dependencies is with the Julia package manager by typing “] add” at the Julia prompt. The repository contains an introduction to the code usage as well as heavily commented example inputs, including those used to produce the figures in this paper.

Author contributions

AL designed the study, wrote the numerical code and drafted the paper. FJGV, DL, FJPI and AS tested the code and suggested improvements. FJGV, DL, AMR, FJPI, SS, OC, MH, TN and NØ helped with the analysis of the ASIM observations. OC, TN, VR and NØ provided access to the ASIM repository and clarified details of the MMIA design. All authors contributed to the final version of the paper.

Competing interests

The authors declare that they have no conflict of interest.


Francisco José Gordillo-Vázquez, Dongshuai Li, Alejandro Malagón-Romero, Sergio Soler and Alejandro Luque acknowledge financial support from the State Agency for Research of the Spanish MCI through the “Center of Excellence Severo Ochoa” award for the Instituto de Astrofísica de Andalucía (grant no. SEV-2017-0709). ASIM is a mission of the ESA SciSpace program for the scientific utilization of the ISS and non-ISS space exploration platforms and space environment analogues. It is funded by ESA and national contributions through contracts with TERMA and the Technical University of Denmark (DTU) in Denmark, the University of Bergen (UB) in Norway and the University of Valencia (UV) in Spain.

Financial support

This research has been supported by the H2020 European Research Council (eLightning; grant no. 681257). Francisco José Gordillo-Vázquez and Sergio Soler received support from the Spanish Ministry of Science and Innovation (MCI, under grant nos. ESP2017-86263-C4-4-R and PID2019-109269RB-C-43).

The publication fee for this article was covered by the CSIC Open Access Publication Support initiative through its Unit of Information Resources for Research (URICI).

Review statement

This paper was edited by Havala Pye and reviewed by two anonymous referees.


Adachi, T., Sato, M., Ushio, T., Yamazaki, A., Suzuki, M., Kikuchi, M., Takahashi, Y., Inan, U. S., Linscott, I., Hobara, Y., Frey, H. U., Mende, S. B., Chen, A. B., Hsu, R.-R., and Kusunoki, K.: Identifying the occurrence of lightning and transient luminous events by nadir spectrophotometric observation, J. Atmos. Solar-Terr. Phy., 145, 85,, 2016. a

Bates, D. R.: Rayleigh scattering by air, Planet. Space Sci., 32, 785,, 1984. a

Berk, A., Anderson, G. P., Acharya, P. K., Bernstein, L. S., Muratov, L., Lee, J., Fox, M., Adler-Golden, S. M., Chetwynd, J. H., Hoke, M. L., Lockwood, R. B., Gardner, J. A., Cooley, T. W., Borel, C. C., and Lewis, P. E.: MODTRAN 5: a reformulated atmospheric band model with auxiliary species and practical multiple scattering options: update, vol. 5806 of Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series, p. 662,, 2005. a

Blakeslee, R. J.: Non-Quality Controlled Lightning Imaging Sensor (LIS) on International Space Station (ISS) Science Data, NASA Earth Data,, 2019. a

Blakeslee, R. J., Christian, H. J., J., Mach, D. M., Buechler, D. E., Koshak, W. J., Walker, T. D., Bateman, M. G., Stewart, M. F., O'Brien, S., Wilson, T. O., Pavelitz, S. D., and Coker, C.: Lightning Imaging Sensor (LIS) on the International Space Station (ISS): Launch, Installation, Activation, and First Results, in: AGU Fall Meeting Abstracts, vol. 2016, p. AE23A, 2016. a

Boccippio, D. J., Koshak, W. J., and Blakeslee, R. J.: Performance Assessment of the Optical Transient Detector and Lightning Imaging Sensor. Part I: Predicted Diurnal Variability, J. Atmos. Ocean. Tech., 19, 1318,<1318:PAOTOT>2.0.CO;2, 2002. a

Bodhaine, B. A., Wood, N. B., Dutton, E. G., and Slusser, J. R.: On Rayleigh Optical Depth Calculations, J. Atmos. Ocean. Tech., 16, 1854,<1854:ORODC>2.0.CO;2, 1999. a, b

Bohren, C. and Huffman, D. R.: Absorption and Scattering of Light by Small Particles, John Wiley & Sons, New York, USA, 1983. a

Bozem, H., Fischer, H., Gurk, C., Schiller, C. L., Parchatka, U., Koenigstedt, R., Stickler, A., Martinez, M., Harder, H., Kubistin, D., Williams, J., Eerdekens, G., and Lelieveld, J.: Influence of corona discharge on the ozone budget in the tropical free troposphere: a case study of deep convection during GABRIEL, Atmos. Chem. Phys., 14, 8917–8931,, 2014. a

Brunner, K. N. and Bitzer, P. M.: A first look at cloud inhomogeneity and its effect on lightning optical emission, Geophys. Res. Lett., 47, e2020GL087094,, 2019. a

Cecil, D. J., Buechler, D. E., and Blakeslee, R. J.: Gridded lightning climatology from TRMM-LIS and OTD: Dataset description, Atmos. Res., 135, 404,, 2014. a

Chanrion, O., Neubert, T., Mogensen, A., Yair, Y., Stendel, M., Singh, R., and Siingh, D.: Profuse activity of blue electrical discharges at the tops of thunderstorms, Geophys. Res. Lett., 44, 496,, 2017. a

Chanrion, O., Neubert, T., Lundgaard Rasmussen, I., Stoltze, C., Tcherniak, D., Jessen, N. C., Polny, J., Brauer, P., Balling, J. E., Savstrup Kristensen, S., Forchhammer, S., Hofmeyer, P., Davidsen, P., Mikkelsen, O., Bo Hansen, D., Bhanderi, D. D. V., Petersen, C. G., and Lorenzen, M.: The Modular Multispectral Imaging Array (MMIA) of the ASIM Payload on the International Space Station, Space Sci. Rev., 215, 28,, 2019. a

Chern, J. L., Hsu, R. R., Su, H. T., Mende, S. B., Fukunishi, H., Takahashi, Y., and Lee, L. C.: Global survey of upper atmospheric transient luminous events on the ROCSAT-2 satellite, J. Atmos. Solar-Terr. Phy., 65, 647,, 2003. a

Christian, H. J., Blakeslee, R. J., Boccippio, D. J., Boeck, W. L., Buechler, D. E., Driscoll, K. T., Goodman, S. J., Hall, J. M., Koshak, W. J., Mach, D. M., and Stewart, M. F.: Global frequency and distribution of lightning as observed from space by the Optical Transient Detector, J. Geophys. Res.-Atmos., 108, 4005,, 2003. a, b

Dwyer, J. R., Smith, D. M., and Cummer, S. A.: High-Energy Atmospheric Physics: Terrestrial Gamma-Ray Flashes and Related Phenomena, Space Sci. Rev., 173, 133,, 2012. a

Ebert, U., Nijdam, S., Li, C., Luque, A., Briels, T., and van Veldhuizen, E.: Review of recent results on streamer discharges and discussion of their relevance for sprites and lightning, J. Geophys. Res.-Space Phys., 115, A00E43,, 2010. a

Goodman, S. J., Blakeslee, R. J., Koshak, W. J., Mach, D., Bailey, J., Buechler, D., Carey, L., Schultz, C., Bateman, M., McCaul, E., and Stano, G.: The GOES-R Geostationary Lightning Mapper (GLM), Atmos. Res., 125, 34,, 2013. a

Hale, G. M. and Querry, M. R.: Optical constants of water in the 200-nm to 200-micrometer wavelength region, Appl. Optics, 12, 555,, 1973. a

Haynes, W.: CRC Handbook of Chemistry and Physics, 97th edition, CRC Press, Boca Raton, FL, USA, 2016. a

Iwabuchi, H.: Efficient Monte Carlo Methods for Radiative Transfer Modeling, J. Atmos. Sci., 63, 2324,, 2006. a, b

Kneizys, F. X., Shettle, E. P., Gallery, W. O., Chetwynd, J. H., J., Abreu, L. W., McClatchey, R. A., Fenn, R. W., and Selby, J. E. A.: Atmospheric transmittance/radiance: Computer code LOWTRAN 5, Unknow, 1980. a

Koshak, W. J., Solakiewicz, R. J., Phanord, D. D., and Blakeslee, R. J.: Diffusion model for lightning radiative transfer, J. Geophys. Res., 99, 14361–14371,, 1994. a, b, c

Krapivsky, P., Redner, S., and Ben-Naim, E.: A Kinetic View of Statistical Physics, Cambridge University Press, Cambridge, UK, 2010. a

Le Vine, D. M.: Sources of the strongest RF radiation from lightning, J. Geophys. Res., 85, 4091,, 1980. a

Light, T. E., Suszcynsky, D. M., Kirkland, M. W., and Jacobson, A. R.: Simulations of lightning optical waveforms as seen through clouds by satellites, J. Geophys. Res., 106, 17103–17114,, 2001. a, b

Liu, F., Zhu, B., Lu, G., Qin, Z., Lei, J., Peng, K.-M., Chen, A. B., Huang, A., Cummer, S. A., Chen, M., Ma, M., Lyu, F., and Zhou, H.: Observations of Blue Discharges Associated With Negative Narrow Bipolar Events in Active Deep Convection, Geophys. Res. Lett., 45, 2842,, 2018. a

Luque, A.: aluque/CloudScat.jl: 1.0 release candidate 1 (Version 1.0rc1), Zenodo,, 2020. a

Mach, D. M., Christian, H. J., Blakeslee, R. J., Boccippio, D. J., Goodman, S. J., and Boeck, W. L.: Performance assessment of the Optical Transient Detector and Lightning Imaging Sensor, J. Geophys. Res.-Atmos., 112, D09210,, 2007. a

Minschwaner, K., Anderson, G. P., Hall, L. A., and Yoshino, K.: Polynomial coefficients for calculating O2 Schumann-Runge cross sections at 0.5 cm−1 resolution, J. Geophys. Res., 97, 10103,, 1992. a

Molina, L. T. and Molina, M. J.: Absolute absorption cross sections of ozone in the 185- to 350-nm wavelength range, J. Geophys. Res., 91, 14501,, 1986. a, b

Neubert, T., Østgaard, N., Reglero, V., Chanrion, O., Heumesser, M., Dimitriadou, K., Christiansen, F., Budtz-Jorgensen, C., Kuvvetli, I., Rasmussen, I. L., Mezentsev, A., Marisaldi, M., Ullaland, K., Genov, G., Yang, S., Kochkin, P., Navarro-Gonzalez, J., Connell, P. H., and Eyles, C. J.: A terrestrial gamma-ray flash and ionospheric ultraviolet emissions powered by lightning, Science, 367, 183,, 2020. a, b, c

Pan, L. L., Homeyer, C. R., Honomichl, S., Ridley, B. A., Weisman, M., Barth, M. C., Hair, J. W., Fenn, M. A., Butler, C., Diskin, G. S., Crawford, J. H., Ryerson, T. B., Pollack, I., Peischl, J., and Huntrieser, H.: Thunderstorms enhance tropospheric ozone by wrapping and shedding stratospheric air, Geophys. Res. Lett., 41, 7785,, 2014. a

Pasko, V. P., Yair, Y., and Kuo, C.-L.: Lightning Related Transient Luminous Events at High Altitude in the Earth's Atmosphere: Phenomenology, Mechanisms and Effects, Space Sci. Rev., 168, 475,, 2012. a

Peck, E. R. and Reeder, K.: Dispersion of Air, J. Opt. Soc. Am., 62, 958–962, 1972. a

Peterson, M.: Using Lightning Flashes to Image Thunderclouds, J. Geophys. Res.-Atmos., 124, 10175–10185,, 2019. a, b, c

Peterson, M. and Rudlosky, S.: The Time Evolution of Optical Lightning Flashes, J. Geophys. Res.-Atmos., 124, 333,, 2019. a

Rison, W., Krehbiel, P. R., Stock, M. G., Edens, H. E., Shao, X.-M., Thomas, R. J., Stanley, M. A., and Zhang, Y.: Observations of narrow bipolar events reveal how lightning is initiated in thunderstorms, Nat. Commun., 7, 10721,, 2016. a, b

Sato, M., Takahashi, Y., Kikuchi, M., Suzuki, M., Yamazaki, A., and Ushio, T.: Lightning and Sprite Imager (LSI) Onboard JEM-GLIMS, IEEJ Transactions on Fundamentals and Materials, 131, 994,, 2011. a

Smith, D. A., Shao, X. M., Holden, D. N., Rhodes, C. T., Brook, M., Krehbiel, P. R., Stanley, M., Rison, W., and Thomas, R. J.: A distinct class of isolated intracloud lightning discharges and their associated radio emissions, J. Geophys. Res., 104, 4189,, 1999. a

Soler, S., Pérez-Invernón, F. J., Gordillo-Vázquez, F. J., Luque, A., Li, D., Malagón-Romero, A., Neubert, T., Chanrion, O., Reglero, V., Navarro-Gonzalez, J., Lu, G., Zhang, H., Huang, A., and Østgaard, N.: Blue Optical Observations of Narrow Bipolar Events by ASIM Suggest Corona Streamer Activity in Thunderstorms, J. Geophys. Res.-Atmos., 125, e2020JD032708,, 2020. a

Thomason, L. W. and Krider, E. P.: The Effects of Clouds on the Light Produced by Lightning., J. Atmos. Sci., 39, 2051,<2051:TEOCOT>2.0.CO;2, 1982. a, b, c

Tilles, J. N., Liu, N., Stanley, M. A., Krehbiel, P. R., Rison, W., Stock, M. G., Dwyer, J. R., Brown, R., and Wilson, J.: Fast negative breakdown in thunderstorms, Nat. Commun., 10, 1648,, 2019.  a

United States Committee on Extension to the Standard Atmosphere: U.S. standard atmosphere, 1976, National Oceanic and Amospheric [sic] Administration: for sale by the Supt. of Docs., U.S. Govt. Print. Off., Washington D.C., USA, 1976. a

van de Hulst, H.: Light Scattering by Small Particles, Dover Books on Physics, Dover Publications, New York, USA, 1981. a

Warren, S. G. and Brandt, R. E.: Optical constants of ice from the ultraviolet to the microwave: A revised compilation, J. Geophys. Res.-Atmos., 113, D14220,, 2008. a

Wilkman, O.: MieScatter, available at: (last access: 10 November 2020), 2013. a

Zhang, D., Cummins, K. L., Bitzer, P., and Koshak, W. J.: Evaluation of the Performance Characteristics of the Lightning Imaging Sensor, J. Atmos. Ocean. Tech., 36, 1015,, 2019. a

Short summary
Lightning flashes are often recorded from space-based platforms. Besides being valuable inputs for weather forecasting, these observations also enable research into fundamental questions regarding lightning physics. To exploit them, it is essential to understand how light propagates from a lightning flash to a space-based observation instrument. Here, we present an open-source software tool to model this process that extends on previous work and overcomes some of the existing limitations.