- Email: [email protected]
Empirical model of the TEC response to the geomagnetic activity over the North American region
Available online at www.sciencedirect.com
Advances in Space Research 48 (2011) 1041–1048 www.elsevier.com/locate/asr
Empirical model of the TEC response to the geomagnetic activity over the North American region B. Andonov, P. Mukhtarov, D. Pancheva ⇑ Geophysical Institute, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Block 3, 1113 Sofia, Bulgaria Received 3 March 2011; received in revised form 2 May 2011; accepted 4 May 2011 Available online 12 May 2011
Abstract The paper presents an empirical model of the total electron content (TEC) response to the geomagnetic activity described by the Kpindex. The model is built on the basis of TEC measurements covering the region of North America (50°W–150°W, 10°N–60°N) for the period of time between October 2004 and December 2009. By using a 2D (latitude-time) cross-correlation analysis it is found that the ionospheric response to the geomagnetic activity over the considered geographic region and at low solar activity revealed both positive and negative phases of response. The both phases of the ionospheric response have different duration and time delay with respect to the geomagnetic storm. It was found that these two parameters of the ionospheric response depend on the season and geographical latitude. The presence of two phases, positive and negative, of the ionospheric response imposed the implementation of two different time delay constants in order to properly describe the two different delayed reactions. The seasonal dependence of the TEC response to geomagnetic storms is characterized by predominantly positive response in winter with a short (usually 5–6 h) time delay as well as mainly negative response in summer with a long (larger than 15 h) time delay. While the TEC response in March and October is more close to the winter one the response in April and September is similar to the summer one. Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Total electron content; Geomagnetic storm; Time delay; Positive and negative response
1. Introduction Solar activity such as flares and coronal mass ejections often produce large variations in the particle and electromagnetic radiation incident upon the Earth. Such variations can lead to a significant perturbation of the “quiet-time” ionosphere due to large variability in the ionospheric density distribution, total electron content (TEC), and the ionospheric current system. The geomagnetic storms have important terrestrial consequences such as disrupting satellite communications or detection and tracking of aircraft, distort navigation and rarely an interruption of the flow of electrical energy over power grids. Thus, it is important to monitor ionosphere particularly during the geomagnetic storms, to study its response to them, and, if ⇑ Corresponding author. Tel.: +359 2 9793308; fax: +359 2 9713005.
possible, to forecast the evolution of the ionospheric variability. The development of the Global Navigation Satellite System (GNSS) during the last decade has provided a number of possibilities for studying the spatial distribution and temporal variability of ionospheric electron density disturbances forced by external or internal sources. The measure of the phase delay of the propagating through the ionosphere high frequency (HF) radio waves is a key parameter for defining the total electron content (TEC) between the satellite and the receiver (Liu et al., 1996). The accurate measurement of the phase delay is a strong requirement for the reliable performance of many applications including HF communications, satellite positioning and navigation applications. The necessity of the data corrections obtained from the navigation satellites imposes the requirement not only the regular variability of the electron density, but also
E-mail address: [email protected] (D. Pancheva). 0273-1177/$36.00 Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2011.05.007
1042
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
the variability related to the geomagnetic storms or other drivers to be taken in mind (Jakowski et al., 2005). Before sufficiently long TEC time series to be accumulated by using the navigation satellites the problem for the TEC prediction had to be solved on the basis of the existing electron density profile models (Radicella and Zhang, 1995). The profiles were constructed by using the main ionospheric parameters determined by the ground based ionosonde stations. The response of the main ionospheric characteristic foF2 (known as the critical frequency of the ionospheric F-region or respectively) to the geomagnetic activity and its dependence on the geographic/geomagnetic coordinates, season and local time (LT) has been investigated by Kutiev and Muhtarov (2003). This empirical model is based on the all available ionosonde measurements. The foF2 prediction according to the level of the geomagnetic activity provides a possibility for an approximate assessment of the TEC variability. This approach however has a significant disadvantage: the ionosonde measurements cannot provide any information about the electron density profile above the F-region maximum. The upper part of the electron density profile however contributes to the TEC and cannot be ignored. The recent investigations of the TEC response to the geomagnetic storms based on the measurements directly obtained from the navigation satellites have reported the main features of the TEC response. Zhao et al. (2007) presented the latitude/longitude distribution of the relative deviation of the TEC (rTEC) during winter, summer and equinoxes, as well as the effect of the LT. Wang et al. (2008) studied the correlation dependences between the low latitude TEC and the equatorial Dst-index. Habarulema et al. (2007) suggested a TEC prediction method by means of the Neural Network (NN) based model. One of the input parameters of this method is the Kp-index. This study shows also the difference between the TEC defined by the navigation satellites and that defined by ionosonde measurements (i.e. only up the F-region maximum). The results indicated that usually the TEC defined by the ionosonde profile is on the average 75% from the entire TEC and that the variability of the both types of TEC data can be significantly different particularly for the considered in the paper examples. Afraimovich et al. (2009) investigated the dependence of the TEC variability on the Kp-index for equatorial, middle and high latitudes in the Northern Hemisphere (NH). The authors obtained linear regressions of the dependence of the relative TEC deviation (rTEC) on the geomagnetic index for the considered latitude range. Only the positive dependences between the rTEC and Kp-index are found; the response at high latitudes (50°N–80°N) is stronger than that at low latitudes. The authors found also significant longitudinal differences (geographic frame is used in the above mentioned study). Liu et al. (2010) studied the time delay of the TEC response by using the cross-correlation analysis. Its global distribution is presented separately for the positive and negative phases depending on the season.
The obtained values for the time delay are between 2 and 17 h, i.e. a range close to that found by Kutiev and Muhtarov (2001), but on the basis of ionosonde data. Stankov et al. (2001) and later Stankov et al. (2004) presented a correlation method for the prediction of TEC depending on the geomagnetic activity for the geographic region (20°W–40°E, 32.5°N–70°N). This method is based on the ideas used for the foF2 prediction and reported by Muhtarov et al. (2001) and Kutiev and Muhtarov (2001). In the above mentioned correlation method a time delay constant of 18 h and a periodic Fourier function accounted for the LT effect were introduced. Very recently Stankov et al. (2010) have applied a common epoch analysis on data representing nearly 300 storm events from the last solar cycle. It was found that the storm-time behavior of TEC shows clear positive and negative phases with amplitudes that tend to increase during the more intense storms. The most pronounced positive phase is observed during winter, while the strongest and yet shortest negative phase is detected during equinox. The TEC response to the geomagnetic activity for the American sector was reported by Araujo-Pradere et al. (2006). The authors found consistent features from storm to storm were found, and these feature became more apparent when the data were separated between the “driven” phase of the storm, when the integral of ap-index is rising, and the “recovery” to the storm, when the integral of ap-index is declining. The existence of the positive phase over the continental USA at the beginning of the storm period appears to be related to the timing of the peak of the perturbation; a positive phase will be observed when the peak of the perturbation occurs near midnight UT. The basic aim of the present paper is to find a mathematical dependence between the most probable values of the rTEC and the geomagnetic activity described by the Kp-index as a function of the calendar month, geographic latitude and LT. The wanted empirical dependence is obtained by using a long time series of TEC measurements over the region covering the North America. A particular attention will be paid to both the inertness of the ionospheric processes and the phase (positive/negative) of the TEC response. A significant improvement of the present study is the inclusion of both positive and negative phases of the TEC response to the geomagnetic storms. While the used TEC measurements (2004–2009) belong mainly to moderate and partly extremely low solar activity conditions, the latter being responsible for anomalous conditions at least in the thermosphere (Emmert et al., 2010). Therefore, the TEC model presented here is valid mainly for low solar conditions. The solution of the above described task will improve the forecast of the evolution of ionospheric irregularities during the geomagnetic storms. 2. Data The time series of vertical TEC (VTEC) for the period of time between October 2004 and December 2009 are used in
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
this study. The VTEC data are hourly data and cover the following geographic region: longitudes from 150°W to 50°W with longitude resolution of 1° and latitudes from 10°N to 60°N with a latitude resolution also of 1°. The VTEC data were downloaded from the NOAA National Geophysical Data Center (NGDC) web site: http:// www.ngdc.noaa.gov/stp/iono/ustec/index.html. The geomagnetic activity is defined in this study by the global Kp-index as it describes well the geomagnetic activity from tropical to middle/high middle latitudes. The Kpindex data are downloaded from the Space Physics Interactive Data Resource (SPIDR), Boulder, Colorado for the considered period of time. In the model the VTEC response to the geomagnetic activity is investigated by the relative deviation of the VTEC defined as: rVTEC = (VTECobs VTECmed)/ VTECmed. In this way the effect of the regular seasonal, diurnal and solar changes is removed from the VTEC variability. The term VTECmed represents the monthly median value. Data are grouped into 12-month bins, as every bin contains all the available hourly data within the respective month of the year.
January
60
0.4 0.36 0.32
50
Lat [deg]
0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
40 35 30 25 20 15 -8
0
8
16
24
32
40
48
56
64
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20
April
-16
-8
0
8
16
24
32
40
48
64
56
45 40 35 30 25 20
45 40 35 30 25 20 15
10 -24
-16
-8
0
8
16
24
32
40
48
56
64
10 -24
72
July
60
Lat [deg]
45 40 35 30 25 20
20
16
24
32
40
48
56
64
72
October
-8
0
8
16
24
32
40
48
56
64
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
45 40 35 30 25 20
45 40 35 30 25 20
-16
-8
0
8
16
24
32
40
48
56
64
72
November
40 35 30 25 20
10 -24
-16
-8
0
8
16
24
32
Time lag [hours]
40
48
56
64
72
64
72
20
-16
-8
0
8
16
24
32
40
48
56
64
72
September
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20
-16
-8
0
8
16
24
32
40
48
56
64
72
December
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
45
10 -24
56
25
60
50
15
48
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
40
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
60
55
55
32
72
August
60
50
24
June
15
10 -24
8
16
30
15
0
8
40
10 -24
-8
0
35
10 -16
-8
45
15 -24
-16
50
15
60
Lat [deg]
-16
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
30 25
10 -24
60
55
35
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
40
60
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
45
10 -24
72
May
60
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
10 -24
72
March
60
15
-16
60
Lat [deg]
We use the planetary Kp-index as an indicator of the geomagnetic activity in this study. This index reflects both types of variability: from the equatorial ring current and the auroral currents. Another important reason for using this index is connected with the fact that there are online models for its prediction. It is well known that during the recovery phase of the ionospheric storms with geomagnetic origin the ionospheric reaction continues some time after the geomagnetic storm attenuation. This phenomenon aggravates the relationship between the Kp-index and the ionospheric anomalies. Later a method will be suggested for avoiding of this shortcoming which is based on the integration of the Kp-index. The effect of geomagnetic activity, described by the Kp-index, on the rVTEC is studied by cross-correlation analysis between the both parameters. A two-dimensional (time-latitude) cross-correlation function is calculated for each month of the year (from the 12month bins) because its seasonal dependence is well known. The calculated cross-correlations functions for all months are shown in Fig. 1. The most important result is that for
55
45
10 -24
3. Cross-correlation analysis between rVTEC and Kp-index
February
60
55
1043
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20 15
-16
-8
0
8
16
24
32
Time lag [hours]
40
48
56
64
72
10 -24
-16
-8
0
8
16
24
32
40
48
56
64
72
Time lag [hours]
Fig. 1. Two-dimensional (latitude-time) cross-correlation functions calculated between the rVTEC and Kp-index for all months of the year; the zero line is shown by dashed white line.
1044
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 0.80 Polynomial Coefficients Degree 0: -6.00067 Degree 1: 3.07795 Degree 2: -0.524284 Degree 3: 0.0298974 Data
rVTEC
0.60
Best Fit 0.40
0.20
0.00 4.0
5.0
6.0
7.0
8.0
9.0
Kp
Fig. 2. The empirical dependence between the Kp and rVTEC.
1.6
2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8
November 2004
December 2005
1.4
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
1.2 1
rVTEC
rVTEC
all months the ionospheric response is composed by two phases, positive and negative, with different duration and time delay. The results from the figure clearly indicate that the cross-correlation depends on both latitude and time lag (which reflects the delayed reaction of rVTEC with respect to Kp changes) and this dependence is different for each month of the year. The cross-correlation function is positive (reaching up to 40% for December) with time lag not longer than 10 h for tropical (10°N–30°N) and high (higher
than 50°N) latitudes, while the cross-correlation function for middle (30°N–50°N) latitudes is negative (up to 30% for April) with time lag up to 6 times longer than that for tropical and high latitudes. Hence, during winter (from November to February) the cross-correlation function for tropical and high latitudes is characterized by high positive values while for middle latitudes it is predominantly negative and with time lag longer than 40 h. For the spring months (March–April) the cross-correlation function has both high positive and increasing negative values but the latter are reached for a shorter than in winter time lag of 20 h. The cross-correlation function for May and June indicates a decrease of the positive values while its negative behavior is similar to that in spring. The cross-correlation function for the months from July to October is similar to that for March–April with an exception for August when the negative cross-correlation is small and confined only to latitudes near 40°N. We note that the result for August could be affected by the decreased quality of the VTEC data for this month of the entire period of time (2004–2009). In summary, it can be concluded that the crosscorrelation function: (i) for almost all months of the year at low and high latitudes is mainly positive with a time lag of up to 10 h, (ii) for winter it has predominantly positive values which move toward the middle latitudes, and
0.8 0.6 0.4 0.2 0 -0.2 -0.4
06
07
08
09
10
11
12
12
13
14
Day
16
17
18
19
Day
9
9
November 2004 8
December 2005
8
KpTs KpTl Kp Data
7
KpTs KpTl Kp Data
7
6
6
5
5
Kp
Kp
15
4
4
3
3
2
2
1
1
0
0
06
07
08
09
Day
10
11
12
12
13
14
15
16
17
18
19
Day
Fig. 3. (Upper row of plots) Comparison between the observed (dash line) and model (solid line) rVTEC at the geographic point (20°N, 130°W) for the geomagnetic storm in November 2004 (left plot and December 2005 (right plot); (bottom row of plots) comparison between the measured Kp-index (tick solid line with dots) and modified K p T s - (dash line) and K p T l - (solid line) indices for November 2004 (left plot) and December 2005 (right plot). The time delay constants for November 2004 geomagnetic storms are: Ts = 5 h; Tl = 15 h, while those for December 2005 are: Ts = 6 h; Tl = 48 h.
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
(iii) for equinoxes and summer it has mainly negative values at middle latitudes with time lag of 20 h; the negative cross-correlation values in equinoxes are larger than those in the summer. The above mentioned results clearly indicated that two different time constants, Ts and Tl, should be introduced in the empirical model, representing the different time delays of the rTEC with respect to the geomagnetic activity described by the Kp-index. 4. Empirical model The investigations on the foF2 response to the geomagnetic activity (Muhtarov et al., 2001; Wang et al., 2008) indicated that this is a delayed response which can be satisfactorily modeled by assuming that the geomagnetic influence is imposed on the inertial system described by an inhomogeneous differential equation from a first order (Kutiev and Muhtarov, 2003). These studies found that the time delay constant is near 18 h. The above presented cross-correlation analysis however indicated that in this case the TEC reaction to the geomagnetic activity has to be presented by a sum of two responses with different time delay constants and with different sign of the cross-correlation function. A possible reason for the observed difference between the geomagnetic influence on the foF2 and TEC is probably related to a difference in the altitude structure of the electron density profile response.
If we assume that the impact of the geomagnetic activity on the TEC is accomplished be two mechanisms with different time delay constants then the variability of rVTEC can be described as follows: rVTECðtÞ ðfT s ðK p T s ðtÞÞ þ fT l ðK p T l ðtÞÞÞf ðLTÞ
ð1Þ
where the function f(LT) represents the dependence of the response on the LT at equal other conditions. The parameters K p T s and K p T l are the modified with time delay constants respectively Ts and Tl values of the Kp-index. These modified parameters are solutions of the equations shown below and are obtained easily by a numerical integration:
Ts Tl
dK p T s ðtÞ þ K p T s ðtÞ ¼ K p ðtÞ dt dK p T l ðtÞ dt
ð2Þ
þ K p T l ðtÞ ¼ K p ðtÞ
ð3Þ
It is convenient the unknown functions fT s and fT l from (1) to be expressed by their Taylor time series expansions while the dependence on the LT to be presented by a Fourier time series as follows:
0.6
0.6
April 2008
April 2005
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
0.4
0.4
0.2
rVTEC
rVTEC
1045
0
0.2
0 -0.2
-0.4
-0.2
01
02
03
04
05
06
07
08
09
01
02
03
04
Day
06
07
08
09
10
Day
8
8 April 2005
7
April 2008 7
KpTs KpTl Kp Data
6
KpTs KpTl Kp Data
6
5
5 Kp
Kp
05
4
4
3
3
2
2
1
1
0
0 01
02
03
04
05
Day
06
07
08
09
01
02
03
04
05
06
07
08
09
10
Day
Fig. 4. The same as Fig. 3 but for the geomagnetic storms in April 2005 (left column of plots) and April 2008 (right column of plots); the time delay constants for both geomagnetic storms in April 2005 and 2008 are: Ts = 6 h; Tl = 36 h.
1046
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 2
3
2
3
fT s ðK p T s Þ ¼ a0s þ a1s K p T s ðtÞ þ a2s K p T s ðtÞ þ a3s K p T s ðtÞ þ fTl ðK p T l Þ ¼ a0l þ a1l K p T l ðtÞ þ a2l K p T l ðtÞ þ a3l K p T l ðtÞ P flt ðLTÞ ¼ b0 þ bi cosði 2p LT ui Þ 24 i
ð4Þ As the basic aim of this study is to build a model which presents the rVTEC changes trough the geomagnetic activity described by the Kp-index, then it means that a functional dependence between the Kp and rVTEC should be found. The most appropriate type of the wanted functional dependence can be defined empirically. For this purpose the mathematical expectations of the rVTEC are calculated at a preliminary grouping of the Kp values into subintervals with a given length. The range of Kp values (the min Kp is the most left value while the max Kp is the most right value on the axis) is divided into subintervals with equal length which are arranged in an ascending line. The mathematical expectation and the standard deviations of the rVTEC are calculated for each subinterval. The empirical dependence between the Kp and rVTEC is presented in Fig. 2. It is evident that the functional dependence between Kp and rVTEC is close to the cubic function. As the dependence between the Kp and rVTEC shown in Fig. 2 is close to the cubic function then in the Taylor expansions (4) only four terms have to be included. The next step is to obtain the most probable values of the
coefficients: ais, bi, ail, Ts and Tl from (4). This is a nonlinear optimizing task that can be solved by applying the “trial-and-error” method in a way that the best approximation in a sense of minimum least squares deviation has to be assured. In order to solve the problem the following steps are made: (i) it is given a range of Ts changes from 0 to 10 h with a time resolution of 1 h and a range of Tl changes from 11 to 60 h with a time resolution again of 1 h; (ii) for each grid point of the built in this way grid the coefficients ais, bi, ail are found by using the method of least squares best fit, and (iii) the coefficients ais, bi, ail, Ts and Tl at which the best approximation (in a sense of minimum least squares deviation) is obtained are accepted as optimal coefficients for the model. 5. Model results The upper row of plots of Fig. 3 show the observed (dash line) and model (solid line) variability of the rVTEC during two strong geomagnetic storms taking place in November 7–9, 2004 (left plot) and December 14–16, 2005 (right plot) for a geographical point (20°N, 130°W). These geomagnetic storms are characterized by a predominantly positive response of the rVTEC. The bottom row of plots of Fig. 3 show the comparison between the Kp-index and the determined modified Kp indices with the respective time constants Ts and Tl. It is evident from the figure that 0.8
0.6
July 2006
July 2006
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
0.6
0.4
rVTEC
rVTEC
0.4 0.2
0.2 0
0 -0.2 -0.2
-0.4
01
02
03
04
05
06
07
25
08
26
27
28
Day
30
31
01
29
30
31
01
8
8
July 2006 7
July 2006
7
KpTs KpTl Kp Data
6
KpTs KpTl Kp Data
6 5
Kp
5
Kp
29
Day
4
4
3
3
2
2
1
1 0
0
01
02
03
04
05
Day
06
07
08
25
26
27
28
Day
Fig. 5. The same as Fig. 3 but at the geographic point (40°N, 130°W) for the geomagnetic storms in beginning of July 2006 (left column of plots) and at the end of July 2006 (right column of plots); the time delay constants for both geomagnetic storms in July 2005 are: Ts = 10 h; Tl = 20 h.
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 0.4
October 2008 rVTEC Data rVTEC Model
rVTEC
0.2
0
-0.2
10
11
12
13
14
15
16
15
16
Day 8
October 2008
7
KpTs KpTl Kp Data
6
Kp
5 4 3 2 1 0
10
11
12
13
14
Day
Fig. 6. The same as Fig. 3 but for the geomagnetic storm in October 2008; the time delay constants for the considered geomagnetic storm are: Ts = 10 h; Tl = 36 h.
the maxima of the considered geomagnetic storms are near midnight for both storms: on November 8 in the first case and on December 15 in the second case. The maximum of the ionospheric response is detected 6 h later in both cases: at 6 UT on November 8 in the first storm and at 6 UT on December 15 for the second storm. A similar result from the cross correlation analysis was shown in Fig. 1 for November and December. From the presented cross correlation analysis shown in Fig. 1 it was found that the rVTEC response to the geomagnetic activity can be both positive and negative but with different time delay constants. It turns out to be that the ionospheric response to one and the same geomagnetic storm can have both positive and negative phases. Such cases are displayed in Fig. 4. The upper left plot of Fig. 4 shows the comparison between the observed (dash line) and model (solid line) rVTEC for the period of time between 1 and 9 of April 2005. The geomagnetic storm starts near noon time on April 4. A positive response is seen during the first 24 h of the geomagnetic storm and a negative one in the following 24–36 h. A similar response is evident in the upper right plot of Fig. 4 which shows the comparison between the observed and model rVTEC for the period of time between 1 and 10 of April 2008. The ionospheric response to the weak geomagnetic storm (Kp is less than 5) during 6–7 of April is described comparatively well.
1047
The suggested in this study method for empirical modeling of the rVTEC response to the geomagnetic activity is able satisfactorily to model the observed alternative responses with similar amplitudes of both phases because of the availability of two delayed mechanisms with different time constants. A significant role in improving the modeling of such cases plays also the introduced dependence on the LT. Fig. 5 displays two examples of the middle latitude (at geographic point (40°N, 130°W)) ionospheric responses to summer geomagnetic storms of medium intensity. In these cases a predominant positive response is seen but also a prolonged negative one. Fig. 6 presents a tropical ionospheric response to a geomagnetic storm in autumn, i.e. in the middle of October 2008. The maximum of the geomagnetic storm is near noon time of October 11 and the ionospheric response is evident again near 6 h later. A predominant positive response is seen from the figure. The above presented ionospheric responses to the geomagnetic storms cover all seasons, winter (Fig. 3), spring (Fig. 4), summer (Fig. 5) and autumn (Fig. 6). It has been already mentioned that the dependence of the rVTEC on the Kp-index, presented by the cross-correlation analysis, experiences latitudinal and seasonal variability. The model coefficients and respectively the values of rVTEC from (1) are calculated for all points of the grid with available data and for all months of the year. The reliability of the all model rVTEC values is practically the same as those shown in the Figs. 3–6. The error of the model varies from 0.05 to 0.15 rVTEC units for different months of the year and for different latitudes and longitudes. 6. Discussion and summary On the basis of the 2D (latitude-time) cross-correlation analysis, described in Part 3, it follows that the ionospheric response (expressed by rVTEC) to the geomagnetic activity (described by Kp-index) experiences both latitude and seasonal dependences. While the cross-correlation function at low (10°N–30°N) and high (higher than 55°N) latitudes is predominantly positive with a short time delay, up to 10 h, that at middle latitudes is predominantly negative with long tide delay, up to 60 h. The seasonal dependence is characterized by predominantly positive response in winter with a short time delay and mainly negative response in summer with a longer than 15 h time delay. The presence of both positive and negative rVTEC responses to the geomagnetic activity imposed the implementation of two different time delay constants Ts and Tl in order to describe the two different delayed reactions. It has been shown that the ionosphere can react to one and the same geomagnetic storm with both positive and negative phases (Fig. 4) which have different time delay constants (Ts is 6 h while Tl is up to 36 h). This example provides evidence that really the two time delay constants which describe the two different reactions have
1048
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
to be included in the regression equation (1). As the model rVTEC values are calculated for each point of the grid – longitude, latitude and LT it is possible the ionospheric response to be tracked out according to latitude, longitude and in time. In other words, the ionospheric response to a given geomagnetic storm can be tracked out along latitudes and longitudes at a fix time, or a sequence of regional maps representing the temporal evolution of the TEC response can be prepared as well. As a next step the presented here empirical model can be used for a short-term prediction of the rVTEC values depending on the LT and geomagnetic activity. This is a possible task because there are available models which predict the geomagnetic activity with a reliable accuracy. An example of such models is a MAK model described by Andonov et al. (2004). This model provides online prediction of the Kp-index and is implemented on the web site: http://www. geophys.bas.bg/NIGGG1/kp_for/kp_mod.php. If the ionospheric response to the geomagnetic activity can be predicted (i.e. rVTEC values can be predicted) then the correct phase delay of the propagating trough the ionosphere radio waves can be predicted as well. Such prediction will improve significantly the accuracy of the geodetic and navigation data which have increasing importance in resolving both scientific and practical tasks. Acknowledgement We are grateful to the NOAA National Geophysical Data Center for access to the VTEC data on http:// www.ngdc.noaa.gov/stp/iono/ustec/index.html. References Afraimovich, E., Lesyuta, O., Ushakov, I., Voeykov, S. Geomagnetic storms and the occurrence of phase slips in the reception of GPS signals. Ann. Geofis. 45 (1), 55–71, 2009. Araujo-Pradere, E.A., Fuller-Rowell, T.J., Spencer, P.S.T. Consistent features of TEC changes during ionospheric storms. J. Atmos. Sol.Terr. Phys. 68 (16), 1834–1842, 2006.
Andonov, B., Muhtarov, P., Kutiev, I. Analogue model, relating Kp index to solar wind parameters. J. Atmos. Sol.-Terr. Phys. 66, 927–932, 2004. Emmert, J.T., Lean, J.L., Picone, J.M. Record-low thermospheric density during the 2008 solar minimum. Geophys. Res. Lett. 37, L12102, doi:10.1029/2010GL043671, 2010. Habarulema, J.B., Mc Kinnel, L.A., Cilliers, P.J. Prediction of GPS total electron content using neural network over South Africa. J. Atmos. Sol.-Terr. Phys. 69 (15), 1842–1850, 2007. Jakowski, N., Stankov, S.M., Klaehn, D. Operational space weather service for GNSS precise positioning. Ann. Geophys. 23 (9), 3071– 3079, 2005. Kutiev, I., Muhtarov, P. Modelling the midlatitude response to geomagnetic activity. J. Geophys. Res. 106, 15501–15510, 2001. Kutiev, I., Muhtarov, P. Empirical modeling of global ionospheric foF2 response to geomagnetic activity. J. Geophys. Res. 108, A11021, 829– 832, 2003. Liu, J., Zhao, B., Liu, L. Time delay and duration of ionospheric total electron content responses to geomagnetic disturbances. Ann. Geophys. 28, 795–805, doi:10.5194/angeo-28-795-2010, 2010. Liu, J.Y., Tsai, H.F., Jung, T.K. Total electron content obtained by using the global positioning system. Terr. Atmos. Ocean. Sci. 7, 107–117, 1996. Muhtarov, P., Kutiev, I., Cander, L. Geomagnetically correlated autoregression model for short-term prediction of ionospheric parameters. Inv. Probl. Phys. 17, 1–17, 2001. Radicella, S., Zhang, M. The improved DGR analytical model of electron density height profile and total electron content in the ionosphere. Ann. Geofis. 38 (1), 35–41, 1995. Stankov, S.M., Kutiev, I., Jakowski, N., Wehrenpfennig, A. A new method for total electron content forecasting using global positioning system measurements, in: Proceedings of the ESTEC Workshop on Space Weather, Noordwijk, The Netherlands, 17–19 December 2001, pp. 169–172, 2001. Stankov, S.M., Kutiev, I.S., Jakowski, N., Wehrenpfennig, A. GPS TEC forecasting based on auto-correlation analysis. Acta Geod. Geophys. Hung. 39 (1), 1–14, 2004. Stankov, S.M., Stegen, K., Warnant, R. Seasonal variations of storm-time TEC at European middle latitudes. Adv. Space Res. 46 (10), 318–1325, 2010. Wang, X., Sun, Q., Eastes, R., Reinisch, B., Valladares, C.E. Short-term relationship of total electron content with geomagnetic activity in equatorial regions. J. Geophys. Res. 113, A11308, doi:10.1029/ 2007JA012960, 2008. Zhao, B., Wan, W., Liu, L., Mao, T. Morphology in the total electron content under geomagnetic disturbed conditions: results from global ionosphere maps. Ann. Geophys. 25, 1555–1568, doi:10.5194/angeo25-1555-2007, 2007.
Advances in Space Research 48 (2011) 1041–1048 www.elsevier.com/locate/asr
Empirical model of the TEC response to the geomagnetic activity over the North American region B. Andonov, P. Mukhtarov, D. Pancheva ⇑ Geophysical Institute, Bulgarian Academy of Sciences, Acad. G. Bonchev Str., Block 3, 1113 Sofia, Bulgaria Received 3 March 2011; received in revised form 2 May 2011; accepted 4 May 2011 Available online 12 May 2011
Abstract The paper presents an empirical model of the total electron content (TEC) response to the geomagnetic activity described by the Kpindex. The model is built on the basis of TEC measurements covering the region of North America (50°W–150°W, 10°N–60°N) for the period of time between October 2004 and December 2009. By using a 2D (latitude-time) cross-correlation analysis it is found that the ionospheric response to the geomagnetic activity over the considered geographic region and at low solar activity revealed both positive and negative phases of response. The both phases of the ionospheric response have different duration and time delay with respect to the geomagnetic storm. It was found that these two parameters of the ionospheric response depend on the season and geographical latitude. The presence of two phases, positive and negative, of the ionospheric response imposed the implementation of two different time delay constants in order to properly describe the two different delayed reactions. The seasonal dependence of the TEC response to geomagnetic storms is characterized by predominantly positive response in winter with a short (usually 5–6 h) time delay as well as mainly negative response in summer with a long (larger than 15 h) time delay. While the TEC response in March and October is more close to the winter one the response in April and September is similar to the summer one. Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. Keywords: Total electron content; Geomagnetic storm; Time delay; Positive and negative response
1. Introduction Solar activity such as flares and coronal mass ejections often produce large variations in the particle and electromagnetic radiation incident upon the Earth. Such variations can lead to a significant perturbation of the “quiet-time” ionosphere due to large variability in the ionospheric density distribution, total electron content (TEC), and the ionospheric current system. The geomagnetic storms have important terrestrial consequences such as disrupting satellite communications or detection and tracking of aircraft, distort navigation and rarely an interruption of the flow of electrical energy over power grids. Thus, it is important to monitor ionosphere particularly during the geomagnetic storms, to study its response to them, and, if ⇑ Corresponding author. Tel.: +359 2 9793308; fax: +359 2 9713005.
possible, to forecast the evolution of the ionospheric variability. The development of the Global Navigation Satellite System (GNSS) during the last decade has provided a number of possibilities for studying the spatial distribution and temporal variability of ionospheric electron density disturbances forced by external or internal sources. The measure of the phase delay of the propagating through the ionosphere high frequency (HF) radio waves is a key parameter for defining the total electron content (TEC) between the satellite and the receiver (Liu et al., 1996). The accurate measurement of the phase delay is a strong requirement for the reliable performance of many applications including HF communications, satellite positioning and navigation applications. The necessity of the data corrections obtained from the navigation satellites imposes the requirement not only the regular variability of the electron density, but also
E-mail address: [email protected] (D. Pancheva). 0273-1177/$36.00 Ó 2011 COSPAR. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.asr.2011.05.007
1042
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
the variability related to the geomagnetic storms or other drivers to be taken in mind (Jakowski et al., 2005). Before sufficiently long TEC time series to be accumulated by using the navigation satellites the problem for the TEC prediction had to be solved on the basis of the existing electron density profile models (Radicella and Zhang, 1995). The profiles were constructed by using the main ionospheric parameters determined by the ground based ionosonde stations. The response of the main ionospheric characteristic foF2 (known as the critical frequency of the ionospheric F-region or respectively) to the geomagnetic activity and its dependence on the geographic/geomagnetic coordinates, season and local time (LT) has been investigated by Kutiev and Muhtarov (2003). This empirical model is based on the all available ionosonde measurements. The foF2 prediction according to the level of the geomagnetic activity provides a possibility for an approximate assessment of the TEC variability. This approach however has a significant disadvantage: the ionosonde measurements cannot provide any information about the electron density profile above the F-region maximum. The upper part of the electron density profile however contributes to the TEC and cannot be ignored. The recent investigations of the TEC response to the geomagnetic storms based on the measurements directly obtained from the navigation satellites have reported the main features of the TEC response. Zhao et al. (2007) presented the latitude/longitude distribution of the relative deviation of the TEC (rTEC) during winter, summer and equinoxes, as well as the effect of the LT. Wang et al. (2008) studied the correlation dependences between the low latitude TEC and the equatorial Dst-index. Habarulema et al. (2007) suggested a TEC prediction method by means of the Neural Network (NN) based model. One of the input parameters of this method is the Kp-index. This study shows also the difference between the TEC defined by the navigation satellites and that defined by ionosonde measurements (i.e. only up the F-region maximum). The results indicated that usually the TEC defined by the ionosonde profile is on the average 75% from the entire TEC and that the variability of the both types of TEC data can be significantly different particularly for the considered in the paper examples. Afraimovich et al. (2009) investigated the dependence of the TEC variability on the Kp-index for equatorial, middle and high latitudes in the Northern Hemisphere (NH). The authors obtained linear regressions of the dependence of the relative TEC deviation (rTEC) on the geomagnetic index for the considered latitude range. Only the positive dependences between the rTEC and Kp-index are found; the response at high latitudes (50°N–80°N) is stronger than that at low latitudes. The authors found also significant longitudinal differences (geographic frame is used in the above mentioned study). Liu et al. (2010) studied the time delay of the TEC response by using the cross-correlation analysis. Its global distribution is presented separately for the positive and negative phases depending on the season.
The obtained values for the time delay are between 2 and 17 h, i.e. a range close to that found by Kutiev and Muhtarov (2001), but on the basis of ionosonde data. Stankov et al. (2001) and later Stankov et al. (2004) presented a correlation method for the prediction of TEC depending on the geomagnetic activity for the geographic region (20°W–40°E, 32.5°N–70°N). This method is based on the ideas used for the foF2 prediction and reported by Muhtarov et al. (2001) and Kutiev and Muhtarov (2001). In the above mentioned correlation method a time delay constant of 18 h and a periodic Fourier function accounted for the LT effect were introduced. Very recently Stankov et al. (2010) have applied a common epoch analysis on data representing nearly 300 storm events from the last solar cycle. It was found that the storm-time behavior of TEC shows clear positive and negative phases with amplitudes that tend to increase during the more intense storms. The most pronounced positive phase is observed during winter, while the strongest and yet shortest negative phase is detected during equinox. The TEC response to the geomagnetic activity for the American sector was reported by Araujo-Pradere et al. (2006). The authors found consistent features from storm to storm were found, and these feature became more apparent when the data were separated between the “driven” phase of the storm, when the integral of ap-index is rising, and the “recovery” to the storm, when the integral of ap-index is declining. The existence of the positive phase over the continental USA at the beginning of the storm period appears to be related to the timing of the peak of the perturbation; a positive phase will be observed when the peak of the perturbation occurs near midnight UT. The basic aim of the present paper is to find a mathematical dependence between the most probable values of the rTEC and the geomagnetic activity described by the Kp-index as a function of the calendar month, geographic latitude and LT. The wanted empirical dependence is obtained by using a long time series of TEC measurements over the region covering the North America. A particular attention will be paid to both the inertness of the ionospheric processes and the phase (positive/negative) of the TEC response. A significant improvement of the present study is the inclusion of both positive and negative phases of the TEC response to the geomagnetic storms. While the used TEC measurements (2004–2009) belong mainly to moderate and partly extremely low solar activity conditions, the latter being responsible for anomalous conditions at least in the thermosphere (Emmert et al., 2010). Therefore, the TEC model presented here is valid mainly for low solar conditions. The solution of the above described task will improve the forecast of the evolution of ionospheric irregularities during the geomagnetic storms. 2. Data The time series of vertical TEC (VTEC) for the period of time between October 2004 and December 2009 are used in
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
this study. The VTEC data are hourly data and cover the following geographic region: longitudes from 150°W to 50°W with longitude resolution of 1° and latitudes from 10°N to 60°N with a latitude resolution also of 1°. The VTEC data were downloaded from the NOAA National Geophysical Data Center (NGDC) web site: http:// www.ngdc.noaa.gov/stp/iono/ustec/index.html. The geomagnetic activity is defined in this study by the global Kp-index as it describes well the geomagnetic activity from tropical to middle/high middle latitudes. The Kpindex data are downloaded from the Space Physics Interactive Data Resource (SPIDR), Boulder, Colorado for the considered period of time. In the model the VTEC response to the geomagnetic activity is investigated by the relative deviation of the VTEC defined as: rVTEC = (VTECobs VTECmed)/ VTECmed. In this way the effect of the regular seasonal, diurnal and solar changes is removed from the VTEC variability. The term VTECmed represents the monthly median value. Data are grouped into 12-month bins, as every bin contains all the available hourly data within the respective month of the year.
January
60
0.4 0.36 0.32
50
Lat [deg]
0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
40 35 30 25 20 15 -8
0
8
16
24
32
40
48
56
64
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20
April
-16
-8
0
8
16
24
32
40
48
64
56
45 40 35 30 25 20
45 40 35 30 25 20 15
10 -24
-16
-8
0
8
16
24
32
40
48
56
64
10 -24
72
July
60
Lat [deg]
45 40 35 30 25 20
20
16
24
32
40
48
56
64
72
October
-8
0
8
16
24
32
40
48
56
64
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
45 40 35 30 25 20
45 40 35 30 25 20
-16
-8
0
8
16
24
32
40
48
56
64
72
November
40 35 30 25 20
10 -24
-16
-8
0
8
16
24
32
Time lag [hours]
40
48
56
64
72
64
72
20
-16
-8
0
8
16
24
32
40
48
56
64
72
September
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20
-16
-8
0
8
16
24
32
40
48
56
64
72
December
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
45
10 -24
56
25
60
50
15
48
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
40
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
60
55
55
32
72
August
60
50
24
June
15
10 -24
8
16
30
15
0
8
40
10 -24
-8
0
35
10 -16
-8
45
15 -24
-16
50
15
60
Lat [deg]
-16
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
30 25
10 -24
60
55
35
55 0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
40
60
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
45
10 -24
72
May
60
55
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50
15
10 -24
72
March
60
15
-16
60
Lat [deg]
We use the planetary Kp-index as an indicator of the geomagnetic activity in this study. This index reflects both types of variability: from the equatorial ring current and the auroral currents. Another important reason for using this index is connected with the fact that there are online models for its prediction. It is well known that during the recovery phase of the ionospheric storms with geomagnetic origin the ionospheric reaction continues some time after the geomagnetic storm attenuation. This phenomenon aggravates the relationship between the Kp-index and the ionospheric anomalies. Later a method will be suggested for avoiding of this shortcoming which is based on the integration of the Kp-index. The effect of geomagnetic activity, described by the Kp-index, on the rVTEC is studied by cross-correlation analysis between the both parameters. A two-dimensional (time-latitude) cross-correlation function is calculated for each month of the year (from the 12month bins) because its seasonal dependence is well known. The calculated cross-correlations functions for all months are shown in Fig. 1. The most important result is that for
55
45
10 -24
3. Cross-correlation analysis between rVTEC and Kp-index
February
60
55
1043
0.4 0.36 0.32 0.28 0.24 0.2 0.16 0.12 0.08 0.04 0 -0.04 -0.08 -0.12 -0.16 -0.2 -0.24 -0.28
50 45 40 35 30 25 20 15
-16
-8
0
8
16
24
32
Time lag [hours]
40
48
56
64
72
10 -24
-16
-8
0
8
16
24
32
40
48
56
64
72
Time lag [hours]
Fig. 1. Two-dimensional (latitude-time) cross-correlation functions calculated between the rVTEC and Kp-index for all months of the year; the zero line is shown by dashed white line.
1044
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 0.80 Polynomial Coefficients Degree 0: -6.00067 Degree 1: 3.07795 Degree 2: -0.524284 Degree 3: 0.0298974 Data
rVTEC
0.60
Best Fit 0.40
0.20
0.00 4.0
5.0
6.0
7.0
8.0
9.0
Kp
Fig. 2. The empirical dependence between the Kp and rVTEC.
1.6
2.8 2.6 2.4 2.2 2 1.8 1.6 1.4 1.2 1 0.8 0.6 0.4 0.2 0 -0.2 -0.4 -0.6 -0.8
November 2004
December 2005
1.4
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
1.2 1
rVTEC
rVTEC
all months the ionospheric response is composed by two phases, positive and negative, with different duration and time delay. The results from the figure clearly indicate that the cross-correlation depends on both latitude and time lag (which reflects the delayed reaction of rVTEC with respect to Kp changes) and this dependence is different for each month of the year. The cross-correlation function is positive (reaching up to 40% for December) with time lag not longer than 10 h for tropical (10°N–30°N) and high (higher
than 50°N) latitudes, while the cross-correlation function for middle (30°N–50°N) latitudes is negative (up to 30% for April) with time lag up to 6 times longer than that for tropical and high latitudes. Hence, during winter (from November to February) the cross-correlation function for tropical and high latitudes is characterized by high positive values while for middle latitudes it is predominantly negative and with time lag longer than 40 h. For the spring months (March–April) the cross-correlation function has both high positive and increasing negative values but the latter are reached for a shorter than in winter time lag of 20 h. The cross-correlation function for May and June indicates a decrease of the positive values while its negative behavior is similar to that in spring. The cross-correlation function for the months from July to October is similar to that for March–April with an exception for August when the negative cross-correlation is small and confined only to latitudes near 40°N. We note that the result for August could be affected by the decreased quality of the VTEC data for this month of the entire period of time (2004–2009). In summary, it can be concluded that the crosscorrelation function: (i) for almost all months of the year at low and high latitudes is mainly positive with a time lag of up to 10 h, (ii) for winter it has predominantly positive values which move toward the middle latitudes, and
0.8 0.6 0.4 0.2 0 -0.2 -0.4
06
07
08
09
10
11
12
12
13
14
Day
16
17
18
19
Day
9
9
November 2004 8
December 2005
8
KpTs KpTl Kp Data
7
KpTs KpTl Kp Data
7
6
6
5
5
Kp
Kp
15
4
4
3
3
2
2
1
1
0
0
06
07
08
09
Day
10
11
12
12
13
14
15
16
17
18
19
Day
Fig. 3. (Upper row of plots) Comparison between the observed (dash line) and model (solid line) rVTEC at the geographic point (20°N, 130°W) for the geomagnetic storm in November 2004 (left plot and December 2005 (right plot); (bottom row of plots) comparison between the measured Kp-index (tick solid line with dots) and modified K p T s - (dash line) and K p T l - (solid line) indices for November 2004 (left plot) and December 2005 (right plot). The time delay constants for November 2004 geomagnetic storms are: Ts = 5 h; Tl = 15 h, while those for December 2005 are: Ts = 6 h; Tl = 48 h.
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
(iii) for equinoxes and summer it has mainly negative values at middle latitudes with time lag of 20 h; the negative cross-correlation values in equinoxes are larger than those in the summer. The above mentioned results clearly indicated that two different time constants, Ts and Tl, should be introduced in the empirical model, representing the different time delays of the rTEC with respect to the geomagnetic activity described by the Kp-index. 4. Empirical model The investigations on the foF2 response to the geomagnetic activity (Muhtarov et al., 2001; Wang et al., 2008) indicated that this is a delayed response which can be satisfactorily modeled by assuming that the geomagnetic influence is imposed on the inertial system described by an inhomogeneous differential equation from a first order (Kutiev and Muhtarov, 2003). These studies found that the time delay constant is near 18 h. The above presented cross-correlation analysis however indicated that in this case the TEC reaction to the geomagnetic activity has to be presented by a sum of two responses with different time delay constants and with different sign of the cross-correlation function. A possible reason for the observed difference between the geomagnetic influence on the foF2 and TEC is probably related to a difference in the altitude structure of the electron density profile response.
If we assume that the impact of the geomagnetic activity on the TEC is accomplished be two mechanisms with different time delay constants then the variability of rVTEC can be described as follows: rVTECðtÞ ðfT s ðK p T s ðtÞÞ þ fT l ðK p T l ðtÞÞÞf ðLTÞ
ð1Þ
where the function f(LT) represents the dependence of the response on the LT at equal other conditions. The parameters K p T s and K p T l are the modified with time delay constants respectively Ts and Tl values of the Kp-index. These modified parameters are solutions of the equations shown below and are obtained easily by a numerical integration:
Ts Tl
dK p T s ðtÞ þ K p T s ðtÞ ¼ K p ðtÞ dt dK p T l ðtÞ dt
ð2Þ
þ K p T l ðtÞ ¼ K p ðtÞ
ð3Þ
It is convenient the unknown functions fT s and fT l from (1) to be expressed by their Taylor time series expansions while the dependence on the LT to be presented by a Fourier time series as follows:
0.6
0.6
April 2008
April 2005
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
0.4
0.4
0.2
rVTEC
rVTEC
1045
0
0.2
0 -0.2
-0.4
-0.2
01
02
03
04
05
06
07
08
09
01
02
03
04
Day
06
07
08
09
10
Day
8
8 April 2005
7
April 2008 7
KpTs KpTl Kp Data
6
KpTs KpTl Kp Data
6
5
5 Kp
Kp
05
4
4
3
3
2
2
1
1
0
0 01
02
03
04
05
Day
06
07
08
09
01
02
03
04
05
06
07
08
09
10
Day
Fig. 4. The same as Fig. 3 but for the geomagnetic storms in April 2005 (left column of plots) and April 2008 (right column of plots); the time delay constants for both geomagnetic storms in April 2005 and 2008 are: Ts = 6 h; Tl = 36 h.
1046
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 2
3
2
3
fT s ðK p T s Þ ¼ a0s þ a1s K p T s ðtÞ þ a2s K p T s ðtÞ þ a3s K p T s ðtÞ þ fTl ðK p T l Þ ¼ a0l þ a1l K p T l ðtÞ þ a2l K p T l ðtÞ þ a3l K p T l ðtÞ P flt ðLTÞ ¼ b0 þ bi cosði 2p LT ui Þ 24 i
ð4Þ As the basic aim of this study is to build a model which presents the rVTEC changes trough the geomagnetic activity described by the Kp-index, then it means that a functional dependence between the Kp and rVTEC should be found. The most appropriate type of the wanted functional dependence can be defined empirically. For this purpose the mathematical expectations of the rVTEC are calculated at a preliminary grouping of the Kp values into subintervals with a given length. The range of Kp values (the min Kp is the most left value while the max Kp is the most right value on the axis) is divided into subintervals with equal length which are arranged in an ascending line. The mathematical expectation and the standard deviations of the rVTEC are calculated for each subinterval. The empirical dependence between the Kp and rVTEC is presented in Fig. 2. It is evident that the functional dependence between Kp and rVTEC is close to the cubic function. As the dependence between the Kp and rVTEC shown in Fig. 2 is close to the cubic function then in the Taylor expansions (4) only four terms have to be included. The next step is to obtain the most probable values of the
coefficients: ais, bi, ail, Ts and Tl from (4). This is a nonlinear optimizing task that can be solved by applying the “trial-and-error” method in a way that the best approximation in a sense of minimum least squares deviation has to be assured. In order to solve the problem the following steps are made: (i) it is given a range of Ts changes from 0 to 10 h with a time resolution of 1 h and a range of Tl changes from 11 to 60 h with a time resolution again of 1 h; (ii) for each grid point of the built in this way grid the coefficients ais, bi, ail are found by using the method of least squares best fit, and (iii) the coefficients ais, bi, ail, Ts and Tl at which the best approximation (in a sense of minimum least squares deviation) is obtained are accepted as optimal coefficients for the model. 5. Model results The upper row of plots of Fig. 3 show the observed (dash line) and model (solid line) variability of the rVTEC during two strong geomagnetic storms taking place in November 7–9, 2004 (left plot) and December 14–16, 2005 (right plot) for a geographical point (20°N, 130°W). These geomagnetic storms are characterized by a predominantly positive response of the rVTEC. The bottom row of plots of Fig. 3 show the comparison between the Kp-index and the determined modified Kp indices with the respective time constants Ts and Tl. It is evident from the figure that 0.8
0.6
July 2006
July 2006
rVTEC Data rVTEC Model
rVTEC Data rVTEC Model
0.6
0.4
rVTEC
rVTEC
0.4 0.2
0.2 0
0 -0.2 -0.2
-0.4
01
02
03
04
05
06
07
25
08
26
27
28
Day
30
31
01
29
30
31
01
8
8
July 2006 7
July 2006
7
KpTs KpTl Kp Data
6
KpTs KpTl Kp Data
6 5
Kp
5
Kp
29
Day
4
4
3
3
2
2
1
1 0
0
01
02
03
04
05
Day
06
07
08
25
26
27
28
Day
Fig. 5. The same as Fig. 3 but at the geographic point (40°N, 130°W) for the geomagnetic storms in beginning of July 2006 (left column of plots) and at the end of July 2006 (right column of plots); the time delay constants for both geomagnetic storms in July 2005 are: Ts = 10 h; Tl = 20 h.
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048 0.4
October 2008 rVTEC Data rVTEC Model
rVTEC
0.2
0
-0.2
10
11
12
13
14
15
16
15
16
Day 8
October 2008
7
KpTs KpTl Kp Data
6
Kp
5 4 3 2 1 0
10
11
12
13
14
Day
Fig. 6. The same as Fig. 3 but for the geomagnetic storm in October 2008; the time delay constants for the considered geomagnetic storm are: Ts = 10 h; Tl = 36 h.
the maxima of the considered geomagnetic storms are near midnight for both storms: on November 8 in the first case and on December 15 in the second case. The maximum of the ionospheric response is detected 6 h later in both cases: at 6 UT on November 8 in the first storm and at 6 UT on December 15 for the second storm. A similar result from the cross correlation analysis was shown in Fig. 1 for November and December. From the presented cross correlation analysis shown in Fig. 1 it was found that the rVTEC response to the geomagnetic activity can be both positive and negative but with different time delay constants. It turns out to be that the ionospheric response to one and the same geomagnetic storm can have both positive and negative phases. Such cases are displayed in Fig. 4. The upper left plot of Fig. 4 shows the comparison between the observed (dash line) and model (solid line) rVTEC for the period of time between 1 and 9 of April 2005. The geomagnetic storm starts near noon time on April 4. A positive response is seen during the first 24 h of the geomagnetic storm and a negative one in the following 24–36 h. A similar response is evident in the upper right plot of Fig. 4 which shows the comparison between the observed and model rVTEC for the period of time between 1 and 10 of April 2008. The ionospheric response to the weak geomagnetic storm (Kp is less than 5) during 6–7 of April is described comparatively well.
1047
The suggested in this study method for empirical modeling of the rVTEC response to the geomagnetic activity is able satisfactorily to model the observed alternative responses with similar amplitudes of both phases because of the availability of two delayed mechanisms with different time constants. A significant role in improving the modeling of such cases plays also the introduced dependence on the LT. Fig. 5 displays two examples of the middle latitude (at geographic point (40°N, 130°W)) ionospheric responses to summer geomagnetic storms of medium intensity. In these cases a predominant positive response is seen but also a prolonged negative one. Fig. 6 presents a tropical ionospheric response to a geomagnetic storm in autumn, i.e. in the middle of October 2008. The maximum of the geomagnetic storm is near noon time of October 11 and the ionospheric response is evident again near 6 h later. A predominant positive response is seen from the figure. The above presented ionospheric responses to the geomagnetic storms cover all seasons, winter (Fig. 3), spring (Fig. 4), summer (Fig. 5) and autumn (Fig. 6). It has been already mentioned that the dependence of the rVTEC on the Kp-index, presented by the cross-correlation analysis, experiences latitudinal and seasonal variability. The model coefficients and respectively the values of rVTEC from (1) are calculated for all points of the grid with available data and for all months of the year. The reliability of the all model rVTEC values is practically the same as those shown in the Figs. 3–6. The error of the model varies from 0.05 to 0.15 rVTEC units for different months of the year and for different latitudes and longitudes. 6. Discussion and summary On the basis of the 2D (latitude-time) cross-correlation analysis, described in Part 3, it follows that the ionospheric response (expressed by rVTEC) to the geomagnetic activity (described by Kp-index) experiences both latitude and seasonal dependences. While the cross-correlation function at low (10°N–30°N) and high (higher than 55°N) latitudes is predominantly positive with a short time delay, up to 10 h, that at middle latitudes is predominantly negative with long tide delay, up to 60 h. The seasonal dependence is characterized by predominantly positive response in winter with a short time delay and mainly negative response in summer with a longer than 15 h time delay. The presence of both positive and negative rVTEC responses to the geomagnetic activity imposed the implementation of two different time delay constants Ts and Tl in order to describe the two different delayed reactions. It has been shown that the ionosphere can react to one and the same geomagnetic storm with both positive and negative phases (Fig. 4) which have different time delay constants (Ts is 6 h while Tl is up to 36 h). This example provides evidence that really the two time delay constants which describe the two different reactions have
1048
B. Andonov et al. / Advances in Space Research 48 (2011) 1041–1048
to be included in the regression equation (1). As the model rVTEC values are calculated for each point of the grid – longitude, latitude and LT it is possible the ionospheric response to be tracked out according to latitude, longitude and in time. In other words, the ionospheric response to a given geomagnetic storm can be tracked out along latitudes and longitudes at a fix time, or a sequence of regional maps representing the temporal evolution of the TEC response can be prepared as well. As a next step the presented here empirical model can be used for a short-term prediction of the rVTEC values depending on the LT and geomagnetic activity. This is a possible task because there are available models which predict the geomagnetic activity with a reliable accuracy. An example of such models is a MAK model described by Andonov et al. (2004). This model provides online prediction of the Kp-index and is implemented on the web site: http://www. geophys.bas.bg/NIGGG1/kp_for/kp_mod.php. If the ionospheric response to the geomagnetic activity can be predicted (i.e. rVTEC values can be predicted) then the correct phase delay of the propagating trough the ionosphere radio waves can be predicted as well. Such prediction will improve significantly the accuracy of the geodetic and navigation data which have increasing importance in resolving both scientific and practical tasks. Acknowledgement We are grateful to the NOAA National Geophysical Data Center for access to the VTEC data on http:// www.ngdc.noaa.gov/stp/iono/ustec/index.html. References Afraimovich, E., Lesyuta, O., Ushakov, I., Voeykov, S. Geomagnetic storms and the occurrence of phase slips in the reception of GPS signals. Ann. Geofis. 45 (1), 55–71, 2009. Araujo-Pradere, E.A., Fuller-Rowell, T.J., Spencer, P.S.T. Consistent features of TEC changes during ionospheric storms. J. Atmos. Sol.Terr. Phys. 68 (16), 1834–1842, 2006.
Andonov, B., Muhtarov, P., Kutiev, I. Analogue model, relating Kp index to solar wind parameters. J. Atmos. Sol.-Terr. Phys. 66, 927–932, 2004. Emmert, J.T., Lean, J.L., Picone, J.M. Record-low thermospheric density during the 2008 solar minimum. Geophys. Res. Lett. 37, L12102, doi:10.1029/2010GL043671, 2010. Habarulema, J.B., Mc Kinnel, L.A., Cilliers, P.J. Prediction of GPS total electron content using neural network over South Africa. J. Atmos. Sol.-Terr. Phys. 69 (15), 1842–1850, 2007. Jakowski, N., Stankov, S.M., Klaehn, D. Operational space weather service for GNSS precise positioning. Ann. Geophys. 23 (9), 3071– 3079, 2005. Kutiev, I., Muhtarov, P. Modelling the midlatitude response to geomagnetic activity. J. Geophys. Res. 106, 15501–15510, 2001. Kutiev, I., Muhtarov, P. Empirical modeling of global ionospheric foF2 response to geomagnetic activity. J. Geophys. Res. 108, A11021, 829– 832, 2003. Liu, J., Zhao, B., Liu, L. Time delay and duration of ionospheric total electron content responses to geomagnetic disturbances. Ann. Geophys. 28, 795–805, doi:10.5194/angeo-28-795-2010, 2010. Liu, J.Y., Tsai, H.F., Jung, T.K. Total electron content obtained by using the global positioning system. Terr. Atmos. Ocean. Sci. 7, 107–117, 1996. Muhtarov, P., Kutiev, I., Cander, L. Geomagnetically correlated autoregression model for short-term prediction of ionospheric parameters. Inv. Probl. Phys. 17, 1–17, 2001. Radicella, S., Zhang, M. The improved DGR analytical model of electron density height profile and total electron content in the ionosphere. Ann. Geofis. 38 (1), 35–41, 1995. Stankov, S.M., Kutiev, I., Jakowski, N., Wehrenpfennig, A. A new method for total electron content forecasting using global positioning system measurements, in: Proceedings of the ESTEC Workshop on Space Weather, Noordwijk, The Netherlands, 17–19 December 2001, pp. 169–172, 2001. Stankov, S.M., Kutiev, I.S., Jakowski, N., Wehrenpfennig, A. GPS TEC forecasting based on auto-correlation analysis. Acta Geod. Geophys. Hung. 39 (1), 1–14, 2004. Stankov, S.M., Stegen, K., Warnant, R. Seasonal variations of storm-time TEC at European middle latitudes. Adv. Space Res. 46 (10), 318–1325, 2010. Wang, X., Sun, Q., Eastes, R., Reinisch, B., Valladares, C.E. Short-term relationship of total electron content with geomagnetic activity in equatorial regions. J. Geophys. Res. 113, A11308, doi:10.1029/ 2007JA012960, 2008. Zhao, B., Wan, W., Liu, L., Mao, T. Morphology in the total electron content under geomagnetic disturbed conditions: results from global ionosphere maps. Ann. Geophys. 25, 1555–1568, doi:10.5194/angeo25-1555-2007, 2007.