Interactive comment on “ Semi-annual oscillation ( SAO ) of the nighttime ionospheric D-region as detected through ground-based VLF receivers ” by I

General remarks This paper describes VLF observations in two mid latitude stations of the semi-annual oscillation in the D region of the lower ionosphere and suggests that NOx molecule transport from the lower atmosphere to the night time D region could be at the origin of this oscillation. Such VLF observations are new and could be worth for publication in ACP journal. However, some parts of the paper could be improved taking into account the following remarks:


Introduction
Earth's middle and upper atmosphere exhibit several dominant large-scale oscillations in many measured parameters.These oscillations can be found at all latitudes, from the Equator to the mid-and high latitudes.One of these oscillations is the semi-annual oscillation (SAO).Among different parameters, the SAO can be detected in neutral atmospheric measures, e.g., the wind regime in the mesosphere-lower thermosphere (MLT) (e.g., Groves, 1972;Gregory and Manson, 1975;Lysenko et al., 1994); MLT temperatures (e.g., Groves, 1972;Takahashi et al., 1995;Taylor et al., 2005;Huang et al., 2006;Shepherd et al., 2006); and concentrations of atmospheric species, such as atomic oxygen at 80-115 km altitudes (e.g, Russell et al., 2004) and excited hydroxyl (OH * ) molecules around 87 km (e.g., Takahashi et al., 1995;Marsh et al., 2006;Shepherd et al., 2006;Gao et al., 2010).In addition, the charged part of the atmosphere (i.e., the ionosphere), experiences the SAO, which was observed and derived in and from measurements of several parameters, e.g., the ionospheric lower transition height at ∼ 180-260 km (a level where atomic and molecular ion concentrations become equal) (e.g., Lei et al., 2004); the electron and plasma density within the daytime D, E, and F regions of the ionosphere (e.g., Lauter and Nitzsche, 1967;Bremer and Singer, 1977;Forbes et al., 2000;Peters and Entzian, 2015); and also ionospheric total electron content (TEC) (e.g., Zhao et al., 2007;Opio et al., 2015).
Measurements of the D region of the ionosphere, which lies at altitudes of ∼ 80-95 km during nighttime and expands downwards to lower altitudes (∼ 60 km) during daytime (mainly due to direct solar extreme ultraviolet (EUV) Published by Copernicus Publications on behalf of the European Geosciences Union.and X-ray radiation) (Brasseur and Solomon, 2005;Inan et al., 2010), are usually made using remote-sensing techniques, because these altitudes are too high for weather balloons and too low for in situ measurements by satellites.In addition, remote-sensing techniques usually do not possess the limited spatial and temporal coverage of rocket lofted experiments (Rodger and McCormick, 2006).One of these remote-sensing techniques involves the use of very low frequency (VLF) radio waves, spanning a frequency range of 3-30 kHz.These waves, which are generated both by natural and man-made sources, propagate thousands of kilometers within the Earth-ionosphere waveguide, reflected off the Earth's surface and inside the ionosphere's D region, while experiencing a very weak attenuation of ∼ 2 dB Mm −1 (Barr et al., 2000;Wait, 1957).Due to the significant difference in the D region's characteristics (electron and ion densities) between day and night (Hargreaves, 1995), the region's conductivity changes dramatically, causing the D region reflection height of VLF signals to change from as low as ∼ 60 km during daytime to ∼ 85 km during nighttime (Hargreaves, 1995;Inan et al., 2010).Thus, as received VLF signals inherently contain information of the ionosphere and its variability within the reflection region (Inan et al., 2010;Rodger et al., 2012), these signals probe different altitudes during day and night.
Because the D region's formation and chemistry are tightly bound to the neutral MLT (Brasseur and Solomon, 2005), it is believed that the D region is affected by the same forcings, experiencing similar oscillations (e.g., Schmitter, 2011;Silber et al., 2013;Marshall and Snively, 2014).As far as the current authors know, there are no previous works showing the SAO dominating the natural long-term oscillations in the nighttime D region apart from the work of Toledo-Redondo et al. (2012), who presented a SAO indication within the equatorial latitudes by using space-based ELF-VLF data from the Detection of Electro-Magnetic Emissions Transmitted from Earthquake Regions (DEMETER) micro-satellite.In this paper, we present evidence of a strong SAO, detected in the low-and midlatitude nighttime D region, through ground-based VLF measurements in both hemispheres.

Instrumentation and methodology
During this study, we used ground-based VLF narrowband (NB) signals, which are generated by VLF transmitters.These man-made transmitters are used nowadays worldwide primarily for communication with military submarines (Clilverd et al., 2009;Rodger and McCormick, 2006).However, they are extremely well suited to long-range remote sensing of the D region, because of their high radiated power, their nearly continuous operation, and their fixed location and frequency band (e.g., Clilverd et al., 2009;Inan et al., 2010).As the VLF signals travel from the transmitters to the receivers along a great circle path (GCP), a time series of their recorded amplitude and phase give an indication on the changes of the D region along the GCP.
The NB signals were recorded at two VLF receiving stations.The first station is located at the Emilio Segre' Observatory of the Israeli Cosmic Ray and Space Weather Center, at Mt. Hermon (MH), in the north part of Israel (33.18 • N, 35.47 • E).This VLF receiving system is part of the Atmospheric Weather Electromagnetic System for Observation, Modeling, and Education (AWESOME) network (Cohen et al., 2010).The antenna is built from two orthogonal triangular loop antennas that measure the two horizontal components of the VLF magnetic field.Each loop has a baseline of 2.6 m and height of 1.3 m, giving an area of approximately 1.69 m 2 for each loop, and has a total number of 12 turns.The loop antenna impedance is 0.85 mH and 1 .
The second VLF receiving station, which is located in Dunedin (DN), New Zealand (45.8 • S, 170.5 • E), is part of the Antarctic-Arctic Radiation-belt (Dynamic) Deposition -VLF Atmospheric Research Konsortium (AARDDVARK) network (Clilverd et al., 2009) and is operated by the University of Otago.Its antenna measures the normal component of the VLF electric field, thus making both of the stations' measurements equivalent.The VLF signal is recorded by an "OmniPAL" narrowband VLF receiver (Dowden et al., 1998).
Global lightning activity responds to the Earth's surface air temperature on both the semiannual and annual timescales, within the tropical belt and the mid-to high latitudes, respectively (Williams, 1994).As most of the electromagnetic (EM) energy generated by lightning discharges (termed "sferics") are radiated within the ELF (extremely low frequency) and VLF bands (peaking between 5-10 kHz) (Cummer, 1997;Rakov and Uman, 2003), and due to the weak attenuation of VLF signals (as mentioned above), sferics originating in the tropics and at mid-to high latitudes can be easily detected at MH or DN.Although lightning pulses are very powerful, their duration is very short (up to the order of 10 −4 s) (Rakov and Uman, 2003).Nevertheless, more than 1000 active thunderstorm are present on average, at any given moment (Mezuman et al., 2014).Altogether, global lightning activity generates significant VLF EM fields that can produce interference with the NB measurements, inducing the lightning activity natural oscillations within the NB data.Therefore, in order to investigate long-term oscillations in NB amplitudes, the background VLF noise, which is created mainly by lightning discharges (Barr et al., 2000), should be removed.
The MH receiving system records the NB signals continuously.However, as discussed above, these signals are potentially biased by the background VLF sferics.Thus, we decided to use the broadband (BB) data, which consist of signals of the whole VLF band and are recorded in a synoptic mode; i.e., each data file is a recording for 1 min every 15 min.For the extraction of a NB signal for a certain VLF transmitter frequency from the BB data, every minute of data was filtered using a Parks-McClellan finite impulse response (FIR) band-pass filter, with a passband width of 300 Hz, where the VLF transmitter's central broadcast frequency lies at the middle of the passband.In order to represent the noise at the transmitter's frequency, the original BB data were filtered again, one time when the middle of the passband was 300 Hz above the VLF transmitter's central broadcast frequency, and another time when the middle of the passband was 300 Hz below the VLF transmitter's central broadcast frequency.The average of the two filtered noise time series was then subtracted from the filtered NB signal time series, assuming that the average noise represents the background noise in the frequency band of the transmitter's broadcast.Thus, we received the NB data with a strongly reduced bias of the background lightning noise.We then average the resultant NB signal every 10 s, thus receiving up to 24 values per hour made up of 10 s averages across the 1-in-15 min observation cadence.
Unlike the MH VLF receiver, the DN VLF receiving station operates the OmniPAL receiver, which also records the NB signals continuously.During the data acquisition process, this receiver's software uses a noise-clipping algorithm that reduces the effect of lightning impulses over the NB data (see Dowden et al., 1998).Thus, there are no additional procedures needed in order to investigate the NB received signals, which are being recorded at 1 Hz.
In this paper, we limit ourselves to observations from two VLF transmitters, with call signs NSY (38.00 • N, 13.50 • E, broadcasting at 45.9 kHz) and NWC (21.82 • S, 114.17 • E, broadcasting at 19.8 kHz).These transmitters were chosen as their operating power was constant along the data sets, their signal-to-noise ratio was relatively high, and their data were fairly complete.The NSY transmitter's signal was recorded at MH during the years 2009-2012, while the NWC signal was recorded at DN during the years 2005-2010 (Fig. 1).We extracted from each time series of transmitter-receiver amplitudes two data sets of the average amplitude during the midday and midnight hours, when the solar elevation angle at the middle of the GCP during equinox was at its maximum and minimum, respectively.It should be mentioned that atmospheric tides, which have time periods that are subharmonics of a solar day (Oberheide et al., 2003) can possibly affect our data.Therefore, a 24 h data average should have been used in order to remove the tidal effect.However, as the ionosphere (and, as a result, the received NB signals) changes significantly between day and night, this was not possible.As a result, we used 1 h mean amplitudes and not the daytime or nighttime amplitudes spanning those entire time periods, in order to reduce tidal averaging effects (assuming they exist) to the minimum possible.In other words, because daytime and nighttime lengths change throughout the year, in the case of complete daytime and nighttime averaging, a different window size of the tidal oscillation would have been averaged every day.

Results
The midday and midnight 1 h mean 30-day running-average time series for the MH-NSY and DN-NWC GCP deviation from the mean amplitude (of the entire time series) are shown in Fig. 2 (black solid curves).As can be seen, all the time series exhibit a strong oscillatory behavior, with higher amplitudes in the midnight data than in the midday data.The midday data in both GCPs show a dominant oscillation with longer time periods than the midnight data in both GCPs.Examination of the apparent time periods of the large oscillations shows that they appear to correspond to the annual oscillation (AO) and SAO.Therefore, we wanted to fit these harmonics to the data and examine the level of agreement with the different time series.
However, in addition to these oscillations, a trend can be seen in the time series: negative in both DN-NWC time series and positive in the MH-NSY midnight data, but hard to determine in the MH-NSY midday data, as the trend seems to shift from negative to positive around mid-2010.Examination of the data set time span shows that the DN-NWC data were acquired during a period when solar activity was dropping towards a minimum, while the MH-NSY data were acquired when the Sun started to become active again, as part of its 11-year cycle.Thomson and Clilverd (2000) showed a positive correlation between VLF amplitudes and solar activity.Therefore, as both of our data sets show a general positive correlation with solar activity, we may conclude that the trend is a result of solar activity.Nevertheless, because we did not have enough data to cover a full 11-year solar cycle, it was problematic to fit the data with an 11-year harmonic.Thus we decided that a linear fit to the data would be best.Therefore, the time series were fitted with curves, described by following equation: where A fit is the fitted curve; t represents the time steps (in days); A 0 is the mean amplitude (which is equal to 0 in this case); S is the linear fit coefficient; A SAO and A AO are the fitted SAO and AO amplitudes, respectively; and t SAO (t AO ) represents the SAO (AO) maximum time of year.Both the linear and the harmonic fits were made using a least-squares method over all of the data points.The fitted curves are shown in Fig. 2 (dashed red).As can be seen, these simple curves follow the VLF amplitude patterns fairly well.Pearson's correlation coefficients between the time series and the fitted curves were calculated and are shown at the bottom right of each panel.As the correlation coefficients span from values of 0.53 up to 0.84 (all statistically significant), we can deduce that the simple curve may explain from 28 up to ∼ 70 % of the midday and midnight long-term variability.
The simple model's parameters described in Eq. ( 1) can be investigated as well.Comparison of the two oscillation am-plitudes (A SAO and A AO ) shows that during midday A SAO is 3 times weaker than A AO , but during midnight it is stronger than A AO by up to ∼ 60 %.Moreover, A SAO appears to have a very strong peak-to-peak amplitude of 3.3 dB in MH-NSY and 4.2 dB in DN-NWC.By running the Long Wave Propagation Capability (LWPC) model (Ferguson, 1998), we found that these strong amplitudes are equivalent to h (ionospheric base) change of no more than 1.8 and 1.3 km in DN-NWC and MH-NSY, respectively, or 0.13 and 0.15 km −1 in β (electron density profile sharpness), respectively.When the standard D region electron number density profile is used (Wait and Spies, 1964), where h is the altitude and n e is the electron number density, it can be shown that the peak-to-peak amplitude change caused by the SAO alone is equivalent to more than doubling of the electron number density at 85 km.The actual solution for the ionospheric profile is probably a combination of changes in h and β, but as we do not have reliable phase measurements the actual solution cannot be calculated at the moment.
When the t SAO values for the midnight data are examined, it is found that the SAO maxima occur up to a month prior to Earth's winter and summer solstices (not shown).This is quite surprising, as we would normally expect the maxima of a SAO-driven forcing to occur around equinox (e.g.Opio et al., 2015;Taylor et al., 2005;Williams, 1994).The fitted curve for MH-NSY GCP data shows the SAO maxima dur-ing mid-November and mid-May, while the oscillation maximizes at the beginning of December and June in the DN-NWC GCP data.Thus, a 16-day phase difference exists between the two data sets.
In order to confirm our findings of the apparent dominating SAO in the nighttime NB measurements (and possibly within the nighttime D region), spectral analysis was performed.Because the data were unevenly sampled (due to transmitter off-times, receiver malfunctions, etc.), it was not possible to use fast Fourier transforms (FFTs), which demand constant time steps between samples.In addition, the FFT-calculated frequencies are directly determined by the data set length and sampling rate; hence there is no option to choose which exact frequencies to inspect.Moreover, a transformation of the FFT frequencies into the oscillations time periods (via T = 1/f ) results in a very high time period resolution at very low values (i.e., high frequencies) and very poor resolution at high values (i.e., low frequencies, more likely to represent the SAO long-term oscillations).Therefore, we have analyzed the midnight 1 h mean data using the Lomb-Scargle (LS) periodogram (Lomb, 1976;Scargle, 1982), which results in spectral power of the data at user-determined frequencies (and, hence, time periods) and allows the spectral analysis of unevenly sampled data (Press and Rybicki, 1989).
Figure 3 shows the LS periodogram of the midnight (unsmoothed) MH-NSY (top panel) and DN-NWC (bottom panel) 1 h mean VLF amplitude anomalies, with arbitrary power units (as a result of the LS periodogram procedure).The dashed red line denotes 95 % confidence, which was calculated using the quantile function (Wilks, 2006).The inspected time periods range from 2 to 730 days, thus spanning from as short as the data set's Nyquist frequency up to 2 years.Examination of the MH-NSY periodogram confirms that the SAO at ∼ 180 days is by far the most dominant and significant oscillation within the data.The second peak of the periodogram is of 343 days.We ascribe this peak to the AO, while we assume that the difference in periodicity from 365 days is attributable to our data set covering only 4 years and thus containing only 4 AOs, which are not fully sampled, due to the data gaps (data coverage of ∼ 57 % along the time span).Secondary peaks, which do not pass the 95 % significance threshold, are seen at time periods of 47, 96, 137, and 212 days.Some of these oscillations might be higher harmonics of the SAO, but it is not possible to explain them at the moment, leaving this topic for future studies.
The SAO at ∼ 180 days appears to be even more pronounced and significant in the DN-NWC periodogram and is the dominant oscillation within the midnight data.As can be seen in the lower panel of Fig. 3, the second-highest peak is of a time period of 241 days (∼ 8 months), an oscillation which is quite unexpected, but does not appear in the MH-NSY data.The probable signature of the AO seen in this periodogram is also statistically significant, peaking at 366 days.Here, the secondary peaks which do not pass the 95 % significance threshold are located at time periods of 152 and toward 730 days (∼ 2 years); the latter might be hinting at a very weak quasi-biennial oscillation (QBO) affect.
As the nighttime measurements of VLF amplitude are highly variable, robustness tests were made to examine our nighttime data and methodology, by removing 20 % of the raw data's measured points.Later, we added Gaussian noise into the raw data's measured points, with a standard deviation equal to that of the raw data.Each of these tests was repeated for 100 000 iterations.The results showed that in both data sets the SAO and the solar cycle trend passed these robustness tests in 100 % of the runs, and they therefore strengthen our findings.The AO passed the robustness tests in 100 % of the runs for the DN-NWC data set but did not prove to be statistically significant for the MH-NSY data set, as only 91 % of the data-removal runs and 58 % of the noise-addition runs kept this oscillation statistically significant.

Discussion
In this study, we analyzed several years of VLF NB data received in both hemispheres, during midday and midnight hours.The analysis shows that the AO dominates midday VLF amplitudes, and the SAO is the strongest oscillation during the 1 h long period centered around midnight.Both the SAO and the essential differences between daytime and nighttime dominating oscillations should be explained.The probable reasons for both of these observations are of chemical and dynamical origin, which take place in the transport of species and tidal forcing.These sources shall now be discussed.

D region ions and dynamical transport of neutral species
When analyzing Earth's ionosphere in general (and the D region in particular), we can assume electrical neutrality (Kelley, 2009): where n + is the positive ion number density, n − is the negative ion number density, and n e is the electron number density.VLF radio signals interact and are reflected off the D region mainly by electrons rather than ions, as a result of their much lower mass (Inan and Inan, 2000).Using Eq. ( 3), we can assume that the electron number density is determined primarily by the dominant ion densities.Therefore, we can investigate the dominant D region's ions, their production, and distribution.
As mentioned in the Introduction section, VLF reflection height ranges from ∼ 60-70 to ∼ 85 km during daytime and nighttime, respectively (Inan et al., 2010;Rodger and Mc-Cormick, 2006).Although both altitudes are part of the D region, their chemical composition and dynamical processes are different and both are very complicated.The lower altitudes of the D region (below ∼ 80 km) are dominated by positive cluster ions (mostly of the type H(H 2 O) + n ), due to a relatively high neutral density in general and water molecules in particular, which enable an effective three-body reaction (with O + 2 and NO + ions), thus creating this type of ions (Glukhov et al., 1992;Goldberg and Aikin, 1971;Mitra, 1981;Narcisi and Bailey, 1965).Cluster ions have a rapid recombination rate, and many of the chemical reactions involving them are strongly temperature sensitive (Kelley, 2009;Pavlov, 2014).Therefore, the ion composition of this region should be quite variable with season and latitude, and sporadic changes associated with local temperature variations should be observable (Brasseur and Solomon, 2005).
The higher altitudes of the D region (above ∼ 80 km) are dominated by NO + ions, mainly as a result of strong ionization by the solar Lyman-α line (121.6 nm).O + 2 ions are also abundant in this region and are created mostly by solar radiation in the 102.7-111.8nm wavelengths.The threebody reactions that form cluster ions are less frequent in this region, due to the lower neutral density and the lack of water vapor, making cluster ions much less abundant (Mitra, 1981;Narcisi and Bailey, 1965).The Lyman-α radiation scattered by the hydrogen geocorona at the uppermost part of the atmosphere is still a major ionization source during nighttime, though 2-3 orders of magnitude smaller than during daytime (Brasseur and Solomon, 2005).Although NO + and O + 2 ion recombination rates are orders of magnitude slower than cluster ions (Kelley, 2009;Pavlov, 2014), the lifetime of ions in the D region is short compared to the transport timescale; hence the ion concentrations are determined by a photochemical equilibrium between production and loss processes (Brasseur and Solomon, 2005).Nevertheless, in addition to the strong dependence on solar radiation variability -i.e., solar activity, solar zenith angle, etc. -the production rates of NO + and O + 2 ions are proportional to NO and O 2 neutral molecules, respectively (Pavlov, 2014).In the neutral atmosphere, as we rise from ∼ 60 to ∼ 90 km, the chemical lifetime of NO x molecules increases, in comparison with the typical constant for vertical exchanges K zz .Thus, as we reach higher altitudes within the D region, which coincides with these altitude levels, the role of local dynamics (vertical exchange) becomes more significant (Solomon et al., 1982a).In addition, the amplitudes of gravity and planetary waves penetrating into the MLT (as a function of season; Lindzen, 1981) grow exponentially with altitude, as the ambient density drops (Ern et al., 2015;Smith, 2012).Therefore, downward transport of NO molecules, which are created mainly in the lower thermosphere (Solomon et al., 1982a, b), can increase the D region NO + concentrations, depending on the vertical wind patterns (e.g., Clilverd et al., 2006).
Altogether, changes in neutral NO occurring by dynamical forcing will affect NO + ion concentrations within the D region, mainly at higher altitudes.At lower altitudes of the D region, where daytime VLF signals are reflected, we would expect chemical processes to dominate over NO dynamics, and therefore a very strong signature of solar insolation changes (which are seen mainly in the AO), together with relatively weak perturbations caused by other forcings and temperature changes, is observed (e.g., Schmitter, 2011;Silber et al., 2013).At higher altitudes within the D region, where nighttime VLF signals are reflected, dynamical processes are much more pronounced; thus oscillations, such as the SAO, which are driven by dynamical transport of important species as well as dynamical forcing (e.g., gravity and planetary waves) are much stronger and thus more easily detected.
There are additional factors that can be expected to increase the detected SAO amplitudes in NB measurements.The first factor is the SAO of atomic oxygen in the MLT (e.g., Russell et al., 2004).These atoms, which can also be transported from higher regions of the atmosphere, are important for molecular ion chemistry, through numerous chemical reactions (Pavlov, 2014).Thus, they might increase the SAO amplitude in the D region, thereby enhancing the measured oscillation in the received VLF amplitudes.An additional affect comes from the MLT temperatures that, as mentioned in the Introduction section, are also experiencing a SAO.MLT oscillatory temperature changes can influence VLF received signals, by modifying the electron collision frequency, recombination, and attachment rates (see the discussion in Sil-

The VLF SAO phase and its comparison with satellite data
All of the above-mentioned phenomena can explain the measured dominating SAO in the received VLF amplitude differences at midnight.However, the phase of the measured SAO may be determined by a number of effects.The observed difference between the SAO phase in DN-NWC and MH-NSY might originate in the phase differences between the phenomena as a function of latitude, or be due to the changes of the phenomena's phases as a function of latitude and altitude (e.g., Groves, 1972;Lysenko et al., 1994;Russell et al., 2004;Taylor et al., 2005).In order to test our hypothesis, we decided to compare the phase of the SAO in the VLF measurements to the SAO phase detected in peak emission values of the OH * 2.0 µm emission band, using the SABER (Sounding of the Atmosphere using Broadband Emission Radiometry) instrument on board the TIMED (Thermosphere Ionosphere Mesosphere Energetics Dynamics) satellite (Mlynczak, 1997;Russell III et al., 1999).The OH * airglow layer peaks around 87 km (Baker and Stair Jr., 1988), in the altitude vicinity of the VLF nighttime reflection height.We used the 2.0 µm data as they are a direct measure of the chemical reaction which creates the OH * (Mlynczak, 1999;Mlynczak et al., 2013) and are tightly linked to atomic oxygen abundance as well as MLT dynamics (Gao et al., 2010;Marsh et al., 2006).A 60-day running-mean of the data was created, in order to obtain good local time coverage (Marsh et al., 2006).Only the SABER nighttime data (in a 1 • zonal-mean resolution) were examined, as the daytime and nighttime OH * airglow layers behave differently (Gao et al., 2015;Marsh et al., 2006;Smith, 2004).An examination of the 11-year (2002-2012) average peak emission is presented in Fig. 4.An apparent SAO seems to dominate the OH * emission pattern, especially at the equatorial and midlatitudes, and a strong phase propagation of the oscillation towards the poles can be observed.It can also be seen that the OH * emission does not show a constant SAO phase as a function of latitude (see also Gao et al., 2010), similar to the analyzed VLF data.We examined the OH * data (1 • width) at the latitudes that match the middle of the two transmitterreceiver GCPs, during the same time period as the VLF data (2005-2010 for DN-NWC, and 2009-2012 for MH-NSY), by using Eq. ( 1), and finding the SAO maximum phase.The results show that the SAO in the OH * leads by 6 days the MH-NSY VLF data, and by 17 days the DN-NWC VLF data.The differences in the oscillation phase might be due to the long GCPs, which (unlike the data used from SABER measurements) are over 1 • of latitude width.In addition, the OH * emission was zonally averaged, and therefore it was less sensitive to local perturbations than the VLF data.Nevertheless, as the maxima are all in the same 30-day window, we note that both parameters are affected by similar dynamics and forcing.Moreover, according to our knowledge, the exact reason for this phase lag as a function of latitude has not been investigated previously.As the dynamics of the MLT will affect multiple wind and wave fields as well as species concentrations, this topic should be investigated in future studies.

Tidal effects
Finally, we should point out one caveat in this study.Marsh et al. (2006) concluded that the observed SAO in OH * emission is a result of the seasonal change in atmospheric tidal amplitudes and is not caused by changes in diffusive transport, as was previously proposed.This might also be a significant effect in the VLF NB measurements, as tidal amplitudes also experience a SAO which maximizes at equinox and is driven by the strong shears in zonal mean zonal winds (McLandress, 2002a, b).Tides tend to break at ∼ 85 km (Lindzen, 1981) and like other large-scale atmospheric waves modify MLT temperatures by ∼ 5 K at midlatitudes (Marsh et al., 2006).In addition, tides are able to generate several kilometers of atmospheric species transport from above and below (Marsh et al., 2006;Smith, 2004).As was mentioned in the methodology section, we averaged a constant (localtime) hour in order to reduce the tidal effect, as it was not possible to average 24 h of VLF data due to the significant changes in the ionosphere between day and night.However, this procedure will not have completely solved the tidal effect problem for two reasons.Firstly, non-migrating tides have phases that are non-Sun-synchronous.As the phases of this type of tide do not follow the Sun's apparent motion in the sky, their total amplitude and phase during the midnight and midday 1 h average could cause leakage of the tidal oscillation into the VLF data's long-term frequencies, known as tidal aliasing (Oberheide et al., 2003).Secondly, the migrating tides, which are Sun-synchronous, should have no effect over the data when using a constant local time at a single station.However, as the transmitter and receiver are always located at different latitudes and longitudes, the migrating tides are in different stages of their phase along the GCP, and they may have an additional influence over the VLF received signals, amplifying the SAO.

Conclusion
A strong SAO was detected in the nighttime D region using the amplitude of ground-based VLF NB signals.This oscillation dominates over all other long-term oscillations.
We conclude that the main source of the SAO is most likely to be due to NO x molecules transported from upper levels of the atmosphere.This transport results in enhanced ionization and the creation of additional free electrons in the nighttime D region, thus inducing the SAO signature on VLF NB amplitude measurements.Nevertheless, further research and analysis should be undertaken in order to confirm our conclusions.A good test would involve the use of both high-end chemistry and general circulation models due to the complexity of the D region, or analysis of NO + measurements from space, which might be acquired in the future using instruments such as NASA's Middle Atmosphere Sounder and Thermal Emission Radiometer (MASTER) (Mlynczak et al., 2014).
In addition, as far as the authors are aware, no current VLF wave propagation model (e.g., Ferguson, 1980) takes into account SAO forcing of the D region and hence the impact on received VLF signals.As we have shown in this paper, the SAO influence over VLF signal attenuation is significant, affecting the received signal amplitudes by several decibels.VLF signal studies are an important tool for understanding the D region of the ionosphere, being low-cost and having high temporal resolution and potentially high spatial resolution (by using numerous receivers at many different locations).Therefore, propagation models should take this oscil-lation into consideration, in order to acquire better and more precise results, particularly over long time periods.

Figure 2 .
Figure 2. Midday (left panels) and midnight (right panels) 1 h mean 30-day running-average time series of MH-NSY (top panels) and DN-NWC (bottom panels) transmitter-receiver GCP deviation from the mean amplitude (solid black curves).The dashed red curves show the combination of the SAO, AO, and linear fit to the data time series (see Eq. 1).A Pearson's correlation coefficient between the red and black curves is shown at the bottom right of each panel.

Figure 3 .
Figure 3. Lomb-Scargle periodogram of the midnight MH-NSY (top) and DN-NWC (bottom) GCP 1 h mean VLF amplitude anomalies in arbitrary power units.The dashed red line denotes the 95 % confidence level.