Precipitable Water Vapor Converted from GNSS-ZTD and ERA5 Datasets for the Monitoring of Tropical Cyclones

Tropical cyclones (TCs) are intense rotating storm systems, characterized by strong winds, heavy humidity and storm surges, and usually cause an abnormal increase in regional water vapor during the period of their occurrence. Global navigation satellite system (GNSS) is well-established and initially mainly for positioning, navigation and timing. However, due to their global coverage and 24/7 operability under all weather conditions, GNSS has been widely used in meteorology such as the retrieval of precipitable water vapor (PWV). Since PWV time series obtained from a regional ground GNSS tracking stations network has high temporal and spatial resolutions, in this research, the temporal and spatial variation trends of regional PWVs during several TCs’ periods in Hong Kong were investigated, and a new method using high-resolution GNSS-PWV time series to estimate a TC’s movement was proposed. Test data were from 15 GNSS tracking stations equipped with meteorological sensors, and three GNSS stations without meteorological sensors but using meteorological data interpolated from ERA5 (fifth-generation reanalysis dataset of the European Centre for Medium-range Weather Forecasting) dataset, and six TC events occurred in 2017 were selected for cases studies. Moreover, for a better accuracy of GNSS-PWV, which is converted from the zenith wet delay of the GNSS signal multiplying a conversion factor that uses a weighted mean temperature ( $T_{m}$ ), a seasonal multi-factor $T_{m}$ model was established and tested. Results showed that the new $T_{m}$ model outperformed the previous regional $T_{m}$ model used in Hong Kong. In addition, based on the temporal variation trend and spatial distribution of the regional PWV, one can determine if a TC approaches the region. A geometric approach to estimating a TC’s movement was developed and its test results agreed well with the truth. These results suggest the potential of using GNSS signals to monitor the movements of TCs for short-time TC forecasting due to their high temporal resolution.


I. INTRODUCTION
Over recent years, increasing tropical cyclone (TC) events occurred frequently in the North Pacific, Indian Ocean, Southwest Pacific and North Atlantic, causing considerable damage in nearby cities [1]- [5]. Many countries, provinces and cities have set up meteorological departments, which take the early warning and forecasting of TCs as a key task. Water vapor, as a main component of the atmosphere, is a key factor The associate editor coordinating the review of this manuscript and approving it for publication was Gerardo Di Martino .
influencing TCs' development [6], [7]. The monitoring of regional water vapor variation is critical for the forecasting of the movement, intensity and rainfall of TCs [8]- [10].
Traditional techniques for the acquisition of atmospheric water vapor, such as infrared satellite remote sensing, radiosondes (RS) and microwave radiometers, have the disadvantages of insufficient temporal-resolution, high cost and susceptibility to prevailing weather conditions [11], [12]. This has limitation for detecting short-time extreme weather events. Since Bevis et al. [13] proposed the concept of global navigation satellite system (GNSS) meteorology, the GNSS has gradually become one of the most important data sources for retrieval of water vapor at a high accuracy, high temporal and spatial resolution under all weather conditions. The most commonly used two GNSS data processing approaches are the precise point positioning (PPP) and double difference network (DDN) to estimate the unknown parameters, which usually include the position and zenith tropospheric delay (ZTD) over a GNSS station. The precipitable water vapor (PWV) are then converted from the ZTD. Compared to the DDN approach, the PPP method does not need additional GNSS stations, which is more effective for the estimation of the PWV time series over a single GNSS station [14]. This is the main reason why the PPP technique is often adopted for retrieving PWV when highly accurate orbits and satellite clock products are available. Many previous studies have shown that there were no significant differences in the quality of the PWVs resulting from these two approaches, and the accuracy of GNSS-PWV is in the range from 1 mm to 3 mm, compared with those obtained from RS and water vapor radiometers [15]- [18]. At present, some studies proved that the variation in GNSS-PWV could reflect a typhoon's life and so GNSS-PWV can be used to forecast short-term rainfall during the typhoon's period [19], [20]. In addition, some researchers proposed to use PWV time series over several GNSS ground tracking stations in a region to estimate the movements of a super typhoon [21], [22]. Many results were based on only one TC case, especially for super typhoons [23]- [25], which might not sufficiently reveal the true characteristics of water vapor or may not be applicable to other cases. To address this problem, in this study six TC events occurred in 2017 in the Hong Kong region were investigated, and the temporal and spatial variation trends of PWV over the region during the TCs' periods were investigated. Moreover, a new method using GNSS-PWV time series at several GNSS tracking stations in a region to estimate the TCs' movement was developed.
From GNSS data processing, the unknown parameter that is relevant to water vapor and directly estimated is the zenith wet delay (ZWD) (or ZTD, which will be explained later) over the GNSS station. Then the ZWD is converted into PWV using a conversion factor that is a function of the variables of ground temperature, pressure measurements, and weighted mean temperature (T m ) over the station. The first two variables' values can be measured by meteorological sensors, if the GNSS station is equipped with such sensors; while T m needs to be calculated based on a temperature measurement profile. However, not all the GNSS tracking stations in the regional GNSS tracking network in Hong Kong were equipped with meteorological sensors. Thus, ground temperature and pressure measurements cannot be acquired at those stations without meteorological sensors equipped. In this case, the values interpolated from the ERA5 (fifth-generation reanalysis dataset of the European Centre for Medium-range Weather Forecasting) dataset were adopted, although the accuracy of these data would not be as desired as observed ones. Similarly, due to the unavailability of temperature profiles over all the GNSS stations, an empirical T m model is most commonly used. For a better accuracy, in this study a seasonal multi-factor T m model was developed based on local measurement data for Hong Kong. The accuracy of the converted GNSS-PWV over the station that was closest to the regional RS station was assessed by comparing the GNSS-PWV against the RS-derived PWV (RS-PWV). The characteristics of the seasonal variation in GNSS-PWV were also analysed, and the temporal and spatial variation trends shown in the GNSS-PWV over the region in the duration of the TCs' approaching were investigated as well. Finally, a method using GNSS-PWV time series at a 1-min temporal resolution over several stations along the direction of the TCs' path to determine the TCs' movement is proposed.

II. MATERIALS AND METHODOLOGIES
A. STUDY AREA AND DATA SOURCE Hong Kong is a coastal city in South China and its GNSS continuously operating reference stations (CORS) network, called Satellite Positioning Reference Station Network, consists of 18 evenly-distributed stations. Of the 18 stations, 15 stations were equipped with meteorological sensors. From the GNSS data processing, the estimated GNSS-ZWD time series over each station was converted into PWV time series and the latter was used to investigate the characteristics of water vapor variation in both temporal and spatial domains. Apart from the GNSS data from the above regional CORS network, a variety of other data were also needed for the case studies of TCs in 2017. More details for the datasets and products used in this study are as follows: (  The datasets (1) and (2) were used to estimate the ZTDs with a 1 min resolution over all the 18 GNSS stations and PWV time series with the same time resolution over the 15 GNSS stations that were equipped with meteorological sensors. The dataset (3) was used to establish a regional T m model in Hong Kong and estimate RS-PWV (with a 12h resolution) for the validation of the GNSS-PWV. The dataset (4) was used for obtaining ground temperature and pressure (with a 1h resolution) over the three GNSS stations that were not equipped with meteorological sensors, which were for the PWVs (with a 1h resolution) converted from the ZTDs for the three stations. The dataset (5) was used for analyzing the correlation between wind speed and PWV temporal variation during a TC's period. The dataset (6) was used to assess the accuracy of the GNSS-ZTD. The geodetic coordinates and elevations of the 18 GNSS stations, one RS station and one anemometer station were listed in Table 1, and Fig. 1 shows their geographical distribution.

B. PWV DERIVED FROM GNSS AND ERA-5 DATASET
When a GNSS signal propagates through the atmosphere, the GNSS signal delays. The so-called atmospheric delay can be divided into two parts: the ionospheric and tropospheric delays. From dual-frequency GNSS observations, the ionospheric delay can be canceled out using an ionosphere free (IF) linear combination of two observations on the two frequencies. Whilst the tropospheric delay cannot be canceled out in the same way. Instead, it is estimated, along with other unknown parameters, from GNSS data processing. For the purpose of reducing the number of unknown parameters to be estimated, all the slant tropospheric delays from different satellites observed are projected into the zenith direction, i.e. the ZTD, using a mapping function such as the Global Mapping Function (GMF), Neill Mapping Function (NWF) and Gridded Vienna Mapping Function 1 (VMF1) models [26]- [28]. Previous studies showed the GMF model was more accurate than the NWF model and to be consistent with the VMF1 [29], so the GMF model was adopted in this study. The GAMP software package and the PPP static positioning module [30] were used to estimate the ZTD. The data sampling interval and cut-off elevation angle were set to 30s and 5 • , respectively; and cycle slip detection and repair, satellite and receiver phase center offset and change, and ocean tide loading corrections were all carried out in the pre-processing stage using the software.
The ZTD can be divided into zenith dry delay (or hydrostatic delay) (ZHD) and ZWD. The ZHD contributes about 90% of the ZTD and it closely correlates with the surface pressure, temperature and refractive index of the troposphere [31]. The ZHD can be calculated accurately using a standard model and the Saastamoinen model is the most commonly used model with an about 1 mm accuracy of ZHD [32], thus this model was used in this research. The Saastamoinen model can be expressed as: where P s (unit: hPa) is the surface pressure at the GNSS antenna height; ϕ 0 is the geographical latitude of the antenna; h 0 (unit: km) is the geodetic height of the antenna. The ZWD can be obtained from the ZTD subtracting the ZHD, and then it can be converted into PWV using the following formula: where is the conversion factor that can be expressed as: where ρ w is the density of liquid water (unit: kg/m 3 ); R w is the specific gas constant for water vapor; T m is the weighted mean temperature (unit: K) over the stations site, as mentioned before; K 2 and K 3 are atmospheric refractivity constants and their value are (17 ± 10) K·mbar −1 and (3.776 ± 0.014)×105 K 2 ·mbar −1 , respectively. The numerical approximation of T m can be expressed as: where e i , T i and h i are the mean partial pressure (unit: hPa), mean absolute temperature (unit: K) and the height difference (unit: m) between the ith and (i + 1)th layers.
Using (4) to obtain T m requires water vapor pressure and temperature profiles. Temperature profiles can be obtained from RS data (note that the location of the RS station is usually not the same as any GNSS stations); e i requires a conversion [33], [34]: where T i and rh i are the mean temperature (unit: • C) and relative humidity between the ith and (i + 1)th layers. The relative humidity profile can be also obtained from RS data, like the temperature profile. (5) was given by the World Meteorological Organization in 2008 [35]. The above two formulas are only applicable for RS stations, and its resolution is only twice daily. These two drawbacks mean they are hardly useful for the GNSS-PWV converted from the GNSS-ZTD over the GNSS station. In this case, an empirical T m regression model based on the region of interest is commonly used. The Bevis T m model [13] is based on the RS data of 27 • ∼ 65 • north latitude areas, and its RMS error is 4.74 K. The model may not perform well for the Hong Kong region. Thus, it is better to use a Hong Kong regional T m model. The regional model used in this study can be seen in section IV.B in detail.
The meteorological data including temperature and pressure for the selected three GNSS stations without equipping with meteorological sensors were obtained from the ERA5 dataset. Since the GNSS height is based on the ellipsoid of WGS-84, while ERA5 is geopotential height, a conversion for the GNSS and ERA5 height systems needs to be performed. Vedel [36] mentioned the following two steps needed to perform for the conversion: using the Earth Gravitational Model 2008 [37] to convert the WGS-84 ellipsoid height to a geoid height first, then the geoid height is converted to geopotential height using: where H g and H s are the geopotential height (unit: m) and geoid height (unit: m), respectively; G ϕ and G 0 are the normal gravity (unit: m/s 2 ) of latitude ϕ and 45 • , respectively; r ϕ is the effective earth radius (unit: m) at latitude ϕ. The value of G 0 is 9.80665m/s 2 , and G ϕ and r ϕ can be defined as: In addition, due to the difference between the position of the ERA5 grid and the antenna of a GNSS receiver in both height and horizontal domains, the cubic spline interpolation was used to obtain the pressure and temperature values for the antenna height first; then the inverse distance weighed method was further used to obtain the interpolated values for the horizontal position of the antenna, and the four (horizontal) distances used in this interpolation were from the GNSS station to the four corner points of the ERA5 grid surrounding the position of the station. The whole procedure for the obtaining of the GNSS-PWV value at an epoch over each of the 18 GNSS stations tested is illustrated in Fig.2

III. RESULTS AND DISCUSSIONS A. ACCURACY OF ZTD, TEMPERATURE AND PRESSURE
(2) manifests that the accuracy of both ZWD and determines the accuracy of PWV. The accuracy of the ZWD is determined by that of the GNSS-ZTD and the ZHD resulting from the Saastamoinen model expressed in (1). The conversion factor is obtained from T m which can be calculated from a regional T m model as one of those expressed by (12) - (15), and these models are a function of ground pressure and temperature over the site. Wang et al. [38] analyzed a certain error in each of P S , T m and ZTD contributing to the error in the PWV and stated that a 4 mm ZTD error and a 1.3 K T m error would lead to an 0.64 mm and 0.08 mm error in PWV, respectively; a 1.65 hPa error in ground pressure would result in a 0.6 mm error in the PWV. Since ground pressure and temperature observed by meteorological sensors can be regarded the truth, the accuracy of the ERA5-T and ERA5-P interpolated for the position of the antenna of the   For the accuracy assessment of the GNSS-ZTD, the IGS provided post-processed ZTD (IGS-ZTD) with the time resolution of 5 min and RMS error of 5 mm [39] were used as the reference. Since only the HKSL and HKWS stations were in the list of the IGS continuous observation and tracking stations, the results at these two stations were tested. Fig. 3 shows the 1 min GNSS-ZTD and 5 min IGS-ZTD time series in 2017: It can be found that both the GNSS-ZTD and IGS-ZTD time series agreed well at both stations, which is also shown in Table 2 for their statistical results, with the Bias, Std and RMS under 0.8 mm, 14.4 mm and 14.5 mm, respectively. This accuracy meets the threshold value of 15 mm required for the ZTDs as the input to numerical weather prediction models [40], [41]. In addition, seasonal results were also compared. It was found that the accuracy in the summer was not as good as the other seasons, due to the larger ZTD in the summer. This is likely to be caused by active water vapor in summer. This result implies that, generally, the accuracy of the GNSS-ZTD is likely to be worse in summer [42]. Tables 3 and 4 show the statistics of the ground temperature and pressure values interpolated from the ERA5-and ERA5-P grids (with the sizes 0.25 • and 0.125 • ), respectively for the antenna positions of the aforementioned 15 GNSS stations that were equipped with meteorological sensors. The accuracy of the interpolated result was evaluated by its difference from the meteorological observations collected at the 15 stations in 2017. The statistical results indicate insignificant difference between the interpolated result and the reference, and also between the two grid sizes. The interpolated temperature had small negative bias values at most stations means they were underestimated; their RMS values at most stations were under 2 K, and the maximum Bias (2.26 K) was at the HKKS site, and the mean RMS at all stations was under 1.6 K. The interpolated pressure had small negative values at all stations except one also means they were underestimated; their RMSs at most stations were under 0.6 hPa; the maximum Bias (0.65 hPa) was at the HKTK site, and the mean RMS at all stations was less than 0.5 hPa. Although accurate T m can be obtained from RS data, the very low spatial and temporal resolution of RS data means that if the T m is used for GNSS stations, which are usually not all close to the RS station, the accuracy of its resultant PWV cannot be as desired. Thus, a regional T m regression model developed based on surface meteorological measurement data is commonly used. At present, the main regional T m models include three types: 1) single meteorological factor linear model (SMLM) like the relationship between surface temperature T s and T m expressed in [43]- [45]; 2) single meteorological factor non-linear model (SMNLM), like the relationship of T s and T m stated in [46]; 3) multiple meteorological factors model (MMFM), like the relationship among T s , surface pressure P s and T m expressed in [47], [48]. The SMLM, SMNLM and MMFM can be expressed as (9), (10) and (11), respectively.
where a 0, a, b, c and d are all unknown coefficients/ parameters, which are estimated using the least squares estimation method based on sample data, usually from RS. For the modelling of the MMFM, it is necessary to analyze the relationship between the three factors, i.e. T s , P s and T m . Fig. 4 is the scatter plots of T m to T s and P s from the Hong Kong RS data in the four years period from 2014 to 2017.
Figs. 4(a) and 4(b) show a positive linear relationship between T m and T s and a negative linear relationship between T m and P s . Their correlation coefficients (R) were also calculated and the results were: the R value between T m and T s was 0.86, which was higher, and this is why most regional   T m models were based on T s . Whilst the negative R value of −0.75 between T m and P s implies that the two factors can be used to model T m. . Thus, in this study RS data from the three years (2014 -2016) were selected to estimate the parameters for the aforementioned three types of T m models, and T m from RS data in 2017 were used as the reference to evaluate the accuracy of these models. The time series and statistics of these models' results are shown in Fig. 5 and Table 5, respectively.
We can see that the T m values resulting from all the three models agreed well with the reference and the accuracy of MMFM-T m is the highest. The Bias, Std and RMS of the MMFM-T m were 0.35K, 1.71 K and 1.74 K, respectively. The accuracies of the other two models were both better than that of the currently used Hong Kong regional models, with the RMSs of 1.91 K and 1.89 K [49]. However, some studies have shown that the accuracies of T m in different seasons are different because the sun radiates differently to the ground in different seasons [47], [50], [51]. Thus, to further improve the accuracy of T m models developed based on RS data from one year or more, four new seasonal models were developed based on RS data from each season during the period 2014 -2016. Their accuracies were measured by the differences 87280 VOLUME 8, 2020  Table 6, respectively. From Fig. 6 and Table 6, the accuracies of T m resulting from the three types of seasonal models were all improved over the yearly model previously developed in this study, the mean of the RMSs of the four seasonal SMLM, SMNLM and MMFM results decreased by 2%, 1% and 6%, respectively. In addition, there were noticeable differences in the RMS of the T m values in different seasons, and the T m accuracy in the summer was the highest. Among these three types of models, the MMFM was superior to the other two, with the seasonal mean Bias, Std and RMS of 0.34 K, 1.57 and 1.63 K, respectively. As a result, the four seasonal MMFM models, as the optimal models, were adopted in this study to obtain T m for the selected GNSS stations, which were used in the conversion from the GNSS-ZTD into GNSS-PWV. The four models were: T m = −116.5794 + 1.0259T s + 0.0964P s where (12) (13) (14) and (15) were for spring, summer, autumn and winter, respectively for the Hong Kong region.

C. ACCURACY OF PWV
Since PWV-RS is more accurate than GNSS-PWV, the former is often used as the reference for the validation of the latter. Due to the fact that only one RS station (NO:45004) was deployed in Hong Kong (it launches sounding balloons at UTC 00:00 and UTC 12:00 each day) and HKSC (see Fig. 2) was the GNSS station that was closest to the RS station (about 4 km horizontal distance and 46 m height difference), the accuracy of the GNSS-PWV over the HKSC station was assessed. The PWVs converted from GNSS-ZTD combining with one of the three data sources for pressure and temperature: the one interpolated from European Centre for Medium-range Weather Forecasting ERA5 of 0.25 • or 0.125 • VOLUME 8, 2020   The PWV values in the summer and autumn were larger than that in the winter and spring, and their monthly RMSs were also larger as well, as can be seen in Fig. 8.
The differences in the RMSs resulting from the same month but three different approaches were all under 0.1 mm; the maximum and minimum monthly RMSs were in July (about 3.8 mm) and January (1.8 mm), respectively. It should be noted that the mean RMS values of the temperature and pressure interpolated from the ERA5 datasets were 1.59 K ( Table 3) and 0.48 hPa (Table 4), respectively, which may introduce a less than 0.5 mm error in PWV. According to the relationship between an error in ZTD and its resultant PWV error [38], the mean 14.3 mm error in the GNSS-ZTD over the two GNSS stations shown in Table 2 may result in an about 2 mm error in PWV. These results indicate that the error in GNSS-PWV is mainly contributed by the error in the GNSS-ZTD. The large PWVs in the summer and autumn were due to the subtropical monsoon climate of Hong Kong,  Fig. 7, and RS-PWV in the month was the reference. and the atmospheric water vapor over the region is generally very active in the seasons. This is why a variety of water vapor related extreme weather events like typhoons, storms, heavy rainfall often occur in these two seasons. Thus, it will result in an additional error for obtaining ZTD when we use a general troposphere delays model to estimate the ZTD, because existing troposphere delays models such as the Saastamoinen and Hopfield models, are not capable to consider complex refractivity fields [52].

D. TEMPERAL VARIATION IN PWV DURING TROPICAL CYCLONES
The coastal area of Hong Kong suffers from invasions of TCs that form over the Pacific and Atlantic Oceans every year. There were seven TCs influencing the Hong Kong in 2017, which can be classified into the following six categories, according to the maximum sustained surface winds near the center of the TCs (Table 7): one super typhoon named 'Hato'; one severe typhoon named ''KHANUN''; three severe tropical storms named 'merbok', 'PAKHAR', 'MAWAR', respectively; one tropical storm named 'ROKE' and one tropical depression. The last one (tropical depression) did not cause any significant damage thus it has not been named by the Hong Kong Observatory (HKO). Therefore, in this study the other six TC events were investigated for the characteristics of water vapor variations in the duration of the TCs. Fig. 9 shows the movement path of these six TCs.
A long period of heavy rain was observed after each TC made landfall. The statistics of rainfall recorded by several rain-gauge stations in Hong Kong is given in Table 8. Long periods of rainfall require an abundant supply of water vapor, 87282 VOLUME 8, 2020  and the interior of TCs contains large amounts of water vapor that provides conditions for sustained rainfall. During the six TCs' periods, the PWV over the Hong Kong region would first increase dramatically, then it began to decrease significantly as the TCs moved away from Hong Kong and carried away a large amount of water vapor. Fig. 10 shows the time series of RS-PWV and GNSS-PWV over the HKSC station during the periods of the six TC events.
According to Fig. 10, both GNSS-PWV and RS-PWV time series most agreed well. The deviations between both may be partly caused by different observation paths: the sounding balloons usually travel for a while and would deviate from the vertical direction over the RS station, especially in the conditions of strong winds from a TC, whist a GNSS observation was from a moment. It is noted that the two rectangular boxes labelled 'signal' and 'closest' in duration of each TC means the time the strong wind signal was issued by the HKO and the time the highest gale or storm signal was issued by the HKO, respectively. The former also means the TC began to approach Hong Kong and the latter indicates the center of the TC moved closest to Hong Kong and would soon start to depart from Hong Kong as well. It can be seen that when a TC was approaching Hong Kong the PWV values kept increasing for several or a dozen hours; then, when the TC started to depart from Hong Kong the PWV values started to decrease. This temporal variation trend in PWV, i.e. an increase as a TC approaches and a quick decrease after the TC passage, was consistent with that of many other TC cases [21]- [23]. The low temporal resolution of RS-PWV (twice daily) is not  TABLE 9. Correlation coefficients between wind speed and PWV time series before and after the TCs made landfall sufficient for the study of the characteristics of the PWV variation, especially for the weather conditions of TCs, not to mention the possible drift of sounding balloons caused by strong winds during the TC occurrence. Wind facilitates water vapor transport during a TC period, as shown by the wind speed and GNSS-PWV time series during the six TC events in Fig. 11: It can be seen that the trends of the temporal variation in both PWV and wind speed were similar, and both reached their maximum when the TCs were closest to Hong Kong. This phenomenon further indicates that the TC wind field contains a large amount of water vapor, which may be transported through the passage. To understand the relationship between them, Table 9 shows their correlation coefficients before and after the six TCs' landfall.
The variation in atmospheric water vapor is mainly dominated by the thermodynamic and dynamic processes. The thermodynamic process is due to the absorption or release of energy by atmospheric water vapor, thus the water morphology varies such as evaporation and liquefaction of water. The dynamic processes of water vapor include regional or global atmospheric circulation. Table 9 shows that the correlation between PWV and wind speed was close to 0 before the arrival of a TC, while the correlation coefficient was significantly large after landfall. A strong correlation was shown in the case of PAHAR, with a 0.79 value. The significant correlation manifests that the main factor influencing water vapor variation in Hong Kong after the arrival of a TC is the dynamic process of water vapor transport related to wind. Fig. 11 and Table 9 have shown a good agreement between the temporal variation in PWV and wind speed during the VOLUME 8, 2020  TCs' periods, i.e. their correlation was positive, while the correlation between PWV and the distance of the center of the TC was negative. This result was consistent with the ones found in several cases [22]- [24]. In addition, the grade of a TC in a region can be determined by wind speed, as shown in Table 7. Fig. 12 shows the buffer zone of a typical TC, where the higher the level of the wind circle, the larger the wind speed; and the smaller the radius of the wind circle, the nearer to the center of the wind circle.

E. SPATIAL VARIATION IN PWV DURING TROPICAL CYCLONES
From section D, the temporal variation in PWV reflected the life cycle of a TC, e.g. the PWV values over the station increased when a TC approached, and the PWV values decreased when the TC departed. However, for the monitoring of a TC's status, especially when a TC approaches a region, the spatial variation of PWV in the region needs to be monitored and investigated. In this study, hourly PWV time series over all the aforementioned 18 GNSS stations were used as the sample data and the kriging interpolation method was to obtain the interpolated PWVs for the whole region. Resultant PWV maps are depicted by the twodimensional shown in Fig. 12. It is noted that the figure only shows four TCs as examples, and for the case when the four TCs were approaching the region from different directions. From Fig. 13, the spatial distribution of PWV was dependent upon the location of the TC, i.e. the area with the maximum PWV was close to the TC. The spatial distribution of PWV might indicate the movement direction of the TC, e.g. the larger PWVs over the east-south area in Fig. 13(e), east-north area in Fig. 13(f), east-south area in Fig. 13(g) and east Fig. 13(h) area were similar to the location of the TCs 'merbok' in Fig. 13(a), 'ROKE' in Fig. 13(b), 'Hato' in Fig. 13(c) and 'KHANUN' in Fig. 13(d), respectively. In addition, there were significant local minimum PWV over the HKST and HKNP stations in Fig. 13(e -h). The two main reasons for this may be (1) the two stations were both located in the mountainous area, meaning their elevations were much higher than the other stations, as shown in Table 1. When a typhoon passed the mountainous area near the station, the typhoon would be weakened due to the blocking effect of the terrain, which would cause the movement of water vapor from higher ground to nearby lower ground [53]; (2) A TC usually causes precipitation events, and the precipitation tends to increase with the increase in elevation due to the orographic effect of mountainous terrain, which causes the air to be lifted vertically, and the condensation occurs due to adiabatic cooling [54].

F. USING GNSS-PWV SIGNALS TO MONITOR TROPICAL CYCLONES
Section D has indicated that the increase and decrease in the temporal PWV over a site were related to the activities of a TC; Section E has showed that the spatial variation of PWV over the region near the center of a TC was larger than the region away from the center when the TC approached. These results suggest that the temporal and spatial variations in PWV may be used to analyze some characteristic of a TC. Since a large PWV value above each GNSS station was mainly brought by a TC, as a signal at the station, the time of PWV arrival (TOPA) may be used to indicate the movement status of the TC. Based on the times of the signals occurred over several stations, the movement of a typhoon can be estimated, which can be found from very recent researches [21], [22]. However, these study cases were only from super typhoons. Generally speaking, the radii of high-level TCs were larger than that of low-level TCs, as shown in Fig. 14 for different wind circle radii of four presupposed TCs in the Hong Kong region when the TCs approach the region.
In Fig. 14, the wind radii of the TC1, 2, 3 and 4 decrease successively. If the radius of a TC (e.g. TC1) is much larger than the distance between two adjacent GNSS stations, the edge of the TC could be considered a line. In the literatures [21], [22], the TCs with smaller radii were not analyzed. Fig. 15 shows a corrected movement model of a TC for which a curve edge was taken into account for the case that the wind radii of a TC is not large.
In Fig. 15, the red area and small circles represent the wind circle of the TC and GNSS stations, respectively. When the edge of the TC arrives at a GNSS station, the PWV value over the station starts to increase (see the 'signal' in Fig. 10). The TOPA over all these stations covered by the TC can be obtained accurately from the high temporal resolution GNSS-PWV time series. Since the TC moves with time, the TOPA values at two GNSS stations must be different. From Fig. 15, an independent plane rectangular coordinate system was established as below: the NO.1 station is assumed the origin over which the PWV brought by the TC first reaches; the X-axis and Y-axis are in the north and east. When the TC with a direction angle of θ starts move from the NO.1 station to the NO.2 station, the TC moves at a distance of d 12 . If the radius of this TC is much larger than this distance, the edge of the TC will be regarded as a straight line, hence d 12 can be easily obtained from a geometric relationship between the two stations. However, in many cases, this condition can be hardly met. It is noted that d 12 would be approximately the distance of the line l 12 between the two stations when the azimuth of l 12 is approximately θ or θ ± 180 • . For the ith station on the edge of the TC and any jth station (i < j), the azimuth of l ij can be calculated using: where θ ij is the azimuth of l ij ; (x i , y i ) and (x j , y j ) are the coordinates of the ith and jth stations, respectively. If the absolute difference of θ ij and θ is close to 0 • , the distance d ij of the jth station and the edge of the TC along the direction of θ can be calculated from: The TOPA difference t ij between the two stations indicates the time difference of the TC arriving the two stations. Then, the movement speed v ij of the TC during the period from the time the TC arrives the ith station to the time the TC arrives the jth station can be simply calculated:   In fact, (18) is only suitable for the mean movement speed of a TC and its prerequisite is that the TC has passed the two stations. Moreover, the movement speed of a TC usually varies, thus an acceleration parameter a can be introduced to correct (18) as below: where v is the initial movement speed of the TC at the NO.1 station; t i1 is the TOPA difference between the 1st and ith stations. In (19), θ, a and v 1 are all unknown parameters to be estimated. Let t ij = t k1 , t i1 = t k2 , θ ij = θ k , l ij = l k and the kth observation equation is: where n is the number of observations used in (19), and its residual r k is: (20) can be linearized using the first-order Taylor expansion: where θ m 0 , a m 0 and v m 0 are the values of the mth iteration of θa and v, respectively. The matrix form of (21) can be expressed as: where A is the Jacobian matrix of L(θ, a, v). For the mth iteration, A m and β m can be expressed as (24), shown at the bottom of the next page.
In (24), θ, a and v are the unknown parameters, and the equation can be solved using the least squares method (under the condition that redundant observations are available):  Let q m be the estimate result from the mth iteration, then update the approximate values for the unknown parametersfor the (m + 1)th iteration by: Substitute the new approximate values into (22) and repeat the above recursive process until the result converges, i.e. the estimated q values are under pre-defined thresholds or the residual r k+1 r k values are neglectable.
The above approach was performed using the TOPAs over the aforementioned 15 GNSS stations when the TCs approached Hong Kong. The movement speeds of the TCs from a post-processed result provided by the China Meteorological Administration (CMA) were used as the reference for the validation of the movement speed estimated by the models used in this study. The estimated movement speeds and the reference were listed in Table 10, which shows a good agreement, with the difference values of all under 6 km. This result suggests that the model adopted in this study can be used to estimate the movement speed of a TC.

IV. CONCLUSION
Water vapor retrieval with high spatiotemporal resolution is crucial for the study of TCs' activities, and high temporal resolution GNSS-PWV time series over the Hong Kong CORS network were used to study TCs weather events during the year of 2017. The seasonal multi-factor T m model developed based on seasonal sample data in this study, as the optimal model, outperformed the existing regional T m models in Hong Kong, in terms of accuracy. The seasonal T m model was used to obtain the T m values for the GNSS stations to convert GNSS-ZTD into GNSS-PWV. Moreover, based on the differences in the TOPAs over different GNSS stations when a TC approached, a model for the estimation of the TC's movement was proposed and its test results agreed well with the reference provided by the CMA. These results suggest that it is promising to apply GNSS derived water vapor products to the monitoring, or even forecasting, of short-term TC events. However, for further improvement to the proposed approach, some other factors may need to be taken into account, e.g. when a TC makes landfall, its radii and movement direction may change, due to the blocking effects of the terrain, especially in mountainous areas. For a long-term motion state of TCs, the variation in the radii and the moving direction of TCs should be considered in the model. This may be one of our future works. In addition, although this study was based on GNSS post-processing results, its methodologies are most similar to operational real-time tracking system, and the main differences include that the former read data from files while the latter reads data from stream. Another difference is in the accuracies of satellite orbit and clock product. For realtime operational cases, ultra-rapid products of satellites are available, the accuracy of which is lower than that of the IGS products used for post-processing like the ones used in this study. This may affect the performance of real-time applications of the approaches proposed in this research. This will be tested in future.