Dart-Leader and K-Leader Velocity From Initiation Site to Termination Time-Resolved with 3D Interferometry

,


Introduction 1.Brief History of VHF Instrumentation for Lightning Studies
Very High Frequency (VHF) radiation has been used to study lightning since the 1960's because of its ability to penetrate clouds, where most lightning activity occurs.This type of instrument was pioneered by Oetzel and Pierce (1969), who describe an instrument with three antennas ∼ 30 m apart arranged in a right triangle.Their design was narrowband but they suggested signal strengths should be sufficient for any band between 30 MHz and 100 MHz.This design bore a striking resemblance to modern interferometers (INTF), but it only measured the azimuthal direction to a flash.The early INTFs of Warwick et al. (1979), Hayenga (1984) and, Rhodes et al. (1994) used analog phase detection (mixers) due to limitations in digital technology.They operated in narrow frequency bands so that the intermediate frequency signal was within band of available signal processing electronics.These INTFs could measure both the azimuth and elevation of sources, but not the range to the source.Narrow-band INTFs can improve their angular resolution by adding antennas on baselines of different lengths to allow the elimination of phase ambiguity (Rhodes et al., 1994;Shao & Krehbiel, 1996).
Beginning in the late 1990s, improvements in available digitizer speeds allowed the development of broadband digital INTFs with short (∼ 1 µs) recording times that could be triggered multiply during a single flash (Shao et al., 1996;Kawasaki et al., 2000).These instruments took the Fourier transform of signals from multiple antennae to digitally recover the relative phase information.Broadband INTFs can use the higher frequencies in their data-stream to achieve the high angular resolution of a longer baseline narrowband INTF, while using the lower frequencies to remove the phase ambiguity of a narrowband INTF.Despite these advantages, the short recording lengths available for the first broadband INTFs made data interpretation challenging; each VHF event lost the "context" of the entire lightning flash.
More recently (Stock et al., 2014) developed a broadband digital INTF with long continuous recording times (∼ 2 s) using cross-correlation to measure time delays between antennas.Long recording times enabled an entire flash to be captured without loss of context caused by gaps in the data set, while cross-correlation removed the "phasewrap" problem by directly measuring the time-differences between antenna waveforms, making the INTF a short baseline time of arrival (TOA) system.This present study builds heavily on the instrumentation used in Stock et al. (2014) and Stock (2014).The station INTF01 in this study used similar active antennas as Stock with some cost-reduction changes, while the second INTF station (INTF02) used a modified antenna design.Both stations used Stock's processing software.
Another common method of studying lightning through VHF emissions is with time of arrival (TOA) instruments, such as those developed by Proctor (1971), Poehler and Lennon (1979), and Rison et al. (1999).TOA instruments measure the arrival time of individual pulses at a number of different stations several tens of kilometers apart, and use the TOA differences to determine the location of the source in 3D space.The Lightning Mapping Array (LMA) (Rison et al., 1999) is a widely used TOA instrument.With a sufficient number of stations and good line of sight the LMA is able to locate sources over the array to within 12 m RM S horizontally and 30 m RM S vertically (Thomas et al., 2004).
Three-dimensional lightning mapping has also been done with interferometry by combining angular measurements from two different stations.Mardiana et al. (2002) created a three-dimensional INTF (3DINTF) which they estimated to locate sources within 600 m, and Liu et al. (2018) estimated their 3DINTF was accurate within 500 m.These 3DINTFs used segmented recording rather than the continuous recording used in the present study.Hare et al. (2018) have also used a VHF radio telescope, the Low Frequency Array (LOFAR), to map lightning with accuracy on the order 1 m horizontally and 10 m vertically (Hare et al., 2020).3D lightning mapping instruments exist at lower frequencies as well, notably the Huntsville Alabama Marx Meter Array (HAMMA) (Bitzer et al., 2013) operating at 1 Hz to 400 kHz, the Fast Antenna Lightning Mapping Array (FALMA) (Wu et al., 2018) operating at 500 Hz to 500 KHz, and the Position by Fast Antenna (PBFA) instruments of Stolzenburg, Marshall and Karunarathna et al. (2017).

New Contributions
In this paper we will go into some detail on the analysis procedure of one of the first continuous 3D interferometers (3DINTF) and verify its accuracy against a collocated LMA.We will also derive time and space-resolved velocity profiles of repeated K leaders and dart leaders and speculate on what they teach us about leader physics.It thus behooves us to also review what is known about dart and K Leaders.

Terminology
Dart Leaders and K leaders are fast lightning leaders that retrace channels previously created by slower leaders in virgin air.Kitagawa (1957) coined the term K change (K for"Kleine", or small) to describe an electric field change signature which he suggested was caused by the same process as a dart leader.(K changes appear as a smaller version of the return stroke field change that occurs when lightning strikes the ground.)Leaders associated with K changes are often called K leaders or K-processes waves (Winn et al., 2011).Despite the recognized equivalence of the physics behind K leaders and dart leaders (Kitagawa, 1957;Shao et al., 1995), a distinction continues to be made in the lightning community with the general understanding that a dart leader progresses to ground while a K leader remains in the clouds.This distinction is somewhat blurred by authors who refer to failed or attempted dart leaders which do not reach the ground (Shao et al., 1995;Rhodes et al., 1994).It seems clear that there should be a common name which encompasses all such events if the physics behind them is believed to be the same.In search of a common name some authors refer to all such activity as recoil leaders/streamers (Akita et al., 2010;Mazur, 2016), or retrograde leaders (Edens et al., 2012).We suggest to make the term dart leader encompass all leaders of this type, since it was the first term used to describe this phenomenon, it is descriptive of their high velocity, it is agnostic of the detailed physical mechanism which is not yet well established, and it can be in-clusive of both retrograde leaders on positive channels as described by Edens et al. (2012), or prograde leaders on negative channels as observed by Shao et al. (1995).(K leaders could perhaps be called IC dart leaders if it is necessary to specify that they are not followed by a return stroke.)However, until the community reaches a new consensus, we will use existing terminology.In this paper, we analyze two K leaders followed by three dart leaders; all using parts of the same channel.

Properties of Dart Leaders and K Leaders
Dart leaders were identified as early as the 1930's by Schonland et al. (1935) using a Boys camera, where two dimensional average velocities were found to range from 1×10 6 m/s to 23×10 6 m/s , an order of magnitude or two higher than stepped leaders.Further studies of dart leaders are rather consistent with these velocities.Table 1 shows dart leader velocity reported in selected papers, along with some K leader velocities, which fall in the same range.Schonland et al. (1935) also reported that slower dart leaders corresponded with longer intervals between return strokes.This observation was corroborated also by Loeb (1966) and Shao et al. (1995), and in laboratory analogues (Winn, 1965).

Paper
Leader Velocities (m/s) leader type (Schonland et al., 1935) 1 × 10 6 to 23 × 10 6 Dart (Loeb, 1966) 2 × 10 6 to 20 × 10 6 Dart (Jordan et al., 1992) 6 × 10 6 to 50 × 10 6 Dart (Shao et al., 1995) 1 × 10 6 to 10 × 10 6 Dart (Stock et al., 2014) 3 × 10 6 to 17 × 10 6 K leader This study 2 × 10 6 to 20 × 10 6 Dart & K Several studies have reported evidence of dart leaders slowing down as they propagate.In the laboratory, analogues of dart leaders exhibited an effect of the ionization waves slowing down as they propagated along the channel (Winn, 1965).In the field, Jordan et al. (1992) used a streak camera to study triggered lightning in New Mexico and Florida and calculated the average velocity of two short segments of the observed channel near the ground.Jordan found that in New Mexico dart leaders tended to slow down as they approached ground but in Florida they tended to speed up.Schonland et al. (1935)  directly into the data stream on a high order bit.

3D Interferometry
Raw data was first processed separately using three antennas at each station to calculate azimuth and elevation angles according to the methods outlined by Stock et al. (2014) and Stock (2014).Cross-correlation is used to measure the time of arrival difference between each pair of antennas.The time of arrival difference between any two antennas determines the source angle as where α is the incident angle relative to the baseline between the two antennas, c is the speed of light, τ d is the difference in time of arrival, and d is the distance between the two antennas.(Equation (1) is only strictly true for a plane wave, but it is a good approximation for spherical wavefronts when the distance to the source is much greater than the baseline length.)For a set of three antennas arranged in a triangle independent angles α and β can be calculated.These uniquely give the direction to a source, but not its range.From α and β azimuth and elevation may be readily calculated.
Azimuth and elevation angles from two different INTF stations can be combined using the triangulation method as detailed by Thyer (1962) and Liu et al. (2018).A similar method was used by Mardiana et al. (2002).The triangulation algorithm, along with additional details for different station configurations, are included in Appendix A. The triangulation method gives a 3D location for a source given an azimuth and elevation angle measured from two different locations.In any two station method, there is an additional challenge in determining which sources on each station have a correspondence.
In the analysis presented here, we found a correspondence for about half of the detected sources at station 1. Station 2 detected about twice as many sources as station 1.We believe that the antenna and pre-amp design differences discussed in section 2.1 should allow INTF02 to have a better signal to noise ratio than INTF01, but there may be other factors contributing to this difference.This is an ongoing area of investigation.
For a source detected by both stations the maximum possible time of arrival difference is determined by where D is the distance between stations and c is the speed of light.For the two INTF stations used, separated by 16 km, that time is 53 µs.We achieved source matching between the stations by calculating the triangulated locations for every possible pair of sources between the two stations that are separated by 53 µs or less.The source time corresponding to each trial location was calculated using equation 3.
In equation 3, t 0 , x 0 , y 0 , and z 0 are the source time and position for the trial location, and t 1 , x 1 , y 1 , and z 1 are the time of arrival and position for one of the antennas (the choice is arbitrary).Antenna positions were determined by GPS and surveying, while time-of arrival is determined by the GPS-disciplined time stamp of each acquired data sample.
We then determined the best match by calculating a goodness-of-fit parameter, where t obs i is the observed time of arrival at the i-th antenna, t f it i is the time of arrival corresponding to the trial location, and N is the number of antennas.The trial location that gave the lowest χ 2 value was chosen as the best match location.Best matches were first identified for every source from INTF01, so that each source from INTF01 was matched to only one source from INTF02.Best matches were then determined for every remaining INTF02 source so that no source from either station would be included in more than one match in the final set.
We then further refined the 3D locations by using a minimization algorithm to calculate the source location corresponding to the minimum χ 2 value as given by Equation ( 4), similar to the method described in Thomas et al. (2004Thomas et al. ( , 2000) ) for the LMA.The results in this paper used the Gauss-Newton algorithm for minimization, but later processing with the Levenberg-Marquardt algorithm (used by Thomas et al.) yielded similar results.In the comparison with the LMA (discussed in Section 2.3.3 ) minimizing χ 2 only improved the accuracy by about 10% versus triangulation alone, but we are retaining this step because modeling suggests that it will lead to further accuracy increases with improvements in calibration.(These improvements include obtaining centimeter accurate relative antenna locations between stations, and using identical GPS time synchronization at both stations.)

Time correction on INTF01
INTF01 data was synchronized with network time, which provided an accuracy of roughly σ = ±4 ms, however 3DINTF requires timing accuracy of order 100 ns or better.A high-precision GPS attached to one of the digitizer channels of the INTF could provide such accuracy, but none was present for this study.Fortunately, through a synchronization of INTF pulse and LMA pulse arrival times, INTF time can be corrected using the LMA as a reference.This same procedure also results in a set of points which can be used to check the accuracy of the 3DINTF.Details of this synchronization scheme will be discussed in section 2.3.2.

INTF to LMA correlation
To verify the validity of the 3DINTF locations and estimate location accuracy we compare 3DINTF and LMA source locations for a set of four flashes.Matching of LMA and INTF data was carried out as follows: The LMA data used had already been processed by well-established code which turns The LMA locations on this filtered list can now be directly compared against the 3DINTF locations.These results are illustrated by Figures 1 and 2.

Determining Accuracy of 3DINTF
In the procedure described above the LMA was used only to correct the station time on INTF01.Fundamentally, the LMA measurement and the 3DINTF measurements are independent.Thus the matched INTF/LMA pulse pairs derived as per section 2.3.2 can be used to verify the accuracy of the 3DINTF method.To better quantify these errors the absolute median deviation was calculated for each distribution in Figure 2. The result is 90 m for Easting, 80 m for Northing, and 130 m in altitude.The absolute median discrepancy is reported rather than a root mean square (RMS) value because the median is less affected by the extremely long tails.These tails may be caused by spurious and improper matches between the data sets.The LMA only had 8 stations operating, with many sources only located by 6 or 7 stations, and the storm was on the outside edge of the array, so we estimate that the LMA also had errors on the order of 100-200 m for these flashes.An LMA with 13 or more stations can have RMS errors as low as 20-30 m for sources over the array (Thomas et al., 2004).The small number of LMA stations in operation and uncertainty in matching between the 3DINTF and LMA means the measured discrepancies only serve as an upper bound on the uncertainty in the 3DINTF source locations, the true median accuracy may be better than 200 m.
We plan to conduct a more precise measurement of 3DINTF accuracy in the future.

Velocity Estimation
The set of sources in the dart/K leaders displayed a change in position that was generally monotonic in time and there was little VHF activity on other channels during their occurrence.This allowed a simple rolling average (boxcar window) to be used to filter that leader's position vs. time.The values for each coordinate were calculated as: where N is the number of points averaged over, x(k) is the k-th data point in the set of x coordinates for sources in the leader, and x 1 (i) is defined as the x-coordinate for the i-th point in the leader, x being a traditional way to denote the average.The y, z, and time coordinates of the leader were smoothed in the same way.In order to estimate the velocity these coordinates were compared to the next N points, with coordinates defined as x 2 (i) = x 1 (i + N ) and the velocity was calculated as Several different values N were tested.N = 40 was found to be a good balance between channel tracking and noise rejection.To calculate the velocity error we first calculated the standard deviation σ p about the mean position for each rolling average window.In order to arrive at a value representative of the position error of the entire data set the mean of the individual σ p (σ p ), was calculated.
The mean time difference between consecutive N point windows, ∆t was also found.
As expected, it turned out to be approximately N • 1.4 µs.Assuming that the deviations for each window are independent the average 1σ deviation from the true mean velocity σV is then approximately given by Equation ( 7). 3 Results
The plot shows the tripolar structure typical of many thunderstorms, as described by Williams (1989) and Marshall et al. (2005).The lower positive region appears to be around 2.5-3.5 km above mean sea level (MSL), the main negative region appears to extend from roughly 4 km to 6 km, and the upper positive region is spread from 6 km to 9 km MSL.
The altitudes of these regions are significantly lower than those observed by Marshall et al. (2005) and Edens et al. (2012), but this is to be expected as their observations were made in July and August, and this paper discusses results from a storm in late October.The charge regions are known to be defined by temperature rather than altitude (Krehbiel, 1986), and with lower temperatures at the ground in late October we would expect the charge regions to be lower in altitude.
There is no reason 3DINTF data couldn't be used to perform a similar charge analysis in the future, but the data is not currently compatible with the existing tools for the LMA. Figure 6 shows a histogram of the altitudes of 3DINTF sources from 10:00:00 UTC

Flash Overview
In all the plots that follow, the origin of the coordinate system (X=0, Y=0) is located at station INTF01.The INTF02 station was located 13 km north and 9 km west from the INTF01.The measurements were performed in mountainous terrain.The altitude of the two interferometers is 3.16 km and 2.05 km respectively.
Figure 7 gives an overview of the entire lightning flash of interest.The data is presented in the style of an LMA Plot, which has become a familiar way of displaying threedimensional data (Thomas et al., 2004).As described in section 2.3.2, an LMA was used to check the validity of the 3DINTF mapping, but all the 3D locations shown here are solely INTF results.Other than the charge structure plot in Figure 15, no further data obtained from the LMA is presented in this paper.
We first present an overview to orient the reader, and then return to discuss the flash in more detail stroke by stroke.(For a more complete picture of the flash development please refer to the animation included in the supplementary material (Jensen et al., 2021).)The flash begins with two cloud-to-ground (CG) strokes, visible as nearly vertical features in 7a (before -150 ms, colored dark blue).The National Lightning Detection Network (NLDN) (Cummins et al., 1998) also identifies negative CGs at this time (indicated by diamonds in the plots).The second of these CG's is most likely preceded by a dart leader since it appears to go to ground near the same place as the first, but the 3DINTF mapping of the channel in this region is inadequate for useful analysis.Thus we will not discuss this dart leader further in our analysis.(Also, we deemed it less confusing to not count this missed dart leader in our numbering scheme.Thus when we say "first dart leader" we really mean "first well-resolved dart leader" in what follows.This "missed leader" does not, in our view, impact the conclusions of this paper.) The negative charge brought to ground in the initial -CGs results in a pocket of

Blooming and Channel Tracking
As this paper gets deeper into detailed discussions of flash development, it is helpful to define a new term, "blooming".It has been our observation that RF source data (both LMA and INTF) sometimes shows a clear channel, and in the case of an INTF sometimes shows an exceedingly clear subsequent channel (which we discuss further below).Sometimes however, the time development of RF appears more as a cloud of points.
While sometimes these are noise, at other times we believe they reflect the development  of a highly branched channel or a channel with many leader tips that are active at almost the same time.We refer to this phenomenon in which the channel grows in many directions at once where they are not clearly resolved by the instrument as "blooming".
It is purely descriptive.We are not attributing new physics to it.Researchers who have looked at fast videos, may have noted that for some flashes only one or two tips on a channel are active, but some branches may have a dozen or more optically active sub-branches or leader tips at once.These would be barely resolved by an INTF and may well be what an INTF sees as "blooming".

K Leaders
Figure 8 shows the development of the flash up to the first K leader, with sources colored by time.After the two initial negative CGs near the origin the positive leader propagates primarily to the northwest (in panel c).This positive leader appears to remain in a single well defined channel until about -40 ms, when the first K leader begins.
(The K-leader is shown in brick red in Figure 8b,c,d.) After the first K leader the positive leader resumes propagation, and begins blooming into many branches as seen in Figure 9.There are clear large branches to the northwest, east, and south-west, as best seen in Figure 9c.The branch to the north-west develops into the second K leader, which reaches roughly the same point along the channel as the first K leader.Again, the K-leader shows up as a well-defined channel of brick red points within the more scattered branches of the positive leader.

Dart Leaders
Figure 10 shows continued blooming in the north-west and south-west branches of the positive leader, leading to the first dart leader just after 300 ms.The NLDN reports a negative CG at this time, which is consistent with this being a negative leader propagating back down a channel initially created by positive breakdown.
Figure 11 shows further blooming of the positive leader, primarily in the south-west branch, which leads to a second dart leader around 470 ms, again the NLDN reports a negative CG at this time.
Figure 12 shows further positive leader growth in the south-west branch, and the third dart leader around 510 ms, following quite quickly after the second dart leader at 470 ms.This dart leader is also coincident with an NLDN negative CG.

Reduced "branching" of dart leaders and K leaders
Having discussed the flash in detail, a feature of the data in Figure 7 and the subsequent dart and K leader figures should be remarked upon.Back at Figure 7, panels b and c clearly show a brick-red dart leader (the final one) overlaying the earlier positive leader points.In fact, that final dart leader overlays the earlier ones so completely that they cannot be seen.This fact will be useful in our forthcoming velocity calculations, but is itself of note.The great deal of "scatter" or "blooming" visible on all the VHF sources preceding the dart leaders is not a result of poor location precision.Rather, it seems to be characteristic of the much more highly branched nature of leaders into virgin air on the 100-1000 m scale, as compared to dart leaders.Hare et al. (2019) have observed that all of the positive leader sources they see with LOFAR are actually "needles", which are negative breakdown propagating away from the positive channel some distance behind the tip.From this it seems likely that the majority of positive leader sources observed by any similar VHF systems, including the 3DINTF and the LMA, are in fact from needles.If the needles do not emit VHF as a dart leader    We want to clarify that it is also characteristic of the in-cloud portion of dart leaders, as well as K leaders.We suggest that this reduced branching tells us something about the physics of a dart leader and that it is somehow a preferred path for the re-ionization wave.Similar observations were reported by Shao et al. (1995) and seen in the results of Stock et al. (2014).

Velocity
Velocities for both K leaders and all three dart leaders are shown in Figure 13, while Figure 14 shows how the channel segments were aligned.The velocity of each dart or K leader was integrated over time to give a distance along the channel at each point, and the initial offsets of these integrated distances were adjusted so that the shared portions of each dart leader channel would align in the original spatial coordinates (altitude, northing, and easting), as shown in Figure 14.The zero point was arbitrarily chosen to be the beginning of the final dart leader.Figure 14 demonstrates that, for most of their length, all the K leaders and dart leaders share the same three-dimensional path.The leaders were aligned in this way in order to show how the velocity at each point along the channel varied between the K leaders and dart leaders.Some obvious trends appear, most notably the large dip in velocity around 11 km, which occurs in every dart and K leader that passed that point.The second K leader, which is already traveling slowly as it passes the location of the dip, is the only exception to this behavior.A smaller dip is also apparent in all three dart leaders at around 6.5 km (Figure 13).While we do not understand what is special about the locations of the speed dips, it seems that there is some reproducible feature (presumably related to overall charge structure of the storm) which causes the dips to occur repeatedly at the same location in the thundercloud.
The calculated velocities range from 2 × 10 6 m/s to 20 × 10 6 m/s.This agrees very well with the range of velocities other researchers have reported, as listed in Table 1.
With the exception of the first K leader, which started very quickly and then dropped off again, the other 4 leaders generally increased in velocity with each subsequent pass along the channel.This is consistent with the channel being increasingly conditioned by previous leaders and return strokes (Behnke et al., 2005;Rakov & Uman, 1990), although the first K leader shows that there must be multiple factors that determine the velocity of dart leaders and K leaders.Bazelyan and Raizer (1997) hints at a possible mechanism for conditioning in equation 2.17 and his statement that negative Oxygen ions require only 0.5 eV for reionization.Our data provide more clear examples of this poorly understood phenomenon.

Why does leader velocity decrease along the channel?
Figure 13 clearly shows velocity decreasing as the leaders progress along the channel, similar to the behavior Behnke et al. (2005) saw in initial leaders, but now in the context of dart and K leaders.We can think of three mechanisms which could lead to the observed velocity decreases.The first mechanism we will discuss is increasing pres-  sure as the leader propagates toward ground.Based on mean-free-path considerations, the leader might be expected to slow because of increased pressure (da Silva & Pasko, 2013).
However there are two reasons to reject the pressure/speed hypothesis.First, it was previously shown that for negative stepped leaders descending from 10 km to ground the step length decreased with altitude, but the average velocity did not decrease.The step rate increased inversely with the step length (Edens et al., 2014).The second reason to reject the pressure/speed hypothesis is that roughly the first 10 km of channel progresses at a relatively constant 4.5 km of altitude (Figure 14a), while the next 5 km increases in altitude to 5 km.Thus, if there is any pressure effect on speed, we might expect a speedup, contrary to observation.This suggests that pressure is not the driver of leader speed change for most of the leader.Figure 15 is an LMA charge analysis for the 20 minute period surrounding the flash of interest (see section 3.1).The Altitude vs. East-West panel of this figure shows that the eastern region of the storm piles VHF sources to a substantially higher altitude than the western region.Also, around zero kilometers Easting the trilayer charge structure becomes very apparent.The non-inductive charging theory for thunderstorms connects strong updraft regions with the charging engine.Absent other data, we would think that the updraft engine is operating fiercely right around zero kilometer E-W and that the ambient fields might be higher in this region.We have overlaid the path of the dart leaders on top of the LMA data points in black.Figure 15 demonstrates that the dart leaders and K leaders progress from the western extremities of the storm into the eastern region of putatively higher fields.We would naively expect leaders to speed up as they entered higher ambient fields.To the extent that our extrapolations are correct, this second hypothesis for speed change of the leader fails as well.
What remains with highest probability is the third explanation; in which we consider a dart leader as a guided nonlinear ionization wave in a decaying plasma channel (See Bazelyan and Raizer (2000), Section 4.8).In this framework, the rate of electronimpact and thermal ionization dictates how quickly a new leader section can be formed (da Silva & Pasko, 2013), and thus the dart leader speed is proportional to the magnitude of the electric field created at the leading crest of this soliton.(Like a soliton, these leader tips are isolated pulses in a highly non-linear system with limited dispersion.)This high electric field is needed to reionize the decaying channel.As energy is spent reionizing the channel, the magnitude of the leading edge electric field decays, and so does its velocity.A similar process happens in streamer discharges (Naidis, 2009).
Please note that this is different than the transmission line interpretation of a lightning channel.In a transmission line, the amplitude of the wave decays as a function of distance due to the existence of a finite resistance.The wave velocity is a function of the inductance, capacitance, resistance, and frequency content of the wave packet (Rakov, 1998).However for a constant resistance (the assumption of transmission line models), there should be a constant velocity, i.e. it does not decrease.
One might wonder about the source and sink of the reionization energy?This is not well understood so we will make a hypothesis.Once initiated with some initial energy (by a process which is still not understood), a dart leader continues to be fed by the ambient electrostatic environment.This results in a concentration of charge at the leader tip which also represents electrostatic energy.Some of this energy is spent in ionization at the leader tip and is seen as charges left behind on the re-ionized channel.If the energy fed by the environment is insufficient for the expenditure on ionization, the leader field would decrease and the leader could slow, as observed.

Systematic Error
For the 3DINTF there is an extreme amount of scattering in the source locations associated with the initial CGs (shown in dark blue in Figures 7 and 8).This is most likely caused by poor matching in this region as individual branches of the downward leader were well resolved by the INTF01 but not by the INTF02, because the flash started essentially directly overhead of the INTF01.The INTF01 antennas also have a null in their sensitivity at 90 • overhead so there were gaps in the detected sources, and more work needs to be done to precisely measure the relative orientation of the two arrays.
It is not clear if the large dip in velocity around 11 km in Figure 13 is a real signal or simply a systematic error.There is no obvious increase in scattering of sources for any of the dart leaders or K leaders that pass through this region to indicate that the dip is caused by increasing location error, but it is suspicious that the dip occurs almost exactly along the baseline between the two INTF stations.Figure 16 highlights the location of the dip along the third CG dart leader by coloring points by velocity rather than the standard time coloring.The location of the velocity dip in Figure 16a can easily be identified as the channel color changes from yellow, to green, to blue, and back to yellow in a short period, which coincides with the point where the line between INTF01 and INTF02 crosses the channel.Location errors for the Triangulation method can be increased along this line as small changes in measured azimuth lead to large differences in the calculated range to the point, but this would be expected to manifest as a significant broadening of the channel along this line, which does not appear to be the case.
The third dart leader is shown in the plot.As we previously pointed out, all 3 dart leaders and the K leader exhibiting this dip in speed have it at the same channel location.
All three dart leaders have a smaller dip in velocity at 6.5 km along the channel (Figure 13), at a point with no obvious geometrical significance.Along with the lack of scattering, this suggests that these variations may be real.Further evidence is seen in the fact that the second K leader does not exhibit a dip in velocity when it passes through the same point in the channel.If the drop in velocity was a geometrical artifact then it would be expected to show up for any leader passing through that point.Thus we conclude that the observed variations likely reflect true changes in velocity at these points, although further observations and modeling of processing errors will help to verify this claim.
There is also a significant decrease in velocity in all three dart leaders around 17 km, and a subsequent increase in velocity between 18 km and 20 km.Since these coincide with the gap in VHF data directly overhead of the INTF01 and the region of high scatter near the station these trends are almost certainly just artifacts.Flashes which go to ground further from the INTF01 station will need to be observed in order to study the characteristics of dart leaders as they approach the ground.

Conclusions
In this paper we have: 1. Documented in an appendix a double theodolite location method for use in lightning interferometry The original work (Thyer, 1962) contained a typographical error which left the algorithm ambiguous.We have also extended it to allow for interferometers of nonzero size (theodolites can be considered to be points!), and for arbitrary configurations of the two stations.
2. Verified the minimum accuracy of a 2 station 3DINTF Using a two station INTF colocated with a 8-station LMA, we showed that 3DINTF can have good location accuracy (200 m median error), and sufficient time-resolution to observe rapid processes like dart leaders and K leaders in detail.In checking the location accuracy of the INTF for sources located by both instruments, the LMA was assumed to be "correct".In fact, for the 8-station LMA used, the LMA errors might also have been in the 100 m range.

Shown that charge identification can be done with an INTF
The charge layers we identified were consistent with those found be an LMA, though less complete at this time.We speculate here that, because of the high spatial and temporal sample density of an INTF, future studies might allow one to observe the "granularity" of charges in a cloud (in other words, recognizing smaller charge pockets in addition to the overall layers in a storm. time-stamped VHF pulses into located data points that are time-stamped with the time of emission at the source location.Since the goal of LMA/INTF correlation is to locate INTF points at which initially only a time of a arrival at an INTF antenna is known, the LMA source times are updated to the times at which each source would arrive at the selected INTF antenna.Raw VHF data from the selected single INTF antenna was processed with a 60-66 MHz forward-backward (zero-phase) Butterworth filter of net order N=2.The filtered INTF signal is now analogous to the raw LMA receiver data.A Hilbert transform was next applied to produce a power envelope.The largest peak power in a fixed 40 µs window is then recorded along with its time stamp.Since the LMA only records peak powers in fixed windows, the INTF and LMA peaks are assumed to correspond to the same source.This is true often enough for the technique to work.The pairwise time difference between all INTF VHF peaks and all LMA arrival times is plotted in a histogram with one microsecond bins.An initial time offset estimate is determined by taking a weighted average of the histogram peak and neighboring bins.Once the first time offset is calculated, a linear regression between LMA and INTF source times is calculated to produce a refined time offset.Having corrected INTF time with the previous processing step, the strongest INTF pulse in a 40 µs window is determined and a time match to an LMA source within one microsecond is sought.Matches are added to the pulse-pair list.Once a pulse-pair list is available, it is filtered to keep only the pairs which were also 3D located by the INTF.
Figure 1 illustrates the (mostly) good agreement between the LMA locations (red) and the independently calculated INTF locations (blue) for the flash to be analyzed in this paper.There is an obvious disparity visible in panel d between 3 and 6 km northing.It shows LMA (red) sources going to ground that are not visible in INTF data.These are likely altitude errors resulting from the limited number of LMA stations.The INTF, in this case, is correct in NOT detecting them.There is also a burst of INTF (blue) points southeast of the origin not shown by the LMA.In this case, we hold the 3DINTF matching algorithm responsible.Picking the wrong match point results in a time error which becomes a range error.Cropped histograms of the discrepancy in each coordinate direction and the overall discrepancy between the 3DINTF and the LMA are shown in Figure2.(The uncropped histograms are included in the supplementary information in order to verify that we are only excluding long thin tails.(Jensen et al., 2021) ) Overall the distances between matching 3DINTF and LMA points were roughly log-normally distributed, with a median error of about 220 m. (A log-normal distribution has a normally distributed logarithm.)

Figure 1 .
Figure 1.Comparison of 3DINTF and LMA for the 10:00:09 flash on October 23rd, 2018.Plots show altitude vs time (a), altitude vs east/west (b), north/south vs east/west (c), and north/south vs altitude (d).3DINTF sources are marked in blue and LMA sources are in red.See section 2.3.3 for additional discussion of discrepancies between these data sets.

Figure 2 .
Figure 2. Histograms of the measured discrepancies between matching 3DINTF and LMA source locations, east-west discrepancies (a), north-south discrepancies (b), altitude discrepancies (c), and overall 3D discrepancies (d).Histograms are cropped to 800 m because this range includes 80% of the data, the remaining 20% are distributed in thin tails out as far as 16 km.

Figure 3 .
Figure 3.Comparison of rolling averages of the second dart leader using different numbers of points, in a plot of altitude vs time, with the original VHF sources (black) and the rolling average points.Averages are performed using 10 points (blue), 40 points (green), and 100 points (red).The gap in data is caused by a null in antenna sensitivity directly over the INTF01.

Figure 3
Figure 3 of altitude vs. time for the third dart leader displayed with various boxcar windows, demonstrates that a smaller window results in truer tracking of the channel.In contrast, Figure 4 shows that too small a window yields an excessively noisy velocity plot.The fixed N averaging used is approximately equivalent to averaging over a fixed ∆t for this data set.The raw INTF processing used windows 1.4 µs long (256 samples at 180 MS/s), and locates, at most, one source per window.During dart leaders and k leaders we typically observed sources in every 1.4 µs window, meaning the length of the rolling averaging windows was approximately N • 1.4 µs.

Figure 5 .
Figure 5. Charge structure as identified by the LMA for flashes happening between 09:50:00 UTC and 10:20:00 UTC near Langmuir Lab on 2018-10-23.Negative charge is shown as blue, positive as orange/red.

Figure 6 .
Figure 6.A histogram of 3DINTF VHF source altitudes for flashes happening between 10:00:00 UTC and 10:20:00 UTC near Langmuir Lab on 2018-10-23 reduced negative charge in the region above the CG grounding locations.It is thus energetically favorable for a positive leader to issue horizontally from this reduced negative charge into the main negative charge, and this is precisely what is shown in the blue/yellow/red data points that move west from the origin in Figure7b.(Many of these points are covered by the brick red points of dart leaders that occur later in the flash).Once this positive leader has reached a point roughly seven km west of the origin, a negative K leader travels back along the horizontal channel.(This K leader begins around t=-40 ms, and should be coloured blue, but is likewise buried beneath the brick red data points from later dart leaders.)The positive leader resumes its extension with blue to green data points (-40 ms through 230 ms), and there is another negative K leader at t=230 ms.Positive leader growth then continues with the yellow and orange points until the first resolved dart leader occurs around t=310 ms.After the first dart leader, the positive leader continues to grow, leading to a second dart leader at 470 ms, and a third one at 510 ms.Both the NLDN and 3DINTF data suggest that all of the CGs in this flash go to ground at approximately the same location, on the mountain top just east of the INTF01, which is at an altitude of 3.16 km.Having understood the big picture, let us look at each of these sections in more detail.

Figure 7 .
Figure 7. Overview of the entire flash from 10:00:09 UTC on 2018-10-23, with sources colored blue to red according to time.(a) altitude vs time, (b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations, and the time and location of NLDN strokes.(In panel c, the NLDN ground-strike points are immediately to the east of INTF01).Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes (on which the NLDN and 3DINTF agree).

Figure 8 .
Figure 8.The beginning of the flash up to the first K leader.Sources colored blue to red according to time.(a)altitude vs time,(b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations, and the time and location of NLDN strokes (yellow diamonds).Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes.

Figure 9 .
Figure 9. Overview of the second K leader and leader growth that precedes it.Sources colored blue to red according to time.(a) altitude vs time, (b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations.Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes.

Figure 10 .
Figure 10.Overview of the first dart leader and the blooming that precedes it, with sources colored blue to red according to time.(a) altitude vs time, (b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations, and the time and location of NLDN strokes.Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes.

Figure 11 .
Figure 11.Overview of the second dart leader and the blooming that precedes it, with sources colored blue to red according to time.(a) altitude vs time,(b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations, and the time and location of NLDN strokes.Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes.

Figure 12 .
Figure 12.Overview of the third dart leader and the blooming that precedes it, with sources colored blue to red according to time.(a) altitude vs time,(b) altitude vs. east/west position, (c) Plan view (north/south vs. east/west), (d) altitude vs. north/south position.The 3DINTF VHF sources are shown along with the location of the INTF01 and INTF02 stations, and the time and location of NLDN strokes.Dashed vertical lines in (a) are 3DINTF IC strokes, solid lines are CG strokes.

Figure 13 .
Figure 13.Plot of velocity versus distance propagated along the channel for all of the dart leaders and K leaders.Showing K leader 1 (cyan), K leader 2 (red), dart leader 1 (black), dart leader 2 (blue), and dart leader 3 (green).Zero propagation distance is arbitrarily set to be the start of the last dart leader.

Figure 14 .
Figure 14.Plot showing how the channel segments were aligned for each dart leader and K leader.Showing K leader 1 (cyan), K leader 2 (red), dart leader 1 (black), dart leader 2 (blue), and dart leader 3 (green).Zero propagation distance is arbitrarily set to be the start of the final dart leader.The gap between 17-18 km is caused by a null in antenna sensitivity directly above the INTF01 station.

Figure 15 .
Figure 15.Charge structure as identified by the LMA for flashes happening between 09:50:00 UTC and 10:20:00 UTC near Langmuir Lab on 2018-10-23.Negative charge is shown as blue, positive as orange/red.The path of the dart leaders in the 10:00:09 flash is overlaid in black.

Figure 16 .
Figure 16.Plot showing the correspondence between the large dip in velocity at about 11 km along the channel and the baseline between the two INTF stations for the second dart leader.North/south and east/west location of sources and stations (a), and velocity vs propagation distance (b).Sources are colored by velocity.

4.
Profiled in time and space the velocities of dart leaders and K leaders in the cloud

Figure A1 .
Figure A1.A diagram inspired by(Liu et al., 2018) for the double theodolite triangulation, showing station 1 (A), station 2 (B), their azimuth and elevation measurements (AZ1, EL1, and AZ2, EL2).Point D is the point of closest approach along the line of sight for station 1, and point C is the point of closest approach along the line of sight from station 2. Point S is the source location.

Table 1 .
The range of maximum and minimum dart (and K) leader velocities reported by selected studies, and compared to this study.
also observed that dart leaders in natural lightning tended to slow down as