Storm induced large scale TIDs observed in GPS derived TEC

Abstract. This work is a first statistical analysis of large scale traveling ionospheric disturbances (LSTID) in Europe using total electron content (TEC) data derived from GNSS measurements. The GNSS receiver network in Europe is dense enough to map the ionospheric perturbation TEC with high horizontal resolution. The derived perturbation TEC maps are analysed studying the effect of space weather events on the ionosphere over Europe. Equatorward propagating storm induced wave packets have been identified during several geomagnetic storms. Characteristic parameters such as velocity, wavelength and direction were estimated from the perturbation TEC maps. Showing a mean wavelength of 2000 km, a mean period of 59 min and a phase speed of 684 ms −1 in average, the perturbations are allocated to LSTID. The comparison to LSTID observed over Japan shows an equal wavelength but a considerably faster phase speed. This might be attributed to the differences in the distance to the auroral region or inclination/declination of the geomagnetic field lines. The observed correlation between the LSTID amplitudes and the Auroral Electrojet (AE) indicates that most of the wave like perturbations are exited by Joule heating. Particle precipitation effects could not be separated.


Introduction
The ionosphere, which can be characterized by its electron content, strongly depends on the influence of the sun.The high correlation between the variations of the solar radiation (EUV, energetic pearticle flux) and the variation of the total Correspondence to: C. Borries (claudia.borries@dlr.de)electron content (TEC) of the ionosphere has been demonstrated in many publications (e.g., Jakowski et al., 1991;Forbes et al., 2000).
When solar storm enhanced radiation and particle streams hit the earth's atmosphere, significant enhancement of the TEC and plasma redistribution processes were observed (e.g., Baran et al., 2001;Jakowski et al., 2006).The energy input during storm time is realized via coupling processes with the auroral magnetosphere and thermosphere.Any major step-like change (Förster and Jakowski, 2000) is supposed to induce atmospheric gravity waves (AGW) at high latitudes which play an important role in the equatorward redistribution of auroral momentum and energy input (Prölss and Jung, 1978).
The effects of these AGWs are visible in TEC as large scale travelling ionospheric disturbances (LSTID).This was already proposed in early TID studies, e.g. by Hines (1960) and Kersley and Hughes (1989).The category of "very large travelling ionospheric disturbances" was described in detail by Kersley and Hughes (1989).They characterized the LSTIDs by quasi-periods of 20 min to 2 h, speeds of 200 to 400 ms −1 and horizontal wavelengths exceeding 1000 km.Kersley and Hughes also found that auroral generation is efficient as a source for LSTIDs with a speed of more than 270 ms −1 propagating predominantly equatorwards.However, the phase speed of LSTIDs originating in the auroral region has been reported to cover a wide spectrum (Afraimovich et al., 2000).A comprehensive statistical analysis of LSTIDs was performed by Tsugawa et al. (2004), who characterized LSTID observed during 1999 to 2002 over Japan.
In this work ionospheric perturbation TEC maps over Europe are used to analyse LSTIDs induced by geomagnetic storms.The TEC maps are derived from GPS measurements following the procedure described by Saito et al. (1998).Dense data coverage is available over wide parts of Europe.The data is mapped into a regular grid with a high horizontal resolution.The extensions of the map enable to characterize waves of very large scales.The main part of this work is used to characterize the waves according to the occurrence, horizontal scale, speed and direction of the propagation.

GNSS satellites generally transmit signals on different frequencies.
Because the ionosphere is a dispersive medium, the ionospheric influence can be separated using the difference of the phases at two frequencies f (GPS: f(L1)=1575.42MHz, f(L2)=1227.60 MHz).The phase difference is an accurate relative measure of the slant TEC (TECs) along the ray path between satellite and receiver, but having an unknown shift regarding the absolute TECs.For a rough calibration, this shift is minimized using a TEC approximation with the empirical Neustrelitz TEC Model NTCM (Jakowski, 1998).To convert the measured TECs at elevation angle into an equivalent vertical TEC value (TECv) the following mapping function is used: Here, the ionosphere is approximated by a single layer located at the height h i =400 km.Respecting the electron content of the higher ionosphere and plasmasphere h i is chosen higher than the usual F2-layer height.R e represents the earth radius.The geo-reference of the resulting TECv is the piercing point of the ray path and the ionospheric layer at h i .The mapping is afflicted with a certain error in the estimation of the geo-reference of TECv, i.e. the accuracy of the conversion decreases with decreasing elevation.Thus, the accuracy of TECv has to be regarded with respect to its elevation.To reduce this transformation error, a cut-off angle of 30 • was chosen.
Finally, a perturbation TEC value can be deduced by computing the difference of the actual TECv from hourly means.This procedure, which was introduced by Saito et al. (1998), is often used in order to analyse meso scale travelling ionospheric disturbances (MSTIDs) and LSTIDs.
The resulting perturbation TEC is mapped into a regular grid covering Europe (10 • W-30 • E, 30 • N-70 • N, with a spatial resolution of 1 • ×1 • grid size).To estimate the perturbation amplitude at each grid point, a distance dependent weighted mean of the measurements is calculated.This procedure follows the same principle as used for constructing TEC maps described by Jakowski (1998).In addition to the distance dependent Gaussian type weight function defined for each measurement around a selected grid point, we take into account the elevation dependent accuracy of TECv estimations by introducing a Gaussian type elevation dependent weight function.The half widths of both weight functions were fixed at 2 • for the distance from grid points and for the elevation weight functions at 15 • from the zenith.
Based on 30 s sampled GPS data, the perturbation TEC maps are created with a temporal resolution of one minute.This is sufficient for analysing storm induced LSTIDs characterized by periods of about one hour.
The European GNSS data is provided by the EUREF Permanent GNSS Network1 (Bruyninx et al., 2001).The number of GNSS ground stations steadily increased over the years.While there were only about 150 stations per day available on average in 2001, there are more than 400 stations available today.As a measure for the data coverage the average sum of weights is calculated for each grid point of the perturbation TEC map.An example is shown in Fig. 1 for the 23 May 2002.It illustrates that the best data coverage is achieved over central Europe, where the grid points are even overestimated as the sum of weights is bigger than 1 (overestimated map points are encapsulated by a black line).While on 23 May 2002 the percentage of the area of overestimation with respect to the area of the whole map is 7.5% it is already 17.5% on 17 December 2007.This underlines the permanent improvement of the data base.

Selection of geomagnetic storm events
Geomagnetic storms occur when there is a large sudden change in the solar wind dynamic pressure at the magnetopause (Schunk and Nagy, 2000).The geomagnetic D st index is an excellent indicator of storm events.The main attribute of a magnetic storm is a clear decrease of the horizontal intensity of the magnetic field (e.g., Prölss, 2001).The onset phase of a storm is often characterized by a short sudden increase of the D st index (Förster and Jakowski, 2000).
In this paper a rapid D st increase in the onset phase of a geomagnetic storm is used to define the beginning of the storm.A well pronounced negative D st period (at least two days) afterwards characterizes the main and subsequently the recovery phase of a storm event.
The selection and timing of the storm events is important in order to develop reliable storm statistics.The time interval from 2001 to 2007 has been chosen to be analysed for geomagnetic storm events over Europe.In the time before 2001 the data coverage is supposed to be too sparse to provide information about the characteristics of the storm induced waves.
During 2002 to 2004 various stronger and weaker storms occurred.Against this, 2006 and 2007 are years in the solar minimum, thus hardly had strong storm events.A list with the dates of the selected storms can be found in the Appendix A in Table A1.This method may not respect all the ionospheric storms that occurred during 2001 to 2007, but most of them.

Analyses and discussion
Examples of perturbation TEC maps obtained with the method described in Sect. 2 are shown in Fig. 2. The wave crests and minima are marked with solid and dashed lines.The wave front of a LSTID can be seen moving from North to South.For estimating the characteristic scale length and propagation direction, the two-dimensional Fourier analysis (Eq. 2) is applied on each map during the storm events. (2) x is the location vector representing the grid point coordinates in latitude and longitude and k the wave vector.
The wave vector keeps the information about the absolute wavelength λ and direction ϕ (defined as angle from North clockwise) of the wave train.A demonstration of the results of the two-dimensional Fourier analysis on the example map of the 23 May 2002 12:22 UT (Fig. 2, lower panel) is shown in Fig. 3.The shading indicates the amplitude of the waves.Figure 3 demonstrates the occurrence of a wave with 2000 km wavelength and southward direction (ϕ=180 • ).Its amplitude is 1 TECU.Typical properties of the storm induced TIDs are estimated by analysing all perturbation maps during the first 12 h after the selected storm events.The wavelength λ and direction ϕ of the strongest wave of each map are plotted as single dot into a scatter plot as that demonstrated in Fig. 4.Only waves with amplitudes stronger than 0.5 TECU are considered.Not each storm excited such strong oscillations.Especially during the solar minimum (2004)(2005)(2006)(2007) most LSTIDs did not exceed amplitudes of 0.5 TECU.The number of wave occurrence per 10 • ×100 km square of the plot is indicated by grey shading.The highest concentration of points denotes the preferential propagation direction and wavelength of the storm induced LSTID.The propagation directions of the LSTID estimated with this method are listed the the appendix Table A1.
Figure 4 joins the results of all storms.The typical LSTID can be found with a wavelength around 2000 km and southward propagation direction ϕ=180 • .205 maps revealed waves with these properties.It seems to be the typical structure of storm induced TIDs at the European region.
However, the points in Fig. 4  direction.Nevertheless, the number of occurrence illustrates a typical equatorward propagation of the LSTIDs over Europe.It might be assumed that the LSTIDs propagate along the magnetic field lines for which in Europe the declination is very small referring to the geographic meridians (e.g Jakowski et al., 2008).
The LSTIDs are supposed to be generated at high latitudes, where the main storm energy is deposited.Schunk and Nagy (2000) describe an expansion of the plasma convection and particle precipitation pattern in the growth phase of a geomagnetic storm.These changes are accompanied by substantial increases in the Joule and particle heating rates.The Joule Heating in the auroral region is indirectly measured by the auroral electrojet (AE) index.The mean AE in the first 3 h after the selected storm times is listed in Table A1.A sophisticated correlation study of AE and the maximum amplitude of the observed LSTIDs is shown in Fig. 5.During solar maximum (Fig. 5, left panel) the LSTID amplitudes observed during day time are well correlated to the AE, with a correlation coefficient of 0.8 (95% significant).The LSTID amplitudes grow up with increasing AE-index.During night the correlation to AE is low.There seem to be competative dynamic processes interacting.Besides, inaccuracy of the analyses method may also cause some outliers.Further investigation is needed to clarify the processes inducing LSTIDs during night.During solar minimum (Fig. 5, right panel) there are not enough measurements up to now to derive significant information about a correlation between AE and wave amplitudes.It seems that there is a similar correlation but with generally smaller amplitudes of the LSTIDs due to the weaker ionisation during solar minimum.Hence, it is concluded that during day time in solar maximum the observed LSTIDs are excited by a sudden heating of the thermosphere which leads to strong changes in the pressure and thermospheric wind circulation systems.These rapid changes may lead to the generation of AGWs.The strength of the wind and pressure changes is reflected in the amplitude of the AGWs which are observed in the ionosphere as LSTIDs.This mechanism initially needs equatorward winds that uplift the ionospheric plasma along the geomagnetic field lines and cause equatorward propagating LSTIDs.
The observations of an uplifted F2-layer during geomagnetic storms support the preference of this mechanism.This is demonstrated in Fig. 6 with the help of electron density profiles which are derived by radio occultation measurements  onboard the mini satellite CHAMP2 .Although the geographic longitude of the measurements do not correspond to the map extensions, the same effect are supposed to be seen, because of the large longitudinal extension of the observed LSTIDs.The electron density profile in Fig. 6, which was measured during the severe geomagnetic storm on the 29 October 2003 shows an uplifting of the height of maximum electron density (about 75 km) compared to the mean profile of 29 September to 27 October 2003.Strong neutral background winds, induced by the thermospheric heating or electric fields of magnetospheric origin, are supposed to cause this uplifting.Other examples for an uplifted F2layer during geomagnetic storms are reported e.g. by Prölss andJung (1978) andJin et al. (2008).However, strong northward winds, which were described by Shiokawa et al. (2002) and Shiokawa et al. (2003) analysing the geomagnetic storm on the 31 March 2001 conflict with this idea indicating that more detailed investigation is necessary.The 2-D-Fourier analyses showed that the LSTID wave fronts induced by ionospheric storms mainly propagate from North to South in Europe.This is consistent with other studies, e.g.Kersley and Hughes (1989) and Tsugawa et al. (2004).Under this precondition Time-Latitude-Plots (TLP) are efficiently used to trace these wave events and estimate their phase speed.The error due to the 10 • deviation of the propagation direction of some LSTID to the meridional direction is 1.5% (below 20 ms −1 ).Thus, it can be neglected.
Because of the high data density over Central Europe (see Fig. 1) we chose the meridians at 8, 10, 12, 15 and 17 • E to prepare the TLPs.The TLPs start one hour before the storm onset and end 12 h later.
As an example, Fig. 7 shows a TLP demonstrating the propagation of TIDs at 12 • E for the 24 August 2005.Black and white shadings indicate positive and negative amplitudes, respectively.In this case 6 tilted signatures of positive amplitudes can be identified.They are indicated with arrows in Fig. 7.The slope of the amplitude pattern indicates the southward propagation of strong TEC perturbations.It can be noticed that the wave activity starts soon after the start of the storm at 6 a.m.The superposition of the TLPs of all storms (not shown here) revealed that the mean absolute perturbation TEC suddenly increases about 1 h after the defined storm start time, probably related to electric field induced uplifting (cf.Arbesser-Rastburg and Jakowski, 2007).Neutral gas composition changes, horizontal transport processes due to winds and waves are assumed to start immediately after this onset as discussed by Jakowski et al. (1999); Förster and Jakowski (2000).The maximum amplitude appears typically between 2 and 6 h after the storm onset.In Fig. 7 the maximum amplitude is about 6 h delayed.
Most of the TLPs of the selected storms clearly reveal similar wave fronts as shown in Fig. 7 propagating from North to South.Some reveal strong irregular disturbances without a clear structure, or a mixture of stochastic and structured perturbations (wave fronts).Stochastic perturbations mainly occur north of 60 • N, as is visible in Fig. 7.They are most frequent at severe storms and indicate the generation region of the LSTIDs.
The amplitudes of the waves vary from event to event.In most cases they exceed 0.5 TECU.Strong events like those observed on the 29 October 2003 reveal amplitudes exceeding 1.5 TECU.
The slope of the parallel signatures in the TLP indicates the phase speed of the wave fronts.To estimate the phase speed based on TLPs the Radon transform is very efficient.The Radon transform is the projection of the TLP amplitudes along a radius vector oriented at a specific angle α to the x axis and having a distance d from the graphs origin.This straight line is defined parametrically by x(t)=t sin α+d cos α, y(t)=−t cos α+d sin α.The Radon transform of a function f on the plane is defined by Because the TLP is not an infinite function like f but an array with discrete values, a discrete version of the Radon transform R TLP (α m , d n ) is employed.The origin of the coordinate axes is in the center of the TLP.Each value in the array R TLP (α m , d n ) is squared and summed over n to give a discrete measure of |R TLP | 2 dα.Due to the fact that the TLP is a rectangular array and thus does not have the same number of pixels in every direction α m , it is necessary to divide n R TLP by a normalized correction factor C(α n ).C(α n ) can be calculated using the Radon transform on a plain image of the same size as the TLP.The result of this method shows in which direction α m the highest energy of the TLP can be found.Finally, α m needs to be converted into the phase speed.
Applying the Radon transform on all TLPs, a mean phase speed of 684 ms −1 (mean of 34 storms) with a standard deviation of 195 ms −1 has been estimated.The phase speed could not reliably be estimated in all cases because of an unstructured disturbed ionosphere.All results of the phase speed analyses are displayed in Table A1 in the Appendix A. Compared to the phase speed observed by Tsugawa et al. (2004) during geomagnetically disturbed days over Japan (487±145 ms −1 for damped LSTID, 561±171 ms −1 for growing LSTID) the phase speed of the LSTIDs over Europe seems to be high.To explain this difference, miscellaneous aspects have to be considered.On the one hand the regional conditions differ quite a lot.Japan is situated in midlatitudes, whereas Europe covers middle and high latitudes including parts of the auroral region.On the other hand the inclination and declination of the geomagnetic field lines differ due to the latitudinal and longitudinal difference between Japan and Europe.While the declination of the geomagnetic field lines is in Europe approximately zero, it is approximately −7 • in Japan.A sophisticated examination of the phase speed in southern Europe would be helpful to discuss the differing phase speeds in Europe and Japan.Unfortunately, the analysis of the smaller map subsets reach the limits of the reliability of the analysis method.
Finally, we want to discuss the periodicity of the wave trains.Time series from single map grid points are analysed according to the wave periodicities.The wavelet analysis is applied respecting the possibility that the period of the oscillation may change.A modified version of the continuous wavelet transform (Eq. 3) provides additional information about the amplitude of the oscillations.
The wavelet transform is the convolution of a signal f with a translated and scaled version of the mother wavelet ψ 0 .The parameter τ describes the translation and the parameter s the scale.The morlet wavelet ψ 0 (η)=π −0.25 e iω 0 η e −η 2 /2 with a center frequency ω 0 =6 is applied for the mother wavelet.
In order to gain a representative image of the oscillations found in the perturbation TEC the data of 9 different map grid points have been chosen for the time series analyses.All points are part of the area of overestimation (see Fig. 1).The average and standard deviation of the observed periods at each storm event can be found in Table A1 in the Appendix A. The mean period of the observed oscillations is 59 min with a standard deviation of 11 min.Compared to the periodicities of 80±29 min observed over Japan (Tsugawa et al., 2004) these periods are low.But this stands in good agreement with the higher phase speeds which were retrieved from the TLPs.

Summary and conclusions
Storm induced LSTIDs over Europe have been analysed in the variation of GNSS derived TEC.The storm start time was defined by the time of the most rapid increase of the D stindex during the onset phase of the geomagnetic storm.A new data base of perturbation TEC maps has proven as an efficient tool for the observation of LSTIDs over Europe.
The application of different analyses, e.g. the Fourier-and wavelet transform, enabled to identify typical properties of the storm induced LSTIDs.These waves are mainly propagating equatorwards.Because plasma motion is coupled to the geomagnetic field lines storm induced LSTID may propagate preferably along geomagnetic field lines, which have only a very small declination in the European region.
It was found that the LSTIDs are characterized by a mean wavelength of 2000 km, a mean period of 59 min and a mean phase speed of 684 ms −1 .The observed LSTIDs have the same horizontal scale as LSTIDs observed over Japan.However, the mean phase speed is faster than those found during geomagnetically disturbed days over Japan.Damping and differences in the geographic/geomagnetic relationship are suggested be the reason for this disparity.
Due to a clear correlation of the amplitudes of the LSTID and the AE-index, the examined LSTIDs are assumed to be mainly excited by Joule heating in the Auroral region.

Fig. 1 .
Fig. 1.Map of the average sum of all affecting weights at each grid point of the TEC perturbation map on the 23 May 2002.The black line encapsulates the areas where the summed weights are bigger than 1.The triangles indicate the position of the available GNSS ground stations.

Fig. 5 .
Fig. 5. Correlation of the auroral electrojet (AE) index with the maximum amplitude of the wave observed during the first 12 h after the storm onset.Left: solar maximum (2001-2004); Right: solar minimum (2005-2007).Night time (21:00-03:00 UT) storms are marked gray.The dashed lines represent the least squares fit, excluding night time storms.