TID characterised using joint effort of incoherent scatter radar and GPS

Travelling Ionospheric Disturbances (TIDs), which are caused by Atmospheric Gravity Waves (AGWs), are detected and characterised by a joint analysis of the results of two measurement techniques: incoherent scatter radar and multiple-receiver GPS measurements. Both techniques to measure TIDs are already well known, but are developed further in this study, and the strengths of the two are combined, in order to obtain semi-automatic tools for objective TID detection. The incoherent scatter radar provides a good vertical range and resolution and the GPS measurements provide a good horizontal range and resolution, while both have a good temporal resolution. Using the combination of the methods, the following parameters of the TID can be determined: the time of day when the TID occurs at one location, the period length (or frequency), the vertical phase velocity, the amplitude spectral density, the vertical wavelength, the azimuth angle of horizontal orientation, the horizontal wavelength, and the horizontal phase velocity. This technique will allow a systematic characterisation of AGWTIDs, which can be useful, among other things, for statistical analyses. The presented technique is demonstrated on data of 20 January 2010 using data from the EISCAT incoherent scatter radar in Tromsø and from the SWEPOS GPS network in Sweden. On this day around 07:00–12:00 UT, a mediumscale TID was observed from both data sets simultaneously. The TID had a period length of around 2 h, and its wave propagated southeastward with a horizontal phase velocity of about 67 m s −1 and a wavelength of about 500 km. The TID had its maximum amplitude in Tromsø at 10:00 UT. The period length detected from the GPS results was twice the main period length detected from the radar, indicating a different harmonic of the same wave. The horizontal wavelength and phase velocity are also estimated from the radar results using Hines’ theory, using the WKB approximation to account for inhomogeneity of the atmosphere. The results of this estimate are higher than those detected from the GPS data. The most likely explanation for this is that Hines’ theory overestimated the values, because the atmosphere was too inhomogeneous even for the WKB approximation to be valid.


Travelling ionospheric disturbances and atmospheric gravity waves
Travelling ionospheric disturbances (TIDs) are waves in the ionosphere with time periods from tens of minutes up to 2-3 h and wavelengths typically longer than 100 km.The waves can travel globally over distances of thousands of km, including travels between high latitudes and the equatorial region.
TIDs were first noted in ionosonde data in Australia in 1937-1939, as periodic disturbances in the reflection height, which were noted to travel between different locations.Munro (1950) reported the first systematic measurement and analysis of the vertical and horizontal wave parameters of these disturbances, using multiple-location ionosondes around Sydney.
TIDs were extensively theoretically analysed and described in the pioneering study by Hines (1960).According to this, TIDs are caused by atmospheric gravity waves (AGWs) propagating in the neutral thermosphere.These waves are generated in the lower atmosphere and travel upwards.In the ionosphere, the resulting wave patterns of the ionised gas are a measurable signature of the neutral wave.
According to the reviews by Hunsucker (1982) and Hocke and Schlegel (1996), TIDs can be grouped into two categories: large-scale TIDs with horizontal wavelength over 1000 km and oscillation periods of 30-180 min, and medium-scale TIDs with horizontal wavelengths of several hundreds of km and periods of 15-80 min.Large-scale TIDs propagate with phase velocities of 200-1000 m s −1 , comparable to the speed of sound in the thermosphere, while medium-scale TIDs typically propagate more slowly, with phase velocities of 50-250 m s −1 .Modern understanding is that medium-scale TIDs are caused by both AGWs and ionospheric processes, while large-scale TIDs result mainly from magnetosphere-ionosphere coupling processes (Hunsucker, 1982;Hocke and Schlegel, 1996).Kirchengast et al. (1996) described the mechanism of the coupling between AGWs in a neutral atmosphere and TIDs in the ionosphere.They confirmed this mechanism experimentally, by modelling the background atmosphere, the AGW, and this coupling, and finding a good agreement when comparing the resulting TIDs to measurements, giving confidence that at least some TIDs are indeed caused by AGWs.Hocke and Schlegel (1996) gave an extensive review of TIDs, presenting the models and mathematical techniques developed and the observations performed at the time.

Measurements using incoherent scatter radar
An incoherent scatter radar (ISR), such as EISCAT, located in northern Fennoscandia (http://www.eiscat.se/)(Folkestad et al., 1983) is a useful tool to detect TIDs.ISRs measure profiles of electron density, ion velocity, ion temperature and electron temperature in the ionosphere.Even though such measurements have limited horizontal resolution, their high vertical and temporal resolution allows a good detection of the wave features of TIDs.Hines (1960) showed from theoretical analysis that the phase and group velocities of AGW-induced TIDs are considerably different.While the energy in the waves is transported obliquely upwards, the phase fronts propagate obliquely downward.As a result, when observed in a vertical profile, the waves are perceived as propagating downward.This has been confirmed in various measurements.Thome (1964) presented the first observations of TIDs using ISR, using electron density profiles measured by the radar site at Arecibo, Puerto Rico.Crowley et al. (1984) presented the first observations using electron density profiles measured by the EISCAT radar in Tromsø, Norway.Schlegel (1986) detected from EISCAT data both gravity waves and tidal waves by their respective characteristics features, and studied their interference.Also many others have measured and analysed TIDs using ISR.An overview of these observations up to the 1980s can be found in the review by Williams (1989).Nygrén et al. (1990) demonstrated the strength of the detailed height-time diagrams which can be provided by ISR data for the detection and characterisation of TIDs.They analysed wave patterns in electron density and ion velocity in a 4.5 h long data set, for their geometric properties and interactions with other observed phenomena, such as tidal waves, electric field fluctuations and a sporadic E layer.
Also Hocke et al. (1996) showed that the downward propagating phase fronts observed in vertical profile measurements are suitable to distinguish AGW-induced TIDs from other quasi-periodic disturbances which do not have such a significant phase variation with height.They detected TIDs from 4 years of EISCAT data, by looking for wave patterns of similar frequencies in all four measured parameters, and over an extended height range in the E and F layer.Based on the 45 observed TIDs, they presented typical phase and amplitude profile characteristics of the four parameters.Vlasov et al. (2011) used this property of TIDs to detect them from a long-term almost-continuous period of EIS-CAT data, measured in Svalbard during the international polar year (IPY) campaign of March 2007-February 2008.The data of all parameters and at altitudes from about 100 to 450 km were band-pass filtered around period lengths varying from 0.2 to 3 h.From the results, the dominant resonant frequencies were searched.Days with significant dominant wave frequencies at many altitudes were selected, and the height-time plots of these band-pass filtered data of selected days were examined by hand, to find the cases of significant downward propagating phase fronts.Results from the study of Vlasov et al. (2011) are summarised as follows: -AGW-TIDs were observed most frequently in spring and summer.
-Most common period lengths are between 0.5 and 1.3 h; both medium-and large-scale TIDs can be included in this.
-Of all parameters measured by EISCAT, TIDs are mostly observed in ion velocity data.
-Amplitude profile results were in agreement with previous studies.
-No significant relation between the occurrence of TIDs and geomagnetic activity was found.
The study in Sect. 2 of the current paper builds further on the work of Vlasov et al. (2011).In the analysis, the day of 20 January 2010 is chosen as an example case.

Measurements using GPS
Following the original work by Saito et al. (1998), there have been a number of studies addressing properties of TIDs based on the data from networks of ground-based GPS receivers located at mid-latitude: e.g. in Japan (Tsugawa et al., 2004), Southern California (Kotake et al., 2007) and a wider GPS network covering mainland United States and southern Canada (Nicolls et al., 2004;Tsugawa et al., 2007).In these studies the TIDs were detected using maps of differential Total Electron Content (TEC) derived from dualfrequency carrier phase measurements of GPS navigation signals.The studies have confirmed that the occurrence of large-scale TIDs correlates with periods of high geomagnetic activity, while medium-scale TIDs tend to be quiettime events.Kotake et al. (2007) separated medium-scale TIDs (MSTID) observed at mid-latitudes into three classes: (1) daytime MSTIDs, being an ionospheric manifestation of atmospheric gravity waves as suggested by Hines (1960); (2) night-time MSTIDs, possibly caused by electrodynamic phenomena, and (3) dusk MSTIDs generated by the moving solar terminator.Nighttime MSTIDs have a maximum occurrence rate during summer while daytime MSTIDs occur more frequently during winter.
Extensive and numerous GPS receiver networks over mainland United States, Japan, Central and South Europe (maintained for geodetic purposes as well as earthquake forecasting) facilitated studies of TID properties at mid-latitudes.However, at high latitudes, in the vicinity of the auroral oval, such studies have not been conducted.This was due mainly to the fact that GPS networks in Arctic regions were sparse or unavailable to the ionospheric research community.In the last few years new extensive networks of GPS receivers have been established in Fennoscandia, including ∼ 90 receivers in Finland (operated by Geotrim Oy), over ∼ 180 receivers in Sweden (operated by SWEPOS) and over 100 receivers in Norway (operated by Statens Kartverk).This opens an opportunity to study TIDs at auroral and sub-auroral latitudes, where the ionosphere is highly dynamic (due e.g. to auroral activity) and where TIDs have often been found to originate from.This also allows more joint analyses together with the many other ionospheric instruments which are located there, including the EISCAT incoherent scatter radar in northern Fennoscandia.
The study in Sect. 3 of this paper will analyse an example of a medium-scale TID detected from data of the Swedish GPS network on the same day as analysed in Sect.2: namely, 20 January 2010.The study demonstrates that TIDs can be detected from GPS data at auroral latitudes equally well as at less turbulent mid-latitude regions.The wave parameters are calculated from the data more systematically than in previous studies.

Three dimensions
Because TIDs propagate in both the horizontal and the vertical dimension, it would be most useful to study them 3dimensionally, to characterise their properties in all dimensions.Unfortunately, most measurement techniques mainly retrieve information in either the horizontal or the vertical.Several actions have however been taken to combine measurements in the different dimensions.Munro (1950) observed TIDs using multiple ionosondes around Sydney.Observing the fluctuations in reflection height, he detected the horizontal motion of a few TIDs, while at the same time noting a vertical downward progression of the waves.Lanchester et al. (1993) analysed a few events of TIDs observed in vertical electron density profiles measured by EISCAT.They compared them with simultaneous measurements using ionospheric sounder measurements of a 200 × 200 km region around Tromsø, from which horizontal wave motions could be detected.Nicolls et al. (2004) combined GPS TEC data over North America and ISR profiles from the Arecibo observatory in Puerto Rico.They were able to detect horizontal and vertical wave patterns from the two data sets.Albeit useful, all the above techniques require the different horizontal and vertical wave parameters to be determined by hand, which makes a systematic or long-term analysis difficult.Ma et al. (1998) employed a single 3-dimensional technique, using incoherent scatter radar measurements alternated between four directions.They were able to retrieve systematically both horizontal and vertical wave numbers from these data.However, since this data set consists in the horizontal dimension only of a 4-point measurement, with maximum distance of 170 km at 300 km height, corresponding to a few tenths of a wavelength for medium-scale TIDs, this makes the accurate retrieval of the horizontal wave parameters still difficult in the case of noisy or distorted signals.
In this paper, the issue is approached in the same way as Nicolls et al. (2004): a combination of GPS and ISR, which ensures a high resolution and large range in both horizontal and vertical dimensions.At the same time, a systematic retrieval of the wave parameters in both dimensions is used.This will enable repeatable and long-term data analyses in the future.
the software, and band-pass filtered data for these dominant frequencies were produced.After this, these filtered data were examined by hand, to detect downward-propagating phase fronts, which is the main criterion to decide the presence of an AGW-TID.
Albeit thorough, such an examination by hand has the disadvantage of not being objective and repeatable.For this reason, in the current analysis this detection of the downwardpropagating phase fronts is automated.At the same time, this means that the selection of the dominant frequency components (which was for Vlasov et al. (2011) mainly meant to pre-select the data that were to be examined by hand) can be skipped, and all data, band-pass filtered for all frequency bands of interest, can be automatically examined.This procedure will also allow us to determine systematically the vertical phase velocities of these waves.
The radar data from the EISCAT radar near Tromsø, Norway (69.59 • lat; 19.23 • long) on 20 January 2010 are processed to analyse AGW-TIDs.Although all measured parameters were analysed, only the electron density is presented in this paper.The raw data have a temporal resolution of 1 min, and an altitude resolution varying from 3 km (below 100 km altitude) to 20 km (at 300 km altitude).
The data set contains measurements taken alternatingly in three different directions.From these, only the measurements are selected which were taken in the vertical direction.This reduces the time resolution of the data to 2 min, which is still sufficient for the analysis of waves with periods of at least 10 min.
The altitudes of the profiles measured by EISCAT are slightly variable, as these are the result of the preprocessing calculations and hence dependent on the state of the ionosphere itself.In the study of this paper, a fixed set of altitudes is defined, and all profiles are linearly interpolated for these fixed altitudes.The set of altitudes ranges from 100 to 260 km with steps of 10 km.This covers the altitude range in which AGW-TIDs show their typical downward propagating phase fronts.
Figure 1 shows the electron density between 07:00 and 15:00 UT after these steps.This exposes between 09:00 and 13:00 some downward propagating features, which are typical for AGW-TIDs.

Filtering
In order to expose the AGW-TIDs of certain frequencies, the times series data, at every altitude separately, are band-pass filtered at many different frequency bands.The frequency bands used in this study are centred around 23 period lengths T i (h) which are logarithmically spaced ranging from 0.2 to 11.5 h.The filter bands for each period length T i are between lower frequency f li and upper frequency f ui , which are given by The filters used are second-order Butterworth filters, similarly as used by Vlasov et al. (2011).Band-pass filtering is achieved by low-pass filtering at two frequencies, and subtracting the results.The practical complication that a Butterworth filter requires a certain "startup-time" is overcome by adding relaxation intervals before and after each data set to be filtered, of a length of 100 samples.These relaxation intervals are filled with sine waves with frequencies equal to the central frequency of the pass band and amplitudes equal to the closest data sample value, and multiplied by cos 2 -shaped windows to make the amplitude decrease toward zero at the outer edges.The phase shift introduced by Butterworth filters is overcome by always applying the same filter twice: once forward and once backward.Some data time series contain gaps, which complicate filtering.These data gaps are linearly interpolated before filtering.Afterwards, the filtered data points corresponding to gaps in the original data are removed again, if these gaps are longer than 1/(4f ui ) (i.e.longer than 1/4 of any period in the pass band).The relaxation intervals before and after the data are also removed after filtering.Figure 2 shows an example of filtered data: the electron density filtered for periods around 1.05 h.This graph also shows clearly some downward-propagating fronts, which demonstrates that this is one of the frequency components present in the features seen in Fig. 1.The further procedure is designed to automatically detect and characterise such inclined structures.

Cross-correlation analysis
For these data, the cross-correlation between the data at various altitudes is examined.The cross-correlation coefficient between data at two adjacent altitudes is calculated as where X(t,h) represents the data, at altitude h, selected for t t − T < t < t t + T , and multiplied by a Hann window to smooth the edges; h is the height difference between two adjacent altitude layers = 10 km, E[ ] denotes an average over t, σ ( ) denotes a standard deviation over t.The Hann window is described as The length of the data section from which the correlation is calculated is equal to twice the period length T of interest, which means that two wave fronts of the supposed period length would fit in the section of data.
The time shift t between the signal X at two altitudes corresponds to a vertical velocity v z (m s −1 ) as From the same data, also the amplitude spectral density of the detected structures is estimated, as follows: The background to this formula is as follows.Of the variances of the signal X at the two altitudes, the geometric mean is calculated as the square root of their product (equivalent to an average on a logarithmic scale).The resulting mean variance, assumed to be equal to the integrated signal power, is divided by the filter bandwidth f ui − f li , to obtain a power spectral density.The factor 8/3 compensates for the loss of signal due to the Hann window.The square root converts the variance to a standard deviation.Finally, the factor √ 2 converts a standard deviation to an amplitude of a sine-wavelike signal.The resulting amplitude spectral density has the dimension of electrons m −3 Hz −1/2 .As a result, there is a value of this amplitude spectral density A connected to every value of a correlation coefficient C as defined above.
The correlation coefficient C and amplitude spectral density A are calculated for: -23 values of the period length T logarithmically spaced from 0.2 to 11.5 h,

TID detection: vertical phase velocity profile
For each value of T and t t , the resulting values of v z are subsequently examined as functions of the altitude h.(Note that for each value of h, T and t t , there can be none, one, or more values of v z .)From these values, a dedicated piece of software looks for patterns where the vertical phase velocity increases with height within certain limits (or is constant with height).A pattern is searched according to the following rules: -the values of v z in the pattern can be approximated as where h R is a reference height taken as 175 km, and b ≥ 0. This assumed profile is chosen because TIDs have often been found to have a vertical phase speed increasing with height (e.g.Ma et al., 1998).
-The procedure starts by taking the value of v z corresponding to the largest amplitude spectral density A as defined above, and proceeds consecutively with those with decreasing amplitudes, to check whether the values match the pattern.-There does not need to be a value of v z matching the pattern at every altitude, but there must be one at least three altitudes.
A curve is fitted to the values of v z found, according to the relation of Eq. ( 6).The goodness of fit of this curve is examined as follows: where The "quality factor" Q can range from 1 to infinity, being higher for a better fit and for a fit with more points.In general, it was found that a reasonable fit has Q > 1.5, and a good one Q > 2.0.
Figure 3 shows two examples of the values of v z vs. h as found from the correlation analysis on the data of Fig. 2 (T = 1.05 h).At 10:00, a pattern was found and a curve was fitted to it, with Q = 2.51.The curve is included in the graph.Note that the fact that all values of d h for which |d h | > d M are replaced by d M , effectively means that all points further away from the curve than a factor exp(1/2) (in this case two points) are treated the same as missing points, and do not affect the position of the curve.At 14:00, no pattern was found, according to the above requirements.
Next, consecutive cases (in time) where a pattern was found are analysed together.For every value of T , consecutive values of t t are searched where a pattern was found with Q ≥ 2.0.It is assumed that waves detected at the same period length, and at consecutive times within a few hours, are actually the same wave, and that the wave properties of this wave do not significantly change within this time interval.Therefore, for each patch of consecutive times, the points of v z as a function of h are considered together to fit a single curve to as many of the points as possible.The resulting curve is considered to be the height profile of the vertical phase velocity of the wave.Figure 4 shows an example for the times 09:00-12:00 and T = 1.05 h.The curve fitted to the points is described as Eq. ( 6) with a = 28.439km h −1 and b = 2.3071, and has Q = 4.60.Just like in Fig. 3, all points further away from the curve than a factor exp(1/2) do not affect the curve's position and are effectively ignored.
The above analysis is performed for every value of the period length T .As a demonstration of the result, the resulting quality factors Q of the detected downward propagating phase fronts which characterise TIDs are shown as a function of the period length T and the hour of the day t t in Fig. 5.It can be seen from this figure that almost all significant waves (with Q > 2) were found between 07:00 and 12:00.The period lengths T of these frequency components were all between 0.7 and 4.5 h, with the dominant component at T = 1.05 h.Of the significant fitted curves in this analysis, the resulting values of a and b are shown in Fig. 6 as functions of T .In this figure, the small dots are the fitted curves with Q ≥ 1.5, and the stars are those with Q ≥ 2. This figure shows that, at least on this day, the slope b of the height profile of phase velocity (equal to the strength of the height dependence) tended to increase with the period length, and the offset a (equal to the vertical phase velocity at 175 km altitude) tended to decrease.
Figure 7 shows a few of the corresponding height profiles of vertical phase velocity v z , as a snapshot at 10:00 UT.In this figure, the curves are plotted only at altitude ranges for which any values contributing to the fitted curve were found (see e.g. in Fig. 4).This figure shows that not only does the vertical phase velocity of the TID increase with height for all frequency components found, but the phase velocity also mostly decreases with increasing period length, at all altitudes.Since these frequency components are nevertheless part of the same wave structure, this dependence demonstrates the dispersive property of the atmosphere for AGWs.
Figure 8 shows two plots of the filtered data of electron density, as in Fig. 2. On top of the colour plots, the downward propagating wave patterns as detected from the data by the procedure described above are plotted for each hour of the day.Note that the curves do not necessarily mark the maxima of the waves.They show the motion of the phase point which at the time of detection is located at h = h R = 175 km.Just like in Fig. 7, the curves are only plotted for altitude ranges where they were detected, with Q ≥ 2.
Because the curves in each of the graphs were found during a consecutive block of time, they have been analysed together.Therefore the coefficients a and b are constant in each of the graphs, so all the curves at the different hours in either graph have the same parameters.(Only the different lengths of the curves give an optical illusion that this is not the case.) To establish an accuracy estimate of the vertical phase velocity, the relative error of the points used in each fit with respect to the fitted curve (as in Fig. 4) was calculated.The rms value of these relative errors was around 19 % for all fits (for different times and period lengths).This can be seen as an error estimate for any of the resulting values of v z .

Amplitude
In this section, the variations of the amplitude of the detected TID are assessed.Since the wave has been detected in electron density data, the amplitude obviously has the dimension of electron density [electrons m −3 ].This amplitude can be derived from the amplitude spectral density A, calculated in Eq. ( 5), by integrating it in the domain of period length T over the frequency band of the TID.However, it is hard to assess the limits of this frequency band reliably.Therefore, and since the main purpose of the amplitude assessment is to know its variation with height and time, it is assumed that the frequency band size does not change much, so that the amplitude's developments are reflected by those of the amplitude spectral density (which has dimension electrons m −3 Hz −1/2 ).
The variations of the amplitude spectral density A of the detected TID are determined as follows.Of all the values of A calculated as in Eq. ( 5), the values are selected which led to a detected pattern in Sect.2.4.As a result, there is an amplitude spectral density value at every altitude, hour and period length for which a TID was detected in the previous section.The upper graph of Fig. 9 shows the profiles of A of all the TIDs observed from 08:00 to 12:00, and for all period lengths T between 0.8 and 3.0 h, all combined in one graph.Without trying to distinguish the different lines in this, it can clearly be seen that at the altitudes of 210-220 km, the amplitude of this TID was maximised, generally for all hours and period lengths.Comparing this result with Fig. 1, it can be seen that these are the same heights where at this time the electron density itself was maximised.
The height of maximum electron density amplitude is also similar to those found by Hocke et al. (1996) and by Ma et al. (1998), both of whom also analysed TIDs in electron density measured at EISCAT Tromsø.
The lower graph of Fig. 9 shows the values of the upper graph at the altitude of 220 km as a function of the hour of day.This figure shows that the TID reached its maximum amplitude at this location at 10:00 UT.This is near local geographic noon in Tromsø, which is at 10:43 UT.

TID detection using GPS
This section describes how TIDs are detected in the horizontal plane, using Total Electron Content ("TEC") measurements from GPS signals by the Swedish GPS network SWE-POS.The measurements measured by 143 SWEPOS GPS receivers from all 32 GPS satellites on 20 January 2010 were used.The receiver stations are distributed over latitudes from 55.3 to 68.4 • , and longitudes from 13.4 to 18.8 • .The map of Fig. 10 shows the distribution of the stations around Sweden, and the position of this area relative to that of the EISCAT radar in Tromsø.

TEC measurement using GPS
The raw data contain the phases of the signals received from the GPS satellites at the two frequencies f 1 = 1575 MHz and f 2 = 1228 MHz, at a time resolution of 1 min.Signals passing through the ionosphere suffer a frequency-dependent phase delay proportional to the electron density integrated along the path.Due to this, the integrated electron density can be calculated from the phase difference between the two received signals as (e.g.van de Kamp et al., 2013 where TEC is "total electron content" i.e. integrated electron density along the path; r e is the radius of an electron

Slant/vertical conversion, offset, and background TEC
Equation ( 8) gives the electron density integrated along the slant path.To get the vertical TEC, i.e. the integrated ionisation along a vertical column, it is assumed that the ionosphere is a thin shell at a known height h I .Often in TEC measurements, the mean height of ionisation or the peak of the F layer is used for h I , which means a height of around 350-400 km.However, in this study, where not the ionisation level but the wave structures therein are studied, it is considered appropriate to assume for h I the altitude at which these waves are strongest.Because from the EISCAT data, it was found that the amplitude of the observed TID on this day was maximised at 220 km (Fig. 9), h I is taken as 220 km.
The vertical TEC is then estimated as where ε I is the elevation angle of the path at the height h I .However, before this conversion can be made, another point needs to be taken care of.
Because of the inevitable 2π-ambiguity of phases, the absolute phase difference ϕ between the two signals always has an unknown offset.Consequently, also TEC in Eq. ( 8) has a unknown offset.This offset value stays constant during any section of continuous phase measurements between a receiver and a satellite, as long as the receiver stays in lock; such a section typically lasts maximum about six hours.(In addition, the receivers and satellites all have their own instrument biases, which remain constant for every receiver and satellite, and can therefore be taken care of along with the offset due to the phase ambiguity.) With the unknown offset included, Eq. ( 9) becomes where TEC m is the measured TEC and TEC o is the unknown offset.
The geographic point where this value of VTEC is valid is called the ionospheric pierce point (IPP), i.e. the point where the satellite-receiver path penetrates the (assumed thin) ionosphere at height h I .The value of VTEC thus measured should contain variations of different timescales: slow variations of background ionisation with solar radiation and due to the movement of the IPP through different latitudes, as well as faster variations due to scintillation and the TIDs.
To predict the slower variations, the International Reference Ionosphere model "IRI-2007" is used (http://iri.gsfc.nasa.gov/)(Bilitza and Reinisch, 2008).This model gives the background ionisation pattern that can be expected according to known theory, depending on latitude, longitude, altitude and time, using input information such as solar radiation.The background VTEC has been calculated according to the IRI model for a grid from 45 to 75 • latitude with 1 • resolution, from −20 to 50 • longitude with 2 • resolution, and for the entire day of interest with a time resolution of 5 minutes.The predicted background VTEC for the current measurements has been linearly interpolated from this grid for the latitudes, longitudes and times of the IPPs.This background value will be referred to as VTEC IRI .In Eq. ( 10), VTEC can then be written as the sum of this term and the faster variations, referred to as VTEC var .This modifies the equation as Using this equation, the correct value of the unknown offset TEC o can be estimated as follows.Different assumed values of TEC o will introduce different increasing or decreasing trends in VTEC var (through the factor sin ε I ).However, since the background ionisation is removed, VTEC var should not contain any long-term trends, but only short-term variations.Therefore, the correct value of TEC o is assumed to be the value which minimises the standard deviation of VTEC var .This can be calculated to be , where E[ ] denotes the average taken over the same periods as the integrals.This expression is evaluated over each section of continuous data between any satellite and receiver, during which the receiver stays in lock continuously and proper preprocessing has ensured the correction of any cycle slips in the phase measurements.After this, Eq. ( 11) is applied to obtain the estimated VTEC var .

Detrending
As the next processing step, the data are "detrended", i.e. any remaining slow variations are removed by high-pass filtering.This is done as follows: where VTEC LF is the VTEC var signal which has been lowpass filtered with a cos 2 -filter (a moving average window multiplied by a cos 2 -shaped function).
The filter should remove any slow variations in VTEC var which still remain after Eq. ( 11) -for instance, deviations of the diurnal VTEC variations from those predicted by the IRI model.After filtering, VTEC HF should contain only the faster fluctuations in the signals, caused by the TIDs as well as e.g.scintillation.
The optimum cut-off frequency of this filter is decided by the following considerations.The diurnal variations, which mainly have a period length of 24 h, should mostly have been removed as VTEC IRI in Eq. ( 11), but some components of them might still be remaining, which can have harmonics with periods down to 6 h.The main TID fluctuations have a period length of 1 h, although components up to 4 h are found (see Fig. 5).
However, note that these period lengths would be seen by a fixed observer, while the measurement point in this project is the IPP.Since the IPPs move along with the satellite, the observed period lengths differ from those for a fixed observer due to the Doppler effect.At 220 km altitude, the IPP can move at velocities varying from 32 m s −1 at high elevation angles, to 106 m s −1 at 30 • elevation (see later why this elevation is mentioned).The IPPs can move in all horizontal directions.To make sure that as many as possible variations due to TIDs are measured, while removing as many components of diurnal variations as possible, it was decided to use in Eq. ( 13) a cos 2 -filter with a window length of 5 h.Effectively, this means that TEC HF is a high-pass filtered TEC signal with a cutoff frequency of around 5.55 × 10 −5 Hz.Note that this means that some contributions of TIDs may nevertheless be removed, particularly when the IPP moves at comparable speed and direction to the TID, which could lead to very long observed period lengths.However, these cases are expected to be a relatively small portion of all variations measured.

Elevation window
The assumption in Eq. ( 9) of a thin ionospheric shell at a fixed height leads to the most accurate IPP localisations when the elevation angle is high, and gives progressively larger uncertainties for lower elevation angles.Because of this, an elevation angle window of 30 • is used, and all data at lower elevation angles are discarded.
This step also removes any edge effects of the high-pass filter applied in the previous section, since the start and the end of a satellite pass are always at low elevation angles, and will therefore be removed at this stage.This additional advantage is the reason to apply this step after the detrending step of the previous section.
It might be argued that an elevation threshold of 30 • is still quite low, as at 30 • elevation, an error of 10 km in the assumed height h I already leads to a horizontal location error of 15 km.However, this wide elevation window is necessary, since the elevation angle of GPS signals at high latitudes is inevitably relatively low.Since GPS satellites reach a maximum latitude of 55 • , they can only just reach the zenith position at the southernmost SWEPOS station, and reach only lower elevations at others, and a maximum elevation of only 74 • at the northernmost one.Due to this, the bulk of all the data generally have elevation angles below about 60 • , and therefore a relatively low elevation threshold is necessary to retain a sufficient amount of data.

Azimuth-dependent preparation
The analysis of the orientation of TIDs is done at a time resolution of 5 min.In preparation of the analysis, a range of azimuth angles is defined as 0, 5, 10, 15, . . .,175 • .This half-circle range is enough when analysing the orientation of elongated structures.(Later on, in the analysis of the propagation direction, a ± sign will provide the added information to cover the full circle of directions.)For each of the azimuth angles in this range, a line through a central point of the area of measurement is defined which extends from −1000 to 1000 km from the central point in the specified azimuth direction; see the example in Fig. 11.For the central point is chosen a location of 60 • lat, 16 • long.For each sample point (each 5 min interval), all of the IPP locations within 500 km distance from this line (for any satellite and any receiver station) are projected onto this line.The line is divided into bins of length of 20 km, and for each bin, the VTEC values of all projected IPPs in the bins are averaged.The result of this is a 1-dimensional cross-section of VTEC, for each azimuth angle and for each point in time.

Wave orientation analysis
Next, the orientations of any spatial variations in (detrended) VTEC are analysed.It is assumed that a TID consists of spatially elongated wavefronts, which propagate in a direction perpendicular to their elongation direction (note however that in this section, the motion of the waves is not yet looked at).When these wavefronts are projected onto lines of different orientations, the apparent wavelength of the projected variations should be shortest and most pronounced on the line perpendicular to the elongation direction.This direction, which is also assumed to be the propagation direction of the wave, will be called the "orientation azimuth angle" α x of the wave.
Because of these considerations, in this section the projection orientation is searched in which the wavelength of the projected variations is shortest.This is done as follows.
For each point in time and for each azimuth angle, the signal resulting from the preparation of the previous section (as a function of distance along the line) is spectrally analysed by FFT as a function of wave number.All local maxima of this spectral density as a function of wave number are selected.The inverse of these wave numbers are the apparent wavelengths λ α of structures at this azimuth angle, provided the spectral density at this wave number is at least 2 times the geometric mean of the entire spectrum, and the wave number is not at one of the edges of its range.
Next, from the entire 2-D range of azimuth α and wavelength λ α , the procedure selects the combination of α and λ α where the spectral component is most standing out from its surrounding spectrum.Starting from this point, other wavelengths λ α at other azimuths are searched to fit the equation to an rms error as small as possible.Here, α is the variable azimuth, and λ x and α x are the parameters which are optimised to make Eq.( 14) fit by least-square-error regression to the measured values of λ α .Initially, the factor before the sin 2 in the above equation was also left adjustable.However, it was found that the typical shape of the measured results always agreed with this function with this factor around 1/2.Because of this, the factor was fixed at 1/2, to reduce the number of unknowns.After this, λ x is the resulting value of the wavelength and α x the orientation azimuth angle of the wave.An example of this analysis is shown in Fig. 12 in a polar plot.The apparent wavelength λ α detected as a function of the azimuth is shown in blue, and the sin 2 -shaped curve of Eq. ( 14) fitted to them is shown in red.The green arrow represents the resulting values from this analysis: the direction of the arrow corresponds to the resulting orientation azimuth α x of the TID (the direction in which the apparent wavelength is shortest), and its length corresponds to the resulting wavelength λ x .
Figure 12.Example of the orientation analysis in a polar plot.The apparent wavelengths λ α as a function of the azimuth α (blue), and a sin 2 -shaped curve fitted to them (red).The direction of the green arrow corresponds to the resulting orientation azimuth α x of the TID (104 • ), and its length corresponds to the resulting wavelength λ x (574 km).Note that in the analysis, only α from 0 • to 175 • are considered, but in this graph, for clarity all blue dots are plotted twice: also for α + 180 • .
It was found that also the wavelengths of other local secondary maxima of the FFT spectrum could be fitted with the function of Eq. ( 14), with usually the same value of α x at the same time sample, but with other (mostly smaller) values of λ x .These other values of λ x can be considered harmonic wavelengths of the same wave.The software is at this point not suited to process these wavelengths separately, hence these secondary wavelengths are not further considered in this paper.
This analysis is performed for every (5 min) sample of the entire day.Resulting cases where the rms error of the natural log of Eq. ( 14) was larger than 0.18 were discarded, for being considered bad fits, as well as results based on fewer than 21 points (of the 36 azimuth values), as these are considered unreliable fits.
Next, the results were analysed per hour.For each hour (0,. . .,24), a time window was defined of 2 h centred on the respective hour, and the results of λ x and α x in this window were averaged in the complex plane, as illustrated in Fig. 13.This results in one value for λ x and α x for each hour.
Figure 14 shows the result of this analysis of α x and λ x for the whole day: both the results for the individual samples and the hourly averages.In these graphs, the circles are values which are considered unreliable, either because they are averages of fewer than 12 data points, or because the rms distance between the points and the average in Fig. 13 is larger than 175 km.The stars, which represent the reliable results, show that significant wave orientation results were found in the period 05:00-15:00.In this period, the orientation azimuth angle α x was around 110 • , and the wavelength λ x was mostly around 500-550 km.Only after 12:00, λ x is smaller: down to 400 km.
It should be noted that this orientation azimuth is assumed to be valid at the altitude of 220 km, which is the altitude of maximum TID amplitude found in Sect. 2 and, because of that, is the assumed altitude of the main wave structure in the VTEC calculation of this section.Ma et al. (1998) have noted that the azimuth angle of the wave vector of a TID can vary with height (from 90 • at 170 km to 180 • at 240 km, in their case), so this azimuth should not be assumed to be height independent.

Wave propagation analysis
Next, the propagation properties of the wave are studied, for the period during which a significant orientation was found.First, an overall average azimuth angle and wavelength for this period were estimated by averaging the 5 min λ x and α x in the complex plane, in the same way as shown in Fig. 13, over this entire period.It was verified that the results over this period were stable enough to be represented by one overall average value, by checking that also now, the rms distance between the points and the average was not larger than 175 km.The average values for this period were found as α x = 108 • and λ x = 500 km.
Next, the VTEC data of this period are used which were projected (as described in Sect.3.1.5)onto the line in the direction closest to this average α x .These data are analysed for their propagation of waves with this wavelength λ x , as follows.The mean of the data is subtracted: where d is the distance from the central point along the projection line.The resulting x is multiplied by a complex wave of the expected wavelength, and the product is averaged over the range of d: where λ x = the average wavelength over the period as described above.(This is essentially the same as a Fourier transform at the wavelength λ x .)The resulting F is a complex number with an angle equal to the phase of the frequency component of wavelength λ x in VTEC(d). Figure 15 shows an example of the values of VTEC between 08:00 and 09:00, as well as cosine waves with the wavelength λ x and the phase of the argument of F .It can be seen that with the progression of time, this wave progresses on average toward the right.At the same time, many smaller-scale fluctuations are also seen; these partly consist of the higher-harmonic parts of the same wave, which led to the smaller wavelengths detected in the procedure described in Sect.3.2.However, also random fluctuations unrelated to the wave are observed.
To detect the wave progression, at each point in time, the progression distance is calculated from the (unwrapped) argument of F as , in which case the amplitude of the wave is considered too small to be significant.The propagation velocity, i.e. the horizontal phase velocity, is the time derivative of s p (with t in hours): To find this velocity, a straight line is fitted to the samples of s p , the gradient of which is the phase velocity.The values of s p used and the fitted line are shown in Fig. 16.Here, in spite of some gaps and some variations, the first-order trend of s p (t) is clearly linear and given by the straight line.It can also be seen that the slope of the line is in agreement with the variation from one sample to the next in many of the sepa- rate patches of the data.The resulting value of the horizontal phase velocity for this period is v x = 67.2m s −1 .
As an accuracy estimate, the rms error of the coefficients resulting from the curve-fitting procedure was calculated.The result for the slope was a rms relative error of 2.5 %, which can be seen as an indication of the rms relative uncertainty of v x resulting from this procedure.

Visualisation
For the sake of a visual examination of the TIDs, the VTEC data as described in Sects.3.1.1-3.1.4,after detrending and applying the elevation window, were also collected in bins of a horizontal grid, with bins sizes of 0.5 • in latitude and 1 • in longitude, for every 5 min.This can help the human eye to recognise wave patterns in the data.However, note that these gridded data were not at all used in the automated TID detection described in the previous sections.Note also that these data have not been filtered in any way and therefore also show other variations, which can obscure the view to the TIDs. Figure 17 shows an example of these gridded data in colour plots, at times during which a TID was detected in the previous sections.In these plots, comparing one picture to the next can show the gradual developments.Note that to see significant developments, it is better to concentrate on the variations around zero (yellow, green and cyan) rather than the extremes (red and dark blue).
The wave properties that were found in the previous section for the period during which a moving wave was observed, are also indicated in these plots.The direction of the dashed lines show the elongation of the observed wavefronts (i.e.perpendicular to α x ), and the distance between the two dashed lines in each plot indicates the wavelength λ x .The arrow points in the direction of α x , and its length is equal to the phase velocity v x multiplied by 1 h.
During these times, a wave front was detected elongated northeast-southwest; this is most clearly seen in the graph at 10:20; in the other graphs the wavefront is seen developing.The procedure assumes the wave to propagate perpendicular to its elongation, and therefore in the azimuth direction 108 • .This propagation is best seen in the northernmost part of the map, where a wave is present in all four of the graphs.Note that during the half hour shown here, the wave propagates only over a distance of half the length of the velocity arrow shown, which is about a quarter of the wavelength.
Meanwhile, since the data shown in these graphs have not been filtered, also many other variations are seen, of timescales both comparable to the TID and shorter, and which develop and move in different directions.Part of this movement can be the group velocity of the TID.The group velocity can differ from the phase velocity in magnitude, but unfortunately, using the current data it is not possible to determine the group velocity of the wave.The small-scale variations seen (similarly as those in Fig. 15) can consist partly of higher harmonics of the observed TIDs.On the other hand, both the smaller-and larger-scale variations also show random variations which are unrelated to the TID, and which make detection of the propagating wavefronts from these graphs somewhat difficult.

Comparisons and discussion
In this section, some results of the TID analysis from the GPS data in Sect. 3 are compared with those obtained from the EISCAT data as presented in Sect. 2. Furthermore, results from this paper are evaluated and compared to those from other TID studies.

Direct comparison: period, horizontal wavelength
From the results of the wavelength λ x and the phase velocity v x obtained in Sect.3, the period length T [h] of the TID can be calculated as λ x /3.6v x (with λ x in km and v x in m s −1 ).From the results in the period where a significant TID was found, this gives a values of 2.07 h.The EISCAT results (Sect.2, e.g.Fig. 5) showed in the morning a range of period lengths from 0.7 to 4.5 h, with the dominant period length at 1.05 h.The factor of 2 between these results suggests that the GPS technique detected a lower harmonic of the same wave pattern than the EISCAT technique.
At this point, it is worth noting from Sect.3.2 that the wavelength of 500 km found from the GPS results corresponded to that of only one (the strongest) spectral component in the data.Other harmonic components were also found, with wavelengths which varied in a similar fashion with respect to the azimuth angle as the main one, and which therefore are part of the same elongated structure.These secondary wavelengths were mostly shorter than the main one, and are likely to include the main 1 h harmonic detected from EISCAT.Unfortunately, the current version of the GPS software does not allow a separate systematic analysis of these secondary harmonic components.
Nevertheless, if it is assumed that all wave components observed in the GPS data have (approximately) the same horizontal phase velocity, i.e. 67 m s −1 , the horizontal wavelengths of these components can be estimated from this velocity and the period lengths measured from the EISCAT data.This way, the main component detected from EISCAT at T = 1.05 h would correspond to λ x = 254 km, and the detected band of components between 0.7 and 4.5 h would correspond to λ x between 170 and 1100 km.

Timing; travel velocity
The time during which the TID was significantly detected from the EISCAT data, between 07:00 and 12:00 (see Fig. 5), is in the middle of the time it was detected from the GPS data, 05:00-15:00 (Fig. 14).From the EISCAT data, the largest amplitude was detected at 10:00 (Fig. 9).Unfortunately, the amplitude (and therefore its time of maximum) could not be reliably detected from the GPS data.
Considering that the location of EISCAT Tromsø is north of the area used in the GPS detection, at 1100 km from the central point mentioned in Sect.3.1.5,it might be expected that the "travelling" property of the TID would have caused a significant time shift between the two measurements.However, seen in the direction of orientation detected, the two points are not so far apart: the projection of Tromsø onto the line with azimuth of 108 • , is only 212 km away from the central point.This means that, assuming the wave fronts are elongated in the perpendicular direction over at least about 1000 km, Tromsø and the point in central Sweden are less than half a wavelength apart.More importantly however, it is hard to estimate what time shift should have been expected between these two measurements based on the displacement of the TID, as it travels at the group velocity, which velocity is unknown.After all, note that the horizontal group velocity of the wave need not be the same as its horizontal phase velocity.
Unfortunately, it is not possible to determine the travel velocity or direction of the TID from the current combination of observations.This is mainly because the GPS analysis does not allow a spatial determination: it detects the phase speed and direction of a wave inside the analysis area (of roughly 1000 km across), but does not determine any specific location of the wave pattern within this area at any point in time.Hence, the exact position of the TID in the area cannot be pinned down at any given time, and hence neither can a travel velocity from Tromsø to this location.

Horizontal parameters via Hines' dispersion theory
The results of λ x and v x obtained in Sect. 3 can also be compared with the horizontal wavelength and phase velocity that can be estimated from the parameters as derived from the EISCAT data in Sect.2, using Hines' theory (Hines, 1960).According to this theory, the relation between the horizontal and vertical wave parameters depends on many factors, due to the many influencing parameters, such as the dispersive nature of the propagation medium, the interaction between gravity and pressure, and the interaction with the neutral wind velocity.
Hines' dispersion equation assumes an atmosphere which is stationary in the absence of waves, and uniform in temperature and composition.The gravity field is also assumed uniform.Second-and higher-order terms in the oscillation amplitude are neglected, and the wave is assumed to occur adiabatically.The wave is considered in a 2-D space, where the z dimension is the vertical and the x dimension is the horizontal direction in which the wave is propagating.The wind speed component in this horizontal direction is also considered.The equation (in the form as formulated by Lanchester et al., 1993) is as follows: where ω r is the intrinsic angular frequency of the wave, i.e. in the frame of reference moving with the medium (by the wind) (rad s −1 ) and is equal to ω − k x u w , where ω is the observed angular frequency of the wave, equal to 2π/3600T (rad s −1 ); ω a is the cutoff frequency for acoustic waves (rad s −1 ), equal to γ g/2c; ω b is the buoyancy frequency (rad s −1 ); γ is the ratio of specific heats; g is the gravity acceleration (m s −2 ); c is the velocity of sound (m s −1 ); u w is the horizontal wind speed in the direction of the wave  (m s −1 ); k x is the horizontal wave number (rad m −1 ); and k z is the vertical wave number (rad m −1 ), equal to ω/v z .This equation can be used to estimate the horizontal wavelength of the wave from the measured period length T and vertical phase velocity v z .From the value of k x resulting from the equation, the horizontal wavelength is found as λ x = 2π/1000k x (km), and horizontal phase speed as v x = 2π/3600k x T (m s −1 ).(Note that v x is defined relative to the ground.) The parameters ω a , ω b and c are height dependent, which means that the propagation medium of the wave is not uniform in the vertical direction, violating an assumption used in deriving Hines' dispersion equation.Einaudi and Hines (1970) dealt with this problem using the WKB approximation (Wentzel-Kramers-Brillouin), where the solution consists of wave functions with parameters varying depending on the input parameters.They showed under which conditions this approximation is justified (see below).Following this approximation, the height-dependent parameter values will be inserted into the equation at different altitudes.
The speed of sound c in a gas is given by where k is Boltzmann's constant (= 1.380648813 × 10 −23 J K −1 ), T k is the temperature (K) and m is the molecular mass (kg).For a compound gas such as dry air, m is the average of the molecular masses of its components, weighted by their respective number densities.Up to about 90 km height c varies effectively only with T k , but above this, the change in air composition also has a strong effect on c.
The buoyancy frequency ω b , or Brunt-Väisälä frequency, is for an isothermal atmosphere given by However, for gravity waves in a non-isothermal atmosphere, Einaudi and Hines (1970) showed that under conditions related to the horizontal phase velocity (this will be described in more detail later), the WKB approximation re-mains valid if ω b is modified as where dH /dh is the height derivative of the atmospheric pressure scale height H , given by The ratio of specific heats γ is the ratio between the heat capacity at constant pressure to the heat capacity at constant volume.For gases, this depends on the kind of gas and on temperature.For air, it can be estimated as the average of the values for the component gases, weighted by their respective mass densities.Lange (1973) gives values of γ for N 2 , O 2 , He and Ar dependent on temperature, while for other gases γ can be estimated as 5/3 for monoatomic gases and 1.4 for diatomic gases.This information allows us to estimate the effective value of γ for thermospheric air from the height profiles of temperature and mass densities of the atmospheric components.
The gravitational acceleration g decreases with altitude according to Newton's law of universal gravitation.
Equation ( 19) was applied to the atmospheric configuration of 20 January 2010, 10:00 UT above Tromsø, at altitudes 160, 180, 200, 220 and 240 km.The values of the height-dependent parameters for these altitudes were derived as described above (with Eq. 22 for ω b ).For this purpose, the profiles of temperature and component densities of the atmosphere for this time and location were obtained from the NRLMSIS-00 atmospheric model (Picone et al., 2002).This model can be run online at the website of CCMC (Community Coordinated Modeling Center; http://ccmc.gsfc.nasa.gov/modelweb/models/nrlmsise00.php).The results of this for all relevant parameters are given in Table 1.
To estimate the horizontal wind speed in the direction of propagation, the model TIE-GCM (Thermosphere Ionosphere Electrodynamics General Circulation Model; http:// www.hao.ucar.edu/modeling/tgcm/)(Richmond et al., 1992) has been applied for the same time, place and heights, also www.ann-geophys.net/32/1511/2014/at the CCMC website (http://ccmc.gsfc.nasa.gov/models/modelinfo.php?model=TIE-GCM).The results for this for the wind speed magnitude v w , wind speed azimuth direction α w , and the wind speed u w in the direction of wave propagation (α x = 108 • ), are included in Table 1.According to these values, the thermospheric wind was (around) northward, as is usually the case at high latitudes around local noon.
The acoustic cutoff frequency ω a is the lower frequency limit for acoustic waves, while the buoyancy frequency ω b is the upper limit for gravity waves.The values of the intrinsic frequency ω r of the TID components observed from the EIS-CAT results vary between 0.0004 and 0.0025 rad s −1 .Since this is consistently below the values of both ω a and ω b resulting in Table 1, it can be concluded that this wave is a gravity wave.
The dispersion equation was applied to the values of T and v z as obtained from the EISCAT results at the reference time 10:00 UT.All period lengths detected at that time (with Q ≥ 2) are considered.The vertical phase velocities v z of these wave components increase with height, and decrease with period length, as e.g. can be seen in Fig. 7.
Equation ( 19) is a fourth-order equation in λ x , with up to four solutions, two positive and two negative.Of the positive solutions, the smallest one becomes zero for a zero wind speed, indicating a standing wave whose propagation seen with respect to a fixed reference frame is due to the wind carrying the oscillations along.The other solution, i.e. the largest value of λ x , is the one sought, indicating a freely propagating wave.
The upper graph of Fig. 18 (lines) shows the horizontal wavelength λ x resulting from the equation.The horizontal phase speed v x is shown in the lower graph.According to Hines' theory, λ x and v x should be independent of altitude.Indeed in Fig. 18, the parameters are not very dependent on altitude, apart from small differences.Apart from measurement inaccuracies, these differences may also be due to the approximation of the equation applied to the nonhomogeneous atmosphere.A similar conclusion was also drawn by Lanchester et al. (1993).
The graphs of Fig. 18 also include the values of λ x and v x as observed from the GPS results in Sect.3, for the period length 2.07 h as observed there (stars).It appears that the values as estimated from the EISCAT results using the dispersion equation are much larger than those: for the same period length, these are around 2000 km for λ x and around 250 m s −1 for v x , while those measured from the GPS data are 500 km and 67 m s −1 respectively.
The relative error of the result of the dispersion equation as a consequence of a relative error of 19 % in v z (see Sect. 2.4) was also calculated: this resulted in a relative error in k x (and λ x and v x ) smaller than 20 %, depending on height and period length.The rms relative error of v x from the GPS procedure was 2.5 % (see Sect. 3.3).The discrepancy between the values in Fig. 18 cannot be explained by these uncertainties.It might be considered that the wind data, resulting from a global model, introduce an extra uncertainty in the current calculation, as the momentary wind may have been different.However, also this cannot explain the discrepancy observed, as the wind speed has only a minor impact on the result.For example, in an assumed situation with a wind speed of double the values of Table 1, at T = 2 h the values for λ x and v x would be only up to 15 % lower than the values obtained here.(Conversely, a lower wind speed would lead to higher λ x and v x , and thus even increase the discrepancy.)Also, considering the other positive solution from the fourth-order Eq. ( 19) does not help, as this solution would give wavelengths varying between 3 and 45 km, and strongly dependent on altitude.This is clearly not a realistic solution.
The explanation for the discrepancy may lie in the condition for validity of the WKB approximation of Hines' dispersion equation for gravity waves in an inhomogeneous medium (as briefly mentioned above).Einaudi and Hines (1970) showed that if the condition is not satisfied, and there is a significant difference between Eqs. ( 21) and ( 22) -i.e. a significant impact of the inhomogeneity of the atmosphere (dH /dh = 0) on ω b -then the WKB approximation will break down.Resulting from the current calculation of Hines' equation, at T = 2.07 h, ω r /k x varies from 200 m s −1 at h = 160 km to 450 m s −1 at h = 240 km.This is smaller than c (see Table 1), but not necessarily much smaller.Furthermore, the relative difference between Eqs. ( 21) and ( 22) is within 1 % at h ≥ 200 km, but is 9.4 % at 160 km, which can be significant.On the other hand, note that the strongest effect of atmospheric inhomogeneity on ω b is for the lowest height, while the highest ratio of ω r /k x to c is for the highest height.It can be concluded that this particular case is in the shady area where the validity of the WKB approximation just starts to breaks down, and the results should therefore be taken with caution.
It seems likely that due to this, Hines' dispersion equation in the WKB approximation overestimates the horizontal wavelength (and consequently the horizontal phase velocity) of this particular TID.This overestimation is a result of ignoring the effects of the inhomogeneity of the atmosphere, such as refraction effects due to the inhomogeneous velocity of sound.Because of this, the values of the horizontal wavelength and phase velocity measured from the GPS results are considered more reliable than those estimated from the EIS-CAT results using this theory.Lanchester et al. (1993) obtained somewhat similar conclusions.They also estimated the horizontal wavelength using Hines' dispersion equation from EISCAT vertical profiles measured in Tromsø.They obtained a strong height dependence of this wavelength, which they state is due to the non-uniform composition of the atmosphere.They point out that because of this, the values estimated using the dispersion equation can hardly give much more than the order of magnitude for the wavelengths.As a consequence, when they compared the results to horizontal wavelengths measured using the ionospheric sounder in a 200 × 200 km region around Tromsø, some estimated wavelength values were close to the measured ones, while others were significantly larger, and the comparison remained inconclusive.
It can be concluded that care must be taken when estimating horizontal wave parameters from vertical ones (or vice versa) using Hines' dispersion equation.For gravity waves, it should be verified that either the condition of Eq. ( 24) is satisfied, or the deviation of dH /dh from zero does not have a strong impact on Eq. ( 22).If neither is the case, measured values of the wave parameters should be considered more reliable than those estimated using this equation.

Comparisons with other studies
The TID observed in this study shows similar characteristics as those found in other studies.
The range of period lengths T observed from the EIS-CAT data, i.e. 0.7-4.5 h, with the dominant component at 1 h, and the one observed from GPS, i.e. 2 h, are close to the period lengths observed by others.Especially, other observations at EISCAT Tromsø yield similar results: Nygrén et al. (1990) observed on 10 July 1987 in the evening: T = 0.4 and 0.9 h, Hocke et al. (1996) found from 45 events in the period November 1987 to December 1991 that T was mostly between 0.5 and 1.8 h, Schlegel (1986) found from fifteen 24 h records in January 1982-August 1984 values of T in the range 0.25-3 h, and Ma et al. (1998) on 14 November 1990 found T = 0.73 and 1.2 h.Furthermore, from ISR data from Arecibo, Thome (1964) found T varying from 0.3 h to several hours.
For the horizontal wavelength λ x , while many smaller values have been observed, some observations were also similar to the 500 km of this study.Bristow et al. (1996) found from horizontal HF radar measurements over northeastern Canada (lat 59-75 • ) λ x between 200 and 450 km.Tsugawa et al. (2007) found from GPS measurements over the USA on 28 November 2006 daytime: λ x = 300-1000 km.Chum et al. (2012) found using multiple Doppler sounders in the Czech Republic values of λ x between 100 and 300 km.
The southeastward wave propagation direction of the TID observed agrees with other observations of horizontal parameters of TIDs in winter/daytime.Munro (1950) found from ionosonde measurements in Australia over 1 year from April 1948, that in local winter α x was between about 0 • and 60 • (i.e.northeast, equivalent to southeast in the Northern Hemisphere).Waldock and Jones (1986) found, using Doppler sounders in the UK, α x varying over the time of day, being near midday between 90 • and 150 • .Afraimovich et al. (1999) found from ionospheric interferometer measurements in Irkutsk in the daytime an average of α x = 160 • .Kotake et al. (2007) found from GPS measurement in South California that in winter/daytime α x = 120 • .Also the wave observed by Tsugawa et al. (2007) over the USA on 28 November 2006 in the daytime with λ x similar to this paper, propagated southeastward.Chum et al. (2012) found that the TIDs over the Czech Republic in the winter propagated with azimuths between 90 and 180 • .
Although the horizontal phase velocity v x has often been observed larger than in this study, similar values have also been found.The northeast propagating TIDs in Australia in winter 1948 found by Munro (1950) had v x = 80-180 m s −1 .The TID observed by Thome (1964) in Arecibo (with similar T as in this paper) propagated with v x = 50-100 m s −1 .And the TIDs observed by Chum et al. (2012) in the Czech Republic (with comparable λ x and α x as in this paper) propagated at v x = 70-170 m s −1 .
The waves observed in the above studies are often labelled as "daytime MSTIDs" (e.g.Tsugawa et al., 2007).For the physical background of these waves, the general consensus in the above references is that the waves are generated in lower atmospheric layers, and travel upwards.The prevailing occurrence in the ionosphere of waves of certain period lengths, wavelengths, directions and phase velocities, is considered to be due to the filtering effect of the neutral wind in the mesosphere.Cowling et al. (1971) showed from an analysis using theoretical wind profiles how, depending on the wind, upward travelling waves of certain properties can either be reflected, experience critical coupling (be absorbed), or penetrate the mesosphere toward the ionosphere.The wind profile therefore acts as a filter of the waves.Unfortunately, the wind profiles used in the abovementioned and other theoretical studies are not representative enough of the true momentary wind at any time to provide accurate predictions for the expected TID characteristics.Generally however, Waldock and Jones (1984) showed that, since the mesospheric wind is generally directed antisunward and therefore rotates over the day (clockwise in the Northern Hemisphere), this directional filter rotates similarly.As a rule of thumb, Waldock and Jones (1986) showed that at any given time, waves propagating horizontally antiparallel to the direction of the mesospheric wind present at the same time and up to a few hours in the past, can be expected to be most prominent in the ionosphere.Table 1 shows that according to the TIE-GCM model, during the TID observed the mesospheric wind direction was around northward, agreeing with the general trend at local noon.This means that TID propagation directions between south and east are most likely at this time, which is consistent with the direction observed.

Conclusions
Two techniques to detect TIDs automatically from measured data have been presented: one to detect them from incoherent scatter radar data, and the other from GPS data.Both techniques are already known but have here been developed further, to derive more parameters, and to do this automatically.The techniques have been demonstrated on data on the same period of time: 20 January 2010.Both techniques observe from these data some TIDs around local noon, at approximately the same time despite the fact that the two sources of data are from slightly different locations (Tromsø in Norway, and an area around Sweden).
Each of the methods has its own strength.The incoherent scatter data method uses the vertical profile of TIDs, detecting the typical downward propagating phase fronts.It is able to determine the period length, vertical phase velocity and amplitude spectral density of the TIDs.The method using GPS data uses the horizontal 2-D pattern of the TIDs, and is able to determine the orientation and wavelength of the elongated structures, and their horizontal propagation direction and phase velocity, as well as the period length.
The fact that one method concentrates on the vertical and the other on the horizontal, makes the combination of the two methods a strong tool for the detection of TIDs and the determination of many of their parameters.In addition, this double technique makes it possible to double-check the period length as resulting from both methods.
Around 07:00-12:00 UT on 20 January 2010, a mediumscale TID was observed whose wave propagated southeastward with a horizontal phase velocity of 67 m s −1 and a wavelength of 500 km.The vertical phase velocity was found to be increasing with height and decreasing with period length.The amplitude of the wave was largest at 220 km altitude, which was also the location of the peak of ionisation at that time.The TID had its maximum amplitude in Tromsø at 10:00 UT, near local noon.It was detected over Sweden at the same time as in Tromsø (as far as can be told from the measurements).The GPS technique (which gives one period length at a time) observed a period length of this wave of 2 h.The EISCAT results showed a range of period lengths from 0.7 to 4.5 h, with the dominant period length at 1 h, indicating a higher harmonic of the same wave.
The horizontal wavelength and phase velocity were also estimated from the EISCAT results using Hines' dispersion equation.The results of this are larger than those detected from the GPS data.The most likely explanation for this is that Hines' theory overestimated the values, even in the WKB approximation to account for inhomogeneity of the atmosphere, since the conditions for validity of the WKB approximation were not well fulfilled in this case.
The TID observed in this study agrees in its parameters with those observed at many other high-and mid-latitude sites round midday in winter time.These TIDs were most likely generated at low atmospheric layers, and travelled upward.The concurrent mesospheric wind profile, which is typical for this time of day, acts as a directional filter, allowing in the ionosphere only waves in the southeast direction and only with certain wavelengths and phase velocities.The observed wave properties are therefore a result of this wind filter, rather than of the characteristics of the generating mechanism.

Possible future work
The techniques described in this paper, which allow a systematic characterisation of the parameters of TIDs, can be very useful to perform long-term analyses of TIDs using ISR and GPS data.This may be performed e.g. using the 1-year continuous EISCAT measurement campaign during the International Polar Year of March 2007-February 2008, and GPS data over the same period.Such a long-term analysis will enable us to verify the consistency of the observations and conclusions obtained in this paper, and establish statistical and season-dependent analyses of the wave parameters of TIDs.In addition, in this way relations between TID parameters and other parameters such as meteorological and space weather parameters can be studied.
The accuracy of the GPS results can be improved if the number of measurements is increased, e.g. by including other GNSS constellations (GLONASS, Galileo), or other GNSS receiver networks in Fennoscandia.The latter option may also increase the horizontal size of the analysis area.If that area becomes sufficiently large, this might also make it possible e.g. to split it into segments and analyse each segment separately, which would give more spatial information.If the technique is refined further, it may also be possible to detect the amplitude (as a function of time and location) from these data, as well as the full range of frequencies (period lengths).

Figure 3 .
Figure 3.The phase velocities for which a local maximum in correlation was found, for the data shown in Fig.2, at 10:00 and 14:00.At 10:00, a pattern showing an increasing vertical phase velocity with altitude was found, and a curve fitted to it is included in the graph.

Figure 4 .
Figure 4.The phase velocities as a function of altitude as in Fig.3for 09:00-12:00 and T = 1.05 h, and a curve fitted to all results together.
ah b (i.e. the difference on a logarithmic scale between the values of v z and the fitted curve) at every altitude, E[ ] = ensemble average over all altitudes, d M = 0.5.The maximum allowed value for |d h | is d M , i.e. every value of d h outside the interval [−d M , d M ], as well as every missing value of d h , is replaced by d M .

Figure 5 .
Figure 5.The quality factor Q resulting from the analysis of electron density over the entire day, as a function of the period length T and the hour of the day t t .

Figure 6 .
Figure 6.The slope b and offset a of the phase velocity profiles as found from the data of the entire day, as functions of the period length T .

Figure 7 .
Figure 7.A few examples of height profiles of vertical phase velocity v z as resulting from the data of electron density on 20 January 2010.

Figure 8 .
Figure 8. Electron density measured on 20 January 2010, bandpass filtered around 0.87 h (upper) and 1.52 h (lower), and the wave patterns detected from the data.

Figure 9 .
Figure 9.The amplitude spectral density A of the TID observed from 08:00 to 12:00, for T between 0.8 and 3.0 h, on 20 January 2010.Upper: as a function of altitude; lower: at h = 220 km as a function of time.

Figure 10 .
Figure 10.Map of the positions of the GPS receivers in the SWEPOS network, and of the EISCAT incoherent scatter radar in Tromsø.

Figure 11 .
Figure 11.Illustration of the projection of VTEC data on a line of a specified azimuth angle.The thick solid line has an azimuth ("α") of 120 • .All measured IPPs between the two dashed lines are projected onto the solid line, and all corresponding VTEC values in the same 20 km bin are averaged.

M
. van de Kamp et al.: TID characterised using joint effort of incoherent scatter radar and GPS

Figure 13 .
Figure13.The results of λ x and α x for the samples in the time window centred around 08:00 (blue), and their complex average (red).

Figure 14 .
Figure14.The horizontal orientation azimuths (upper) and the wavelengths (lower) of the spatial variations, as found from the analysis of this section.Small dots: the results for each 5 min sample.Stars: hourly averages.Circles: hourly averages which are considered unreliable (see text).

Figure 15 .
Figure 15.VTEC(d, t)  between 08:00 and 09:00, as functions of the distance d along the projection line (solid lines), and waves detected from these, i.e. cosines with wavelength λ x (in this case 500 km) and phases equal to arg(F ) (the amplitudes are fixed to 0.1 × 10 16 electrons m −2 ).

Figure 16 .
Figure 16.s p calculated from the values of arg(F ), and the fitted straight line, the gradient of which gives the resulting horizontal phase velocity v x .

Figure 17 .
Figure 17.The TEC derived from the GPS data in a 2-D horizontal grid, and the wave properties of the detected TIDs superimposed on the plots, at four times on 20 January 2010.The parameters of the detected wave are α x = 108 • , λ x = 500 km, v x = 67.2m s −1 .

Figure 18 .
Figure18.Horizontal wavelengths (upper graph) and horizontal phase speeds (lower graph).Lines: parameters calculated from the EISCAT results as functions of period length and altitude, using Hines' dispersion equation.Stars: parameters detected from the GPS results.Note that all values are with respect to a reference frame fixed to the Earth.

Table 1 .
The input parameters of Hines' dispersion equation on 20 January 2010, 10:00 UT, above Tromsø, for the specified heights.