Lower Ionosphere Effects on Narrowband Very Low Frequency Transmission Propagation: Fast Variabilities and Frequency Dependence

The man‐made narrowband very low frequency transmission from Rhauderfehn, Germany, is studied using high accuracy, microsecond time resolution measurements in Bath, UK. The high time resolution enables a novel comparison of the measurements with a detailed simulation of the transmitted signal. It is found that the wave propagation frequency response exhibits nonlinear amplitude and phase changes with frequency over the narrow transmission bandwidth during the time of the measurements (13 May 2011 15:00:03–15:00:09 UTC). The high time resolution also enables measurements of fast variabilities in the wave propagation <5 ms. The fast propagation variabilities are likely to originate from integrated ionospheric variability over the propagation path or fast ionospheric processes. The wave propagation frequency response measurement has potential benefits in the study of the lower ionosphere, in particular during highly variable perturbations such as those caused by lightning.


Introduction
The lower ionosphere D and lower E regions from ∼50 to 90 km altitude exhibit significant spatial and temporal variations driven by a number of phenomena such as lightning discharges, neutral atmosphere and solar variabilities, and space weather. Such events are difficult to observe with sufficient coverage due to the altitude range which can only be reached from the ground by rockets but is too low for satellites. Remote sensing methods can address this challenge, and the propagation of very low frequency (VLF) radio waves (3-30 kHz) in the Earth-ionosphere waveguide provides a convenient way to observe conductivity perturbations in the lower ionosphere (e.g., Barr et al., 2000;Inan et al., 2010;Rodger, 2003;Silber & Price, 2017, and references therein).
In particular, thunderstorms produce localized conductivity perturbations via a variety of mechanisms. Lightning-induced electron precipitation events have been identified from disturbances on man-made narrowband VLF transmissions (e.g., Helliwell et al., 1973). Another class of disturbances was identified with similar VLF signatures but without observable delay from the causative lightning discharge-suggesting more direct lightning-ionosphere coupling (e.g., Armstrong, 1983;Inan et al., 1988). These "early' disturbances sometimes occur simultaneously with transient luminous events such as sprites and elves, and it is discussed whether the two are causally linked (e.g., Haldoupis et al., 2013;Inan et al., 2010;Rodger, 2003). Candidate mechanisms for these early disturbances are, for example, quasi-electrostatic effects from lightning discharge-induced charge moment changes and lightning electromagnetic pulse. Typically, it takes ∼5-200 s for early disturbances to recover to undisturbed levels, but early disturbances with longer recoveries have recently been identified (Cotts & Inan, 2007;Haldoupis et al., 2012) leading to various studies on potential mechanisms for persistent ionospheric perturbation (Gordillo-Vázquez et al., 2016;Haldoupis et al., 2013;NaitAmor et al., 2013;Salut et al., 2012). Additionally, thundercloud electrostatic fields have an effect on the lower ionospheric conductivity profile (Lay et al., 2014;Pasko et al., 1998;Salem et al., 2016). The monitoring of man-made narrowband VLF transmissions is thus relevant to understand thunderstorm impact on the lower ionosphere and the different physical mechanisms. The complexity of the possible shape and uncertainty about the size and duration of the perturbations and their causative mechanisms drive the need for more detailed measurements and modeling. For example, holographic strip imaging of perturbations using a number of receivers arranged in a line perpendicular Radio Science 10.1002 to the propagation path was attempted to determine the extent of lightning ionospheric perturbations (Chen et al., 1996;Johnson et al., 1999;Moore et al., 2003) and various studies have considered Gaussian (Moore et al., 2003;NaitAmor et al., 2017) and column (Hardman et al., 1998;Marshall et al., 2006) shapes.
Existing observations of VLF transmission disturbances are mostly single receiver measurements with ∼20 ms resolution. The frequency modulation of the VLF transmissions mean that the transmission switches between two sideband frequencies around the carrier frequency (described further in section 3)-a characteristic that has been used to make dual frequency observations (Dowden & Adams, 1989). However, most observations approximate the VLF transmissions to a monotone at the carrier frequency. Recently, high accuracy, microsecond time resolution small aperture array measurements reveal that the propagation of the VLF transmission signal exhibits variabilities on time scales ≪20 ms , thus sparking the motivation to investigate the cause of the propagation variabilities. These propagation effects should be observable as VLF transmission disturbances on a single receiver. This work presents a method to derive the VLF transmission disturbances on short time scales which elucidates further details of the multipath propagation and its fast variabilities. The purpose is to push existing techniques for lower ionospheric observation by showing that more physical detail can be inferred from detailed measurement.
The method is described in section 3 and applied to an equivalent single receiver recording derived from the array data so that the propagation variabilities may be investigated further. The method is shown to be able to extract the variabilities from the equivalent single receiver data in section 4 and that some of the variability is closely associated with the changing frequency content of the transmission. The array processing of multiple high bit depth, high sampling rate recordings to observe the wavenumber vectors is computationally expensive and so is currently impractical for long-term monitoring. Gathering such data is also logistically difficult; array baselines need to be in the kilometer range to make sensible observations of radio waves with wavelengths that are sensitive to lower ionospheric conductivities. While these obstacles are not insurmountable, it is more convenient to use a single receiver. The data gathered with the array were during an arbitrary period when ionospheric conditions were stable, so it is unknown how the propagation variabilities will change with the ionosphere. Thus, analysis of single receiver data (a different data set to the array) when the terminator is moving along the propagation path was performed with the method. The results are shown in section 5 and illustrate the effect of the diurnal change in the ionosphere on the frequency dependence of propagation.

Experimental Data
In this study, we use a recording of the VLF transmission from Rhauderfehn, Germany, at 23.4 kHz (denoted DHO), averaged from a 10-receiver array situated within a square kilometer area near Bath, UK (Mezentsev & Füllekrug, 2013). The reason for using the array data set is to investigate the propagation variabilities concluded from the array analysis by . This data set covers a short arbitrary period when ionospheric conditions were stable.
Each receiver is a wideband vertical electric field sensor that samples every 1 μs, with an amplitude resolution of ∼35 μV/m and timing accuracy of ∼20 ns (Füllekrug, 2010). The accuracy of the array measurements has been shown to be close to instrumental limits (Füllekrug et al., 2014, section 6). The time of the recording analyzed in this study is a 6-s period from 15:00:03 to 15:00:09 UTC on 13 May 2011. The size of the array (∼1 km 2 , with the largest distance between a pair of receivers being ∼1,260 m) is unique in the study of subionospheric propagation-it is much smaller than the distance from either the transmitter (∼700 km) or the wavelength (∼13 km). This means that the wavefield traveling over the array can be approximated with a plane wave. The recordings can be aligned in time by multiplying with a reference plane wave, e ikr l , where k is the wavenumber vector and r l is the position vector of each receiver with index l in the array.
Additionally, each recording is frequency shifted so that the angular transmission carrier frequency c is at baseband after multiplying with a phasor e −i c and subsequently a low-pass filter is applied to isolate the transmission of interest. This produces narrowband complex analytic signals, y l (t), of the baseband transmission from which the amplitude and phase can be computed for each individual sample. These signals are finally averaged across all receivers in the array. (1)

10.1002/2017RS006456
The proximity between the receivers in the array (<1,260 m) means that the signal sampled by each receiver is from approximately the same propagation path. Thus, the result can be considered as a recording from an equivalent single receiver.
Man-made sources of interference present in the data (e.g., the sources between ∼20 and 250 kHz documented in  are already effectively removed by the narrowband filter applied. Natural sources of interference in the data are weak lightning signals from far away storms (>900 km away; Füllekrug et al., 2016). Most of the electromagnetic energy radiated by lightning is centered around 10 kHz. Additionally, higher frequencies from sferics attenuate faster with distance so that the far away sferics in the filtered data are of small amplitude. This noise is further reduced by the averaging process, although the effect is small-all the results obtained with the equivalent single receiver data are similar to results computed with only the data from one receiver in the array. This may not always be the case in general, and there may be periods when there is strong sferic interference. The effect of strong sferic interference is considered in section 3.
The derivation of ionospheric perturbation signals from disturbances on VLF transmissions requires an understanding of the source field-the detailed model of the VLF transmission used in this work to extract the transmission disturbances is described in section 3. A key assumption of this method is that the transmission at source is exactly as described by the model. The assumption is largely valid as transmissions follow a known modulation scheme. In practice, VLF transmissions such as DHO are operated for broadcast communication purposes and only to the accuracy required by their specific purpose and thus have unknown imperfections.
Existing long-term observations note that VLF transmissions have good amplitude and phase stability, although this is not always the case. There may be periods when there are random phase jumps in the transmission, or the transmission clock is not locked to a frequency standard and the transmission experiences a frequency drift . For example, the frequency of DHO has been reported to vary by a small fraction (<100 μHz) of its published center frequency of 23.4 kHz (Thomson et al., 2017). In this case, the phase disturbance calculated (see section 3) would ramp in time with a gradient equivalent to Δ , the angular frequency difference between the assumed center frequency and the actual transmitted center frequency. Over a 1-s processing period, its effect is <0.036 ∘ , that is, around a magnitude smaller than the variabilities seen in the results and so is not significant in this work. There are also no amplitude or phase jumps in the data used.

Measurement of VLF Transmission Disturbances
A modulation scheme common among VLF transmissions is Gaussian Minimum Shift Keying (GMSK)-a spectrally efficient variant of Minimum Shift Keying (MSK). MSK transmissions can be described in complex notation as Ae i (t) e i c t , a signal with time-varying phase (t) (which contains the message) but constant amplitude A. The amplitude of the received transmission is understood as an indicator of the ionospheric disturbance signal, assuming that the variations of the source amplitude of the transmitted signal are small.
The time-varying phase from the transmitted message has to be removed to get the phase disturbance.
where m(t) is a sequence of symbols represented by either positive or negative unity values. The length of each symbol in the sequence is T s = 1∕f s , where f s = s ∕2 is the symbol rate. The angular baseband frequency can be derived using where can only be one of two values since m can only be either positive or negative unity, that is, The frequency separation between the sideband frequencies can be seen to be f s ∕2 = 100 Hz. The maximum phase difference between the start of one symbol and the next is ∼ ±90 ∘ or ± ∕2.
Therefore, a MSK transmission can be considered to be a special case of frequency shift keying where the message is transmitted by the pattern of switching between two specific sideband frequencies, ± s 4 , spaced s 2 apart. The two sideband frequencies are symmetric around the carrier frequency so that at baseband they have opposite polarities. A Gaussian filter is applied to smooth the symbol transitions to improve spectral efficiency. Specifically, equation (2) becomes where m g (t) is the Gaussian filtered m(t). This variant of MSK is known as GMSK.
The instantaneous frequency of the received signal, f i = i ∕2 , can be calculated from the received phase r = arg(y r ) = arctan( ℑ(y r ) ℜ(y r ) ), given i = d dt r (t). The symbol sequence m(t) can be inferred from the polarity of i as per equation (4), which is essentially a demodulation operation. Once the symbol sequence has been determined, it can be modulated using the GMSK scheme as per equation (5) to estimate the transmitted phase (t).
The phase of this estimated transmission for a 300-ms period (15:00:05.485-15:00:05.785 UTC) is shown in Figure 1, left panel, in blue and red solid lines. Its corresponding instantaneous frequency is shown in the right panel, also in blue and red solid lines. The blue and red colors represent the periods when a "+" symbol or a "−" symbol, respectively, is transmitted and the vertical grid lines delineate each symbol. During a + symbol, the phase increases by a rate determined by the positive sideband frequency. For DHO, f s = 200 Hz, so the positive sideband frequency is 50 Hz as per equation (4)-this is shown in the instantaneous frequency in blue. During a − symbol, the phase decreases at an equal rate but of opposite polarity to that of a + symbol as per the negative sideband frequency (−50 Hz, shown in the instantaneous frequency in red). The phase and instantaneous frequency of the received transmission (shown in green dashed lines on the same plots) are similar to those of the estimated transmission with small deviations which are more visible in instantaneous frequency (the phase of the estimated and received transmissions is shown with an arbitrary offset for clarity).
The received phase includes a time delay, p , introduced by the propagation from transmitter to receiver. Consider a period t 1 − p to t 1 when the transmission is transmitting a symbol m 1 than the transmission has a frequency 1 in this period. This signal propagates at a phase velocity which has a magnitude at some proportion of the speed of light; the exact value is dependent on parameters of the Earth-ionosphere  (Right). Phase disturbance during the positive (blue lines) and negative (red lines) symbols and the transitions between them (gray dashed lines) when the instantaneous frequency is not constant. The vertical grid lines separate each symbol of period 5 ms. The received amplitude and phase disturbance exhibit a frequency dependence (compare blue to red lines) that is more pronounced in the amplitude than in the phase. They also exhibit fast temporal variabilities <5 ms.
waveguide at VLF, so that the absolute delay introduced by the propagation is unknown. The phase of the received transmission at t 1 is as follows: This propagation delay can be considered to be composed of a constant part that results from the average geometry of the Earth-ionosphere waveguide and a time variable part produced by changes in the ionospheric conductivity profile.
In practice, the absolute start phase of the transmission is also unknown (equivalent to the integration constant of equation (5) and the transmitted phase is calculated from the received signal which contains the unknown propagation delay p . This estimated transmitted phase is thus e = (t 1 ) − 1 e − C, where e is approximately the propagation delay and C is the constant offset introduced by the unknown transmission start phase. The actual propagation delay can vary with both time and frequency, that is, p (t, ), whereas the delay in the estimated transmitted phase ( e ) is taken to be a single time delay value that can be interpreted as the constant part of the propagation delay. This value ( e ) is referred to thereafter as the delay estimate. The calculated phase disturbance is thus Hence, d is a relative measure with offset 1 e + C.
The received amplitude and phase disturbance can be calculated per sample of the recording, resulting in a very high time resolution (1μs) observation of ionospheric perturbations via the VLF transmission. This is shown in Figure 2 for the same 300-ms period as in Figure 1 and using the same format-blue and red solid colors represent the periods when a + symbol or a − symbol, respectively, is transmitted and the vertical grid lines delineate each symbol. During transitions between the two symbols, the instantaneous frequency is not constant-the samples around such transitions are plotted in gray dashed lines.
During periods when there is strong interference, certain symbols in m(t) from the demodulation can be missed. In the case of a single missing symbol, d will be offset by ≫ 2 after the missing symbol, which makes the detection of such errors relatively straightforward. For example, consider a symbol +, using equation (2), the transmission phase at the beginning (t = 0) and end of the symbol (t = T s ) is (0) = 0 and (T s ) = 2 . If the symbol is demodulated as −, (T s ) = − 2 . The calculated phase disturbance at the start of the symbol is d (0) = 0, and at the end of the symbol is d (T s ) = . Taking the Gaussian filter of GMSK into account, the offset will be slightly less than but ≫ 2 as stated. Such offsets are monitored in the analysis and removed by flipping the polarity of the demodulated symbol at the offset. It follows that disturbances that produce phase offsets of ≥ are indeterminate in this analysis. However, phase changes ≥ over time scales of one symbol period (5 ms for DHO) are very rare in the data set analyzed in this study.

10.1002/2017RS006456
In the case of interference lasting over consecutive symbols, the correct symbols cannot be determined. A step change of value can still be avoided by making an approximation of the symbols, even though d is ambiguous over the period of the interference. No such interference is present in the data we examine in this study.

Frequency Dependence and Ionospheric Variability
It is important to distinguish the features of the received amplitude and phase disturbance that are VLF radio propagation effects as opposed to local noise or transmitter effects to determine the ionospheric electrical properties over the propagation path. There are two key propagation effects on the signal that have been identified from the array results of : (1) there are "two distinct source locations in the wave number domain. The source location associated with a positive apparent frequency occurs at a smaller wave number such that it would appear at a larger elevation angle in the radio sky when compared to the source location associated with a negative apparent frequency," and (2) accounting for the distinct source locations, there are faster variabilities <5 ms in the wavenumber vectors.
As described in section 3 and shown in Figure 2, right panel, the changing apparent frequency of the received transmission corresponds with the frequency modulation of the transmission. Therefore, the wavenumber vector differences with apparent frequency observed by  appear to be a frequency dependence of the propagation. The propagation effect experienced by a signal is a function of the properties of the propagation path. Frequency-dependent propagation paths mean that different frequency components of the transmission experience different propagation effects. This should be observed on a single receiver recording as a modulation of the frequency spectrum of the transmission by the wave propagation frequency response. Using the equivalent single receiver data produced from the array data set, the received amplitude over a 300-ms period is seen to be different between the two sideband frequencies of the transmission as shown in Figure 2, left panel, where the blue lines (when the + symbol or positive sideband frequency is being transmitted) are consistently higher than the red lines (when the − symbol or negative sideband frequency is being transmitted).
Further, the spectrum of the received signal can be compared to that of the transmission (estimated from the received transmission as per section 3) by taking their transfer function to achieve better frequency resolution of the wave propagation frequency response. The frequency spectrum of the received transmission can be described as Y r ( ) = Y( )h( ) in the frequency domain, where Y( ) is the frequency spectrum of the transmission and h( ) is the frequency response of the wave propagation. In practice, h( ) also contains the relatively constant instrument response. Assuming the estimated transmission to be closely equivalent to the source transmission, then the wave propagation frequency response can be found by taking the transfer func- where Y e ( ) is the frequency spectrum of the estimated transmission and both Y r ( ) and Y e ( ) can be found by taking the Fourier transform of the received and estimated transmissions, respectively.
As the absolute transmission amplitude is unknown, a unity value is used for the estimated transmission amplitude, y e (t) = e i e (t) e i c t , that is, the transfer function calculation compares the received signal spectrum to the spectrum of an estimated transmission with unity amplitude. Therefore, the amplitude of the wave propagation frequency response can only be determined relatively with all frequency components scaled by the unknown transmission amplitude A.
The amplitude of the wave propagation response |h( )| calculated over the whole 6 s period of the data is shown in Figure 3, left panel. At the positive sideband frequency (denoted by the blue vertical line), the amplitude is higher than that at the negative sideband frequency (red vertical line)-consistent with the received amplitude in Figure 2, left panel. The sudden changes in the amplitude response with frequency (e.g., the prominent "spike" at ∼23.36 kHz) are likely to result from noise of unknown origin. However, the amplitude also exhibits a general nonlinear trend over the bandwidth of the transmission. The nonlinearity can be approximated to two discontinuous linear sections shown in the figure by the magenta dashed lines-the amplitude varies with a positive gradient up to around the center frequency where the trend turns into a lower negative gradient.
The phase of the wave propagation response is arg(h( )) = Φ r ( ) − Φ e ( ), where Φ r ( ) and Φ e ( ) are the received and estimated transmitted phases calculated for each frequency component, respectively. When the transmission is at one of the sideband frequencies, 1 , then Φ r ( 1 ) and Φ e ( 1 ) are equivalent to r and e from equation (7) and arg(h( 1 )) = 1 ( e − p ) + C; that is, only a relative measurement is possible due to the unknown transmission start phase and propagation delay as discussed in section 3.
Previous observations at the two sideband frequencies (e.g., Adams & Dowden, 1990) posit that changes in the phase difference between each sideband frequency of the transmission are due to changes in group delay. In this work, any group delay effect results from a difference between the delay estimate e and the actual propagation delay p . For example, if arg(h( )) is calculated for a length of transmission containing multiple frequencies (i.e., when there are transitions between symbols), then the phase of the propagation response would vary with frequency with a constant gradient equal to e − p . This phase-frequency gradient produced by group delay is used in many applications to estimate signal arrival times; for example, Dowden et al. (2002) use the group delay of sferics in a network of receivers to accurately compute lightning location. The value of e is a single value for all frequencies chosen to be as close to p as possible to minimize the phase-frequency gradient in the results.
The phase of the wave propagation response arg(h( )) calculated over the whole 6-s period of the data is shown in Figure 3, right panel. As for the amplitude response, sudden changes in the phase cannot be interpreted as physical as they are likely to be due to noise of unknown origin. However, the phase also exhibits a general nonlinear variation. The nonlinearity can be approximated to three discontinuous linear sections shown in the figure by the magenta dashed lines-the phase increases up to a maximum at ∼23.36 kHz then it decreases to a minimum at ∼24.44 kHz before increasing. This nonlinear trend suggests that the propagation delay is frequency dependent, that is, p ( ), even over the narrow 140-Hz bandwidth shown in the figure. Figure 4 shows the statistical mode of received amplitude (in the left panel) and phase disturbance (in the right panel) every second over the entire 6 s period of the data. The distributions are broken down in time according to which sideband frequency is being transmitted and their statistical modes displayed in blue and red solid lines as per Figure 2. The received amplitude and phase disturbance statistical modes for each sideband frequency are shown to be different in Figure 4, again suggesting frequency dependence of the wave propagation.
This frequency dependence can originate from the transmitter similar to the stability issues discussed in section 2. For example, frequency dependence can arise from the transmitter antenna being a more effective radiator at one sideband frequency than the other. Such imperfections in the transmission cannot be observed by existing methods unlike the stability issues discussed in section 2. In such scenarios where the imperfection is systematic, the frequency dependence seen should be constant, but the statistical modes for each sideband frequency vary independently in Figure 4 over the 6-s period of the data. Moreover, the array analysis carried out previously on the same data by  establishes that the wave propagation is dependent on the sideband frequency of the transmission. Therefore, it is likely that at least part of the observed frequency dependence is a result of wave propagation. The 1-s time period for the distribution calculation of Figure 4 is over several symbols (≫5 ms) so that both sideband frequencies are transmitted in the period and enables the distribution for both sideband frequencies to be calculated for the same second. The relatively long 1-s period is chosen for convenience, and shorter processing periods are possible. The samples around symbol transitions are removed from the analysis as the instantaneous transmission frequency during transitions is not constant (these are the samples shown in gray in Figure 2). The statistical mode of the distribution is used instead of the mean as it is considered more robust to outlier values, for example, lightning interference. Although weak lightning signals from far away storms can be picked up by the array around the time of recording (Füllekrug et al., 2016), the distributions of the received amplitude and phase disturbance are seen to be roughly normal in shape. An example of the distribution is shown in Figure 5 for 1 s of data from 15:00:04 UTC. Since the distribution is roughly normal, the standard deviation can be computed from the probability density at the mode of the distribution, a, so The standard deviation of the received amplitude and phase disturbance distribution is a measure of their variability in the processing period. The shaded areas in Figure 4 represent the standard deviation of each distribution around the statistical mode calculated for each second. The phase disturbance variability in the recording ranges from 0.43 ∘ to 0.80 ∘ , equivalent to a delay of 51-95 ns. Since the clock accuracy of the instrument is within 20 ns, ∼ 70% of the variability cannot be attributed to instrumental noise. This is consistent with the array analysis results of  which estimates a similar proportion of the wavenumber vector variability calculated from the received phase originates from wave propagation. Further, the recording is averaged over several receivers which reduces the contribution from uncorrelated sources (e.g., the instrumental clock accuracy) so that the phase disturbance variability that cannot be attributed to instrumental noise is actually > 70%. Given that wave propagation can produce phase disturbance variability of such magnitude, the effects of wave propagation on VLF transmission phase will still be observable with a true single receiver recording under similar conditions. The received amplitude variability (Figure 4, left panel shaded area) in the recording ranges from 11.7 to 15.7 μVm −1 which is representative of the instrumental accuracy so is likely to be mostly instrumental. Transmitter imperfections can also produce fast variabilities. For example, the unknown clock accuracy of the transmitter can produce fast variabilities in the transmitter phase. However, transmitter imperfections will not produce wavenumber scatter so the origin of the fast variability observed here is likely to be wave propagation.

Diurnal Frequency Dependence and Fast Variability
The sidebands of the DHO transmission are separated only by 100 Hz. Over such a narrow frequency separation, it is unknown how the frequency dependence will change with the lower ionosphere. The D region ionosphere exhibits a predictable, repeated change with the movement of the terminator in the diurnal cycle. Therefore, a separate observation is presented here of the same DHO transmission when the night to day terminator was traveling along the propagation path from 27 August 2016 03:30-07:00 UTC.
During this time, the receiver was operating in a monitoring mode when the data are split into roughly hourly files by automatically restarting the recording process every hour. The restart normally produces a small gap in the data of <30 s, for example, 04:28:18-04:28:42 UTC and 06:28:18-06:28:42 UTC. One of the hourly restarts is during the longer ∼12-min gap between 05:20:05 and 05:32:37 UTC-the extended duration is a result of the automatic processing algorithm disregarding samples with amplitude below a set threshold as unreliable.
On the date of the recording, sunrise time at the DHO transmitter site is 04:31 and 05:14 UTC at the receiver in Bath, UK. Between those times, the terminator was moving from the DHO transmitter site to the receiver. The receiver was located in the University of Bath, UK, so that the propagation path from DHO to the receiver is almost the same as to the array. Figure 6, top row, shows the statistical modes of the received amplitudes (left) and d (right) at each sideband frequency. The received amplitudes at each sideband are seen to vary by up to 0.5 mVm −1 during nighttime conditions before the terminator reaches the DHO transmitter site at 04:31 UTC. After this time, there is a steep monotonic decrease starting ∼04:30 UTC that is attributable to the crossing of the terminator. The received amplitudes then stabilizes with a small steady increasing trend ∼05:30 UTC. It is evident over these few hours that there is a general correlation between the received amplitudes at each sideband, but there is a varying difference between them that can be up to 0.6 mVm −1 or ∼20% of the average value of the received amplitudes. The inset figure in Figure 6 top left panel shows the differences between the sidebands which are greater and vary more during nighttime conditions.
Interestingly, the start of the d response to the crossing of the terminator appears to lead that of the received amplitude by ∼30 min, starting around ∼04:00 UTC. The phase disturbance d increases steadily (faster in the first 15 min) until just past 05:00 UTC, again, before the stabilizing of the received amplitude by ∼30 min. Larger variations are also seen in d during nighttime compared to daytime conditions. The general correlation between the sideband frequencies is also seen in d with the propagation delay (approximately proportional to the difference between the d of each sideband frequency) varying more during nighttime. . The ∼10-min gap around 05:30 UTC is due to the recording process being restarted (see section 5). There is a general correlation between the received amplitudes at each sideband, but there is a varying difference between them. The differences between the sidebands are greater and vary more during nighttime conditions.
The bottom panels of Figure 6 show the standard deviations of the received amplitudes (left) and d (right) which represent the fast variability of the signals. During nighttime conditions, the received amplitude variability is larger (∼40-60 μVm −1 ) but during daytime conditions the variability is smaller so that the effect is likely to be instrumental in nature (<35 μVm −1 , the instrumental noise floor). It is possible that this change in variability is due to the attenuation of variability from the transmitter as the variability change appears to follow the change in received amplitude. However, the variability on both sideband frequencies is similar while there is a difference between the received amplitudes. The variability of d does drop below 50 ns, very similar to the variability seen in the array data set (day time). There is a period ∼04:45-05:30 UTC where the variability shows an increase up to 250 ns before it decreases again. Again, the variability is seen to be similar between the sideband frequencies.

Discussion
The VLF transmission propagates in the Earth-ionosphere waveguide via different wave propagation modes which superpose so that the radio wave amplitude and phase varies spatially when observed on the ground,

Radio Science
10.1002/2017RS006456 that is, they form a modal interference pattern. When the waveguide changes, the radio propagation and hence the modal interference pattern are altered correspondingly. Since the Earth's surface profile and electrical properties can be taken as constant in time, changes in the modal interference pattern are attributed to ionospheric perturbations. For example, measurement of nulls in the interference pattern of VLF transmissions has been used to infer a long-term lowering in the altitude of the D region ionosphere linked to an increase in greenhouse gases . A radio receiver samples the modal interference pattern at a single point; however, the interference pattern can be influenced by ionospheric perturbations anywhere, with the largest influence coming from perturbations along and near the transmitter to receiver great circle path. Therefore, the remote sensing method using VLF wave propagation is effectively an integrated measurement along the transmitter to receiver great circle path. For a sideband frequency of the transmission, that is, a single frequency component, the propagation variability inferred from the microsecond time resolution received transmission amplitude and phase disturbance may be due to the spatial inhomogeneity of the ionosphere over the path. For example, if the ionospheric variation in time at one point is different from that at another point, the integrated effects over the path are variabilities in the received amplitude and phase disturbance that can be faster than the ionospheric variations at each isolated point. The propagation variability can also be an effect of a fast-varying ionosphere as suggested by . Further observations are needed to study the physical meaning of the propagation variabilities. In particular, it is unclear how these variabilities would change with changes in ionospheric forcing over a longer period of time.
From the measured radio signal, the profile of the lower ionospheric electrical properties can be found by solving the inverse radio propagation problem. The radio propagation modeling is usually performed with different electrical property profiles to find a solution that best fits the measured data. To make the process tractable, the lower ionosphere is commonly simplified to a laterally homogeneous, horizontally stratified structure where electron density and collision frequency increase exponentially with altitude z. With this simplification, the ionospheric profile can be described with only two Wait parameters of reference height h ′ and sharpness (Wait & Spies, 1964). This approach has been used to estimate stable D region ionospheric conditions using existing lower time resolution, single frequency measurements of VLF transmissions (e.g., Thomson, 2010;Thomson et al., 2017). Consider a Wait ionosphere in the absence of a geomagnetic field. The ionospheric plasma has a complex permittivity that is a function of its electron density N e (z), collision frequency (z), and the angular wave frequency , = 1 − N e (z)e 2 m e 0 2 (1 − i (z)∕ ) (adapted from Budden, 1988, chapter 3), where e and m e are the electric charge and mass of an electron and 0 is the free space permittivity. It follows that the reflection coefficient for each propagation mode and hence the modal interference pattern are dependent on frequency.
Theory and empirical evidence show that over large frequency ranges, the propagation frequency dependence is nonlinear (e.g., Chapman & Macario, 1956;Dowden et al., 2002). However, the DHO transmission has a very narrow bandwidth-the separation between the two sideband frequencies is only 100 Hz. This may be representative of a more complex ionospheric profile; for example, rocket measurements reveal a vertical ionospheric profile (Friedrich & Torkar, 2001, and references therein) with much more detail than a Wait ionospheric profile, and the lateral homogeneity assumption may not always be valid. The ionospheric profile changes depending on the various different forcing effects, for example, the modulation of solar radiation by time producing diurnal changes in the ionospheric profile. This has a significant effect on the frequency dependence of propagation as shown using the single receiver data set presented in section 5. Ionospheric profiles different to that during the recordings presented here will likely produce different wave propagation responses.
The multifrequency measurement motivates more detailed propagation modeling in the determination of ionospheric profiles. This can potentially improve resolution of finer features or provide higher certainty or accuracy in the inversion. Such improvements are particularly useful in observing complex perturbations such as those produced by lightning in thunderstorms. Another approach to obtain multifrequency measurements is to use averaged broadband sferics (Cummer et al., 1998;Lay & Shao, 2011), which allow measurements over much larger frequency bandwidths than VLF transmissions. However, the source signal and location of VLF transmissions are more certain than those of sferics. Moreover, the multifrequency transmission disturbance analysis can achieve better time resolution (1 s) and is able to provide continuous phase measurements when compared to the intermittent lightning.

Summary
Observing disturbances on VLF transmissions is a commonly used method to probe conductivity perturbations of the lower ionosphere. Most existing analysis provides single frequency, ∼20-ms time resolution observations of the disturbances. This work presents a method to make VLF transmission amplitude and phase disturbance measurements per sample from a 1-μs time resolution record.
The fast measurement is able to resolve the frequency modulation of the VLF transmission and hence observe frequency-dependent propagation effects. The multifrequency measurement can potentially improve reconstructions of the lower ionospheric profile, especially in the study of complex perturbations such as those caused by lightning. Fast variability of the propagation at each sideband frequency of the transmission (<5 ms) can also be observed with the method which may be due to ionospheric variability over the propagation path or fast processes in the ionosphere.