Spatiotemporal characteristics of GNSS-derived precipitable water vapor during heavy rainfall events in Guilin, China

Precipitable Water Vapor (PWV), as an important indicator of atmospheric water vapor, can be derived from Global Navigation Satellite System (GNSS) observations with the advantages of high precision and all-weather capacity. GNSS-derived PWV with a high spatiotemporal resolution has become an important source of observations in meteorology, particularly for severe weather conditions, for water vapor is not well sampled in the current meteorological observing systems. In this study, an empirical atmospheric weighted mean temperature (Tm) model for Guilin is established using the radiosonde data from 2012 to 2017. Then, the observations at 11 GNSS stations in Guilin are used to investigate the spatiotemporal features of GNSS-derived PWV under the heavy rainfalls from June to July 2017. The results show that the new Tm model in Guilin has better performance with the mean bias and Root Mean Square (RMS) of − 0.51 and 2.12 K, respectively, compared with other widely used models. Moreover, the GNSS PWV estimates are validated with the data at Guilin radiosonde station. Good agreements are found between GNSS-derived PWV and radiosonde-derived PWV with the mean bias and RMS of − 0.9 and 3.53 mm, respectively. Finally, an investigation on the spatiotemporal characteristics of GNSS PWV during heavy rainfalls in Guilin is performed. It is shown that variations of PWV retrieved from GNSS have a direct relationship with the in situ rainfall measurements, and the PWV increases sharply before the arrival of a heavy rainfall and decreases to a stable state after the cease of the rainfall. It also reveals the moisture variation in several regions of Guilin during a heavy rainfall, which is significant for the monitoring of rainfalls and weather forecast.


Introduction
Atmospheric water vapor, as a key factor in regulating Earth's climate, plays an important role in global atmospheric radiation, water cycling, and energy balance, and its variation is the main driving force of climate change (Wang et al., 2007). Precipitable Water Vapor (PWV) refers to the total water vapor content of a vertical column per unit area in atmosphere (King et al., 1992). Understanding the spatiotemporal variation of water vapor is of vital significance for the evolution of weather and climate prediction. Various techniques have been developed to acquire the PWV content over the past few decades. Traditional methods mainly include radiosondes, microwave radiometers, and satellite remote sensing. Among them, radiosonde is the primary in situ measurement to obtain PWV because of its high accuracy. However, it is hard to meet the demands of monitoring and forecasting extreme weather at the mesoscale and microscale due to its high cost and low

Open Access
Satellite Navigation https://satellite-navigation.springeropen.com/ *Correspondence: zxmo@glut.edu.cn; xieshaofeng@glut.edu.cn 1 College of Geomatics and Geoinformation, Guilin University of Technology, Guilin 541004, China Full list of author information is available at the end of the article spatiotemporal resolution. With the development of the Global Positioning System (GPS), GLObal NAvigation Satellite System (GLONASS), Galileo navigation satellite system (Galileo), BeiDou Navigation Satellite System (BDS), and other satellite navigation systems, a reliable technique to retrieve PWV with the Global Navigation Satellite System (GNSS), called ground-based GNSS, has been proposed and is fully operational (Bevis et al., 1992;Hein, 2020;Li et al., 2015a;Yang et al., 2020). Groundbased GNSS can detect satellite signal propagation delays in atmosphere and further retrieve the PWV, which has attracted a widespread attention due to its advantages of high temporal resolution, high precision, low-cost, and resistance to all weather conditions (Elgered et al., 1997;Emardson et al., 1998;Vaquero-Martínez et al., 2017).
In the ground-based GNSS technique, the retrieval of PWV from GNSS observations entails the acquisition of the Zenith Wet Delay (ZWD) and the conversion coefficient (Π). GNSS signals are affected by troposphere refraction when propagating through the neutral atmosphere, thus causing a tropospheric delay, which is expressed as the Zenith Total Delay (ZTD). The ZTD consists of two components: the Zenith Hydrostatic Delay (ZHD), which is caused by the dry gases of the troposphere, and the ZWD, which stems from the water vapor. The ZTD can be accurately calculated from GNSS observations using high-precision GNSS processing software. The ZHD can be estimated with tropospheric empirical models. Then, the ZWD can be obtained by subtracting ZHD from the ZTD. To convert the ZWD to the PWV, Bevis et al. (1992) proposed a formula to calculate the atmospheric weighted mean temperature (T m ). T m , as one parameter to compute Π, can be exactly calculated from the vertical profiles of atmospheric temperature and humidity based on numerical integration or estimated from an empirical T m model with surface temperature (T s ) (Bevis et al., 1994;Davis et al., 1985). Many studies have mentioned that the relative error of Π basically comes from T m , and the accuracy of T m is regarded as one of the largest error sources in PWV calculation (Huang et al., 2019a;Liu et al., 2012;Wang et al., 2005). Therefore, precisely estimating T m is of great significance to improve the accuracy of the GNSS-based PWV retrieval. Generally, it is not easy to access the T m in some regions because of the lack of atmospheric profiles. Thus, an accurate empirical T m model is needed to satisfy such demands and provide convenience for users. Bevis et al. (1994) found that T m and T s have a strong linear correlation and proposed an empirical formula (T m = 70.2 + 0.72T s ) using the 8718 radiosonde profiles in North America that is commonly used to estimate T m . Some regional empirical T m models with higher accuracy than the Bevis formula have been formulated. Li et al. (1999) constructed a linear regression equation of T m and T s in eastern China (T m = 44.05 + 0.81T s ) by using the radiosonde data of eastern China, which has a good performance in that region. When obtaining an accurate T m parameter, the error of GNSS-derived PWV can be reduced. In some previous studies, the relative accuracy of the GNSS-derived PWV has been proven to be reliable with the root-mean-square values of 1-3 mm (Gui et al., 2017;Zhao et al., 2019a). Therefore, the GNSS-based PWV has been extensively used in meteorology applications, especially in monitoring and forecasting extreme weather conditions (Calori et al., 2016;Li et al., 2015b;Liu et al., 2019b;Wang et al., 2013;Zhao et al., 2018).
In recent years, a great number of studies have been conducted in the applications of the GNSS retrieving water vapor in the analysis of extreme rainfall conditions. Zhang et al. (2015) investigated the signature of GPSderived PWV with two severe weather case studies. Shi et al. (2015) analyzed a series of rainfall events in Wuhan, China using real-time GPS Precise Point Positioning (PPP) based PWV estimation. Chen et al. (2017) used the Hong Kong GPS network to study the water vapor variability during three heavy rainfall events in Hong Kong. Yao et al. (2017) showed that GNSS-derived PWV has a correlation with rainfall. They also investigated three different rainfall events in Wuhan and Zhejiang regions. Manandhar et al. (2018) studied the trend of PWV values derived at a GPS station in Singapore with diurnal and seasonal variations by separating several rainy and clear days. Barindelli et al. (2018) analyzed the temporal variations in the PWV derived from a regional GNSS network during two heavy rainfall events. Moreover, they found the characteristics corresponding to the intense precipitation events. Chen et al. (2018) constructed a PWV map with the data at Hunan GNSS stations and synoptic observations to reveal the water vapor advection, transportation, and convergence during a heavy rainfall. These studies indicate the feasibility and practicality of GNSSderived PWV to monitor rainfall events. Although the analysis of PWV characteristics during the rainfall is performed, most studies are limited in providing the detailed spatiotemporal evolution with GNSS-derived PWV during heavy rainfall events.
Guilin, which has numerous mountains and rivers, is a well-developed tourist city and has been awarded the only international tourism city in China located in a rain-prone area in low latitudes with an average annual precipitation of more than 1500 mm. In June and July 2017, several heavy rainfalls occurred in Guilin, which induced continuous floods, landslides, and other disasters. Moreover, only one radiosonde station is in Guilin. Thus, the weather changes of the whole area are difficult to monitor. Using the GNSS technique for capturing the signature of PWV with high spatiotemporal resolution is of great significance. The extreme weather in Guilin should be evaluated and forecast by studying the spatial and temporal distribution and evolution of water vapor. In this work, we establish an empirical T m model using radiosonde records for GNSS PWV calculation. Then, the empirical T m model is compared with the other three commonly used T m models to verify its reliability and availability. Later, the hourly GNSS PWV is obtained using the observations at 11 GNSS stations. We evaluate its accuracy by comparing it with radiosonde-derived PWV. Finally, the relationship between GNSS PWV and rainfall is analyzed, and the spatiotemporal characteristics of GNSS PWV under heavy rainfalls in Guilin are investigated.

Study area
Guilin is situated in southern China with a territory of approximately 27,800 km 2 . It has a subtropical humid monsoon climate with four distinctive seasons. Moreover, Guilin has plenty of rainfall. The average annual rainfall varies from 1900 to 2000 mm, which is mainly concentrated from April to July, and the longest continuous precipitation period is 30 days. Heavy showers frequently occur in summer, thus leading to serious threats to the safety of people and property. Using GNSS to monitor the highly variable characteristics of water vapor has a great potential in improving the ability of extreme weather forecasting for Guilin. In this study, the data from the Guilin radiosonde station, 12 ground meteorological stations, and 11 GNSS stations and the ERA5 reanalysis dataset in 2017 are used to analyze the spatiotemporal characteristics of GNSS-derived water vapor during heavy rainfall events in Guilin. The corresponding stations are presented in Fig. 1. The details of GNSS and radiosonde stations, including the station name, geodetic coordinates, elevations, and distance between each GNSS site and the radiosonde station are shown in Table 1.

Radiosonde profiles
The radiosonde is one of the most commonly instruments to measure high-quality meteorological parameters. In this study, radiosonde observations are used to establish the empirical T m model for Guilin and serve as the reference values to evaluate the GNSS PWV results. Radiosonde profiles contain the atmospheric parameters in vertical direction collected with a sounding balloon at 00:00 and 12:00 Coordinated Universal Time (UTC) every day, including surface parameters, such as surface temperature T s and pressure (P s ), and pressure level parameters, such as absolute temperature (T), relative humidity (RH) and geopotential height (H) at every pressure level. From the meteorological parameters at different levels, T m can be calculated using the numerical integration of radiosonde measurements, which can be expressed as follows (Bolton, 1980;Wang et al., 2016): where Δh i represents the height difference of ith layer (unit: m), n represents the number of layers, e i and T i are the average water vapor pressure (unit: hPa) and temperature (unit: K) at the ith layer of the atmosphere, respectively, e s denotes the saturated vapor pressure (unit: hPa) and T d the atmospheric temperature in Celsius (T = T d + 273.15).

ERA5 reanalysis dataset
ERA5 is a newly released fifth-generation European Centre for Medium-Range Weather Forecasts (ECMWF) global climate reanalysis dataset. The ERA5 dataset is produced by using the 4D-Var data assimilation scheme in the CY41R2 model of the ECMWF's Integrated Forecast System (IFS), which is expected to gradually replace the prior ERA-Interim dataset. ERA5 provides hourly RH · e s 100 (3) e s = 6.112 × 10 7.5×T d T d +237.3

Fig. 1
Geographic distribution of the 11 GNSS sites (red triangles), 1 radiosonde station (yellow pentagram), and 12 ground meteorological (MET) stations (blue circles) in Guilin, China atmospheric reanalysis data at 37 pressure levels from 1000 to 1 hPa, with a horizontal spatial resolution of 0.25° × 0.25°. Surface parameter products are also available (such as T s , P s and PWV). In addition, ERA5 has been demonstrated to have a good performance in PWV estimates (Wang et al., 2020;Zhang et al., 2019b). In this study, the gridded T m series derived from ERA5 profiles (37 levels) are used for the evaluation of T m models over Guilin, and the PWV series derived from ERA5 surface products are used to compare with the spatial evolution of GNSS PWV.

Surface meteorological observations
The hourly rainfall data, surface pressure, and surface temperature observations are provided by the China Meteorological Administration (CMA) at 12 MET stations from June to July 2017 over Guilin. The meteorological sensors at these stations are well maintained and regularly adjusted every year by the CMA. As the GNSS and meteorological stations are not collocated ( Fig. 1) and most of the GNSS stations are not equipped with meteorological sensors, the T s and P s at a GNSS site are calculated with their corresponding values at the nearest MET station for GNSS PWV retrieval using the equations below (Suparta & Rahman, 2016): where T MSL and P MSL indicate the mean sea level temperature (unit: °C) and mean sea level pressure (unit: hPa), respectively, T MET and P MET are the temperature and surface pressure at MET station, respectively, T GNSS and P GNSS are the temperature and surface pressure at GNSS station, respectively, and h MET and h GNSS represent the height of the MET and GNSS stations, respectively (unit: km).

GNSS observations
The GNSS observations are collected from the Guilin GNSS network from June to July 2017 ( Fig. 1). Currently, the International GNSS Service (IGS) center provides the ZTD product with a precision of better than 5 mm (Byun & Bar-Sever, 2009;Huang et al., 2019b), which can be used as the reference values for assessing other ZTD results. In this study, the correlation between the tropospheric parameters across the Guilin GNSS network is reduced by introducing six IGS stations (BJFS, SHAO, CHAN, LHAZ, TWTF, and URUM), and the GNSS data are processed with the GAMIT/GLOBK software (a high-precision GNSS data processing software). The ZTD is estimated with a half-hour interval. To confirm the reliability of ZTDs, the ZTDs obtained from the IGS center (IGS-ZTD) are used to evaluate the ZTDs derived from the GNSS observations (GNSS-ZTD) at stations BJFS and TWTF from June to July 2017, which is from Day Of Year (DOY) 152-213. The ZTD differences with their deviations from the mean value larger than three times the STandard Deviation (STD) are removed as gross errors. The results are shown in Fig. 2 From Fig. 2, one can see that the IGS-ZTD and the GNSS-ZTD show good consistency at stations BJFS and TWTF. Compared with the IGS-ZTD, the correlation coefficient (R) and the Root Mean Square (RMS) error of the GNSS-ZTD are 0.998 and 5.3 mm at station BJFS and 0.990 and 5.6 mm at station TWTF, respectively. The results indicate that the ZTD derived from the GNSS observations are reliable and accurate for GNSS PWV retrieval.

Retrieval of PWV from GNSS observations
The basic formula for retrieving PWV from GNSS observations is as follows (Askne & Nordius, 1987): In this study, the ZTD is derived from the Guilin GNSS observations processed with the GAMIT/ GLOBK software. After retrieving the ZTD, the ZWD is calculated by subtracting the ZHD from the ZTD, where the ZHD can be estimated with the Saastamoinen model (Saastamoinen, 1972): where φ is the latitude of station (unit: radians), h o is the height of the station above sea level (unit: km), and P s is the surface pressure of the station (unit: hPa), interpolated from the nearby MET stations. With T m , Π can be calculated as follows: where ρ w is the density of liquid water (1 × 10 3 kg/ m 3 ), R v is the specific gas constant for water vapor (461.495 J kg −1 k −1 ), and k' 2 (22.13 ± 2.2 K/hPa) and k 3 ((3.739 ± 0.012) × 10 5 K/hPa) are the atmospheric physical constants. Given that the humidity and temperature profiles are usually hard to obtain, a regression analysis method, based on the linear correlation between T m and T s , is often adopted to determine the T m :  where a and b indicate the coefficients of regression equation. These coefficients are estimated using the least squares adjustment with the years of radiosonde data, and the formulas are as follows: where In addition, a widely used empirical tropospheric delay model, the Global Pressure and Temperature 2 wet (GPT2w) model, can also provide the T m product (Böhm et al., 2015). The GPT2w model has two horizontal resolutions of 5° × 5° (GPT2w-5) and 1° × 1° (GPT2w-1), and has excellent performance in T m estimation against other models (Huang et al., 2019b).
The summary of GNSS PWV processing is shown in Fig. 3. First, we use the radiosonde profiles to establish a regional empirical T m model. Then, the hourly P s and T s from ground meteorological stations are used to interpolate those for the near GNSS site. The interpolated P s is used to calculate the ZHD by Eq. (9). Then, the interpolated T s is used to obtain the T m with the regional empirical T m model. Subsequently, the ZWD is calculated by subtracting the ZHD from the ZTD, which is obtained from GNSS observations. The conversion coefficient Π is calculated by Eq. (10). Finally, the hourly GNSS PWV is obtained by multiplying the ZWD by Π.

Establishment of T m model for Guilin
Due to the systematic errors of the global T m models, which cannot provide superior performance, a regional optimized T m model needs developing (Huang et al., 2019a). Thus, to reduce the T m calculation error and improve the reliability of GNSS PWV, radiosonde profiles are used to establish the local empirical T m model for Guilin. Before establishing the T m model, the time series and relationship between T m and T s are presented in Fig. 4 using Guilin radiosonde profiles from 2012 to 2017. As seen in Fig. 4a, the T m and T s have a similar trend over time. The T s varies between 270 and 310 K, while the range of T m is smaller than T s , ranging from 260 to 295 K. Figure 4b shows that the R between T m and T s is 0.91, which indicates the T m and T s have a strong linear correlation in Guilin. Thus, the regression equation coefficients of a and b are calculated by Eqs. (12) and (13) with the radiosonde profiles from 2012 to 2017. Then, the regional empirical T m model for Guilin (Guilin model) is established as follows:

Evaluation of T m model
To evaluate the performance of the Guilin T m model, the gridded T m series derived from the ERA5 dataset in 2017, which are different from that used in building Guilin model, are used as references. The widely used Bevis model (Bevis et al., 1994), GPT2w-1 model and the Li model (Li et al., 1999) Table 2. According to Table 2, the Li model has the largest mean bias and RMS, whereas the Guilin model has the smallest RMS and STD values of 2.12 and 2.04 K, respectively, and the mean bias of − 0.51. The Guilin model improves the accuracy (in RMS) of T m estimation by 35%, 40% and 22% with respect to the Bevis, the Li, and GPT2w-1 models, respectively. Thus, the accuracy of the Guilin model is evidently higher than other three models over Guilin.
According to several studies (He et al., 2017;Huang et al., 2019aHuang et al., , 2019bWang et al., 2005Wang et al., , 2016, the impact of the T m values calculated by the four models on their resultant GNSS PWV is analyzed in terms of theoretical function. In this study, a similar method is utilized to analyze the impact of T m values on GNSS PWV. The relationship of the RMS between T m and PWV can be calculated by the following formula: where RMS pwv is the RMS of PWV, RMS T m is the RMS of T m , PWV and T m are given by the radiosonde profiles and the T m with models, respectively. RMS pwv /PWV is defined as the relative error of PWV. RMS pwv and RMS pwv /PWV are used to assess the impact of the errors in T m on its resultant GNSS PWV results. The theoretical results of RMS pwv and RMS pwv /PWV in different T m models are shown in Table 3. Table 3 shows that the mean RMS pwv and RMS pwv / PWV values of the Guilin model are the smallest among the four models. In terms of RMS pwv , the RMS pwv values of the Guilin model are all less than 0.60 mm, and the mean RMS pwv is 0.29 mm. In terms of RMS pwv /PWV, the mean RMS pwv /PWV of the Guilin model is 0.74%,   ranging from 0.71 to 0.77%, which is more stable than other models. From the results, the impact of the T m values calculated from the Guilin model on the GNSS PWV is smaller than those in the Bevis, Li and GPT2w-1 models. The Bevis model and the Li model were established by using the various radiosonde stations in mainland United States (27-65° N area) and eastern China (20-50° N, 100-130° E area), respectively. The empirical model GPT2w-1 was developed based on global ERA-Interim dataset but no input of meteorological parameters. As a result, the systematic errors will be produced when the three models are used in Guilin. Moreover, GPT2w-1 provides one T m value a day, which cannot meet the demand for obtaining hourly T m values. Thus, the T m values derived from other three models are not optimal. Compared with the other three models, the model established in local areas has smaller errors and greater applicability. Hence, the Guilin model is the preferred model, and can be applied to the T m estimation for hourly GNSS PWV retrieval in Guilin.

Comparison of GNSS-derived PWV with radiosonde-derived PWV
After retrieving the PWV, the accuracy of GNSS-derived PWV (GNSS-PWV) needs to be evaluated. In this study, the PWV retrieved from the data at the Guilin radiosonde station at 00:00 and 12:00 UTC per day from June to July 2017 is compared with the GNSS-PWV. The GNSS station JZ87 was selected because it is nearest to the Guilin radiosonde station. The GNSS-PWV values at the corresponding time were extracted for comparison.
The RadioSonde-derived PWV (RS-PWV) is based on the geopotential height system, while the GNSS-PWV adopts the geodetic height system. Thus, it is necessary to unify the different height systems when the GNSS-PWV are compared with the RS-PWV. In this study, we follow the method described in Wang et al. (2016) to unify the different height systems. PWV is greatly affected in the vertical direction, and the PWVs at different heights must be unified to reduce the impact of height differences on the PWV comparisons. The empirical correction model of PWV is used as follows Zhao et al., 2019b): where PWV h 1 and PWV h 2 denote the PWVs corresponding to the heights of h 1 and h 2 , respectively. After the GNSS-PWV and RS-PWV are referred to the same height, they can be compared. The PWV comparison form June to July 2017 is shown in Fig. 5.
As presented in Fig. 5, the GNSS-PWV agrees well with the RS-PWV except some small differences. Moreover, the RS-PWV is generally greater than the GNSS-PWV. The R of the two kinds of PWVs is 0.79, which shows a strong correlation. The mean bias of the PWV differences is − 0.90 mm, and the RMS and STD are 3.53 and 3.43 mm, respectively. Notably, the accuracy of GNSS-PWV meets the requirements of atmospheric research (Wang et al., 2013). The GNSS-PWV and the RS-PWV have a good agreement, and the PWV derived from the GNSS observations is reliable, which can be used to monitor water vapor variation and analyze the PWV spatiotemporal characteristics.

PWV time characteristics
Rainfall intensity is the average amount of rainfall over a certain period. It is also an important indicator to describe the characteristics of a heavy rainfall. The increase in the intensity entails heavier rainfall. In China's meteorological department, the classification of rainfall intensity has generally adopted a standard, as shown in   Table 4. To improve the description of rainfall, we use this standard in the following discussion. From June to July 2017, several continuous rainfalls occurred in many areas of Guilin. The daily precipitation observed by 12 MET stations in Guilin is shown in Fig. 6.
As shown in Fig. 6, several days of rain was recorded from June to July 2017 in Guilin. Among them, especially around 1 July (DOY 182), the precipitation observed by the 12 MET stations in Guilin was almost more than 50 mm, reaching the level of heavy rain. One of these rainfalls had more than 250 mm of precipitation in Yongfu County, reaching the highest level in rainfall intensity. On the basis of the weather condition and the available data, the observations at seven GNSS stations in Guilin from 27 June to 2 July, 2017 and the related rainfall data from MET stations near the GNSS sites are selected as examples to analyze the PWV time-varying characteristics and the relationship between GNSS-PWV and the actual precipitation during the rainfall events, as shown in Fig. 7. Figure 7 shows that the PWVs retrieved from seven GNSS stations (JZ84, JZ85, JZ87, JZ89, JZ92, JZ94, and JZ95) in six days are evidently correlated with the actual precipitation, and have a similar trend. The PWV greatly changes, especially in the case of heavy rain, and most of the precipitation occurs near the peak values of PWV. The PWV usually increases significantly before a heavy rainfall and changes slightly but remains stable during the heavy rainfall. After the heavy rainfall, the PWV drops sharply and reaches the trough. The rainfall that occurred on 1 July 2017 (DOY 182) at station JZ84 is used as an example to analyze the PWV variations in a rainfall process (see Fig. 7a). Before the heavy rainfall, the PWV values fluctuated around 55-60 mm before 28 June (DOY 179). On 28 June, the PWV fell to 50 mm. Then, it began to rise until 1 July (DOY 182) to the peak value of 65.6 mm. When the heavy rainfall stopped, the PWV decreased rapidly and down to the minimum of 44.2 mm on 2 July (DOY 183). Similar variations of PWV and precipitation can also be found in other stations during the heavy rainfall.
To further investigate the relationship between time varying PWV and rainfall, seven rainfall events at different GNSS stations are selected for analysis by grouping them into the before (red circles, see Fig. 7) and after (black circles, see Fig. 7) the rainfall events. The concept of PWV variation and the rate of time varying PWV are defined as follows (Yao et al., 2017): ΔPWV (PWV variation) = Max PWV (the maximum PWV before a rainfall) − Min PWV (the minimum PWV before/after the adjacent maximum PWV); Interval time refers to the period between the Max PWV and the Min PWV; the rate of change of PWV = ΔPWV/Interval time. The results are represented in Tables 5 and 6, respectively.
In Table 5, before rainfall the maximum PWV values reach high levels that are all more than 60 mm, with one value even reaching 64.4 mm. The PWV time variation is also large and can reach 4.3 mm within 1 h before the rainfall. In addition, most of the rates of PWV time variation are greater than 1 mm/h. The results indicate that, as Hourly PWV retrieved from seven GNSS sites and actual precipitation from 27 June to 2 July 2017 (the red circle represents PWV variation before rainfall events, and the black circle represents PWV variation after rainfall events) one of the main preconditions for rainfall, GNSS-PWV quickly increases in a few hours before rainfall. In Table 6, after rainfall the minimum PWV values reach low levels that are all near 50 mm. Moreover, the rate of PWV descent can reach 6.5 mm/h after the rainfall. The results also prove that GNSS-PWV rapidly decreases for a few hours after a rainfall.
To further analyze the time varying characteristics of GNSS-PWV during the rainfall, the hourly PWV at 11 GNSS stations from 27 June to 6 July 2017 in Guilin are selected, as shown in Fig. 8. Figure 8 shows the GNSS-PWV time variation in the selected period. The PWVs retrieved at 11 GNSS stations generally have the same change process. The GNSS can effectively capture the significant variation of PWV over time in Guilin. Before the heavy rainfall, the PWV fluctuated slightly from 27 to 30 June  and remained stable at approximately 60 mm in the high water vapor content state. However, the PWV generally began to show a significant upward trend after 30 June and reached its peak value on 1 July, when the heavy rainfall occurred. After the heavy rainfall, the PWV dropped sharply to its minimum value on 2 July. Then, it gradually increased in the following days. It is finally stabilized at approximately 60 mm. The significant variation of PWV is mainly concentrated in the three days of 30 June to 2 July (DOY 181-183), which is the period of torrential rain. In addition, the PWV derived from station JZ88 is significantly lower than that of other GNSS stations in Fig. 8. Among the 11 GNSS stations, the elevation of station JZ88 is much higher than other stations (see Table 1). This discovery is consistent with the previous studies that the GNSS-PWV tends to decrease with the increasing elevation Onn & Zebker, 2006). From the above analysis, there is a strong correlation between GNSS-derived PWV and rainfall. To some extent the rapid rise of PWV may be a warning of a coming rainfall event since the occurrence of a rain requires sufficient water vapor supply. Though not all high PWV values indicate the occurrence of rainfalls (Yao et al., 2017), the increase of water vapor is one of the necessary pre-conditions for the occurrence of a rain. When rainfall begins to weaken, the PWV drops sharply. Then, the rain finally ceases. Afterward, the consumed water vapor is replenished, and PWV rises gradually until it is stable at a certain value.

PWV space characteristics
The heavy rainfall in the most areas of Guilin is mainly concentrated on 1 July 2017 from the meteorological observations. It is a long process and has a wide range and great intensity. Hence, in this study, we selected the hourly GNSS-PWVs at 10 GNSS stations during the three days from 30 June to 2 July 2017 (the observations of station JZ86 were unavailable in this period). Spatial interpolation is carried out to construct the distribution maps of PWV and to retrieve the corresponding actual precipitation at different times during the three-day Table 5 Statistical PWV variation before several rainfall events (they are also shown by red circles in Fig. 7  rainfall events. Figure 9 exhibits the geographic distribution of the daily accumulated precipitation over Guilin for 30 June to 2 July 2017. The precipitation data are received from the 12 MET stations. Figure 10 shows the spatial evolution of PWV derived from the observations at 10 GNSS stations during the same period, and the spatial evolution of ERA5 PWV is also presented in Fig. 11 for comparison. The PWV distribution maps are hourly, and we present the PWV maps with a time interval of 6 h because of space constraints. In addition, the UTC is adopted in the precipitation maps and PWV maps. As shown in Fig. 9, several rainfall events were observed over the most parts of Guilin. On 30 June, the precipitation increased from about 5 mm at the southeast to 65 mm at the northwest (see Fig. 9a). Then, the rainfall was getting heavy and developed into torrential rain on 1 July over the most parts of Guilin, particularly in western Guilin, even reaching its highest precipitation of 250 mm (see Fig. 9b). The accumulated precipitation from 30 June to 1 July was increasing. On 2 July, the precipitation was getting light and even ceased on the most area except for several rainfall events in the south Guilin (see Fig. 9c).
On 30 June, the GNSS-PWV experienced a significant increase from the southwest Guilin to the central part from 00:00 to 24:00 UTC, and it spread upward from southwest to north of Guilin until it covered the whole Guilin (see Fig. 10a-e), indicating that a large amount of water vapor from the southwest flowed into Guilin. The GNSS-PWV overall increased from 60 to 65 mm. This is consistent with the precipitation pattern displayed in Fig. 9a, wherein the west Guilin has more precipitation than the east. The spatial variations of the GNSS-PWV agree well with the ERA5 PWV displayed in Fig. 11a-e in that the PWV gradually increased from southwest to the central part of Guilin. The increase of water vapor transport from southwest offered a favorable condition for the heavy rainfall. High moisture conditions portended that the heavy rainfall was coming soon. In addition, the southwest jet stream provides sufficient moisture and dynamic conditions for the heavy rainfalls in Guilin according to the report from Guilin Meteorological Bureau. Similarly, the PWV in the southwest increased relatively earlier than the other parts of Guilin, as presented in Figs. 10 and 11. On 1 July, Fig. 10e-f shows that Guilin has high PWV values from 0:00 to 6:00 UTC, and the GNSS-PWV value was approximately 65 mm overall. The water vapor around Guilin gradually gathered and then covered the entire Guilin. In this circumstance, torrential rains occurred in the most parts of Guilin and reached the maximum precipitation on this day (see Fig. 9b). The GNSS-PWV values had small fluctuations but did not change much during this period. From 6:00 to 24:00 UTC, as shown in Fig. 10f-i, the GNSS-PWV in southern and northern Guilin gradually showed a trend of polarization, that is, the significant GNSS-PWV decreases were observed in the northern, while the southern experienced a slight increase in GNSS-PWV. Finally, four gradient forms of GNSS-PWV values emerged in Guilin at 24:00 UTC, which are approximately 35-45 mm, 45-55 mm, Fig. 9 Evolution of actual precipitation in Guilin from June 30 to July 2, 2017 55-65 mm, and 65-70 mm from northwest to southeast, respectively. The phenomenon is consistent with the precipitation pattern displayed in Fig. 9b, wherein the rainfall gradually decreased from west to east. Different intensities of rainfall in different regions caused the multilevel differentiation of PWV values over Guilin. Meanwhile, the GNSS-PWV maps shown in Fig. 10e-i also agree well with the ERA5 PWV variations, but GNSS-PWV seems a little larger than ERA5 PWV due to the lower spatial resolution and the different heights in PWV.
On 2 July, as shown in Fig. 10i-l, the GNSS-PWV differentiation was gradually weakened until the overall PWV values became stable at approximately 55 mm. The rainfalls gradually weakened and ceased in the most parts of Guilin, whereas the south experienced several rainfall events (see Fig. 9c). Unlike GNSS-PWV maps, the ERA5 PWV maps show the PWV in the south Guilin is significantly greater than that in the north from 12:00 to 18:00 UTC in Fig. 11k-l, which is more consistent with the actual precipitation as shown in Fig. 9c. This is benefited from the fact that ERA5 has a higher spatial resolution than GNSS, which can reduce the error caused by spatial interpolation.
As stated above, we can find the corresponding relationship between the accumulation and release of PWV and the actual precipitation process. In general, the accumulation of moisture increased before the torrential rain, and the PWV also increased largely, reaching the peak value afterward. At the peak value of the PWV, the rainfall arrived after a short period. As shown in Figs. 10 and 11, the high level of PWV was sustained for nearly a day. Then, Guilin experienced its heaviest rainfall during the three-day rainfall events. Therefore, the longer the high level PWV was sustained, the heavier the rainfall will be. During the torrential rain, the PWV values had slight fluctuations. The torrential rain lasted for a while, and continuous heavy rainfall caused a large amount of moisture loss. We can see that the PWV decreased rapidly and dropped to the valley value after a period of heavy rainfall. Finally, the PWV increased slightly and was gradually stabilized to a fixed value.
Moreover, the PWV varied from place to place. Guilin has a subtropical monsoon climate due to the difference in land and sea thermal properties, the prevailing wind of warm moist air mostly comes from the southwest and southeast in summer, and the precipitation is degressive from the southeast to northwest. The PWV in northwest Guilin is approximately 5 mm smaller than southeast Guilin because of its proximity to the inland. The increase of PWV mainly starts from south Guilin, which may be influenced by the subtropical monsoon climate and geographical location. In addition, we mentioned that GNSS stations in higher elevations tend to have smaller GNSS-PWV in the previous section. The GNSS-PWV maps also show that the PWV near the station JZ88 is usually lower than that in the surrounding areas in Guilin, as shown in Fig. 10. The spatial evolution of GNSS-PWV and ERA5 PWV is generally consistent except several small discrepancies, indicating the ERA5 PWV can also reveal the moisture variations during the rainfall events in Guilin, which can be combined with GNSS-PWV to analyze the moisture variations with higher spatial resolution in future work.

Conclusions
The distribution condition of moisture must be obtained, and the impacts of the spatiotemporal characteristics of moisture on the evolution of severe weathers and global climate change must be understood. Monitoring the climate change in Guilin is limited because of insufficient observations. As the ground-based GNSS rapidly developed over the past few decades, such measurements have become crucial for measuring highly dynamic water vapor. In this research, GNSS-derived PWV from the observations at 11 GNSS stations in Guilin are used to analyze the spatial and temporal characteristics of GNSS water vapor during heavy rainfalls from June to July 2017 in Guilin. PWV is the indicator used for monitoring heavy rainfall events.
To improve the estimation accuracy of the GNSS PWV retrieval, a new T m model is established for Guilin using the radiosonde profiles from 2012 to 2017. The results show that the new T m model has better performance than the widely used Bevis, Li, and GPT2w-1 T m models over Guilin with the mean bias and RMS of − 0.51 and 2.12 K, respectively. Additionally, the impact of different T m models on the resultant GNSS PWV is investigated, showing that the mean RMS pwv and RMS pwv /PWV are 0.29 mm and 0.74%, respectively.
To assess the reliability of GNSS-PWV, the GNSS-PWV retrieved from GNSS station JZ87 is compared with the RS-PWV, and a good agreement with them is obtained. The correlation coefficient between the two kinds of PWVs is high of 0.79, and the mean bias and RMS of GNSS-PWV are − 0.9 and 3.53 mm, respectively. According to the analysis of the spatial and temporal characteristics of water vapor during the heavy rainfall events that occurred in Guilin from June 27 to July 2, 2017, the variation of PWV derived from GNSS had a direct relationship with the rainfall events. PWV increases before the arrival of a heavy rainfall and decreases to a stable value after its cease. The GNSS-PWV can effectively reveal the significant variation of water vapor and the actual rainfall in several regions of Guilin. Furthermore, given the subtropical monsoon climate and geographical location, the PWV in Guilin mainly starts to increase from the south before a rainfall in summer.
This research demonstrates the possibility of revealing the moisture transportation and variation during heavy rainfalls using GNSS-derived PWV, which is of great significance and an indication for the prediction of rainfall and severe weather. In future, we will mainly focus on the following four issues: (1) considering additional factors that affect moisture, such as geography, pressure, temperature, and humidity; (2) using GNSS tomographic technique to investigate the relationship between the three-dimensional PWV and rainfall; (3) combining multisource data in the construction of PWV maps and improving its reliability; and (4) retrieving PWV with higher spatial and temporal resolutions.