Journal Pre-proofs A global empirical orthogonal function model of plasmaspheric electron content Peng Chen, Lixia Liu, Yibin Yao, Wanqiang Yao PII: DOI: Reference:
S0273-1177(19)30717-3 https://doi.org/10.1016/j.asr.2019.09.039 JASR 14464
To appear in:
Advances in Space Research
Received Date: Revised Date: Accepted Date:
15 January 2019 12 September 2019 23 September 2019
Please cite this article as: Chen, P., Liu, L., Yao, Y., Yao, W., A global empirical orthogonal function model of plasmaspheric electron content, Advances in Space Research (2019), doi: https://doi.org/10.1016/j.asr. 2019.09.039
This is a PDF file of an article that has undergone enhancements after acceptance, such as the addition of a cover page and metadata, and formatting for readability, but it is not yet the definitive version of record. This version will undergo additional copyediting, typesetting and review before it is published in its final form, but we are providing this version to give early visibility of the article. Please note that, during the production process, errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
© 2019 COSPAR. Published by Elsevier Ltd. All rights reserved.
A global empirical orthogonal function model of plasmaspheric electron content Peng Chen1, Lixia Liu1, Yibin Yao2, 3 and Wanqiang Yao1 1. College of Geomatics, Xi’an university of Science and Technology, Xi’an 710054, China 2. School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China 3. Key Laboratory of Geospace Environment and Geodesy, Ministry of Education, Wuhan University, Wuhan 430079, China Abstract: Based on the method for establishing a global plasmaspheric model using observations from COSMIC and MetOp-A orbit determination GNSS receivers, Chen et al. (2017) obtained a global plasmaspheric total electron content product with a spatial resolution of 2.5°×5° and a time resolution of 4h. In this paper, we use those global plasmaspheric electron content product in 2008, 2010, 2011, 2012, 2014 for 1446 days to establish a global plasmaspheric empirical model based on empirical orthogonal function (EOF). The model can well characterize the spatiotemporal variation of plasmaspheric electron content (PEC) and the influence of solar radiation on it. Only the first four orders of EOF sequences can characterize the 98.43% features of the original PEC dataset. The principal component coefficient Pk is decomposed twice during modeling, and the combination of trigonometric function and linear function is used to model Pk to characterize the solar cycle, annual cycle, semi-annual cycle and quarter-cycle variation. We compare the PEC model values with the actual observation data, the results show that the empirical PEC model values are highly correlated with the actual observations. The correlation between the two is above 0.96, and the RMS maximum of the difference between the PEC model values and the observed values are 0.70 TECU, and the average of the difference between the PEC model values and the observed values are -0.18 TECU, respectively. In addition, we validate the reliability of the global plasmaspheric model established by two empirical orthogonal function decomposition method using actual observation data, according to the global distribution of the differences between the PEC model values and the observed values in low solar activity and high solar activity, it can be seen that under low solar activity and high solar activity conditions, the model has good adaptability. Key Words: PEC; Empirical Orthogonal Function; EOF model; accuracy validation 1
1. Introduction The plasmasphere, which is located above the ionosphere, contains low energy plasma and can extend to an altitude 3-5 times the Earth’s radius. The eiectron density in the plasmasphere is several orders of magnitude less than that of the ionosphere, however, under some conditions (for example, at night and during low solar activity), the relative contribution of the plasmaspheric electron content (PEC) to the total electron content (TEC) can account for a large proportion. As an indispensable and important part of the earth’s space environment, they has an important impact on satellite navigation, positioning, radio communication, solar space environmental monitoring and human space activities. Moreover, there is a mechanism of material exchange and coupling between the plasmasphere and the ionosphere (Geiss et al., 1981; Horwitz J L et al., 1990; Khazanov et al., 1994; Belehaki et al., 2004; Zhang et al., 2013; Chugunin et al., 2017). Therefore, it is of great significance for the research and modeling of plasmasphere. The Earth’s plasmasphere contributes essentially to total electron content (TEC) measurements from ground or satellite platforms. (Jakowski et al., 2018). And the plasmaspheric electron content is an important physical parameter to describe the plasma morphology and structure. Kilmenko et al.,(2015) studied the contribution of the global plasmaspheric and ionospheric electron content (PEC and IEC) into total electron content (TEC), and the experimental PEC was estimated by comparison of GPS TEC observations and FORMOSAT-3/COSMIC radio occultation TEC measurements. The empirical model, based on statistical analysis of a large number of datasets, uses appropriate functions to represent the intrinsic characteristics of the plasmasphere and its variations, and has become means to study the plasmaspheric electron content, density distribution characteristics and temporal and spatial variation of different geomagnetic locations, local times, solar and geomagnetic activity conditions. In view of the application requirements of the plasmasphere in the delay error correction of satellite navigation and precise positioning systems, it is necessary to reveal the electron content and density variation of the plasmasphere in different geophysical conditions and the mechanism of mass exchange and coupling with ionosphere. Therefore, for decades, many scholars have done a lot of research on this. KI.Gringauz et al., (1960) discovered a sharp drop in density at the top of plasmasphere, and then D.Carpenter et al., (1963) discovered the same phenomenon. Webb and Essex (2000) established a global ionospheric/plasmaspheric electron 2
density model on the basis of IRI (International Reference Ionosphere) model, but the accuracy and reliability of plasmaspheric electron density obtained by this model are low. O'Brien et al. (2003) established a model of plasmapause location using the database of over 900 plasmapause crossings from in situ CRRES (Combined Release & Radiation Effect Satellite) electron density observations, and showed that, in addition to Kp, other quantities, such as Dst and AE can be used to build a quantitative model of the plasmapause location. D.L. Carpenter and R.R. Anderson (1992) developed an empirical model of equatorial electron density in the magnetosphere. Berube et al., (2004) used the ground-based Ultra Low Frequency (ULF) wave diagnostic data to establish the first empirical model of plasmasphere mass density in the equatorial region and used RPI (Radio Plasma Imager) electron density data to study changes of the average particle quality in plasmasphere under different geomagnetic intensities. Kotova et al., (2015) proposed Physics-based model of the plasmaspheric density. Cherniak et al.,(2018) examined performance of two empirical profile-based ionospheric models, and found that the major source of the model-data discrepancies in the EC/TEC domain comes from the topside ionosphere/plasmasphere system. EOF analysis, also known as Principal Component Analysis (PCA), was first invented by Pearson (1900). The method can compress the original dataset into several main orthogonal basis functions and related variables, and is sufficient to characterize the intrinsic characteristics of the data and most of the changes, so that the characteristic root sequences converge quickly. Based on advantages, the empirical orthogonal function (EOF) has become a powerful tool for extracting the main components of meteorological factors and temporal and spatial variations in the field of Meteorology (Lorenz, 1956). It is generally used as a key algorithm to precisely construct the ionospheric parameter (Zhao et al., 2005; Ercha et al., 2012). The method is originally introduced into meteorology as a method for extracting the dominant modes of spatial variability in meteorological fields (Mao et al., 2008). It is also applied to data analysis and empirical modeling of ionosphere and thermosphere (Daniell et al., 1995; Mao et al., 2005; Liu et al., 2008; A et al.,2011; Ercha et al., 2012; Chen et al., 2015). Mao et al., (2008) used the total electron content (TEC) data of International GNSS Service (IGS) and Crustal Movement Observation Network of China (CMONOC) during 1996-2004 to establish an empirical orthogonal function model of TEC over China in space weather and verified the accuracy and reliability of the model. And through the 3D data assimilation 3
experiments of EOF model and IRI model, it is proved that the EOF empirical model has better proximity forecasting ability and can provide more realistic background information. Chen et al., (2015) used the method of EOF decomposition using the GPS TEC data from 2001 to 2012 to establish the regional empirical model of the total electron content of the ionosphere in North America, and the result shows that the temporal and spatial variation characteristics of ionospheric total electron content (iTEC) can be characterized by using different EOF basis functions and related principal component coefficients of less order. Mukhtarov P. et al., (2013) used the Centre for Orbit Determination in Europe (CODE) TEC data from 1999 to 2011 to establish a global TEC model capable of responding to geomagnetic activity. This model can describe the spatial distribution and time variation of most TEC anomalies caused by geomagnetism. Ercha et al., (2012) used the global ionospheric maps in 1999-2009 provided by the Jet Propulsion Laboratory (JPL) to establish a global ionospheric total electron content map based on empirical orthogonal function analysis. Among them, the EOF basis function and correlation coefficient can well characterize the change of TEC and the influence of solar and geomagnetic activities on TEC. It is shown that 99.4% of the global ionospheric TEC can be represented by the first fourth-order EOF sequence. Based on the advantages of EOF analysis modeling, we extend this method to the plasmasphere to establish a global plasmaspheric electron content model under different solar conditions. The modeling dataset and modeling approach are described in the next section. The method of establishing global PEC empirical model based on EOF decomposition is introduced in detail. The distribution and variation characteristics of orthogonal basis functions of each order are explained. The periodic variation characteristics of correlation coefficient and its correlation with the solar activity index are also analyzed. The third section uses the global PEC data in 2009 and 2013 to verify the accuracy and reliability of the global plasmaspheric electron content model based on EOF decomposition. We draw our discussions and summaries in Section 4. 2. PEC dataset and modeling technology 2.1 PEC dataset The TEC along the signal propagation path from the GNSS satellite to receiver consists of two components: ionospheric electron content (IEC) and PEC. The altitude of COSMIC and MetOp-A is over 800km above Earth’s surface, and therefore, the 4
TEC obtained by observations of COSMIC and MetOp-A can be regarded as approximations of the PEC. MetOp-A is a joint mission of the ESA (European Space Agency) and EUMESAT (European Organization for the Exploitation of Meteorological Satellites), and it was launched into a Sun-synchronous circular polar orbit with 98.7°. The COSMIC is a six-satellite radio occulation mission. Their final orbits are at ~800 km, with a separation angle between neighbouring orbital planes of 30° longitude. Chen et al. (2017) proposed a method to build a global plasmaspheric model using observations from COSMIC and MetOp-A orbit determination GNSS receivers. A spherical harmonic function with 9 × 9 order is used. The noise of the pseudorange observation is reduced by the method of phase smoothing pseudorange. The cutoff elevation angle of projection function is 15°. The differential code bias (DCB) of the satellites and receivers is estimated as an unknown parameter and solved along with the model parameters. The 4-hour observation data is used as a time period to solve a set of model parameters. The spatial resolution of the global plasmaspheric total electron content grid product is 2.5°×5° and the time resolution is 4h. This paper uses those global plasmaspheric total electron content product in 2008, 2010, 2011, 2012 and 2014 for 1446 days to model. 2.2 PEC Data Decomposition and Model Construction The variation of plasmaspheric electron content is closely related to other factors such as solar radiation, geomagnetic latitude and local time. It should be emphasized that our experiments were carried out under the assumption of quiet geomagnetic activity, so we did not consider the geomagnetic activity factor, which means that the model we built is not suitable for geomagnetic storm activity and its vicinity. Therefore, the following global PEC datasets in 2008, 2010, 2011, 2012, and 2014 were modeled based on the EOF decomposition principle. The EOF analysis is a mathematical method that transform a data set composed of a number of multiple inter-correlated variables into a new dataset generally composed of small number of uncorrelated variables that still contains most of the information in the original dataset (Dunteman, 1989; Jolliffe, 2002; Ercha et al., 2012). The main goals of the EOF analysis are the reduction of dimensionality of a dataset, the discovery of hidden structures in a dataset, and their classification according to how much percentage of variance in the dataset they account for (Dunteman,1989). The orthogonal basis mode, which is also referred to as EOF, is defined as a function of geomagnetic longtitude and latitude to represent the spatial pattern in this study. The amplitude coefficient 5
contains a time-varying scalar amplitude of the spatial pattern. The temporal resolution of the amplitude coefficients is one day. By preserving the spatial patterns that capture most of the observed variance, the EOF analysis enables dimensionality reduction of the original dataset (Hannachi et al., 2010; Preisendorfer, 1988). The spatial resolution of the PEC dataset is 2.5° × 5°, and the time resolution is 4 hours. Once the PEC dataset are known, the original PEC dataset can be decomposed in terms of EOF base functions Fi and associated coefficient Pi . The decomposition process of global PEC dataset is as follows: N
PEC ( LT , lat , UT , d ) Fi ( LT , lat ) Pi (UT , d )
(2)
i 1
where, Fi ( LT , lat ) is an empirical orthogonal basis function that varies with local time (LT) and geomagnetic latitude lat. This basis function can not only characterize the spatial distribution of PEC but also show the variation characteristics of PEC with local time. Pi (UT, d) is the correlation coefficient corresponding to the EOF basis function, which is related to universal time (UT) and day of year (DOY), characterize the long-period temporal variation of the original PEC dataset. N represents the order of basis functions and correlation coefficients. Table 1 lists the variable contribution rates of the first eight order EOF basis function sequences to the original dataset. We calculated the percentage of the EOF series of the first EOF decomposition to the original dataset participating in the modeling. From the table, we can find that the contribution of the first four order of empirical orthogonal basis functions to the original dataset is 95.91%, 1.76%, 0.53% and 0.23%, respectively, and their cumulative contribution accounts for 98.43% of the original dataset. This phenomenon indicates that the original PEC dataset can be represented only by the first four order empirical orthogonal functions and their corresponding principal components, thereby greatly reducing the modeling parameters while ensuring high model accuracy and reliability. Table 1. Summary of the variances through PEC decomposition for the first eight EOF series EOF series
variances (%)
cumulative variances (%)
F1×P1
95.910
95.910
F2×P2
1.760
97.670
F3×P3
0.530
98.200 6
F4×P4
0.230
98.430
F5×P5
0.140
98.570
F6×P6
0.130
98.700
F7×P7
0.110
98.810
F8×P8
0.100
98.910
Figure 1 is the distribution of the first four orders orthogonal basis function Fk after the orthogonal decomposition of the PEC. It is apparent from figure 1 that the first-order basis function F1 exhibits hemispherical symmetry. It is worth noting that the first-order EOF basis function F1 reflects the global average composition of PEC, representing the overall spatiotemporal variation of PEC. At the same time, it can be seen that PEC exhibits north-south symmetry on a global scale, and it is mainly distributed in low latitude region (40°S-40°N). With the increase of latitude, the PEC value decreases, and shows a rapid drop at around ±40°latitudes. Moreover, PEC changes with local time, it reaches the peak at 12:00-15:00 LT. The most significant feature of the second-order EOF basis function F2 is the hemispheric asymmetry, which is positive in the northern hemisphere while negative in the southern hemisphere. This phenomenon is consistent with the law of the EOF decomposition of the ionosphere by Mao et al. (2008), the reason is that the annual variation of solar radiation from summer to winter is uneven. The main feature of the third-order EOF basis function F3 is that there are significant LT changes in the middle latitudes. The EOF sequence values of 10:00-20:00 LT are all negative values, and during other periods are mainly positive. At high latitudes, the value of the third-order EOF is basically negative. The fourth-order orthogonal basis function F4 shows that PEC exists semi-diurnal cycle difference in the northern and southern hemisphere. b) F2 0.025
60 40
0.02
20 0
0.015
-20 -40
0.01
-60 -80 0
5
10
15
20
Geomagnetic Latitude(Degree)
Geomagnetic Latitude(Degree)
a) F1 80
80
0.015
60
0.01
40
0.005
20
0 -0.005
0
-0.01
-20
-0.015
-40
-0.02
-60
-0.025
-80 0
Local Time(Hour)
5
10
15
Local Time(Hour)
7
20
d) F4 0.03
60 40
0.02
20
0.01
0
0
-20
-0.01
-40 -60
-0.02
-80 0
5
10
15
0.04
80
Geomagnetic Latitude(Degree)
Geomagnetic Latitude(Degree)
c) F3 80
60
0.03
40
0.02
20
0.01
0
0
-20
-0.01
-40
-0.02
-60
-0.03
-80
20
0
5
Local Time(Hour)
10
15
20
Local Time(Hour)
Figure 1. The distribution of the first four EOF base functions through EOF decomposition of PEC as a function of geomagnetic latitude and local time. Figure 2 shows the distribution of the first four principal component coefficients of the PEC after EOF decomposition. It can be seen from the figure that the coefficient corresponding to the first-order EOF exhibits a semi-annual cycle variation, and the coefficient corresponding to the second-order EOF shows annual cycle and semi-annual cycle variation, and the coefficient value is much smaller than the coefficient value corresponding to the first-order EOF, the coefficient corresponding to the third-order EOF exhibits an annual cycle variation, and the coefficient corresponding to the fourth-order EOF exhibits a semi-annual cycle and a seasonal cycle variation. a) P1
600 550
15
500 450
10
400 350
5
300
20 Universal Time(Hour)
Universal Time(Hour)
20
2010
2011
2012
2014
2015
50
15
0
10
-50
5
-100
250
0 2008
b) P2
650
0
-150
2008
2010
Universal Time(Year)
2011
2012
2014
2015
Universal Time(Year)
c) P3
d) P4 40
20
20
60
15
0 -20
10
-40
5
Universal Time(Hour)
Universal Time(Hour)
20 40
15
20
10
0 -20
5
-60 -40
0 2008
0
-80
2010
2011
2012
2014
2015
2008
Universal Time(Year)
2010
2011
2012
Universal Time(Year)
8
2014
2015
Figure 2. The distribution of the first four EOF coefficients through EOF decomposition of PEC To better represent the variation of solar activity, the Mg II index is used (http://www.iup.uni-bremen.de/gome/gomemgii.html). The emission core of the Mg II doublet (280 nm) exhibits the largest natural solar irradiance variability above 240 nm. It is frequently used as a proxy for spectral solar irradiance variability from the UV to EUV associated with the 11-yr solar cycle (22-yr magnetic cycle) and solar rotation (27d). Since the solar radio flux at 10.7 cm ( F10.7 ) might not represent the solar EUV radiation well during solar minimum (Viereck et al., 2001). Figure 3 shows the long-period variation of composite MgII index and Mg II smooth during 2008-2014. The Mg II index is used to represent the variation of solar activity. Mg II smooth is the Mg II mean with a smooth window of 81 days. It can be seen from the comparison between figure 2 and figure 3 that the coefficient of the first-order EOF exhibits a significant solar cycle change, which is consistent with the variation trend of Mg II index characterizing the solar activity intensity. For example, during the quiet period of solar activity from 2008 to 2011, Mg II index is basically near 0.15 and change little. At this time, the coefficient value of the first-order EOF function is also very small around 250. After 2011, Mg II index shows an increase, and in 2011 and 2014, the peaks of these years appear, and the coefficients of the first-order EOF function show a consistent phenomenon. 0.17
180 Mg II index
Mg II smooth
160
0.165
140
0.16
120
0.155
100
0.15
80
0.145
60 2008
2009
2010
2011
2012
2013
2014
Mg II index
F10.7
F10.7
0.14 2015
Universal Time (Year)
Figure 3. Long-term variation of Mg II index (orange line) 、 Mg II smooth (purple line) and F10.7 (blue line) during 2008–2014. Table 2 lists the correlation coefficients between the first fourth-order EOF series and the solar activity indices. It can be seen from Table 2 that the first-order EOF 9
sequence has a strong positive correlation with the solar activity index Mg II ,and its correlation coefficient reaches 0.877. This phenomenon indicates that the overall change of the PEC is mainly caused by the solar extreme ultraviolet radiation (solar extreme ultraviolet radiation is also the main driving force for the overall change of the electron content in plasmasphere). Another significant phenomenon is that there is a negative correlation between the EOF sequences in Table 2 and solar activity index. This may be because the response of the electron content in the plasmasphere to the solar activity is not strictly synchronized. That is, there may be a time delay between the driving force and the response of the PEC, so there is a negative correlation between them. The time delay in the response of electron content to the solar activity also occurs in the ionosphere (Mao et al., 2008). Table 2. Correlation coefficients of the first four orders of the EOF series with the solar indices Mg II index
Mg II smooth
1st
0.877
0.903
2nd
-0.481
-0.404
3rd
-0.570
-0.677
4th
0.330
0.343
The analysis above shows that the EOF decomposition method can decompose the overall changes of the PEC dataset into different components and their relative contributions generated by the corresponding driving forces and mechanisms. However, it should be noted that the experiment was carried out under the assumption of quiet geomagnetic conditions, the EOF decomposition method used in the experiment is only applicable to the quiet period of geomagnetic activity. Pk UT , d is obtained by the first EOF decomposition. Here, the feature vector Pk UT , d is decomposed a second time by the EOF decomposition method. The decomposition process is as follows: M
Pk (UT , d ) BFkn (UT ) Akn (d )
(3)
n 1
where, BFkn UT is the n-order basis function obtained by EOF decomposition of
Pk , and varies with UT.
Akn d
is the n-order correlation coefficient 10
corresponding to BFkn UT . Akn d denotes the long term variation of Pk (UT , d ) .
M represents the order of basis functions and correlation coefficients.
2000
20
1000
0
2010
2011
2012
0
2010
500
50
0
2014
-100 2015
2012
2014
200
150 A41
0
2011
2012
A33
-50 2015
Coefficients A41
A32
Coefficients A32,A33
Coefficients A31
A31
2010
2011
Universal Time (Year)
Universal Time (Year)
-500 2008
A23
0
-500 2008
-20 2015
2014
A22
A42
100
100
0
50
0 2008
A43
2010
2011
2012
2014
Coefficients A42,A43
0 2008
100 A21
A13
Coefficients A21
A12
Coefficients A12,A13
Coefficients A11
500
40 A11
Coefficients A22,A23
3000
-100 2015
Universal Time (Year)
Universal Time (Year)
Figure 4. Distribution of the associated coefficients acquired from the second EOF decomposition of PEC Figure 4 shows the distribution of the first four orders correlation coefficient
Akn d
obtained by the second EOF decomposition. It can be seen from the figure
that the coefficient
Akn d
exhibits solar activity period, annual period and
semi-annual period, which is also related to solar activity. Moreover, we use the Fast Fourier transform (FFT) to extract the period of the correlation coefficient
Akn d
obtained by the second EOF decomposition. The result shows that the correlation coefficient
Akn d
exhibits obvious annual, semi-annual and seasonal periodicity.
Consequently, according to the analysis above, the correlation coefficient
Akn d
can be fitted by periodic function to represent its periodic variations, and the
solar activity index Mg II and its 81 days smooth value are employed as model parameters to reflect the impact of solar activity on the PEC. The periodic function is 11
as follows: Akn d Vkn1 d Vkn 2 d Vkn3 d Vkn 4 d Vkn5 d Vkn 6 d Vkn 7 d
V
d k n1
V
d kn 2
V
d kn3
V
d kn 4
V
d kn5
V
d kn6
(4)
M g II d c M g II _ s m o o th d a b k n1 k n1 k n1 d M g II d * M g II _ s m o o th ( d ) k n1
M g II d c M g II _ s m o o th d b a kn 2 kn 2 kn 2 s in d M g II d * M g II _ s m o o th ( d ) kn 2
M g II d c M g II _ s m o o th d b a kn3 kn3 kn3 cos d M g II d * M g II _ s m o o th ( d ) kn3
b M g II ( d ) c M g II _ s m o o th d a kn 4 kn 4 kn 4 s in 2 d M g II d * M g II _ s m o o th ( d ) kn 4
M g II ( d ) c M g II _ s m o o th d b a kn5 kn5 kn5 c o s 2 d M g II d * M g II _ s m o o th ( d ) kn5 M g II ( d ) c M g II _ s m o o th d b a kn6 kn6 kn6 s in 4 d M g II d * M g II _ s m o o th ( d ) kn6
M g II ( d ) c M g II _ s m o o th d b a kn7 kn7 c o s 4 V d kn7 kn7 d M g II d * M g II _ s m o o th ( d ) kn7
where, Vkn1 d corresponds to the constant term, correspond to the annual cycle, semi-annual cycle,
Vkn 6 d
Vkn4 d
and
Vkn 5 d
Vkn 2 d
and
(5)
Vkn 3 d
correspond to the
and Vkn 7 d correspond to the 1/4-year period, ε refers
doy to residual error of the model, 2 365.25 , MgII represents solar radiation
level, and Mg II smooth is the Mg II mean with a smooth window of 81 days. The parameters akn , bkn , ckn , d kn
are all calculated by using the least squares
parameter estimation using equations (4).
The Mg II index is based on a spectral
structure near 280 nm and can be used to measure solar radiation levels. After the model is established, the coefficient
of the secondary EOF decomposition
Akn d
can be obtained from MgII, MgII smooth index at that moment. Combining
BFkn U T
with equation (2), the first EOF decomposition coefficient Pk can be
determined. Finally, the first decomposed EOF basis function is used to determine the global PEC value at the moment under the condition of solar activity. Figure 5 shows the comparisons between the observed and fitted coefficients of the secondary EOF decomposition. It can be seen from the figure that the fitting 12
performance is good for A11, A21, A31, A33, this suggests that the variation of the A11, A21, A31, A33 are regular and stable with time and solar activity. In addition, the more periodic of the associated coefficients acquired from the second EOF decomposition of PEC, the better the model fitting effect is. 2500
1000 500
A21
A11
2000 1500
0
1000 500 2008
2010
2011
2012
2014
2015
-500 2008
Universal Time (Year)
2011
2012
2014
2015
Universal Time (Year) 10
200
5
A31
A33
0 -200 -400 2008
2010
0 -5
2010
2011
2012
2014
-10 2008
2015
2010
2011
2012
2014
2015
Universal Time (Year)
Universal Time (Year)
Figure 5. The time series of the observed (blue dots) and fitting (red line) coefficients of the secondary EOF decomposition 3. Model Validation To analyze the effectiveness and reliability of the model and verify its accuracy and stability, two kinds of data sources are used here: the first is observational data which is used in modeling, and the other is observational data which is not used in modeling. To test the effectiveness of the EOF model, we examined the relation between the modeled values and observational data (these data are involved in modeling). Fig.6 shows a comparative scatter plot between modeled PEC and observed PEC in 2008, 2010, 2011, 2012, 2014. It can be seen from the figure that the correlation coefficient between modeled PEC and observed PEC in 2008, 2010, 2011, 2012, 2014 are all above 0.97, that is, modeled PEC has strong correlation with observed PEC. In addition, as can be seen from the figure, the RMS maximum of the difference between modeled PEC and observed PEC in 2008, 2010, 2011, 2012, 2014 is 0.61 TECU, the RMS minimum of the difference is 0.37 TECU, and the average values of the differences are basically 0.03 TECU, indicating that the established PECmodel has high precision and reliability. 13
Figure 6. Scatter plots of observed PEC against the corresponding modeled PEC in 2008, 2010, 2011, 2012 and 2014 To verify the accuracy and reliability of the global plasmaspheric electron content model constructed by using twice empirical orthogonal function decomposition method, the observational data in 2009 and 2013 which is not involved in modeling is used to verify the accuracy and reliability of the model. Figure 7 shows a comparative scatter plot between modeled PEC and observed PEC in 2009 and 2013. Where the observed PEC is the abscissa, the modeled PEC is the ordinate, and the slope of the red line is equal to 1. It can be seen from the figure that the correlation coefficient between modeled PEC and observed PEC in 2009 and 2013 are all above 0.96, that is, modeled PEC has strong correlation with observed PEC. At the same time, it can be seen that the correlation coefficient R is lower in the high solar activity year 2013 than in the low solar activity 2009, not only because the 14
modeling data we used did not cover the entire solar cycle for 11 years and most of them were data of low solar activity, but also the influences of the different geomagnetic activity. In 2013, the solar activity was intense, so the simulation effect of the high solar activity year was not as good as the low solar activity. In addition, as can be seen from the figure, the RMS of the difference between modeled PEC and observed PEC in 2009 is 0.45 TECU, while the RMS of the difference between modeled PEC and observed PEC in 2013 is 0.70 TECU. A possible reason for this phenomenon is that in the high solar activity year 2013, the solar activity index fluctuates greatly with respect to other years which can be seen from Fig.3, resulting in a low degree of model fit and greater RMS in the high solar activity year.
Figure 7. Scatter plots of observed PEC against the corresponding modeled PEC for 2009 and 2013 In order to better illustrate the temporal and spatial variation characteristics of the EOF model, the global distribution of PEC observations and simulations in spring equinox, summer solstice, autumn equinox and winter solstice of low solar activity year 2009 and high solar activity year 2013 are given in figure 8 - 9, respectively. It can be seen from figure that the EOF model can well characterize the change of PEC whether under low solar activity or high solar activity conditions, which indicates that
7
50
6 5
0
4 3
-50
2 1
0
5
10
15
b) Vernal Equinox
7
50
6 5
0
4 3
-50
20
2 1
0
Local Time (Hour)
5
10
15
Local Time (Hour)
15
20
PEC (TECU)
8
Geomagnetic Latitude (Degree)
a) Vernal Equinox
PEC (TECU)
Geomagnetic Latitude (Degree)
the model is reasonable.
5 4
0
3 2
-50
1
10
15
5 4
0 3 2
-50
1
0
20
5
6
50
5 4
0
3 2
-50
1
5
10
15
20
Geomagnetic Latitude (Degree)
7
PEC (TECU)
Geomagnetic Latitude (Degree)
e) Autumnal Equinox
0
6 5
0
4 3
-50
2 1
0
5
0 4 3
-50
2 1
Geomagnetic Latitude (Degree)
6
PEC (TECU)
Geomagnetic Latitude (Degree)
7
15
7
50
0
8
50
10
20
5
10
15
20
Local Time (Hour)
g) Winter Solstice
5
15
f) Autumnal Equinox
Local Time (Hour)
0
10
Local Time (Hour)
Local Time (Hour)
PEC (TECU)
5
6
50
8
h) Winter Solstice
7
50
6 5
0
4
PEC (TECU)
0
d) Summer Solstice
PEC (TECU)
50
Geomagnetic Latitude (Degree)
6
PEC (TECU)
Geomagnetic Latitude (Degree)
7
c) Summer Solstice
3
-50
2 1
0
20
5
10
15
20
Local Time (Hour)
Local Time (Hour)
Figure 8. Global PEC maps calculated by EOF model (right) and observed (left) at 12:00 UT for vernal equinox, summer solstice, autumnal equinox, and winter solstice in 2009
16
8 7
0
6 5 4
-50
3 2
b) Vernal Equinox 10
50
8
0 6
-50
PEC (TECU)
9
50
Geomagnetic Latitude (Degree)
10
PEC (TECU)
Geomagnetic Latitude (Degree)
a) Vernal Equinox
4
2
15
20
0
5
15
20
c) Summer Solstice
10 9
50
8 7
0
6 5 4
-50
3
Geomagnetic Latitude (Degree)
Local Time (Hour)
PEC (TECU)
Geomagnetic Latitude (Degree)
Local Time (Hour)
10
10
d) Summer Solstice
9
50
8 7 6
0
5 4
-50
3 2
2
0
5
10
15
20
0
5
10 9
50
8 7
0
6 5 4
-50
3 2
10
15
Geomagnetic Latitude (Degree)
e) Autumnal Equinox
5
15
20
f) Autumnal Equinox 10
50
8
0 6
4
-50
2
0
20
5
10
15
20
Local Time (Hour)
g) Winter Solstice
14
50
10
0
8 6
-50
4
PEC (TECU)
12
Geomagnetic Latitude (Degree)
Local Time (Hour)
Geomagnetic Latitude (Degree)
10
Local Time (Hour)
PEC (TECU)
Geomagnetic Latitude (Degree)
Local Time (Hour)
0
h) Winter Solstice 14
50
12 10
0
8 6
-50
4 2
2
0
5
10
15
PEC (TECU)
10
PEC (TECU)
5
0
20
PEC (TECU)
0
5
10
15
20
Local Time (Hour)
Local Time (Hour)
Figure 9. Global PEC maps calculated by EOF model (right) and observed (left) at 12:00 UT for vernal equinox, summer solstice, autumnal equinox, and winter solstice 17
in 2013 To further illustrate the global modeling effect of the EOF model, the global distribution of the differences between the PEC model value and the observed value for spring equinox, summer solstice, autumn equinox and winter solstice in 2009 and 2013 are given in figure 10 and 11, respectively. It can be seen from Figure 10 that in the low solar activity year 2009, the differences between the PEC model value and the observed value are kept within ±0.5 TECU in most regions of the world, and the overall difference does not exceed ±2 TECU. It can be seen from Figure 11 that in the high solar activity year 2013 , the differences between the PEC model value and the observed value are maintained at ±0.5 TECU in most regions of the world, and the overall difference does not exceed ±2.5 TECU. It is further illustrated that the model changes
of 1
a) Vernal Equinox 50
0.5
0
0
-0.5
-50
PEC
on
a
global
scale.
b) Summer Solstice 0.5
50 0 -0.5
0
-1
-50
PEC (TECU)
the
Geomagnetic Latitude (Degree)
simulate
PEC (TECU)
Geomagnetic Latitude (Degree)
can
-1.5
-1
5
10
15
20
0
5
1
50
0.5 0
0 -0.5
-50
-1 -1.5
5
10
15
Geomagnetic Latitude (Degree)
c) Autumnal Equinox
0
10
15
20
Local Time (Hour)
PEC (TECU)
Geomagnetic Latitude (Degree)
Local Time (Hour)
d) Winter Solstice
1
50 0.5
0
0
-0.5
-50
20
-1
0
Local Time (Hour)
5
10
15
20
Local Time (Hour)
Figure 10. Difference between modeled and observed PEC values at 12:00 UT for vernal equinox, summer solstice, autumnal equinox, and winter solstice in 2009
18
PEC (TECU)
0
0.5
0
0 -0.5
-50
-1 -1.5
5
10
15
1.5
50
1 0.5
0
0 -0.5
-50 -1 -1.5
20
0
5
50
1
0.5
0
0
-50 -0.5
10
15
Geomagnetic Latitude (Degree)
1.5
c) Autumnal Equinox
5
10
15
20
Local Time (Hour)
PEC (TECU)
Geomagnetic Latitude (Degree)
Local Time (Hour)
0
2
d) Winter Solstice
1.5 1
50
0.5 0
0
-0.5
-50
-1 -1.5
0
20
PEC (TECU)
0
b) Summer Solstice
PEC (TECU)
1
50
Geomagnetic Latitude (Degree)
1.5
PEC (TECU)
Geomagnetic Latitude (Degree)
a) Vernal Equinox
5
10
15
20
Local Time (Hour)
Local Time (Hour)
Figure 11. Difference between modeled and observed PEC values at 12:00 UT for vernal equinox, summer solstice, autumnal equinox, and winter solstice in 2013. 4. Conclusions This paper uses global plasmaspheric total electron content product for 1446 days in 2008, 2010, 2011, 2012, 2014 and EOF decomposition to establish a global PEC model. The main findings and conclusions are summarized as follows: (1). Modeling with EOF decomposition, using only the first few EOF basis function sequences can characterize the overall variation characteristics of PEC. This method has the characteristics of simple method and fast convergence. At the same time, the EOF decomposition method can separate the global PEC from space and time into different parts according to different driving forces or mechanisms, so each basis function and its related coefficients have certain physical meanings. (2). Exploring the relationship between PEC and solar activity, it is found that solar activity is the main driving force of PEC. (3). Using the observation data to verify the model, it is found that the model has global adaptability, under both high solar activity and low solar activity conditions, which is in good agreement with the observation data, indicating that modeling based on EOF is reasonable. 19
(4). In addition, it should be noted that the accuracy of the global plasmaspheric empirical model established by EOF in this paper can be further improvement. One of the reasons is that the original global PEC product used for modeling has limited accuracy and time resolution (only COSMIC and MetOp-A satellites are used and time resolution is only 4 hours). With the launch of more low-orbit satellites, such as China's Hongyan and Hongyun satellite systems, which will make more than 100 low-orbit satellites in operation, allowing for more abundant GNSS observations to establish a global plasmaspheric empirical model with higher precision and higher time resolution. Acknowledgments: We gratefully acknowledge http://www.iup.uni-bremen.de/gome/gomemgii.html for providing composite Mg II index for analyazing the relationship between various EOF sequences and solar activity index. This study was funded by the National Natural Science Foundation of China (41404031) and Outstanding Youth Science Fund of Xi’an University of Science and Technology (2018YQ2-10). References: A, E., Zhang, D.H., Xiao, Z., Hao, Y.Q., Ridley, A.J., Moldwin, M., 2011. Modeling ionospheric for f2 by using empirical orthogonal function analysis. Ann Geophys. 29 (8), 1501-1515. Belehaki, A., Jakowski, N., Reinisch, B.W., 2004. Plasmaspheric electron content derived from GPS TEC and digisonde ionograms. Adv. Space Res. 33 (6), 833-837. Berube, D., Moldwin, M.B., Fung, S.F., Green, J.L., 2004. A plasmaspheric mass density model and constraints on its heavy ion concentration. J. Geophys. Res. Space Phys. 110 (A4). Carpenter, D.L., 1963. Whistler measurements of electron density and magnetic field strength in the remote magnetosphere. J. Geophys. Res. 68 (12), 3727-3730. Carpenter, D.L., Anderson, R.R., 1992. An ISEE/Whistler model of equatorial electron density in the magnetosphere. J. Geophys. Res. Space Phys. 97 (A2), 1097-1108. Chen, Z., Zhang, S.R., Coster, A.J., Fang, G., 2015. EOF analysis and modeling of GPS TEC climatology over North America. J. Geophys. Res. Space Phys. 120 (4), 3118-3129. Chen, P., Yao, Y., Li, Q., Yao, W., 2017. Modeling the plasmasphere based on LEO 20
satellites onboard GPS measurements. J. Geophys. Res. Space Phys. 122 (1), 1221-1233. Chugunin, D.V., Kotova, G.A., Klimenko, M.V., Klimenko, V.V., 2017. Longitudinal dependence of the H+ concentration distribution in the plasmasphere according to INTERBALL-1 satellite data. Cosmic Res. 55 (6), 457-463. Cherniak, I., Zakharenkova, I., 2018. Evaluation of the IRI-2016 and NeQuick electron content specification by COSMIC GPS radio occultation, ground-based GPS and Jason-2 joint altimeter/GPS observations. Adv. Space Res. 63 (2019), 1845-1859. Dunteman, G.H., 1989. Principal components analysis. Sage. 930(2), 527-547. Daniell, R.E., Brown, L.D., Anderson, D.N., Fox, M.W., Doherty, P.H., Decker, D.T., et al., 1995. Parameterized ionospheric model: A global ionospheric parameterization based on first principles models. Radio Sci. 30 (5), 1499-1510. Ercha, A., Zhang, D., Ridley, A.J., Xiao, Z., Hao, Y., 2012. A global model: Empirical orthogonal function analysis of total electron content 1999-2009 data. J. Geophys. Res. Space Phys. 117 (A3). Gringauz, K.I., Bezrokikh, V.V., Ozerov, V.D., Rybchinskii, R.E., 1960. A Study of the Interplanetary Ionized Gas, High-Energy Electrons and Corpuscular Radiation from the Sun by Means of the Three-Electrode Trap for Charged Particles on the Second Soviet Cosmic Rocket. Soviet Physics Doklady, 5, 361. Geiss, J., Young, D.T., 1981. Production and transport of O++ in the ionosphere and plasmasphere. J. Geophys. Res. Space Phys. 86 (A6), 4739-4750. Horwitz, J.L., Comfort, R.H., Richards, P.G., Chandler, M.O., Chappell, C.R., Anderson, P., et al., 1990. Plasmasphere-ionosphere coupling: 2. Ion composition measurements at plasmaspheric and ionospheric altitudes and comparison with modeling results. J. Geophys. Res. 95 (A6), 7949. Hannachi, A., Jolliffe, I.T., Stephenson, D.B., 2010. Empirical orthogonal functions and related techniques in atmospheric science: A review. Int J Climatol. 27 (9), 1119-1152. Jolliffe, I.T., 2002. Principal Component Analysis. J Marketing Res. 87 (100), 513. Jakowski, N., Hoque, M.M., 2018. A new electron density model of the plasmasphere for operational applications and services. J Space Weather Spac. 8, A16. Pearson, K., 1900. On Lines and Planes of Closest Fit to Points in Space. Philos Mag. 2(11), 559-572. Khazanov, G.V., Rasmussen, C.E., Konikov, Y.V., Gombosi, T.I., Nagy, A.F., 1994. 21
Effect of magnetospheric convection on thermal plasma in the inner magnetosphere. J. Geophys. Res. Space Phys. 99 (A4), 5923-5934. Kotova, G.A., Verigin, M.I., Bezrukikh, V.V., 2015. Physics‐based reconstruction of the 3‐D density distribution in the entire quiet time plasmasphere from measurements along a single pass of an orbiter. J. Geophys. Res. Space Phys. 120 (9), 7512-7521. Klimenko, M.V., Klimenko, V.V., Zakharenkova, I.E., Cherniak, I.V., 2015. The global morphology of the plasmaspheric electron content during Northern winter 2009 based on GPS/COSMIC observation and GSM TIP model results. Adv. Space Res. 55 (8), 2077-2085. Lorenz, E.N., 1956. Empirical orthogonal functions and statistical weather prediction. Liu, C., Zhang, M.L., Wan, W., Liu, L., Ning, B., 2008. Modeling M (3000) F2 based on empirical orthogonal function analysis method. Radio Sci. 43 (01), 1-8. Mao, T., Wan, W.X., Liu, L.B., 2005. An EOF Based Empirical Model Of Tec Over Wuhan. Chinese J Geophys-Ch. 48 (4), 827-834. Mao, T., Wan, W., Yue, X., Sun, L., Zhao, B., Guo, J., 2008. An empirical orthogonal function model of total electron content over China. Radio Sci. 43 (2). Mukhtarov, P., Andonov, B., Pancheva, D., 2013. Global empirical model of TEC response to geomagnetic activity. J. Geophys. Res. Space Phys. 118 (10), 6666-6685. O'Brien, T.P., Moldwin, M.B., 2003. Empirical plasmapause models from magnetic indices. Geophys Res Lett. 30 (30), 1-1. Preisendorfer RW, Mobley CD., 1988. Principal component analysis in meteorology and oceanography. Viereck, R., Puga, L., Mcmullin, D., Judge, D., Weber, M., Tobiska, W.K., 2001. The Mg II index: A proxy for solar EUV. Geophys Res Lett. 28 (7), 1343–1346. Webb, P.A., Essex, E.A., 2000. An ionosphere-plasmasphere global electron density model. Phys Chem Earth. Part C: Solar, Terrestrial & Planetary Science, 25 (4), 301-306. Zhao, B., Wan, W., Liu, L., Yue, X., Venkatraman, S., 2005. Statistical characteristics of the total ion density in the topside ionosphere during the period 1996-2004 using empirical orthogonal function (EOF) analysis. Ann Geophys. 23 (12), 3615-3631. Zhang, Q.M., Wang, C., Li, H., Li, C.Q., 2013. Plasma transport between the ionosphere and plasmasphere in response to solar wind dynamic pressure 22
pulsation. Chinese Sci Bull. 58 (33), 4126-4132.
23