*Article* **Long-Term Variations of Plasmaspheric Total Electron Content from Topside GPS Observations on LEO Satellites**

**Shuanggen Jin 1,2,\*, Chao Gao 1,3, Liangliang Yuan 1, Peng Guo 1, Andres Calabia 2, Haibing Ruan <sup>2</sup> and Peng Luo 1,4**

	- Technology, Nanjing 210044, China; andres@calabia.com (A.C.); rhb@nuist.edu.cn (H.R.)

**Abstract:** The plasmasphere is located above the ionosphere with low-energy plasma, which is an important component of the solar-terrestrial space environment. As the link between the ionosphere and the magnetosphere, the plasmasphere plays an important role in the coupling process. Therefore, it is of great significance to study the electron content variation of the plasmasphere for the solarterrestrial space environment. Nowadays, the topside global positioning system (GPS) observations on Low Earth Orbit (LEO) satellites provide a unique opportunity to estimate and study variations in the plasmasphere. In this paper, the plasmaspheric total electron content (PTEC) is estimated, and its long-term variations are studied from topside GPS observations onboard the Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC). The PTEC in the daytime is higher than that in the nighttime, with the peak between 14:00 and 17:00 in the magnetic local time, while the minimum value of PTEC in the belt appears between 3:00 and 6:00 in the magnetic local time before sunrise. For seasonal variations, the PTEC is the highest in spring of the northern hemisphere and the lowest in summer of the northern hemisphere regardless of the state of the solar activity. The long-term variation in PTEC is further analyzed using 11-year COSMIC GPS observation data from 2007 to 2017. A high correlation between PTEC and the F10.7 indices is found. Particularly in the geomagnetic high-latitude region during the daytime, the correlation coefficient reaches 0.93. The worst case occurs during the nighttime in the geomagnetic middle-latitude region, but the correlation coefficient is still higher than 0.88. The long-term variations of plasmaspheric TEC are mainly related to the solar activity.

**Keywords:** plasmasphere; PTEC; GPS; GCPM; F10.7 index

#### **1. Introduction**

With the continuous exploration into deep space and the increasing variety of electromagnetic applications, such as communication and navigation, monitoring and understanding of the solar-terrestrial space environment have become a hot field, including the Earth's neutral atmosphere, ionosphere, plasmasphere, magnetosphere, and so on [1]. The plasmasphere is a part of magnetosphere, also called the inner magnetosphere [2], which starts from the top of the ionosphere and ends at the plasmapause. The plasmasphere is a donut-shaped region surrounding the Earth, containing the coldest plasma of the magnetosphere [3]. It is currently believed that the charged particles in the plasmasphere mainly come from escape of the ionosphere and capture from the solar wind [4,5]. Richards et al. [6] examined the relative importance of ionospheric and thermospheric densities and temperatures in producing the annual variation of the plasmaspheric electron density. Lee et al. [7] compared the global plasmaspheric total electron content (TEC) with

**Citation:** Jin, S.; Gao, C.; Yuan, L.; Guo, P.; Calabia, A.; Ruan, H.; Luo, P. Long-Term Variations of Plasmaspheric Total Electron Content from Topside GPS Observations on LEO Satellites. *Remote Sens.* **2021**, *13*, 545. https://doi.org/10.3390/ rs13040545

Academic Editor: Roberta Giuliani Received: 7 January 2021 Accepted: 29 January 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the ionospheric TEC simultaneously measured by Jason-1 satellite during the declining phase of solar cycle 23, and the results showed that the plasmaspheric density structures fundamentally followed the ionosphere, but there were still significant differences.

Radio signals are refracted by the charged particles, which affects satellite navigation, positioning and microwave remote sensing. When the navigation signals of Global Navigation Satellite System (GNSS) pass through the Earth's ionosphere and plasmasphere, they are delayed due to the refraction. The magnitude of the impact is related to the total electron content (TEC) of the signal path [8]. Although the electron density of the plasmasphere is much lower than that of the ionosphere, the region covered by the plasmasphere is dozens of times larger than that covered by the ionosphere. Therefore, the electron content of the plasmasphere accounts for a considerable proportion of the total electron content, which is usually about 10% in the daytime, but can reach 60% in the nighttime [9,10].

In some practical applications, for example, when using a single frequency GPS receiver for navigation and positioning, it is impossible to eliminate the effects of charged particles by ionosphere-free combined observations, while ionospheric empirical models are not precise enough to eliminate the error caused by such delay. Current main ionospheric models, such as the International Reference Ionosphere (IRI) model and the NeQuick model, can only estimate the electron content of the ionosphere, but ignore the plasmaspheric TEC (PTEC). Therefore, the corresponding delay effect cannot be estimated or corrected precisely, which has an impact on the final positioning results [11,12]. Therefore, it is important to estimate the PTEC for the delay correction. Furthermore, the coupling processes between the plasmasphere and the ionosphere are very complex. The studies on the plasmaspheric electron content variations and dynamic coupling processes between the plasmasphere and the ionosphere are crucial for understanding the plasmasphere [13].

Before the advent of GNSS radio occultation technology, the PTEC was generally difficult to measure directly. The TEC acquired by ground GNSS receivers is the total electron content of the ionosphere and the plasmasphere. The ionospheric TEC (ITEC) can be obtained from the ionosonde or incoherent scattering radar (ISR), and then PTEC is indirectly calculated by subtracting ITEC from the total TEC [14,15]. However, there are a series of problems with these methods. Firstly, normally the electron density profile below the peak value of F2 layer can be obtained by the ionosonde, and the electron density profile above the peak value is extrapolated by the Chapman function [16]. Secondly, the number of observation stations of ionosonde and ISR is relatively small, and the distribution is very sparse. In addition, the cost is a little high, which leads to limited coverages in the global ionospheric observations. Furthermore, there are likely systematic deviations between different observation techniques and methods, which will be involved into the PTEC estimation. Therefore, it is challenging to accurately estimate PTEC and establish the plasmaspheric model.

Nowadays, with the increasing number of GNSS Radio Occultation observations on Low Earth Orbit (LEO) satellites, the topside GNSS observations on LEO satellites provide a unique opportunity to directly estimate PTEC and study its variations in the plasmasphere, particularly Constellation Observing System for Meteorology, Ionosphere, and Climate (COSMIC) with six LEO satellites. The COSMIC can provide more than 2500 occultation events per day during the normal operation period of six LEO satellites [17,18]. In this paper, the PTEC is estimated, and its long-term variations are studied based on topside GPS observations on COSMIC with providing the slant TEC (sTEC) of the signal path. The sTEC is transformed into vertical TEC (vTEC) by a mapping function, and then the grid model of PTEC is established. The spatial and temporal distribution characteristics of PTEC are analyzed as well as its long-time variation characteristics from January 2007 to December 2017. In Section 2, the data and method are introduced, the results and analysis are presented in Section 3, and finally, conclusions are given in Section 4.

#### **2. Data and Method**

The data used in this paper are the precise podTec provided by COSMIC from January 2007 to December 2017, which can be obtained from the COSMIC data analysis and archiving center (https://www.cosmoc.ucar.edu/cdaac/). The PodTec files provide UTC time, three-dimensional coordinates of LEO and GPS satellites, observation elevation of the GPS-LEO observation link at LEO satellite, and the slant TEC on the signal path. It is worth noting that the hardware delays of the transmitters on GPS satellites and the receivers on COSMIC satellites have been deducted from the TEC, and the sampling rate of the observations is 1 s [19]. Since the volume of observation data after 2017 is too small, we only use the observation data from 2007 to 2017 in this paper. In addition, to ensure the consistency of the altitude region covered by the observations, the observation data before LEO satellites lifting their orbits are also eliminated.

To estimate the plasmaspheric TEC, it is necessary to set a thin shell, and the electron content is assumed to be concentrated on the shell. The method of the gridded plasmaspheric TEC model is basically the same as that of the Global Ionosphere Model (GIM) [20,21]. We tested the effects of the thin shell heights on PTEC results, and found a small difference and little effect on the temporal and spatial distribution of PTEC. Therefore, we set the altitude of the thin shell at 1400 km. Detailed plasmaspheric TEC modeling is shown in Figure 1.

**Figure 1.** Diagram of topside GPS observations of LEO satellite and the single layer plasmasphere hypothesis.

Firstly, the baseline between the LEO and GPS satellites is transformed to the stationcentered coordinate system, and then the azimuth angle *A* and the elevation angle *h* are calculated as:

$$X\_{NEII} = \begin{bmatrix} -\sin B \cos L & -\sin B \sin L & \cos B \\ -\sin L & \cos L & 0 \\ \cos B \cos L & \cos B \sin L & \sin B \end{bmatrix} X\_{XYZ} \tag{1}$$

$$A = \arctan(\frac{X\_E}{X\_N}) \tag{2}$$

$$h = \arctan(\frac{X\_{ll}}{\sqrt{X\_N^2 + X\_E^2}}) \tag{3}$$

where *B* and *L* are the geographic latitude and longitude of the LEO satellite, respectively, and *XXYZ* and *XNEU* are the coordinates of the LEO-GPS baseline in Earth-centered Earthfixed (ECEF) coordinate system and the station-centered coordinate system, respectively.

To reduce the mapping errors and multipath effects, only the observations with an elevation angle greater than 40◦ are used to establish the plasmaspheric TEC grid model. The sTEC observations with negative values or over 100 TECU, which can be considered to be unreasonable observations, are also removed. Then, we calculate the obliquity factor *z* between LEO and GPS satellites and the vertical TEC:

$$z = \arcsin\left[\left(R\_r \cosh\right) / \left(R\_E + H\_I\right)\right] \tag{4}$$

$$vTEC = sTEC \cdot \cos(z) \tag{5}$$

where *Rr* is the distance between the receiver on the LEO satellite and the Earth's center, *RE* is the radius of the Earth, *HI* is the altitude of the single layer plasmasphere (here we set it as 1400 km). Then, the geographic longitude and latitude of the plasmaspheric piercing point can be calculated as follows:

$$
\Psi\_I = \pi/2 - h - z \tag{6}
$$

$$\varphi\_I = \arcsin(\sin B \cos \Psi\_I + \cos B \sin \Psi\_I \cos A) \tag{7}$$

$$
\lambda\_I = L + \arcsin(\sin \Psi\_I \sin A / \cos \varphi\_I) \tag{8}
$$

where Ψ*<sup>I</sup>* is the geocentric angle between LEO satellite and the piercing point, and *ϕ<sup>I</sup>* and *λ<sup>I</sup>* are the geographic latitude and longitude of the piercing point, respectively.

The geographic longitude and latitude are converted to geomagnetic longitude and latitude, and the magnetic local time is calculated as follows:

$$m\varphi\_I = \arcsin(\sin(\varphi\_I)\sin(b\_0) + \cos(\varphi\_I)\cos(b\_0)\cos(l\_0 - 1))\tag{9}$$

$$m\lambda\_I = \arctan\left(\frac{\cos(\varphi\_I)\sin(l\_0 - 1)/\cos(m\varphi\_I)}{(\sin(b\_0)\sin(m\varphi\_I) - \sin(\varphi\_I))/(\cos(b\_0)\cos(m\varphi\_I))}\right) \tag{10}$$

$$mLT\_I = UT\_I + (m\lambda\_I - l\_0) / \left(15^{\circ} \times \pi / 180^{\circ}\right) \tag{11}$$

where *mϕ<sup>I</sup>* and *mλ<sup>I</sup>* are the geomagnetic latitude and longitude of the piercing point, respectively, *b*<sup>0</sup> and *l*<sup>0</sup> are the geographic latitude and longitude of the geomagnetic north pole, respectively, and *b*<sup>0</sup> = 80.0◦ , and *<sup>l</sup>*<sup>0</sup> = −72.2◦ , *UTI* and *mLTI* are the universal time and magnetic local time of the observation, respectively.

Finally, we divide all the observations in a month or a season into groups with a latitudinal resolution of 2.5◦ and a temporal resolution of 20 min, and the observations in each group are weighted and averaged according to the value of 1 + *cos*2*h* <sup>−</sup><sup>1</sup> as the PTEC of the corresponding grid point. In all the formulas above, angles are in radians and distances are in kilometers.

#### **3. Results and Analysis**

The COSMIC constellation consists of six LEO satellites, which provide dense global coverage plamaspheric observations. Figure 2 shows the geographical distribution of the piercing points on the 1400 km thin shell on 2 January 2008. Due to the inclination of the satellite orbits and the lowest observation elevation angle of 40◦, the distribution of the piercing points has gaps at the north and south poles. However, in the region between 72◦S and 72◦N, the topside observations of LEO satellites are well-distributed, and all observations can be regarded as observations from GPS satellite altitude (about 20,200 km) to COSMIC satellite altitude (about 800 km). Although this altitude range is not exactly consistent with the real plasmasphere, the observations are homogeneous in the detection altitude, and basically contain most of the charged particles of the plasmasphere, so they

can be regarded as the observations for the plasmasphere. These conditions are quite favorable for plasmaspheric modeling.

**Figure 2.** The distribution of COSMIC piercing points at the 1400 km thin shell on 2 January 2008

#### *3.1. PTEC Estimation from COSMIC*

According to the definition of plasmasphere, the distribution of charged particles is dominated by the Earth's magnetic field. Therefore, the geomagnetic coordinate system is used in modeling. The magnitude of PTEC is closely related to the position of the sun, so the geomagnetic longitude is converted to the magnetic local time. On the basis of topside GPS observations of COSMIC, the plasmaspheric TEC grid model is estimated. This paper mainly analyzes the temporal and spatial distribution and long-term variation of PTEC, considering the volume of data and the convenience of analysis, so we divide the observations separately by month and model. In this way, the effects of long-period variations like solar activity on PTEC are preserved, while the influence of geomagnetic activities, such as magnetic storms and sub-magnetic storms, which are relatively shortlived (from a few hours to one or two days), is averaged over one month's observations, with a minimal impact on modeling. When analyzing the seasonal variation of PTEC, we combine the observations of three months together and re-weight the observations to calculate the PTEC.

The Global Core Plasma Model (GCPM) is the first real global plasmaspheric model, and was established by Gallagher et al. [22] by integrating density distribution models of different regions. It covers the ionosphere, plasmasphere, plasmapause, plasmaspheric poles, and so on. In the ionosphere, GCPM adopts the international reference ionosphere model IRI. The plasmaspheric region is based on the density distribution model of H+ established according to the observations of DE-1 satellite by Gallagher et al. [22]. The Persoon model is adopted in the polar regions [23]. These regional models are integrated by mathematical fitting to create a static three-dimensional plasmaspheric model GCPM, which can extend from the ionosphere to 8 to 9 radii of the Earth.

Figure 3 shows the comparison of PTEC from the topside GPS observations of COSMIC in January 2008 with the GCPM and the differences in the bottom panel. The blank regions in the top and bottom panels are caused by the fact that the observations of COSMIC cannot completely cover the polar region. In January 2008, the state of solar activity was quiet, and the F10.7 indices exhibit little change, at less than 80 sfu. This that the observations of PTEC from COSMIC and GCPM are almost consistent with respect to the overall characteristics. There is a significant belt with higher values of PTEC at geomagnetic latitudes between −45◦ and 45◦, and PTEC values in the daytime are higher than those in the nighttime. In addition, the PTEC values in the top panel decrease slowly from the geomagnetic equatorial region to the geomagnetic polar regions, while in the middle panel, the PTEC values decrease rapidly in the geomagnetic middle-latitude regions. As a result, the differences show significant zonal belt distribution in the bottom panel, and in the middle and high geomagnetic latitudinal regions.

**Figure 3.** Comparison of PTEC from topside GPS observations of COSMIC and GCPM in January 2008 (The altitude range is from about 800 km to 20,200 km): (**a**) PTEC from COSMIC, (**b**) PTEC from GCPM, and (**c**) difference of PTEC.

The left panel in Figure 4 shows the comparison of the PTEC values from topside GPS observations of COSMIC and GCPM in January 2008 at corresponding positions. The correlation coefficient of PTEC values is 0.85, which indicates a good consistency between the PTEC from COSMIC observations and GCPM. The right panel shows the distribution histogram of PTEC differences. Almost all the differences are within ±4 TECU, and the numbers of differences over ±3 TECU are less than 5% of the total statistics, which also shows that the two PTEC results are in good agreement with each other.

In general, it has a relatively high correlation of PTEC between GCPM and the COS-MIC observations, and the correlation is higher in quiet period of solar activity. However, since GCPM is a model fitted by mathematical formulas, the result is very smooth in numerical value and therefore the details cannot be seen from the model. The PTEC estimation from topside GPS observations of COSMIC is based on the actual observations, which contains more rich details.

**Figure 4.** Comparison and difference distribution of PTEC from topside GPS observations of COSMIC and GCPM in January 2008.

#### *3.2. Temporal-Spatial Distribution Characteristics of PTEC*

The long-term variations of plasmaspheric TEC and its temporal-spatial distribution characteristics are analyzed using COSMIC-derived PTEC. According to the solar activity, we selected 2008 and 2014 as the representatives of low and high solar activity years, respectively, and established the PTEC gridded model using COSMIC observations by month and season, and the monthly and seasonal variations characteristics of PTEC under different states of solar activity are analyzed.

Figures 5 and 6 show monthly and seasonal variations of plasmaspheric TEC in 2008, respectively. Observations of January and February 2009 were also used in the subgraph in the bottom right panel of Figure 6. The PTEC values in different months or seasons have the same following basic characteristics: PTEC values at daytime are higher than those at nighttime; PTEC values in lower geomagnetic latitudinal regions are higher than those in higher geomagnetic latitudinal regions; and there are obvious zonal belts with higher PTEC values within the ±45◦ geomagnetic latitudinal region. This is because the solar incidence angle in the geomagnetic low-latitude region is the greatest in the daytime, where the plasmasphere captures the most energy, and thus generates more charged particles through ionization. These phenomena can also prove a close relationship between the plasmaspheric TEC and solar activity, which will be analyzed in the next section.

In the high PTEC value belts within ±45◦ geomagnetic latitude, the peak values of PTEC in monthly and seasonal models all appear in the geomagnetic equatorial region between 14:00 to 17:00 o'clock in the magnetic local time, while the minimum values of PTEC appear between 3:00 and 6:00 o'clock in the magnetic local time. This can be explained by the coupling process between the plasmasphere and the ionosphere. The charged particles in the ionosphere drift upward along the Earth's magnetic field lines to the plasmasphere in the daytime, while the charged particles in the plasmasphere will return to the ionosphere to maintain the electron density of the F layer in the nighttime, so the electron content of plasmasphere will reach the minimum value before sunrise [24]. In Figure 5, we can see that values of plasmaspheric TEC in June, July and August 2008 are significantly lower than those in other months, and values of plasmaspheric TEC in March and November are the highest. At the seasonal scale, the plasmaspheric TEC is the highest in the northern hemisphere spring and the lowest in northern summer, which are related to the variation of the vertical radiation region of the sun. The seasonal variation of the plasmaspheric TEC is consistent with the variation of the ionospheric TEC due to the strong coupling interaction [25].

**Figure 5.** Monthly variation of plasmaspheric TEC in 2008.

**Figure 6.** Seasonal variation of plasmaspheric TEC in 2008.

To analyze the influence of magnetic local time and geomagnetic latitude on plasmaspheric TEC, we divided the monthly PTEC into six parts. There are daytime and nighttime regions: the magnetic local time from 6:00 to 18:00 o'clock is the daytime, and the nighttime is from 18:00 to the second day's 6:00 o'clock. As for geomagnetic latitude, it

is divided into low-latitude, mid-latitude and high-latitude, with boundaries of ±30◦ and ±60◦. The average PTEC of each region was calculated, and is shown in Figure 7. Since the solar and geomagnetic activities were quite calm in 2008, the variations of plasmaspheric TEC in different regions were also small and gentle, with a maximum variation range of about 1.2 TECU. Apparently, the maximum values of plasmaspheric TEC always appeared in the low-latitude region in the daytime, and the maximum value was in November 2008, reaching 6.2 TECU. The sub-maximum values were in the nighttime low-latitude region, with a maximum value of 5 TECU in November. The differences between daytime and nighttime PTEC in the low-latitude region are about 1 TECU. In mid-latitude region, the plasmaspheric TEC is around 3 TECU, and there is a small difference between daytime and nighttime. However, the relative differences of plasmaspheric TEC between daytime and nighttime are quite large for the small value of PTEC in high-latitude region. In general, the effect of geomagnetic latitude on plasmaspheric TEC is more obvious.

**Figure 7.** Monthly mean daytime and nighttime PTEC in different geomagnetic latitudes in 2008.

Figures 8 and 9 show the monthly and seasonal variations of plasmaspheric TEC in 2014, respectively. Observations of January and February 2015 are also used in in the subgraph at bottom right of Figure 9. The basic distribution characteristics of plasmaspheric TEC mentioned above still exist. The most obvious difference is that the belts with higher PTEC values are wider during the solar active period, especially in the daytime, which indicates that during the solar active period, a larger region of the plasmasphere can receive strong solar radiation, thus ionizing to generate more charged particles. In terms of numerical value, the peak values of PTEC are significantly higher than those in 2008, which are almost double in some months. In June, July and August 2014, the values of plasmaspheric TEC are smaller than those in other months, which is the same as 2008. This phenomenon is also reflected on the seasonal scale, whereby the values of plasmaspheric TEC in northern summer are much lower than those in other seasons. An interesting phenomenon is that the local maximum values of PTEC appear around 12:00 o'clock in the region of geomagnetic latitude -80◦ in spring, autumn and winter of the northern hemisphere, which indicates that the charged particles of the plasmasphere will accumulate in this region at noon, and its physical mechanism needs to be further studied.

**Figure 8.** Monthly variation of plasmaspheric TEC in 2014.

**Figure 9.** Seasonal variation of plasmaspheric TEC in 2014.

In the same way, we also analyzed the influence of magnetic local time and geomagnetic latitude on plasmaspheric TEC in 2014, during which year the solar activity was very active. An obvious feature is that the maximum values of plasmaspheric TEC in different regions all appeared in March, and the maximum reached 13.7 TECU in the

daytime low-latitude region (Figure 10). The variations of plasmaspheric TEC in different regions were relatively large, with a maximum variation range of about 5.6 TECU, and the relative variations were also greater than those in 2008. The same characteristics as 2008 were not repeated here. An obvious distinction is that values of plasmaspheric TEC in the daytime high-latitude region were higher than those in the nighttime mid-latitude region, except for January, June and July, which shows that the high-latitude region was more affected in intense periods of solar activity.

**Figure 10.** Monthly mean daytime and nighttime PTEC in different geomagnetic latitudes in 2014.

#### *3.3. Correlation with Solar Activity*

The plasmasphere is located above the Earth's ionosphere with much lower particle density, but it is greatly affected by the solar radiation, which leads to the complex characteristics of plasmaspheric electron density variation. To study the relationship between plasmaspheric TEC and solar activity, we chose the F10.7 index as the reference indicator. There is a strong correlation between the 10.7 cm radio flux and the number or area of sunspots. At present, as one of the most important indices of solar activity, the F10.7 index has been widely used in space weather research and related studies of ionosphere and magnetosphere [26].

Figure 11 shows the monthly mean F10.7 indices and monthly mean PTEC from 2007 to 2017, where the red line represents the monthly mean F10.7 indices, corresponding to the left ordinate, and the blue line represents the plasmaspheric TEC, corresponding to the right ordinate. Before 2011, the monthly mean F10.7 indices were relatively small, and the variation was relatively gentle. During this period, the solar activity was relatively calm. After January 2011, the monthly mean F10.7 index rose sharply, and the variation was very intense, indicating that the solar activity was in an active period. Then, after 2015, the F10.7 index began to decline, and the solar activity decreased accordingly. The variation of monthly mean PTEC (the arithmetic average of PTEC at each grid point of the monthly plasmaspheric model) from 2007 to 2017 is shown in Figure 11, and the standard deviation of each monthly mean PTEC is also calculated and shown with an error bar. Most of the standard deviations are less than 2.5 TECU, which indicates that the monthly mean PTEC is of significance. The values of monthly mean PTEC were also low before 2011, basically no more than 5 TECU, and the variation was relatively gentle. After 2011, the monthly mean PTEC also began to rise, and then began to decline after 2015, which was basically consistent with the variation of the monthly mean F10.7 indices.

**Figure 11.** Monthly mean F10.7 index and monthly mean PTEC from January 2007 to December 2017.

To determine the correlation between PTEC and solar activity in the daytime and nighttime, we divided the plasmaspheric TEC into the daytime region and nighttime region, and the division standard was the same as in the previous section. Similarly, the values of PTEC in the daytime and nighttime were averaged by month, respectively. The mean values and the standard deviations are shown in Figure 12. The variations of PTEC in the daytime and nighttime are basically synchronous with differences less than 2 TECU, and both show consistency with the variation of the monthly mean F10.7 indices.

**Figure 12.** Monthly mean F10.7 index and monthly mean daytime and nighttime PTEC from January 2007 to December 2017.

To further estimate the correlation between PTEC and the F10.7 index, we divided the monthly plasmaspheric TEC into latitudinal regions according to the geomagnetic latitude and the magnetic local time, and the criteria of division are the same as mentioned before. The monthly mean PTEC for each geomagnetic latitudinal region was calculated for daytime and nighttime. Then, they were counted with the corresponding monthly mean F10.7 indices in different subgraphs in Figure 13. There are strong correlations between the monthly mean F10.7 indices and the monthly mean daytime and nighttime PTEC in each latitudinal region. In the same geomagnetic latitudinal region, the correlation between the monthly mean PTEC and the monthly mean F10.7 indices in the daytime is higher than that in the nighttime. This is because the plasmasphere in the daytime is directly exposed to the sun and can receive solar radiation and the energy directly, while in the nighttime, the plasmasphere is more affected by the ionosphere, which reduces the correlation with the solar activity. In the term of geomagnetic latitudinal region, the correlation in the geomagnetic high-latitude region is also the highest, while the geomagnetic middle-latitude region has the lowest correlation. However, among these correlation coefficients, the lowest one is still more than 0.88, indicating a very strong relationship between plasmaspheric TEC and the solar activity.

**Figure 13.** Correlations between monthly mean F10.7 index and monthly mean daytime and nighttime PTEC in different geomagnetic latitude regions

#### **4. Conclusions**

Using the topside GPS observations on COSMIC, the long-term plasmaspheric total electron content was obtained. By comparison with the GCPM, the plasmaspheric TEC from the topside observations of COSMIC was verified. Using the observation data in the solar minimum year 2008 and the solar maximum year 2014, PTEC was estimated at the monthly and seasonal scales, respectively, and its temporal-spatial distribution characteristics under different states of the solar activity were analyzed. The PTEC was mainly distributed in a belt region around the Earth within ±45◦ of the geomagnetic latitude. The plasmaspheric TEC in the daytime is higher than that in the nighttime, which reaches a peak between 14:00 and 17:00 in magnetic local time, while the minimum value of PTEC in the belt appears between 3:00 and 6:00 in magnetic local time before sunrise. For seasonal variations, the plasmaspheric TEC is the highest in the spring of the northern hemisphere and the lowest in the summer of the northern hemisphere, regardless of the state of solar activity. The variations of monthly mean PTEC in different regions were quite gentle during the solar minimum year 2008, while dramatical changes were found during the solar maximum year 2014. Furthermore, the long-term variations of the 11 year plasmaspheric TEC were analyzed with the F10.7 index as the reference indicator. The results showed a strong correlation between plasmaspheric TEC and solar activity.

**Author Contributions:** Conceptualization, S.J., C.G. and L.Y.; methodology, S.J., and C.G; software, C.G.; validation, S.J., C.G. and L.Y.; formal analysis, C.G.; investigation, C.G., and L.Y.; data curation, C.G., and P.L.; writing—original draft preparation, S.J., C.G., and P.L.; writing—review and editing, A.C., H.R., and P.G.; visualization, S.J.; supervision, S.J.; project administration, S.J.; funding acquisition, S.J. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (No. 12073012 and No. 41761134092), Shanghai Leading Talent Project (Grant No. E056061), Jiangsu Province Distinguished Professor Project (Grant No. R2018T20), and Startup Foundation for Introducing Talent of NUIST (Grant No.2243141801036).

**Data Availability Statement:** The data presented in this study are available from the corresponding website.

**Acknowledgments:** We thank the CDAAC (COSMIC Data Analysis and Archive Center) for providing the podTec observation data and SEPC (Space Environment Prediction Center) for providing the F10.7 indices.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


**Mingliang Wu 1,2, Shuanggen Jin 2,3, Zhicai Li 4,\*, Yunchang Cao 5, Fan Ping <sup>6</sup> and Xu Tang <sup>3</sup>**

	- Technology, Nanjing 210044, China; xu.tang@nuist.edu.cn

**Abstract:** The Global Navigation Satellite System (GNSS) plays an important role in retrieving high temporal–spatial resolution precipitable water vapor (PWV) and its applications. The weighted mean temperature (*Tm*) is a key parameter for the GNSS PWV estimation, which acts as the conversion factor from the zenith wet delay (ZWD) to the PWV. The *Tm* is determined by the air pressure and water vapor pressure, while it is not available nearby most GNSS stations. The empirical formular is often applied for the GNSS station surface temperature (*Ts*) but has a lower accuracy. In this paper, the temporal and spatial distribution characteristics of the coefficients of the linear *Tm*-*Ts* model are analyzed, and then a piecewise-linear *Tm*-*Ts* relationship is established for each GPS station using radiosonde data collected from 2011 to 2019. The *Tm* accuracy was increased by more than 10% and 20% for 86 and 52 radiosonde stations, respectively. The PWV time series at 377 GNSS stations from the infrastructure construction of national geodetic datum modernization and Crustal Movement Observation Network of China (CMONC) were further obtained from the GPS observations and meteorological data from 2011 to 2019. The PWV accuracy was improved when compared with the Bevis model. Furthermore, the daily and monthly average values, long-term trend, and its change characteristics of the PWV were analyzed using the high-precision inversion model. The results showed that the averaged PWV was higher in Central-Eastern China and Southern China and lower in Northwest China, Northeast China, and North China. The PWV is increasing in most parts of China, while the some PWVs in North China show a downward trend.

**Keywords:** GPS meteorology; weighted mean temperature; precipitable water vapor; radiosonde

#### **1. Introduction**

Water vapor is an important part of the Earth's hydrosphere and plays a key role in the energy exchange and water cycle in nature. The atmospheric water vapor content is limited by local temperature and pressure and is closely related to the formation of various precipitation, such as clouds, rain, and snow [1]. Accurate measurements of water vapor and its distribution changes have become one of the basic problems in synoptics, weather forecasting, and climate research [2–5]. One of the indicators to measure the amount of atmospheric water vapor is the precipitable water vapor (PWV), which represents a certain height of the water column produced by the condensation of all tropospheric atmospheric water vapor in the column per unit bottom area at any time into liquid water. Since the demand for real-time and accurate weather services is becoming more and more urgent, the traditional detection technologies such as radiosondes, water vapor radiometers, and

**Citation:** Wu, M.; Jin, S.; Li, Z.; Cao, Y.; Ping, F.; Tang, X. High-Precision GNSS PWV and Its Variation Characteristics in China Based on Individual Station Meteorological Data. *Remote Sens.* **2021**, *13*, 1296. https://doi.org/10.3390/rs13071296

Academic Editor: Stefania Bonafoni

Received: 24 February 2021 Accepted: 25 March 2021 Published: 29 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

solar photometers cannot meet the application requirements for continuous high-precision, high-temporal resolution monitoring of water vapor.

Nowadays, Global Navigation Satellite System (GNSS) technology has become an important means to obtain precipitable water vapor with high temporal and spatial resolution, and many scholars have used GNSS technology to study global or regional climate change and characteristics [6–12]. Bevis et al. first proposed the concept of GPS meteorology and obtained the global surface temperature and weighted mean temperature from 8718 radiosonde stations in North America and explained the linear relationship coefficient and specific process of ground-based GPS inversion of water vapor [13]. Ross and Rosenfeld used 23 years of profile data from 53 radiosonde stations from the National Center for Atmospheric Research to calculate the global *Tm*-*Ts* linear coefficient and demonstrated that it was related to the station and the season except the equatorial region [14]. Several *Tm*-*Ts* conversion models were established for different regions and different seasons using limited meteorological observation data [15–18]. Yao et al. established the nonlinear transformation relationship by combining mathematical statistics and derivation, which can improve the accuracy of fitting in China [19]. However, many GPS stations are not equipped with temperature and pressure sensors due to cost and other reasons. Hence, it is impossible to use transformation coefficient from real observations to calculate *Tm*.

Many researchers used the National Centers for Environmental Prediction (NCEP), and the European Center for Medium-term Weather Prediction (ECMWF) reanalyzed data and interpolated data from the ground meteorological observation station to obtain ground temperature and air pressure for ground-based GNSS PWV inversion [20–25]. Yao et al. established a global *Tm* model and a tropospheric delay model using spherical harmonic functions with considering annual, semi-annual, and diurnal changes, which greatly improved the accuracy of the PWV estimation [26,27]. However, the empirical models have relatively low accuracy. The research on *Tm*-*Ts* models in the region of China was mostly focused on climate zones and seasonal divisions. Since many GPS stations in China lacked meteorological data in the past, they cannot obtain precise PWV, as well as investigate their long-term variation characteristics of GPS PWV.

In this study, we collected more radiosonde stations observations at co-located GNSS stations from the infrastructure construction of national geodetic datum modernization and Crustal Movement Observation Network of China (CMONC). The temporal and spatial distribution characteristics of the *Tm*-*Ts* model were analyzed from radiosonde data, and a piecewise-linear model was established at each radiosonde station in China. With this model, the nine-year GPS PWV time series at GNSS stations was obtained from CMONC GNSS observations and meteorological data. The precision of the PWV was evaluated, and the distribution characteristics and their changes were investigated. Section 2 shows the data and methods, evaluation and comparison are presented in Section 3, variations in the characteristics of GNSS PWV are presented in Section 4, and finally, the conclusions are given in Section 5.

#### **2. Data and Methods**

#### *2.1. Observation Data*

The infrastructure construction of national geodetic datum modernization in China was launched in 2012 and completed in 2017, which contained 210 GPS stations. A highprecision, dynamic, and unified modern surveying and mapping datum system was established to provide coordinate frame services, including data products, real-time positioning, processing and analysis, and other services [28]. Crustal Movement Observation Network of China (CMONC) was established in 2006 and completed in 2012, which contained 360 GPS stations [29,30]. The National Geomatics Center of China (NGCC) provided the hourly ZTD data and 1-h measured temperature and atmospheric pressure data at these GPS stations from 2011 to 2019. We used BERNESE 5.2 software [31] to process the raw data, including the implementation of daily solutions and adjustments [32].

The Integrated Global Radiosonde Archive (IGRA) has released radiosonde and pilot balloon observations with more than 2700 stations around the world since 1905, including air pressure, temperature, geopotential height, and relative humidity (ftp://ftp.ncdc.noaa. gov/pub/data/igra, accessed on 24 January 2021), whose temporal resolution is 12 h. We used IGRA-released radiosonde profiles in China from 2011 to 2019 to calculate the *Tm* at each radiosonde station. Figure 1 shows all the used GPS stations and radiosonde stations.

**Figure 1.** Distribution of GPS stations from the infrastructure construction of national geodetic datum modernization and Crustal Movement Observation Network of China (CMONC) in mainland China. The green circle is the radiosonde station, and the red triangle is the GPS station.

#### *2.2. Establishment of Site-Specific Piecewise-Linear Tm-Ts Relationship*

*Tm* is related to the temperature and vapor pressure at different altitudes in the atmosphere, which can be obtained from the IGRA. The method to calculating the 12-h *Tm* can be expressed as [13]:

$$T\_m = \frac{\int\_{z\_0}^{+\infty} \left(\frac{c}{T}\right) dz}{\int\_{z\_0}^{+\infty} \left(\frac{c}{T^2}\right) dz} \approx \frac{\sum\_{i=1}^N \left(\frac{c\_i}{T\_i}\right) \Delta z\_i}{\sum\_{i=1}^N \left(\frac{c\_i}{T\_i^2}\right) \Delta z\_i} \tag{1}$$

where *e* (hPa) refers to the water vapor pressure, *T* (K) is the corresponding temperature, and *Z* (m) denotes the starting height of integration. Δ*zi* (m) denotes the altitude of the *i*th atmospheric layer, *N* denotes the number of the atmospheric layer, *ei* (hPa) denotes the water vapor pressure of *i*th atmospheric layer, and *Ti* (K) denotes the temperature of *i*th atmospheric layer. We converted the geopotential height into the geoid height in the calculation process.

The empirical formula based on the long-term radiosonde data in the study area and the linear relationship between surface temperature *Ts* and *Tm* can be established by a regression analysis as:

$$T\_m = a \cdot T\_s + b \tag{2}$$

where *a* and *b* are the linear regression equation parameters.

The *Tm* variation has a significant annual variation. In some studies, the half-year variations and the daily variations are also considered when building the model. In order to get the time-varying characteristics of the *Tm*-*Ts* coefficient, a piecewise-linear least squares fitting is performed for each radiosonde station for one month.

#### *2.3. PWV from Site-Specific Piecewise-Linear Tm-Ts Relationship*

Three hundred and seventy-seven GPS stations were selected for estimating the PWV time series. Since the cubic spline method has a higher accuracy than the nearest-neighbor interpolation method and linear interpolation method, we used the cubic spline method to obtain the *Tm*-*Ts* model coefficients *a* and *b* at the GPS stations in the corresponding time period.

ZTD is the sum of Zenith Hydrostatic Delay (ZHD) and Zenith Wet Delay (ZWD), as

$$\text{ZTD} = \text{ZHD} + \text{ZWD} \tag{3}$$

The Saastamoinent model has been widely used for ZHD computation [33], as

$$\text{ZHD} = 0.002277 \cdot \frac{P}{1 - 0.0026 \cdot \cos\left(2\phi\right) - 0.00028 \cdot h\_0} \tag{4}$$

where *P* (hPa) denotes the ground pressure of the GPS station, *φ* denotes the latitude of the GPS station, and *h*<sup>0</sup> (m) denotes the elevation of the GPS station.

The ZWD is caused by water vapor in the atmosphere under nonstatic equilibrium. Generally, empirical models and meteorological parameters at GPS station are used to obtain the ZHD, and then, the ZHD is deducted from the ZTD to obtain the ZWD.

$$\text{ZWD} = \text{ZTD} - \text{ZHD} \tag{5}$$

The linear relationship between the ZWD and PWV can be expressed as [13]

$$\text{PWV} = \Pi \cdot \text{ZWD} \tag{6}$$

$$\Pi = \frac{10^6}{\rho\_w \cdot \frac{R}{m\_w} \cdot \left[\frac{k\_3}{T\_m} + k\_2 - \frac{m\_w}{m\_d} \cdot k\_1\right]} \tag{7}$$

where Π is the conversion factors between the ZWD and PWV; Π is a function of *Tm*; *ρ<sup>w</sup>* represents the density of the liquid water; *R* is the universal gas constant and *R* = 8314 Pa·m3·K−1·kmol<sup>−</sup>1; *mw* represents the molar mass of water vapor and *mw* = 18.02 kg·kmol<sup>−</sup>1; *md* represents the molar mass of the dry atmosphere and *md* = 28.96 kg·kmol<sup>−</sup>1; and *<sup>k</sup>*1, *<sup>k</sup>*2, and *k*<sup>3</sup> are constants (*k*<sup>1</sup> = 77.604 ± 0.014 K/hPa, *k*<sup>2</sup> = 70. 4 K/hPa, and *k*<sup>3</sup> = (3.776 ± 0.014) × 105 K2/hPa) [13].

In general, approximately 6.7 mm of ZTD error will cause a PWV error of 1 mm. Hence, the ZTD with a Root Mean Square Error (RMSE) of greater than 6.7 mm are eliminated. In addition, ZTD data and meteorological parameters are missing at some stations in a few time periods. For sites with missing data for more than one year, they will not be used when analyzing the long-term changes. Then, Equations (3)–(6) were used for high-precision PWV estimations based on the site-specific piecewise-linear *Tm*-*Ts* relationship in China. Finally, the PWV time series derived from 377 GPS stations in China from 2011 to 2019 were obtained.

#### *2.4. PWV from Radiosonde*

The PWV at the radiosonde station is calculated as follows:

$$\text{PWV} = \int\_0^{p\_0} \frac{q}{\rho\_{w\mathcal{G}}} d\_{\mathcal{P}} \tag{8}$$

where *p* (hPa) denotes the atmospheric pressure, *p*<sup>0</sup> (hPa) is the ground pressure at the GPS station, *<sup>q</sup>* (g·kg−1) denotes the specific humidity, *<sup>g</sup>* (m·s−2) is the acceleration of gravity, and *<sup>ρ</sup><sup>w</sup>* (g·cm−3) refers to the density of the liquid water. By discretizing Equation (8), the integral equation of PWV from the ground to the top of the atmosphere is obtained as follows:

$$\text{PWV} = -\frac{1}{\mathcal{S}} \sum\_{p=0}^{p} q \times p \tag{9}$$

#### *2.5. Fitting Function of the PWV Time Series*

Jin and Luo [34] analyzed the PWV series derived from 155 globally distributed GPS sites observations and found that most of the periods of the PWV series were 1 year, 0.5 years, 1 day, and 0.5 days. To fit the PWV time series of 377 GPS stations in China, we established the following equation:

$$\begin{split} \text{PWV} &= k\_0 + k\_1 \cdot \cos\left(\frac{DOY - c\_1}{365.25} \cdot 2\pi\right) + k\_2 \cdot \cos\left(\frac{DOY - c\_2}{365.25} \cdot 4\pi\right) + k\_3 \cdot \cos\left(\frac{HOD - c\_3}{24} \cdot 2\pi\right) \\ &+ k\_4 \cdot \cos\left(\frac{HOD - c\_4}{24} \cdot 4\pi\right) + \varepsilon \end{split} \tag{10}$$

where *k*<sup>0</sup> is a constant term; *k*1, *k*2, *k*3, *k*4, *c*1, *c*2, *c*3, and *c*<sup>4</sup> are the amplitude and phase at the period (1 year, 0.5 years, 1 day, and 0.5 days); DOY is the day of year; HOD is the hour of day; and *ε* is the residual. The least square method was used to determine the unknown parameters in Equation (10) with the PWV time series.

#### **3. Evaluation and Comparison**

#### *3.1. Spatial Distribution and Time-Varying Characteristics of the Tm-Ts Coefficient*

Figure 2 shows the spatial distribution of the *Tm*-*Ts* relationship coefficients. The slope coefficient ranges from 0.5 to 0.7 in South China, around 0.7 in Northwest China, and about 0.8 in Central and Northeast China. As shown in the right panel, the intercept coefficient is from 80 to 140 in Southern China, about 60 in Central and Northwestern China, and approximately 20 in Northeastern China.

**Figure 2.** Distribution of the *Tm*-*Ts* fitting coefficients a and b at each station by *Tm* = *a* ∗ *Ts* + *b*. (**a**) The slope coefficient *a* and (**b**) the intercept coefficient *b*.

Figure 3 shows the distribution diagram of the slope coefficient with the elevation, latitude, and longitude. No evident correlation is found with the elevation and longitude, whereas there is strong positive correlation with the latitude, especially in low-latitude areas. Some studies have shown that when the same *Tm*-*Ts* coefficient is used at a global scale, it will cause different errors in the *Tm* of different latitudes [35,36]. Wang et al. found that the *Tm* derived from the Bevis *Tm*-*Ts* relationship has a cold bias in the tropics and subtropics and a warm bias in middle and high latitudes, and furthermore, the RMS was dominated by the mean bias rather than the random error. Therefore, the distribution of the *Tm*-*Ts* coefficients is mostly related to the latitude. The slope coefficient *a* is generally

less than 0.6 in low-latitude areas, which is consistent with the previous results of some global *Tm*-*Ts* models.

**Figure 3.** Distribution of the *Tm*-*Ts* model slope coefficient with the elevation (**a**), latitude (**b**), and longitude (**c**).

The monthly coefficients of the *Tm*-*Ts* model are obtained. The slope coefficients have obvious annual cycle at most stations. The time series of slope coefficients at four radiosonde stations are randomly selected and shown in Figure 4. The information of the four stations is shown in Table 1. The slope coefficients of the *Tm*-*Ts* model vary greatly with the time, from about 0.5–1. The regression slope and its changing tendency are smaller at the SIMAO station, which is due to the lower latitude.

**Figure 4.** Monthly slope coefficient of the *Tm*-*Ts* model from 2011 to 2019.

**Table 1.** Information about the 4 stations used in Figure 4.


#### *3.2. Comparison with Bevis Tm-Ts Relationship*

We calculated the root mean square error (RMSE) (in unit of K) and accuracy improvement (%) of *Tm* for 94 radiosonde station, which are shown in Figure 5. Compared with the Bevis model, the site-specific piecewise-linear model has a significant improvement in the regression accuracy at most stations due to the consideration of the temporal and spatial distributions of the *Tm*-*Ts* conversion coefficient. According to the statistics, the *Tm* accuracy with 86 radiosonde stations is increased by more than 10% and by more than 20% with 52 radiosonde stations, and the lowest is increased by 6% (station 58457, 30.23◦ N, 120.17◦ E), and the highest is increased by 69% (station 56691, 26.87◦ N, 104.28◦ E).

**Figure 5.** RMSE (K) (**a**) and accuracy improvement of the *Tm* (**b**) calculated by the site-specific piecewise-linear and Bevis *Tm*-*Ts* relationship.

Figure 6 shows the RMSE distribution of the Bevis model and the site-specific piecewiselinear model. As we can see, compared with the Bevis model, the single-station piecewiselinear model has a greater accuracy improvement in the southern, southwest, and northeastern regions by more than 30%. The increase in the central and northwestern regions is relatively small and approximately 15%.

**Figure 6.** RMSE (K) distribution of the *Tm* calculated by the site-specific piecewise-linear relationship (**a**), Bevis *Tm*-*Ts* relationship (**b**), and reduction of the RMSE (%) by the site-specific piecewise-linear *Tm*-*Ts* relationship (**c**).

#### *3.3. Comparison with GPS-Derived PWV and Radiosonde PWV*

The selected GPS stations are closer to the radiosonde, in which the horizontal distance is less than 10 km and the elevation difference is less than 100 m. The water vapor obtained by the integral of the radiosonde profile data was used to evaluate the ground-based GPS PWV based on the Bevis model and the site-specific piecewise-linear *Tm*-*Ts* model. We calculated the deviation bias (in unit of mm) and relative errors. For example, Figure 7 shows the results at GXHC station in 2018. The accuracy of the PWV is better based on the site-specific piecewise-linear *Tm*-*Ts* relationship when compared to the Bevis's model.

**Figure 7.** Bias (mm) (**a**) and relative error (%) (**b**) of the precipitable water vapor (PWV) calculated based on different *Tm*-*Ts* models at GXHC station in 2018. The blue dots are the PWV based on the Bevis model, and the red dots are the PWV based on the site-specific piecewise-linear model.

The distribution of water vapor in the atmosphere is not uniform. Thus, the groundbased GPS water vapor can get the average distribution of water vapor in each satellite signal direction. In addition, the layered meteorological parameter data of the radiosonde is not strictly vertical. These factors will bring errors into the PWV results. However, it still can be seen that the accuracy of the PWV is significantly improved based on site-specific piecewise-linear *Tm*-*Ts* relationship, especially when there are more atmospheric water vapors in the summer.

#### **4. Variations Characteristics of GNSS PWV**

#### *4.1. Spatial Distribution of PWV in China*

The annual averaged PWV at all GPS stations from 2011 to 2019 in China range from 0 to 48 mm. As shown in Figure 8, the annual averaged PWV in Central-Eastern China and Southern China are relatively high, reaching above 25 mm, while the annual averaged PWV in Northwest, Northeast, and Northern China regions are lower, below 15 mm. Southeast China has a low latitude and is close to the East China Sea and the South China Sea, which are subject to subtropical monsoons. Monsoons transport water vapor from the sea to these areas, resulting in the higher annual averaged PWV [37]. Northwest China has a relatively high latitude and is an inland region with a temperate continental climate, so the annual averaged PWV is lower.

The GPS stations are divided into four regions: the eastern central region (24.5◦–37◦ N, 105◦–123◦ E), the southern region (19◦–24. 5◦ N, 105◦–120◦ E), the northwestern region (42◦–49◦ N, 80◦–90◦ E), and the northeastern region (37◦–50◦ N, 110◦–130◦ E). Figure 9 shows the nine-year daily averaged PWV of these four regions. The variation characteristics of the PWV in each region are consistent throughout the year, with the peaks in June and

July and the droughts in January and December. Among them, the peak in Southern China appeared the earliest, which was affected by the south-to-north monsoon. Atmospheric circulation transported the water vapor to the Yangtze River Basin and then continued to transport it to other regions [38].

**Figure 8.** Annual averaged PWV (mm) in China from 2011 to 2019.

**Figure 9.** Nine-year averaged daily PWV (mm) in four regions of China from 2011 to 2019.

#### *4.2. Seasonal Variations of PWV in China*

Figures 10 and 11 show the annual and semiannual PWV variation amplitudes at 377 GPS sites. The spatial distribution of the annual PWV variations is similar to the annual average distribution of the PWV. Central China, Southern China, and the southeast coastal areas have higher annual PWV variation amplitudes, reaching about 15 mm, while the annual cycle amplitude in Northwest China is lower, below 10 mm. The semiannual PWV variation amplitudes are relatively small, about 3–9 mm. Among them, the semiannual PWV variation amplitudes are the highest in Central China and Southwest China, reaching

above 6 mm, while the semiannual PWV variation amplitudes in Southern China and Northwestern China are lower, below 5 mm.

**Figure 10.** Annual PWV variation amplitudes (mm) at 377 GPS sites.

**Figure 11.** Semiannual PWV variation amplitudes (mm) at 377 GPS sites.

#### *4.3. Long-Term Variation Trend of PWV in China*

Figure 12 shows the long-term variation trends of the PWV (mm/year) of all GPS stations from 2011 to 2019. It can be seen that the PWV has been increasing in most parts of China, while the PWV in Northeast China shows a downward trend. To understand the variations of the PWV in more detail, a monthly anomaly of the PWV is obtained by subtracting the monthly averaged PWV of every month from the mean value of the monthly averaged nine-year PWV in each region, which can be seen in Figure 13. The monthly anomaly of the PWV was mainly negative before 2015 and then turned negative. The PWV in Northwest China, Northeast China, and Northern China were relatively stable, and the monthly mean value changed little, which was related to the dry inland climate.

**Figure 12.** Long-term variation trend of the PWV. The red upward arrows (**a**) stand for the increase of the PWV variation trend (mm/year), and the green downward arrows (**b**) represent the decrease of the PWV variation trend (mm/year).

**Figure 13.** Monthly anomaly (mm) of the PWV in four regions of China.

#### **5. Discussion**

Most of the published *Tm* models that considered the temporal and spatial distributions of the relationship between the *Tm* and meteorological measurements have been empirical models on a global scale. Their accuracy remains to be verified in a specific area and a period of time. Here, we estimated the monthly coefficients at each station, and Table 2 shows the comparison between the piecewise-linear model and two recently published representative models: the time-varying global-gridded *Ts*–*Tm* model (TVGG) and neural network-based *Tm* model (NN). We tested our model with radiosonde data from 2011 to 2019, and the results showed that our model has the best accuracy.

**Table 2.** Statistics of the *Tm* estimates for different models.


The spatial distribution characteristics of the PWV are consistent with other studies in different years [39,40]. There are bimodal characteristics of the PWV in Southern China, and the formation mechanism of the bimodal characteristics remains to be further studied. The past results of some studies in China indicated that the PWV showed a downward trend from 1995 to 2012 [41–43], while the trends of this article are upward from 2011 to 2019. Most PWV in Central and Eastern China and Southern China show an upward trend and some PWV in North China is downward from 2011 to 2019. The monthly averaged PWV increased significantly in 2015, which was affected by the El Niño.

#### **6. Conclusions**

In this study, we analyzed the temporal and spatial distribution characteristics of the coefficients of the linear *Tm*-*Ts* model and showed that the distribution of *Tm*-*Ts* coefficients is mainly related to the latitude. The *Tm*-*Ts* conversion coefficient changes with the time and has an obvious annual cycle. Based on the spatial distribution and time-varying characteristics of the *Tm*-*Ts* coefficients, a site-specific piecewise-linear model was established. Compared with the Bevis model, this model reduced the *Tm* RMS by more than 20% for the most of the tested radiosondes. The accuracy of GPS PWV is better based on the site-specific piecewise-linear *Tm*-*Ts* relationship when compared to the Bevis model. Furthermore, the PWV time series at 377 GNSS stations were further obtained and analyzed from the GPS observations and meteorological data from 2011 to 2019. The results showed that the average PWV in Central and Eastern China and Southern China is higher, reaching more than 25 mm, while the average value is lower and below 15 mm in Northwest China, Northeast China, and North China. The PWV is increasing in most parts of China, while some PWV in North China show a downward trend.

**Author Contributions:** Conceptualization, S.J. and Z.L.; methodology, S.J. and M.W.; software, M.W.; validation, S.J., Z.L. and F.P.; formal analysis, M.W.; investigation, M.W. and S.J.; data curation, Z.L. and S.J.; writing—original draft preparation, M.W. and S.J.; writing—review and editing, F.P., Y.C. and X.T.; visualization, M.W.; supervision, S.J.; project administration, S.J.; and funding acquisition, S.J. and Z.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Strategic Priority Research Program Project of the Chinese Academy of Sciences (Grant No. XDA23040100), National Natural Science Foundation of China (NSFC) Project (Grant No. 12073012), and National Key Research and Development Program of China (No. 2016YFB0501405).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Radiosonde data are available from the Integrated Global Radiosonde Archive (IGRA) (ftp://ftp.ncdc.noaa.gov/pub/data/igra, accessed on 24 January 2021).

**Acknowledgments:** The authors are grateful to the infrastructure construction of the national geodetic datum modernization in China and Crustal Movement Observation Network of China (CMONC) for providing the GPS observation data.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

