B-Spline Model of Ionospheric Scintillation For – High Latitude Using In-situ Satellite Data

1. Introduction

Ionospheric scintillation is a popular phenomenon among space scientists and GNSS users. It has been widely discussed and studied in past but still difficult to model and predict on large scales. Ionospheric scintillations are caused by rapid random variations of the phase and amplitude of the radio waves passing through the ionosphere. As the signal propagation continues after passing through the region of irregularities in the ionosphere, phase and amplitude scintillation develops through interference of multiple scattered waves. After propagation to a receiver, the irregular phase may combine either constructively or destructively to increase or decrease the wave amplitude. Another possibility is that the cause of either increased or decreased phase velocity may be refractive when an electromagnetic wave enters a medium [8].

The first empirical model of scintillation was proposed by Fremouw and Rino in 1973 [6]. This model could estimate the scintillation index S₄ on VHF/UHF, under weak scatter conditions. Weak scatter condition is often violated near the equatorial anomaly and auroral regions. This model led the foundation of more advanced model “WBMOD”. Aarons developed analytic model in 1985 [1] using 15-min peak to peak scintillation indices (not S₄) taken over 5 years at Huancayo, Peru using LES 6 satellite transmitted at 254 MHz. Next comes India model by Iyer and his group in 2006 [7]. They used cubic-B spline technique to develop an empirical model of magnetic quiet time scintillation occurrence at Indian equatorial and low latitudes. 250 MHz signal from FLEETSAT satellite was measured for 2 years at Trivandrum, near magnetic equator and Rajkot at the crest of equatorial anomaly. To describe the structure and extent of the radio scintillation generated by turbulence around and within the equatorial plumes a physical model has been developed by J. M. Retterer [10].

The first climatological model WBMOD has been developed by Northwest Research Associates, Inc. in which the user can specify his operating scenario. As the output the model returns: the phase scintillation spectral index p, the spectral strength parameter T, S₄, and phase scintillation index σ_ϕ. GISM has been described by Beniguel and Buonomo in year 1999 [4]. The model consists of two parts the NeQuick model and the scintillation model based on multiple phase screen algorithm. 2nd part of the model needs statistical information about irregularity as input. The algorithm is used to calculate the scintillation index at the receiver.

Basu and his group used first time satellite in situ data in scintillation modeling in 1976 [2]. They assumed a 3D power law irregularity spectrum with a constant spectral index of 4. They prepared another high latitude scintillation model in 1981, 1988 [3] using Atmospheric Explorer D data. Due to limited availability of data the model was suitable for northern winter under sunspot minimum condition.

Wernik et al. [12] used the Dynamics Explorer B data to estimate the irregularity spectral index and turbulence strength parameter, the factors that are required to calculate the scintillation index [11]. Their approach has been extended by Liu et al. [16] by introducing the finite outer scale.

Present model makes use of Dynamic Explorer 2 plasma density data covering period of August 1981 to February 1983. This period was near to maximum solar activity. In this model we are using the turbulence strength parameter C_s and the spectral index derived from Wernik et al. [12]. Simplest phase screen model described by Rino has been used to derive S₄ index. The parameters derived from Dynamic Explorer 2 satellite data are used with IRI model [5]. For comparing present model to the WAM model we produce maps in magnetic local time (MLT) and invariant latitude.

We present a spline model for the high latitude ionospheric scintillation using satellite in situ measurements made by the Dynamic Explorer 2 (DE 2) satellite. This analytical model is based on products of cubic B-splines and coefficients determined by least squares fit to the binned data. This product is constrained to make the fit periodic in 24 hours of geomagnetic local time, periodic in 360 degree of invariant longitude, in geomagnetic indices and solar radio flux. Discussion of our results clearly shows the seasonal and diurnal behavior of ionospheric parameters important in scintillation modeling for different geophysical and solar activity conditions. We also show that results obtained from our analytical model match observations obtained from in situ measurements. DE 2 satellite measurements give observations only along satellite orbit but our interpolation model fills the gaps between the satellite orbits.

2. Data preparation

The input data to our scintillation model are DE 2 retarding potential analyser (RPA) measurements of the ion density, equivalent by the charge neutrality to the electron density N_e. The altitude of its orbit was between about 300 and 1000 km (perigee: 309 km, apogee 1012 km, inclination: 89.99°, period: 98 min). The satellite was on a nearly polar orbit. The sampling frequency of RPA was 64 Hz, corresponding to every 120 m along the satellite orbit. These measurements were grouped over 8 s (512 samples) long segments [12].

The parameters derived from Dynamic Explorer 2 satellite (DE 2) data have been grouped separately into seasonal bins and for specific duration of case study (for studying geophysical events). These binned data have been appended from the few parameters obtained from IRI model and NASA’s Goddard Space Flight Centre (http://omniweb.gsfc.nasa.gov/html/omni2_doc.html).

NASA’s data server was used to get K_p index, F10.7 cm solar radio flux, Geo-magnetic field and sun spot number data. IRI model was used to append the binned data set with mean electron density, electron density peak height, and calculated irregularity layer thickness. The amplitude scintillation index S₄ was calculated using Rino’s [11] weak scatter phase screen formula. In this calculation we took only those values of one dimensional spectral index p, which was less than 4.

3. B-spline model derivation

The reason of using DeBoor B-spline function [15] is one of the most famous property of B-spline functions that they has minimal support to a given degree of freedom, smoothness and domain partition. B-spline models are best because they provide similar results, even when using low-degree splines to the models produced using higher degree polynomials while avoiding instability at the edges of an interval (Runge’s phenomenon).

Parameter's derived from DE 2 satellite data as a function of local time, day/season/month, geographic coordinates, K_p index and solar flux value F10.7, is expressed as simultaneous product of univariate normalized B-splines as given below

$$\begin{array}{c}Parameter\left(t,d,geo.coord.,\mathrm{ }kp,\mathrm{ }F10.7\right)=\\ \sum _{i=1}^{24}\sum _{j=1}^{365}\sum _{k=1}^9\sum _{l=1}^6\begin{array}{c}{a}_{i,j,k,l}*N_{i,4}\left(t\right)N_{j,2}\left(d\right)N_{k,2}\left(Kp\right)N_{l,2}\left(F10.7\right)\end{array}\end{array}\mathrm{ }$$E1

$$\begin{array}{c}{S}_4\left(t,d,geo.coord.,\mathrm{ }kp,\mathrm{ }F10.7\right)=\mathrm{ }\\ \sum _{i=1}^{24}\sum _{j=1}^{365}\sum _{k=1}^9\sum _{l=1}^6\begin{array}{c}{a}_{i,j,k,l}*N_{i,4}\left(t\right)N_{j,2}\left(d\right)N_{k,2}\left(Kp\right)N_{l,2}\left(F10.7\right)\end{array}\end{array}$$E2

where a_i,j,k,l are monthly mean of amplitude scintillation index and/or parameters derived from RPA measurements for each interval of magnetic local time, invariant coordinate, K_p index and solar radio flux F10.7 cm. N_i,4 is a b-spline basis function of degree 4 and other b-spline basis function are of degree 2.

These all 4 B-spline basis functions are non-vanishing over limited intervals. They all add up to one at all magnetic local time, season, K_p index and F10.7 solar radio flux interval. For more details about properties of basis function one should refer to DeBoor [15]

We have used 3-hourly time nodes for magnetic local time. They vary between 0 to 23 hour as 3, 6, 9, 12, 15, 18, 21 and 24. For seasonal maps we have used 3-hourly magnetic local times and 10 day (the maps resolution is 3 h X 10 days) median of binned data. For K_p index 9 nodes were chosen which vary between 1 to 9. 6 individual nodes have been chosen for F10.7 cm radio flux value, they are 80, 130, 180, 230, 280 and 330 respectively.

The number and placement of magnetic local time nodes for each season and solar flux interval were individually chosen to account for large variability in amplitude scintillation index and other modeled parameters. It is tricky to cleverly observe rapid changes in amplitude scintillation index and parameters. Consequently more basis functions are needed to account for these rapid changes. Therefore, placement and number of magnetic local time nodes are different for different seasons and geophysical cases.

We could have used the higher density mesh of basis functions for all geophysical case studies. This gives freedom to the programmer and one can approach closer to the real observational results. Which simply means one should see the actual behavior of modeled parameter and cleverly chose the number and placement of more or less basis function in order to derive same parameter from the model.

The coefficient a_i,j,k,l were determined by least square fit to the binned data and constrained to make the fit periodic in 24 hour and 360 geomagnetic longitude. The local time and diurnal B-spline functions are shown in Fig. 1. As we have already discussed for quiet geophysical conditions smaller number of basis functions are sufficient which should be

Figure 1.

Basis function for the local time (upper figure) and month of the year (lower figure).

Considered as a quality of this consolidated model. This is because of less need of the modelling coefficients (e.g. basis functions) for modelling quiet geophysical condition (e.g., low solar flux). Using more B-spline function and their placement at right positions we can upgrade our model, which makes us enable of modelling disturbed geophysical conditions. This freedom outweighs our consolidated model. For comparing our model with observation we have prepared contour maps for real observation which use “contour” MatLab subroutine which uses linear interpolation method for plotting. Our modelled contours use B-spline interpolation method. From in situ measurements we derived the turbulence strength parameter C_s and the spectral index using the method discussed in WAM model. With the IRI model [5] the C_s parameter was rescaled to get its value at the height of the maximum electron density. The IRI model was also used to estimate the irregularity layer thickness. To convert the parameters derived from RPA measurements to the equivalent scintillation index one should rely on the scintillation theory [12].

4. Result and Discussion

From the weak scatter phase screen model introduced by C. L. Rino [11] amplitude scintillation index S₄ can be expressed as as

$S_4^2={\left(r_e\lambda \right)}^{2}L\text{csc}E{C}_S{Z}^{p/2}\frac{\Gamma \left(1-p/4\right)}{{\pi }^{1/2}\Gamma \left(p/4+0.5\right)p}F\left(a,b\right)$E3

where r_e is classical electron radius, λ is wavelength of the signal, L is the irregular layer thickness, θ is the zenith angle, C_s is turbulence strength parameter and Z is the Fresnel zone parameter and F is Fresnel filter factor.

The one dimensional spectral index p is an important parameter which determines the scintillation level [12]. Therefore, variation of spectral index with ionospheric changes becomes significantly important for studying the scintillation effect on trans-ionospheric communication links. Fig 2 shows the behaviour of spectral index in equinox for geomagnetic quiet conditions.

Figure 2.

Spectral index for equinox for K_p ≤ 3. Left contour is built from real data. Right picture is spectral index modelled using B-spline technique.

It is evident from the above figure that in quiet geomagnetic condition one dimensional spectral index intensifies near magnetic noon. For invariant latitude >70 degree, spectral index is more intense which in start of auroral zone. It seems expanding from auroral boundary to polar cusp region in night time. One can easily observe that in geomagnetic quiet condition modelled map is in good agreement with the one produced from real observation. Fig. 3 below shows the behaviour of spectral index for geomagnetic disturb days in equinox months.

Figure 3.

Spectral index for equinox for K_p > 3. Left contour is built from real data. Right one is spectral index modeled using B-spline technique.

It is evident from the figure that spectral index seems expanding from the auroral boundary towards the equator. The spectral index level is high in geomagnetic disturbed conditions of the equinox than that of the geomagnetic quiet days of the equinox months. Though the maxima is visible near the magnetic noon. Modelled results are in good agreement with the observations.

Figure 4.

Spectral index for K_p ≤ 3. Left contour is for summer. Right one is spectral index for winter. Both are modelled using B-spline technique.

Fig. 4 shows the modelled contours for winter and summer in quiet geomagnetic condition. During summer and low geomagnetic conditions largest mean value of spectral index is observed at high latitude (> 70 degree invariant latitude). At low latitudes, in summer the mean spectral index is smaller but larger than that in winters. Fig. 5 represents geomagnetic disturbed behavior of spectral index for summer and winter. One can easily observe that in summer at high latitude greater than 70 degree spectral index is independent of geomagnetic effect. Nevertheless, spectral index at low latitudes (i.e. less than 70 degree) increases with geomagnetic activity both in summer and winter.

Figure 5.

Spectral index for K_p > 3. Left contour is for summer. Right one is spectral index for winter. Both are modelled using B-spline technique.

Seasonal behaviour of amplitude scintillation index and other ionospheric parameter have been modelled for geomagnetic quiet and disturb conditions. Fig. 6 shows the turbulence strength parameter for equinox when the K_p≤3. As we know invariant latitude is a parameter that describes where a particular magnetic field line touches the Earth. On the Earth’s surface, the invariant latitude is equal to the geomagnetic latitude. Fig. 6 shows that C_s maximizes in auroral region near and after magnetic midnight. The left figure is created from the real observation and the right contour map is prepared using our B-spline model. As it is evident from the figure above, the model and observations are is good agreement. Fig. 7 below shows the variation of C_s with invariant latitude and magnetic local time for equinox month and in geomagnetic disturbed condition.

Figure 6.

Turbulence strength parameter C_s of plasma density fluctuation for k_p ≤ 3. Left map is from the real observation and right map is produced using B-spline model.

Invariant latitude >70 degree seems dependent of geomagnetic activity. Near magnetic noon a significant increase is visible at auroral boundary which is expanding to the polar cusp near the midnight.

Figure 7.

Turbulence strength parameter C_s of plasma density fluctuation for K_p > 3. Left contour is for summer derived from observation, right one is for winter which is modelled using B-spline technique.

Figure 8.

Turbulence strength parameter C_s of plasma density fluctuation. Left contour is for K_p ≤ 3. Right one is spectral index for K_p > 3. Both are modelled using B-spline technique.

Fig. 8 above shows modelled turbulence strength parameter C_s for K_p≤ 3 and K_p >3. In both the high and low magnetic conditions we observe a considerable enhancement in C_s for invariant latitude >70 degree. This maximum seems very much coherent with the maximum of spectral index near polar cusp region and noon sector of the polar cap. The only difference is the time of enhancement. For low magnetic activity condition maximum is near magnetic noon. But, for strong geomagnetic activity maxima is visible near magnetic dusk. Fig. 9 is modelled turbulence strength parameter for winter. Left contour is for weak magnetic activity while right one is for strong magnetic activity condition. The value of C_s is one order of magnitude higher than that in summer. Here our results are fully consistent with the WAM model results [12]. In winter C_s is independent of magnetic local time. In winter low latitude C_s is always smaller as compared to the high latitude C_s. During disturbed geomagnetic conditions a significant enhancement in C_s is visible in the dawn and dusk time for invariant latitude > 70 degree. Up to now our model is in excellent coherence with the observational results in both the magnetic weak as well as strong magnetic activity conditions.

Figure 9.

Turbulence strength parameter C_s of plasma density fluctuation. Left contour is for K_p ≤ 3 and right one is for K_p > 3. Both are modelled using B-spline technique.

Figure 10.

Amplitude scintillation index during equinox K_p ≤ 3. Left contour is prepared from calculated S₄ but, right contour is modelled using B-spline technique.

As we know from Rino’s work on phase screen model, if we know turbulence strength parameter C_s and one dimensional spectral index p then it is possible to calculate equivalent amplitude scintillation index. While calculating equivalent amplitude scintillation index S₄ we have considered ionospheric irregularity to be isotropic for which Fresnel’s filter factor F=1(see equation number 3). Fig. 10 shows the computed and modelled behaviour of amplitude scintillation index for equinox during magnetic weak condition.

For high invariant latitude maxima is visible near polar cusp. It is evident from the figure that scintillation index is varying with magnetic local time but, enhancement is more clear near mid-night. Modelled contour is in good agreement with the computed S₄ index map using DE 2 retarding potential (RPA) measurements.

Figure 11.

Amplitude scintillation index during summer Kp ≤ 3. Left contour is prepared from calculated S₄ but, right contour is modelled using B-spline technique.

Fig. 11, 12, 13 and 14 are contour maps for summer and winter months for high and low magnetic activity conditions. These maps show that in summer the scintillation is much weaker than that in winter. During summer and low magnetic activity, the strongest scintillation is noted around magnetic noon and midnight at latitudes corresponding to the polar cusp and auroral zone, respectively.

During winter and low magnetic activity, strongest scintillation is observed in the polar cap. The maps show that with increasing magnetic activity the regions of the most intense scintillation expand equatorward. In winter the expansion of the scintillation zone is less defined in the mid-night sector, and the polar cap scintillation intensity weakens. Spline model magnetic activity variations of the high latitude scintillation zone are consistent with those found in the scintillation and other ionospheric irregularity-sensitive measurements [12, 13 and references therein].

Figure 12.

Amplitude scintillation index during winter k_p ≤ 3. Left contour is prepared from calculated S₄ but, right contour is modelled using B-spline technique.

Figure 13.

Amplitude scintillation index during summer k_p > 3. Left contour is prepared from calculated S₄ but, right contour is modelled using B-spline technique.

Figure 14.

Amplitude scintillation index during winter k_p > 3. Left contour is prepared from calculated S₄ but, right contour is modelled using B-spline technique.

5. Summary and conclusion

Our scintillation model makes use of the DE 2 retarding potential analyzer plasma density data [14] covering the period from August 1981 to February 1983, near to the solar maximum activity. DE2 in-situ measurements of plasma density fluctuations provide direct information of structure and morphology of irregularity that are responsible for scintillation of radio waves on trans-ionospheric links.

Described model is for northern hemisphere high latitude ionosphere which uses DE 2 RPA measurements. For geomagnetic activity dependence of scintillation there is good agreement between model and measurements. Spline model are the best because they provide similar results, even when we use low-degree splines, to the models produced using higher degree polynomials while avoiding instability at the edges of an interval (Runge’s phenomenon). This provides a reasonable realistic description of scintillation index and other ionospheric parameters.

Like any other model, our model also has certain limitations. Since it is an empirical model therefore we derive model from real observations. In any case our model will give suitable and convincing results for the geophysical condition which would be closely similar to the duration in which the data is recorded. Nevertheless, present model gives an average behavior of ionospheric parameters during different geophysical conditions. It can be compared with the observations performed in different solar activity condition and we strongly believe that the comparison would be convincing. DE 2 was working during moderate solar activity period when the sun spot number was between 80-140. Therefore our model is only valid for moderate solar activity conditions. We are using IRI model for the irregularity slab thickness and the height of peak electron density. IRI model often fails to give real behavior of ionospheric parameters for high latitude. There is possibility of erroneous calculation in our model similar to WAM model [12]. Third and most serious limitation of our model is placement and number of B-spline basis function. The number of data points for high geomagnetic activity condition is always less than that of weak geomagnetic condition. Sometimes it seems that in weak geomagnetic activity conditions the contour maps are smoother than that in high geomagnetic condition. This may also be considered as a limitation which can’t be overcome since it is natural. As we have already discussed that for individual geophysical situations we choose different number of basis function and keep on experimenting with the placement of basis function in order to get more convincing results. It is always possible that some one can use different set of B-spline basis functions and could be able to produce better modelling than we did here. But, we take this limitation positively. It is because that we feel confident that there is always possibility of upgrading in our model.

Acknowledgments

This research work is undertaken in the scope of the TRANSMIT ITN (www.transmit-ionosphere.net), funded by the Research Executive Agency within the 7th Framework Program of the European Commission, People Program, Initial Training Network, Marie Curie Actions – GA No. 264476.