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

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 destruc‐ tively 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].


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 4 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 4 ) 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 4 , 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 4 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.

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 4 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.

B-spline model derivation
The reason of using DeBoor B-spline function [15] is one of the most famous property of Bspline 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 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 Bspline 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 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].

Result and Discussion
From the weak scatter phase screen model introduced by C. L. Rino [11] amplitude scintillation index S 4 can be expressed as as 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.
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. 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.    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.
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.    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.    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].

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.