www.ann-geophys.net/25/2175/2007/

We present refinements of a method of iono- spheric D-region sounding that makes opportunistic use of powerful (10 9 -10 11 W) broadband lightning radio emissions in the low-frequency (LF; 30-300 kHz) band. Such emis- sions are from "Narrow Bipolar Event" (NBE) lightning, and they are characterized by a narrow (10-µs), simple emis- sion waveform. These pulses can be used to perform time- delay reflectometry (or "sounding") of the D-region under- side, at an effective LF radiated power exceeding by orders- of-magnitude that from man-made sounders. We use this op- portunistic sounder to retrieve instantaneous LF ionospheric- reflection height whenever a suitable lightning radio pulse from a located NBE is recorded. We show how to correct for three sources of "regular" variability, namely solar zenith angle, radio-propagation range, and radio-propagation az- imuth. The residual median magnitude of the noise in re- flection height, after applying the regression corrections for the three regular variabilities, is on the order of 1 km. This noise level allows us to retrieve the D-region-reflector-height variation with solar X-ray flux density for intensity levels at and above an M-1 flare. The instantaneous time response is limited by the occurrence rate of NBEs, and the noise level in the height determination is typically in the range ±1 km.


Introduction and background
Structure and variability of the lower ionosphere ("D"-layer, altitude <90 km) are difficult to monitor, for two main reasons: First, the low plasma frequency at D-layer heights necessitates use of LF (low-frequency; 30-300 kHz) or VLF (very-Correspondence to: A. R. Jacobson (abramj@u.washington.edu) low-frequency; 3-30 kHz) techniques. For these frequencies, suitable radio-sounder antenna dimensions (on the order of the half-wavelength, or λ/2=5 km for f =30 kHz) cannot be realized in practice, so that one must make do with inefficient, low-gain antennas.
Second, the D-layer's high electron-neutral collision rate causes high signal loss in D-layer radio sounding, further worsening the signal-to-noise problem already implicit in the low-gain antennas.
A feasible alternative to wideband-pulse radio sounding of the D-layer is opportunistic reception of extremely powerful, narrow-band VLF transmissions that are transmitted for other purposes (see, e.g., Piggott et al., 1965, and references therein; Thomson and Clilverd, 2001; Thomson and Rodger, 2005). Thomson's and his colleagues' use of this approach over long paths in the day-lit hemisphere has provided a remarkably accurate proxy for instantaneous solar X-ray flux density. Reception of extremely powerful, narrow-band VLF transmissions has served also in the study of local ionospheric perturbations, e.g. "Sprites" (Dowden et al., 1996), D-region heating (Cho and Rycroft, 1998), and the optical signature of D-region heating, known as "elves" (Fukunishi et al., 1996).
Another approach retains wideband features, but uses opportunistic reception of the natural VLF emissions ("sferics") from lightning strokes (Cummer, 1997;Cummer et al., 1998). This exploits the large effective antenna dimensions (>1 km) and extremely high peak powers (∼GW) of natural lightning-discharge currents. Neither these large antenna dimensions, nor these high peak powers, could be straightforwardly achieved in active radio-sounders ("ionosondes") at VLF. Using opportunistic reception of wideband VLF sferics, the mean D-layer electron-density profile parameters have been obtained (Cummer, 1997;Cummer et al., 1998), and transient D-layer responses to both lightninginduced heating/ionization (Cheng and Cummer, 2005) and energetic-particle precipitation have been observed (Cheng et Published by Copernicus Publications on behalf of the European Geosciences Union. 2176 A. R. Jacobson et al.: Low-frequency ionospheric sounding with NBE lightning radio emissions al., 2006). The method uses averaging (over several serics) of the mean sferic spectrum, from which structural parameters of the D-layer can be retrieved.
Our present work is intended to be an extension of the general approach of Cummer and of Cheng, that is, using opportunistic reception of powerful lightning sferics. Our present work will depart from the earlier work in three respects: First, we will use the LF spectrum and a class of temporally narrow lightning strokes, for higher time resolution. Second, we will emphasize lightning-to-receiver paths that tend to be shorter (200-800 km), so that the ionospheric firsthop echo and the ground wave are both received, but are well time-resolved from each other. Third, by recording both the ground wave and the resolved ionospheric echo, we will attempt to retrieve D-layer structural parameters from a single sferic, potentially allowing us to track the time-development of fast perturbations.
A key goal of our present work is to develop a technique suitable for diagnosis of local (100s of km horizontal), rapidly-varying (seconds to hours) ionospheric perturbations. Examples of such transient and local electron-density perturbations include those caused by energetic-particle precipitation  and by powerful lightning discharges's radiated fields. Local (horizontal scale <1000 km) disturbances of these sorts have been inferred from longrange VLF-beacon-signal amplitude and phase perturbations in work pioneered by the Stanford group (see references in, e.g., Lev-Tov et al., 1995). Their work strongly suggests that a steep-incidence LF sounder operating underneath a localized disturbance would observe transient changes in D-layer reflection height. This motivates the present, preliminary article, which characterizes the noise and variability in the opportunistic LF sounding method.

General approach
Our technical approach is based on the recording of narrow LF signals from the unique lightning stroke known as an "NBE", or Narrow Bipolar Event (Jacobson, 2003;Jacobson and Light, 2003;Smith et al., 1999). NBEs are received by LF recorders of the transient vertical electric field at ground level, e.g. LASA, or Los Alamos Sferic Array (Smith et al., 2002). The pertinent facts about NBEs for this article include: (a) The NBE main pulse is fast (∼10-µs width) and powerful (peak effective radiated power in range 10 9 -10 11 W). This provides an opportunity for high-time-resolution timedelayed-reflectometry (TDR) from the ionosphere.
(b) NBEs are a common (but not ubiquitous) intracloud discharge in many electrified storms (Jacobson et al., 2007;Jacobson and Heavner, 2005;Suszcynsky and Heavner, 2003;Wiens et al., 2007). The incidence of NBEs has been studied for storms both in Florida (Jacobson et al., 2007;Jacobson and Heavner, 2005;Suszcynsky and Heavner, 2003) and in the United States Great Plains (Wiens et al., 2007). NBEs do not appear to constitute a consistent proportion of total lightning in a storm. NBEs' incidence is, instead, extremely inconsistent, although there is a tendency to be most common in extremely intense thunderstorms (Wiens et al., 2007). Some storms have no observable NBEs, while other storms have as many as 20% NBE incidence. The pattern is extremely inconsistent (Jacobson et al., 2007;Jacobson and Heavner, 2005;Suszcynsky and Heavner, 2003;Wiens et al., 2007).
(c) Recording of the NBE ground wave, followed by both the ionospheric and the ground-then-ionospheric echoes, allows retrieval of the source and "ionospheric reflector" heights (Smith et al., 2004), for lightning-to-receiver ranges 200-800 km.
Using NBEs as a transmitter-of-opportunity for measuring reflective height of the D-layer, we will attempt to answer the following questions: (a) What are the major sources of reflection-height variability, and can this variability be determined and corrected?
(b) Can the residual variabilities after correction be interpreted in terms of solar-flare activity during the last solar maximum?
(c) What is the sensitivity of the method? We will address these three questions within the constraint of a sharp-boundary (i.e., mirror-like) reflection model. This model obviously lacks the richness of information provided by the full-wave, diffuse-profile VLF model used in the works of Thomson et al., of Cummer, and of Cheng et al. They are able to estimate the ionospheric gradient length, while we cannot do so under the constraints of our simple model. A future publication will apply a full-wave ionospheric-reflection model to our LF data recordings.

Reflection model and height retrieval
Typical LASA multi-station recordings of a single NBE near Florida are shown in Fig. 1. Each panel contains the recording from a different station (e.g., "ta") at a different range (e.g., 300 km) from the NBE location. The NBE location itself has been determined by time-difference of arrival (TDOA) by LASA (Smith et al., 2002). The time scale is zeroed for each station at that station's trigger time. The two vertical dashed lines indicate the expected echo-arrival times for, respectively, the ionospheric and the ground-then-ionospheric echoes, based on a model which assumes a discrete mirror at an ionospheric height H i =61.2 km and a source height H s =13.1 km. The standard LASA height-retrieval (Smith et al., 2004) uses this approach and has obtained serviceable estimates of H s for lightning research. The standard retrieval idealizes the ionosphere as a sharp boundary between vacuum (underneath) and a finiteconductivity (e.g., 2.2 micro-mhos/m, as recommended by VLF experience (Wait and Spies, 1964)), isotropic conductor. Inadequacies of the sharp-boundary reflection model cause some inaccuracy in retrievals of H i but only small errors in H s , which is of primary interest for lightning research.
In the present approach we will modify our retrieval of ionospheric height and source height in three key ways: Modification 1: We will let each qualifying station provide a separate retrieval for the two heights. It turns out that to the accuracy we presently seek in H i , the diffuse (as opposed to sharp-boundary) ionospheric profile results in a systematic and monotone decrease of retrieved H i with range. We will continue to model a sharp-boundary, isotropic-medium reflection, but will let each range make its own estimate. We will then find an empirical regression of H i with respect to range, and correct all H i retrievals to "short" (200-km) range. Figure 2 shows a close-in view of the two echoes for the data in the second through fifth panels of Fig. 1. In each panel in Fig. 2, the light curve is the data for that station, while the heavy curve is a model of the data, based on a retrieval (see legend in each panel) for that station only. Modification 2: Our model imposes a fixed low-pass filter on the input (groundwave) waveform. Rather than generate an echo model using a partially-conducting reflector (Smith et al., 2004), our present approach is to impose a frequency roll-off that is found to work better than the approach of Smith et al at mimicking the data, over the majority of events analyzed. By this, we mean that the modified approach results in a modeled ionospheric echo that more closely resembles the observed ionospheric echo, relative to that of an echo model using a partially-conducting reflector (Smith et al., 2004). Each echo waveform in Fig. 2 is derived from the groundwave waveform E gnd by a constant-phase filter F (f ) in the frequency domain: followed by return to the time domain and finally imposition of a time delay for (first echo) the lightning-ionosphere-receiver path, or (second echo) the lightning-ground-ionosphere-receiver path. The filter is of the form . (2) This is a low-pass filter with a knee at 60 kHz, and a half-width of 30 kHz. This filter suppresses the higher frequencies ( LASA model (given by reflection off a conductor with conductivity 2.2X micro-mhos/m).
Modification 3: The standard LASA retrieval of heights (Smith et al., 2004) performs a correlation between the modeled and observed echo waveforms, to determine the "best" echo delays. In the present work, which seeks to minimize systematic errors in H i , we avoid correlation between two waveforms that are intrinsically somewhat dissimilar, and instead compare derived peak times from the observations and from the model. Given the observed and modeled double echo, we search jointly in H i and H s to minimize the sum of the square of the residuals between modeled and observed peak times of the each of the two echoes. This requires some care, because the data echo is corrupted by additive noise; see, e.g., the lower two traces (c and d) in Fig. 2. This noise makes it imperative not simply to identify the echo peak with the highest point in the noisy data, but rather to fit a smooth function to the echo and then find the peak of the best-fit smooth function. For each echo peak, we fit the upper half of the echo with a fourth-order polynomial, and we then find the peak of the best-fit polynomial. This results in an estimate of echo peak that is more robust against noise than would be the case if we simply found the raw echo data's highest point.
The modified height-determination method just described incurs errors due to the imprecision of the fitted peak delay for the ionospheric echoes. As part of the automated calculation, we tabulate the estimated error in the peak location, based on the goodness of fit to the smooth polynomial model. Typical estimated errors in the fit are in the range 2-6 µs. This causes reflection-height propagated errors in the range 0.5-1.5 km, due to uncertainties in the smooth-model fit.

Building a statistical database of H i retrievals
We stored all LASA (Smith et al., 2002) NBE data from three years (2000,2001,2002) near the most recent solar maximum. This data was for the Florida subarray of LASA, centered on 28 deg N and −81.5 E. We restricted NBEs to those lying within 400 km from this array center. We then applied automatic height retrievals to all of those data, along the lines of Sect. 2.2 above. Several obvious quality factors were then applied to select against noisy data or against echo waveforms that were poorly matched by the model. After this automated, criteria-driven selection, we are left with 65 000 height retrievals, peaked in the year 2002, with fewer in 2001 and still fewer in 2000. (This year-to-year variation was wholly due to our reducing the trigger threshold as the system became more capable of handling greater volume of data; the year-to-year variation was not geophysical.) There are on average about two accepted height determinations per NBE, as opposed to the example in Fig. 2, which provides four height determinations. Figure 3 shows the distribution of our 65 000 height retrievals over each of three types of parameters. Figure 3a shows the distribution of height determinations over great-circle-path range from the receiver to the NBE location. The strong rolloff at range >300 km is due in part to the two selection criteria we enforce: First, we select only for data lying within 400 km from the array center (28.5 deg N, −81.5 deg E). Second, we perform height retrievals only in the distance interval 200 km<range<800 km from the recording station (Smith et al., 2004). (The rationale for this distance interval is as follows: For shorter ranges than ∼200 km, the ionospheric echo signal is suppressed due to antenna nulls at zenith of both the vertical-monopole-on-groundplane receiver antenna, and the vertical-dipole NBE transmitter (Smith et al., 2002(Smith et al., , 2004  1999). On the other hand, for longer ranges than ∼800 km, the ionospheric echo becomes poorly resolved (in time) from the ground-wave signal (Smith et al., 2004)). Together these two selection criteria contribute to the behavior of Fig. 3a.

Geophysical variabilities
Because data were gathered for several months each year (late Spring through early Fall), there is not a unique relationship between local solar time and solar zenith angle. We shall see if zenith angle, without regard for day of year, can account for local-time variability in H i . Figure 3b shows the distribution of zenith angle (at the ground). The obvious control by zenith angle derives from the main source of ionospheric electron density, namely, photoionization by solar extreme ultra-violet radiation.
We use data from the National Oceanographic and Atmospheric Administration (NOAA) GOES-10 satellite on solar X-ray flux density, in 5-min averages. The X-ray flux density (W/m 2 ) is monitored by GOES-10's SEM payload both in a "longwave" channel (1-8 µm) and in a "shortwave" channel (0.5-4.0 µm), called X l and X s , respectively. Each height retrieval is associated with a pair X l and X s from the 5-min period in which the NBE occurred. Figure 3c shows the distribution of X-ray flux density tabulated for those 5-min periods. Also shown are the conventional "C-flare", "M-flare", and "X-flare" thresholds for X l . Figure 4 is a scatter-plot of the retrieved H i for all 65 000 height retrievals as a function of instantaneous solar zenith angle (at the ground). Obviously zenith angle accounts for the majority of natural variability of H i . We define the minimum range (200 km) as R 0 . In order to separate the H i variability due to zenith angle, from H i variability due to range, we take all the points in Fig. 4, select only those shortrange retrievals within an arbitrarily selected range <250 km from the recording station, and fit the zenith dependence using only these short-range height retrievals "H i0 ". Recall that we do not perform any height retrievals for range <200 km, so the bottom 50 km of range lies in the domain 200<range<250 km. This is shown in Fig. 5. The fit of H i0 (Z) to zenith angle Z is shown as a line, and is based on a model

Empirical correction for main variabilities in H i
where the best-fit parameters for amplitude, center, width, and offset are: A=8.9 km, Z cent =78 deg, Z=22 deg, and H=79 km. We can now apply the regression on Z to the ionospheric height, then use all the data together, regardless of zenith, to regress the range dependence of retrieved height. Let us define "scaled height" and "scaled range" as normalized by the best-fit ionospheric height at that solar zenith angle but at shortest range (200 km): scaled range = range/H i0 (Z) (4b) Figure 6 shows the scatterplot of scaled height versus scaled range. The decrease of scaled height with increasing range is due to the departure of the D-layer reflector from an ideal sharp boundary; it is, in reality, diffuse. As range increases, then so does angle of incidence, and thus the penetration of the radio wave into the medium decreases.
Range is obviously correlated with ionospheric angle-ofincidence: If the range is zero, the angle-of-incidence is zero, and as the range increases, the angle-of-incidence also monotonically increases. The D-layer underside is not a sharpboundary reflector but rather is a diffuse profile, with electron www.ann-geophys.net/25/2175/2007/ density tending to increase smoothly versus altitude in the range 70-90 km. This implies deeper penetration for LF radiation at smaller incidence angles, and shallower penetration for grazing incidence (i.e., larger incidence angles). We capture that trend via the sorting parameter of range, but it could equally be sorted via the parameter of angle-of-incidence. We fit scaled height to a quadratic in scaled range; the 2ndorder-fit coefficient vector is: (1.04, −0.014, 5.6 × 10 −5 ) This best-fit ionospheric height will henceforth be called "H i,fit ", and we will deal with a dimensionless height (difference) corrected for solar zenith angle and for range: This dimensionless height difference is plotted as a scatterplot versus azimuth (deg clockwise from N, looking from the lightning location to the receiver) in Fig. 7. The reason for plotting the data versus azimuth is that the ionosphere is an anisotropic medium due to the geomagnetic field (Cummer, 1997;Cummer et al., 1998;Pitteway, 1965), and it is in principle possible that the observed reflection heights might have some systematic, detectable dependency on azimuth. Figure 7 reveals a slight azimuth dependence, subtle compared to the dependencies on either zenith or range. Nonetheless we empirically fit the data in Fig. 7 to a function of the form δ(dimensionless height) = A a X cos(az − az 0 ) The best-fit amplitude and center are, respectively, A a =0.012 and az 0 =66 deg. The fit is shown in Fig. 7 as a solid curve. From now on we include this correction δ(dimensionless height) into the dimensionless height (see Eq. 5).  . 7. Scatter-plot of dimensionless height corrected for zenith and range (first differenced from range-controlled fit, then divided by the modeled height at the appropriate zenith angle but at short range) versus radio-propagation azimuth. The solid curve is a fit to a diurnal cosine, parametrized by amplitude and phase (see text).
The median absolute value of the final residual to the fit is about 0.015. That is, after removal of the "predictable" variabilities due to zenith, range, and azimuth, the residual variations in dimensionless height difference may be described by a median magnitude of 0.015, corresponding to ∼1 km unexplained median variability in daytime (when H i0 (Z=0) ∼70 km).

Simple regression
The 5-min-averaged GOES-10 long-and short-wavelength X-ray flux densities (see Fig. 3) will now be compared with their contemporaneous LASA height retrievals. Figure 8 shows scatter plots of dimensionless height difference (see Eq. 5) versus instantaneous long-wavelength X-ray flux density, for daytime (Z<70 deg, Fig. 8a) and for nighttime (Z>120 deg, Fig. 8b). The points where the X-ray flux density attains flare level "M" or "X" are marked with larger square symbols. Figure 8 excludes the dawn/dusk region (70<Z<120 deg), because the D-region is in rapid photochemical transition during dusk or dawn and may not be suitable to demonstrate a stationary height/X-ray-flux relationship. From Fig. 8 we observe: (a) The solar-X-ray influence is much more apparent during sunlight (Fig. 8a) than during darkness (Fig. 8b). This is logical, since only during daylight can the direct X-ray photo-ionization process occur. If there is a nocturnal height effect correlated with solar X-rays, the effect could be only an indirect one, communicated to the nighttime hemisphere by magnetosperic substorm dynamics.
(b) Recently long-path VLF propagation phase has been shown to follow an exceptionally clear and reproducible proportionality to the logarithm of X-ray flux density (Thomson and Clilverd, 2001;Thomson and Rodger, 2005;Thomson et al., 2004). The main path studied was the Seattle (USA) to Dunedin (NZ) propagation of the 24.8 kHz radio beacon "NLK" transmitted near Seattle. During times when the path was sun-lit, the observed propagation phase and amplitude responses to solar-X-ray flux density were analyzed in terms of an exponential model of the ionospheric boundary, with two best-fit parameters: a height H i and a logarithmic gradient β. In the propagation model employed by Thomson and colleagues, the path attenuation is primarily influenced by β, while the path phase is mainly influenced by H i . The sensitivity of H i to the logarithm of long-wave X-ray flux density can be obtained by scaling from Fig. 11 in Thomson and Rodger (2005): dH i /d(log 10 (X l )) = −5.1 km (7) That result is for extremely long-range propagation (Seattle to Dunedin) at lower frequency (24.8 kHz) than our LF measurements. Nonetheless, let us compare our data to the scaling of Eq. (7). For a mean daytime, short-range height H i0 (Z=0) of 70 km (see our Fig. 5 above), the scaling of Thomson et al would imply that our dimensionless height difference (see Eq. 5 above) should vary with X-ray flux density according to d((H i − H i,fit )/H i0 (Z))/d(log 10 (X l )) = −0.073 at least on the high-flux tail where the radiative effects exceed our system noise and other unexplained variabilities. This implied dependency is shown as the solid line in Fig. 8a. Remarkably, our daytime data, which is both noiser than the long-range VLF beacon method and involves a completely different approach to height retrieval, appears to be well described by the VLF (Thomson and Rodger, 2005) scaling. Figure 9 shows similar scatter plots as does Fig. 8, but for the short-wave X-ray flux density. The symbol size is determined by long-wave X-ray level, but otherwise the abscissa value is for short-wave X-rays. Both Figs. 8 and 9 show dimensionless-height data only when the GOES X-rayflux-density data is available. Occasionally the GOES data is indicated to be "bad" and is easily recognizable because in the archive it is pinned at a non-physical level. For example, the three low points near (1-2)×10 −6 (W/m 2 ) in Fig. 9b are lacking in Fig. 8b. That is because the X l values are indicated to be "bad" in the NOAA archive.

Short-term clusters of data
Frequently a lightning storm in this Florida array area will last for up to a few hours without migrating more than about 50 km. This allows us to identify clusters of events from the same receiver station and from the same lightning storm. Using such a cluster allows us to hold constant the observational geometry, and thus to monitor possible ionospheric height changes with minimal systematic variation. If the corrections for range and azimuth (see Sect. 3 above) were perfect, then there would be no added value in looking at clusters of data from the same storm and the same station. However, we do not believe that those corrections are such as to allow us to skip this prudent step.
We divide all time of these observations, starting with the first observations in 2000, and finishing with the last observations in 2002, into time windows of 6-h duration each, advanced by 3 h (for 50% overlap). Consider the horizontal vector leading from a lightning stroke to the receiver station. There is one such vector for each height determination. During each time window, we group all lightning-to-receiver vectors (km Eastward, by km Northward) into square bins of 50-km×50-km size, and these in turn are spaced at 25km intervals along both the East-West and North-South axes, for 50% overlap along each axis. We then iteratively search for the bin containing the most vectors. When we find that bin, we label its occupant height determinations as a new cluster of data. We then take those lightning events off the list of candidates for the next iteration, and repeat the process. The search is iterated until we can no longer find any 50-km×50-km bin containing at least 2 not-yet-used lightning events vectors. At this point we have collected all the useable clusters of lightning storm/receiver vectors. Within each of the vector bins, the observing geometry is almost constant. Within each vector bin, we subtract the first retrieval of dimensionless-H i (see Eqs. 5 and 6 above) from all subsequent dimensionless-H i retrievals within that bin, in order to null-out any residual bias from the cluster's observing geometry. We refer to this product as "relative dimensionless height". In this manner we find 13 257 clusters of data. Due to the temporal overlapping of the technique, this total count of clusters turns out to be about twice the number of truly independent storm/baseline combinations. In order not to miss vector clusters at the edge of regular time windows, we overlap the windows, and this results in a tendency toward double-counting.
Sometimes a data cluster's time window includes part or all of a solar-flare X-ray transient. Figure 10 shows such a view of a flare, in this case an X-1 flare on 26 July 2002. The height determinations span about two hours and are irregularly spaced, due to the irregular occurrence of NBE lightning events. Figure 10a shows the 5-min-averaged X l flux density from GOES-10 at the times of the lightning events. Figure 10b shows the relative dimensionless height, that is, the dimensionless height subtracted from the first dimensionless height in the cluster. Thus the first "relative" dimensionless height in the cluster is, by definition, zero. Figure 10c shows the solar zenith angle at the ground.
The earlier finding that the D-layer "height" in VLF waveguide propagation scales linearly with (minus) the logarithm of X l (Thomson and Rodger, 2005;Thomson et al., 2004) suggests that we define a cluster-relative flux density as the logarithm of X l minus the logarithm of X l of the cluster's first point. We then can plot all daytime cluster-relative dimensionless heights versus the corresponding cluster-relative X-ray flux densities, as in Fig. 11. Most of the clusters do not contain a dramatic X-ray transient, and this majority of points constitutes the ball of data near the origin. In order to highlight cluster-relative X-ray flux-density magnitudes exceeding 0.5, we use a larger square symbol for those points. We perform a regression on only these highlighted points; the best-fit slope is −0.079. This is to be compared with the slope of −0.073 implicit in Fig. 11 from Thomson and Rodger (2005). Our slope determination is certainly much noisier than that of the earlier long-path VLF work, as is expected from a localized measurement such as steep-incidence sounding. Nonetheless we find it noteworthy that the LASA statistical regression of LF, steep-incidence height sensitivity to X-ray flux density agrees so well with the more precise regression from long-path VLF.

Discussion
As mentioned in Sect. 2.1 above, we posed three questions to be answered in this article. The next three subsections deal in sequence with those questions.

Major sources of regular reflection-height variability, and their correction
By "regular" variability, we mean a variability of the reflection height straightforwardly related to, and predictable by, "independent" geophysical parameters. Of these regular variabilities in reflection height, we find the most significant to be due to solar zenith angle. Over several months, the D-layer insolation (at a set local time) is not constant, and zenith angle-which alone determines the insolation geometry -is preferable to local time as a proxy for insolation. The next most significant regular variability is controlled by lightning-to-receiver range. Finally, the weakest regular variability is controlled by the propagation azimuth. All three of these regular variabilities are appropriately best-fit and then corrected in the data. We combine all corrected data in a "dimensionless height" that is scaled by the best-fit short-range height for the particular solar zenith angle.

Residual variabilities and solar-flare transients
We find that the response of dimensionless height to solar Xray flux-density transients is noisy but consistent with the far less noisy predictions based on long-range VLF propagation in the Earth-ionosphere waveguide (Thomson and Rodger, 2005;Thomson et al., 2004).

Implications for the sensitivity of the method
After correction for the "regular" variabilities, and for variability controlled by solar X-ray intensity, we are left with a Fig. 11. Scatter plot showing dimensionless height for all heights in all clusters (after subtracting the dimensionless height of the first point in each cluster), versus the logarithm of the long-wave X-ray flux density relative to the first point in each cluster. Points whose logarithm of the flux density differs from its cluster's first point's flux density by more than ±0.5 are marked by square symbols (see text).
residual median dimensionless-height variability of around 0.015. For a mean daytime short-range height of 70 km, this implies a residual median noise magnitude of ∼1 km. This is on the order implied by the estimated statistical errors from the smooth-polynomial fit to the echo peak delay (see Sect. 2.2): The range of statistical errors propagates to a height error of 0.5-1.5 km.
We were able to estimate the sensitivity in daytime because of the availability of solar-flares as a calibration source. During nightime, we do not have such a well-established and hemispherically homogeneous calibration source. However, we suspect that the nighttime sensitivity should be no worse, and possibly better than, the daytime sensitivity, due to a simple propagation effect: During nighttime, the LF propagation is less lossy, and the ionospheric reflection coefficient is enhanced (Pitteway, 1965) relative to daylit condition. This is the same reason why short-range VLF studies of ionospheric transients have been conducted during nighttime (Cheng and Cummer, 2005;Cheng et al., 2006).