Improved Drought Monitoring Index Using GNSS-Derived Precipitable Water Vapor over the Loess Plateau Area

Standardized precipitation evapotranspiration index (SPEI) is an acknowledged drought monitoring index, and the evapotranspiration (ET) used to calculated SPEI is obtained based on the Thornthwaite (TH) model. However, the SPEI calculated based on the TH model is overestimated globally, whereas the more accurate ET derived from the Penman–Monteith (PM) model recommended by the Food and Agriculture Organization of the United Nations is unavailable due to the lack of a large amount of meteorological data at most places. Therefore, how to improve the accuracy of ET calculated by the TH model becomes the focus of this study. Here, a revised TH (RTH) model is proposed using the temperature (T) and precipitable water vapor (PWV) data. The T and PWV data are derived from the reanalysis data and the global navigation satellite system (GNSS) observation, respectively. The initial value of ET for the RTH model is calculated based on the TH model, and the time series of ET residual between the TH and PM models is then obtained. Analyzed results reveal that ET residual is highly correlated with PWV and T, and the correlate coefficient between PWV and ET is −0.66, while that between T and ET for cases of T larger or less than 0 °C are −0.54 and 0.59, respectively. Therefore, a linear model between ET residual and PWV/T is established, and the ET value of the RTH model can be obtained by combining the TH-derived ET and estimated ET residual. Finally, the SPEI calculated based on the RTH model can be obtained and compared with that derived using PM and TH models. Result in the Loess Plateau (LP) region reveals the good performance of the RTH-based SPEI when compared with the TH-based SPEI over the period of 1979–2016. A case analysis in April 2013 over the LP region also indicates the superiority of the RTH-based SPEI at 88 meteorological and 31 GNSS stations when the PM-based SPEI is considered as the reference.


Introduction
Climate change in drought/humidity is closely related to human activity and production and has an important impact on socioeconomic and regional ecological development [1,2]. Some drought indexes, such as Palmer drought severity index (PDSI), standardized precipitation index (SPI), and standardized precipitation evapotranspiration index (SPEI), have been proposed to monitor and quantify drought [1,[3][4][5]. Drought is a water deficit phenomenon caused by low precipitation in a given period of time [4]. Under the background of global change characterized by warming, changes in temperature and precipitation and other factors affect the balance of surface water revenue and expenditure, which ultimately changes the surface drought/humidity situation [6,7].
From northwest to southeast, the maximum height difference is approximately 2000 m [32]. The climate of the LP belongs to the typical continental monsoon climate, which is mainly affected by the polar dry and cold airflow in winter and spring, and the climate is cold and windy. In summer and autumn, the climate is mainly affected by the high western Pacific subtropics and the low pressure in the Indian Ocean, which is hot and rainy [33]. The annual average precipitation is approximately 466 mm, with evident spatial variation. The annual precipitation increases from northwest to southeast, with rainfall of 200-750 mm, and 65% of the precipitation distributes from June to September. The estimated annual potential ET is considerably higher than that of precipitation, ranging from 865 mm to 1274 mm. The northwest region is a typical semi-arid region, which makes it a typical region for drought monitoring [33].

Retrieval of GNSS and ERA-Interim PWV
PWV can be obtained from some techniques, such as radio sounding, remote sensing, and lidar but with corresponding shortcomings such as low temporal-spatial resolutions, high cost, and low precision [34]. GNSS technique can also retrieve PWV with the advantages of low cost, all-weather conditions, and high precision. However, zenith total delay (ZTD) is first obtained before calculating the PWV using the GNSS technique. The ZTD is an average value by projecting some values of slant path with different azimuths and elevation angles into a vertical direction, and the accuracy of GNSS-derived ZTD is approximately 4 mm [35]. When the radio signal crosses the troposphere, the corresponding delay will be generated, which is influenced by the atmospheric refraction effect [36]. In this study, the ZTD data is calculated using the GNSS observations derived from Crustal Movement Observation Network of China (CMONOC) based on GAMIT/GLOBK software (http://www-gpsg. mit.edu/~{}simon/gtgk/down.htm). ZTD consists of zenith hydrostatic delay (ZHD) and zenith wet delay (ZWD), where the ZHD can be precisely calculated based on an empirical model, such as the Saastamoinen model [37]. Therefore, the ZWD parameter, which is used to calculate PWV, can be obtained by extracting ZHD from ZTD. The specific formula for calculating PWV from ZWD is as follows [38]: where T m is the water-vapor-weighted atmospheric mean temperature, which can be calculated by the empirical model based on the observed temperature [39,40]. k 2 , k 3 and R are the physical constants with values of 16.25 K/hPa, 3.776 × 10 5 K 2 /hPa, and 461.495 J/K/kg, respectively. The PWV value of the meteorological stations is interpolated using the PWV data of the four surrounding grid points. This data is from the fourth-generation atmospheric reanalysis of the global climate (ECMWF ERA-Interim), which covers rich data from 1979 to the present. ERA-Interim provides PWV, pressure, temperature, and other meteorological variables at grid points with different horizontal resolution (0.125 • × 0.125 • to 3 • × 3 • ) globally. In this study, the monthly PWV values from the ERA-Interim reanalysis data with the spatial resolution of 0.25 • × 0.25 • were selected. The PWV values at meteorological stations are obtained based on the bilinear interpolation method using the above data.

Meteorological Data
Daily meteorological data of 88 stations that are evenly distributed in the LP region over the period of 1979-2016 were collected (Figure 1), which include highest and lowest temperature, precipitation, relative humidity, sunshine duration, and 2 m wind speed. These data are downloaded from the China Meteorological Sharing Network (http://data.cma.cn/) and have been rigorously checked by CMA's staff. The work includes checking and marking the wrong and missing data, as well as supplementing those data before release. The monthly corresponding data were obtained using the daily meteorological data. The monthly temperature and precipitation of the Yang Kun dataset derived from the Qinghai-Tibet Plateau Scientific Data Center (http://www.tpedatabase.cn/portal/MetaDataList.jsp) were also used to interpolate the corresponding meteorological data at GNSS stations [41]. The specific process of generating such a dataset can be found in [41] and [42]. The temperature and precipitation are interpolated using the corresponding data of the four surrounding grid points derived from the Yang Kun dataset to calculate the SPEI value at selected GNSS stations. from the China Meteorological Sharing Network (http://data.cma.cn/) and have been rigorously checked by CMA's staff. The work includes checking and marking the wrong and missing data, as well as supplementing those data before release. The monthly corresponding data were obtained using the daily meteorological data. The monthly temperature and precipitation of the Yang Kun dataset derived from the Qinghai-Tibet Plateau Scientific Data Center (http://www.tpedatabase.cn/portal/MetaDataList.jsp) were also used to interpolate the corresponding meteorological data at GNSS stations [41]. The specific process of generating such a dataset can be found in [41] and [42]. The temperature and precipitation are interpolated using the corresponding data of the four surrounding grid points derived from the Yang Kun dataset to calculate the SPEI value at selected GNSS stations.

Theory of ET and SPEI Calculation
ET is calculated based on the Thornthwaite (TH) and Penman-Monteith (PM) models, respectively, in this paper. TH model is widely used because it is easy to use and only needs monthly average temperature and latitude [43,44]. However, this model does not consider the influence of wind speed, air humidity, and other factors, and the value of ET is usually underestimated. Based on the energy balance equation and water vapor diffusion theory, the PM model not only considers the physiological characteristics of crops but also considers the changes in aerodynamic parameters. The calculated ET value based on the PM model is accurate. However, this model requires many meteorological parameters, such as the highest and lowest daily temperature, relative humidity, sunshine hours, and 2 m wind [17]. The specific processes of calculation ET based on TH and PM model have been given in Appendix A.
SPEI is an acknowledged good index for drought monitoring because it considers the superiority of PDSI and SPI. This index includes the temperature variability of PDSI and the multiscalar characteristic of SPI. Therefore, this index is suitable for drought monitoring with the rapidly changing climate scenarios. Before the calculation of SPEI, the difference between monthly precipitation and ST is first obtained, and then accumulating this difference in different time scales. Finally, the SPEI value can be calculated by standardized the accumulated difference. The specific calculation process has been given in Appendix B.

RTH Model
In this paper, the revised TH (RTH) model can be established following the three key steps: 1. Calculating the ET residual between TH and PM model:

Theory of ET and SPEI Calculation
ET is calculated based on the Thornthwaite (TH) and Penman-Monteith (PM) models, respectively, in this paper. TH model is widely used because it is easy to use and only needs monthly average temperature and latitude [43,44]. However, this model does not consider the influence of wind speed, air humidity, and other factors, and the value of ET is usually underestimated. Based on the energy balance equation and water vapor diffusion theory, the PM model not only considers the physiological characteristics of crops but also considers the changes in aerodynamic parameters. The calculated ET value based on the PM model is accurate. However, this model requires many meteorological parameters, such as the highest and lowest daily temperature, relative humidity, sunshine hours, and 2 m wind [17]. The specific processes of calculation ET based on TH and PM model have been given in Appendix A.
SPEI is an acknowledged good index for drought monitoring because it considers the superiority of PDSI and SPI. This index includes the temperature variability of PDSI and the multiscalar characteristic of SPI. Therefore, this index is suitable for drought monitoring with the rapidly changing climate scenarios. Before the calculation of SPEI, the difference between monthly precipitation and ST is first obtained, and then accumulating this difference in different time scales. Finally, the SPEI value can be calculated by standardized the accumulated difference. The specific calculation process has been given in Appendix B.

RTH Model
In this paper, the revised TH (RTH) model can be established following the three key steps:

1.
Calculating the ET residual between TH and PM model: 2.
Analyzing the time series of ET residual and fitting the ET residual model using the GNSS-derived PWV and temperature:

3.
Obtaining the initial ET value ET 0 using TH model and establishing the ETH model using the ET residual:

Missing Data Interpolation
The GNSS stations from the CMONOC are selected, where the GNSS-derived PWV is obtained following the method proposed by [45]. A total of 31 GNSS stations exist over the LP region (Figure 1), and the average data missing rate at those stations is 13.8%. Therefore, an interpolated method referred to as SSAM (singular spectrum analysis for missing data) proposed by Wang et al. was used to interpolate the missing data at each GNSS station, and this method works well in interpolating the missing data of long-term PWV time series [46]. 2. Analyzing the time series of ET residual and fitting the ET residual model using the GNSS-derived PWV and temperature: 3. Obtaining the initial ET value 0 ET using TH model and establishing the ETH model using the ET residual:

Missing Data Interpolation
The GNSS stations from the CMONOC are selected, where the GNSS-derived PWV is obtained following the method proposed by [45]. A total of 31 GNSS stations exist over the LP region (Figure 1), and the average data missing rate at those stations is 13.8%. Therefore, an interpolated method referred to as SSAM (singular spectrum analysis for missing data) proposed by Wang et al. was used to interpolate the missing data at each GNSS station, and this method works well in interpolating the missing data of long-term PWV time series [46].   Figure 3a,b presents the calculated average ET values at 88 meteorological stations using TH and PM models over the period of 1979-2016. It can be observed that the accuracy of TH-derived ET values is poor at some stations when the PM-derived ET values are regarded as the reference. In addition, the ET differences of selected stations between the PM and TH models have been presented in the LP area. It can be observed that the maximum ET difference is almost 45 mm while the minimum ET difference is approximately 27 mm. Such a phenomenon demonstrates the necessity of establishing the improved model to correct the ET value.  Figure 3a,b presents the calculated average ET values at 88 meteorological stations using TH and PM models over the period of 1979-2016. It can be observed that the accuracy of TH-derived ET values is poor at some stations when the PM-derived ET values are regarded as the reference. In addition, the ET differences of selected stations between the PM and TH models have been presented in the LP area. It can be observed that the maximum ET difference is almost 45 mm while the minimum ET difference is approximately 27 mm. Such a phenomenon demonstrates the necessity of establishing the improved model to correct the ET value. Here, the ET residual between PM and TH models are obtained at 83 meteorological stations, and the relationships between ET residual and PWV/T are presented in Figure 4. It can be concluded that the ET residual has an evident negative linear relationship with PWV and T when T > 0 °C, whereas a positive linear relationship exists between ET residual and T when T < 0 °C. The different Here, the ET residual between PM and TH models are obtained at 83 meteorological stations, and the relationships between ET residual and PWV/T are presented in Figure 4. It can be concluded that the ET residual has an evident negative linear relationship with PWV and T when T > 0 • C, whereas a positive linear relationship exists between ET residual and T when T < 0 • C. The different phenomena presented between T and ET residual for cases of T > 0 • C and T < 0 • C is mainly related to the 0 values of ET when T < 0 • C (Appendix A). The relationship between ET residual and altitude of meteorological stations is also analyzed, and no evident relationship exists. Therefore, a linear model between ET residual and PWV/T is established and expressed as follows: Here, the ET residual between PM and TH models are obtained at 83 meteorological stations, and the relationships between ET residual and PWV/T are presented in Figure 4. It can be concluded that the ET residual has an evident negative linear relationship with PWV and T when T > 0 °C, whereas a positive linear relationship exists between ET residual and T when T < 0 °C. The different phenomena presented between T and ET residual for cases of T > 0 °C and T < 0 °C is mainly related to the 0 values of ET when T < 0 °C (Appendix A). The relationship between ET residual and altitude of meteorological stations is also analyzed, and no evident relationship exists. Therefore, a linear model between ET residual and PWV/T is established and expressed as follows:

Validation of the RTH Model
To validate the proposed RTH model using the PWV and T over the LP region, only the meteorological data over the period of 1979-2014 were used to fit the coefficients of the ET residual model, whereas the data over the period of 2015-2016 were used to perform the validation experiment at 88 meteorological stations. Here, the ET values calculated by the PM model are regarded as reference. Figure 5 presents the comparisons of RMS and mean absolute error (MAE) of the ET residual between PM and TH/RTH models at 88 stations. It can be found that the values of RMS and MAE of the RTH model are smaller than those from the TH model at all stations when the values derived from the PM model are considered as reference. The statistical result reveals that the RMS and MAE derived from the RTH model are 12.3 and 10.2 mm, respectively, which are smaller than that from the TH model (34.2 and 29.2 mm, respectively). The RMS improvement rate of the RTH model is also calculated when compared with the TH model, as shown in Figure 6. It can be

Validation of the RTH Model
To validate the proposed RTH model using the PWV and T over the LP region, only the meteorological data over the period of 1979-2014 were used to fit the coefficients of the ET residual model, whereas the data over the period of 2015-2016 were used to perform the validation experiment at 88 meteorological stations. Here, the ET values calculated by the PM model are regarded as reference. Figure 5 presents the comparisons of RMS and mean absolute error (MAE) of the ET residual between PM and TH/RTH models at 88 stations. It can be found that the values of RMS and MAE of the RTH model are smaller than those from the TH model at all stations when the values derived from the PM model are considered as reference. The statistical result reveals that the RMS and MAE derived from the RTH model are 12.3 and 10.2 mm, respectively, which are smaller than that from the TH model (34.2 and 29.2 mm, respectively). The RMS improvement rate of the RTH model is also calculated when compared with the TH model, as shown in Figure 6. It can be observed that the RMS improvement rate at different stations ranges from 22% to 75%, and the statistical result reveals that the average improvement rate of the RTH model is approximately 61.6%. observed that the RMS improvement rate at different stations ranges from 22% to 75%, and the statistical result reveals that the average improvement rate of the RTH model is approximately 61.6%.  To further evaluate the complete accuracy of the established RTH model over the LP region, the scatter plot of monthly ET values calculated by TH, RTH, and PM models is shown in Figure 7. It can be observed that the RTH-derived ET has a better agreement with PM-derived ET when compared with that from TH-derived ET. The correlation coefficient between PM-and RTH-derived ET is 0.98, observed that the RMS improvement rate at different stations ranges from 22% to 75%, and the statistical result reveals that the average improvement rate of the RTH model is approximately 61.6%.  To further evaluate the complete accuracy of the established RTH model over the LP region, the scatter plot of monthly ET values calculated by TH, RTH, and PM models is shown in Figure 7. It can be observed that the RTH-derived ET has a better agreement with PM-derived ET when compared with that from TH-derived ET. The correlation coefficient between PM-and RTH-derived ET is 0.98, To further evaluate the complete accuracy of the established RTH model over the LP region, the scatter plot of monthly ET values calculated by TH, RTH, and PM models is shown in Figure 7. It can be observed that the RTH-derived ET has a better agreement with PM-derived ET when compared with that from TH-derived ET. The correlation coefficient between PM-and RTH-derived ET is 0.98, whereas the value between PM-and PH-derived ET is 0.93. The long-term time series of average ET calculated  Figure 8. The time series of ET derived from the RTH model agrees well with that derived from the PM model, whereas the value of the TH model is relatively poor, especially for those times with low ET values, where the ET differences between PM and TH models are large. However, the RTH model established in this study offers an improvement for those values well, which is also the advantage of the RTH model. whereas the value between PM-and PH-derived ET is 0.93. The long-term time series of average ET calculated using different methods over the period of 1979-2016 in the LP region is shown in Figure 8. The time series of ET derived from the RTH model agrees well with that derived from the PM model, whereas the value of the TH model is relatively poor, especially for those times with low ET values, where the ET differences between PM and TH models are large. However, the RTH model established in this study offers an improvement for those values well, which is also the advantage of the RTH model.

Evaluation of RTH-Based SPEI at Meteorological Stations
After the ET value is obtained using the RTH model, the final SPEI can be calculated, which is called RTH-based SPEI in this study. In addition, the ET value is also calculated based on the PM and TH models. Therefore, the SPEI can also be calculated by such models, which is called PM-and THbased SPEI. The quality of the RTH-based SPEI is the key to determining whether it can be further used. Therefore, the accuracy of RTH-based SPEI is initially evaluated and compared with TH-based SPEI when the PM-based SPEI is considered as reference. whereas the value between PM-and PH-derived ET is 0.93. The long-term time series of average ET calculated using different methods over the period of 1979-2016 in the LP region is shown in Figure 8. The time series of ET derived from the RTH model agrees well with that derived from the PM model, whereas the value of the TH model is relatively poor, especially for those times with low ET values, where the ET differences between PM and TH models are large. However, the RTH model established in this study offers an improvement for those values well, which is also the advantage of the RTH model.

Evaluation of RTH-Based SPEI at Meteorological Stations
After the ET value is obtained using the RTH model, the final SPEI can be calculated, which is called RTH-based SPEI in this study. In addition, the ET value is also calculated based on the PM and TH models. Therefore, the SPEI can also be calculated by such models, which is called PM-and THbased SPEI. The quality of the RTH-based SPEI is the key to determining whether it can be further used. Therefore, the accuracy of RTH-based SPEI is initially evaluated and compared with TH-based SPEI when the PM-based SPEI is considered as reference.

Evaluation of RTH-Based SPEI at Meteorological Stations
After the ET value is obtained using the RTH model, the final SPEI can be calculated, which is called RTH-based SPEI in this study. In addition, the ET value is also calculated based on the PM and TH models. Therefore, the SPEI can also be calculated by such models, which is called PM-and TH-based SPEI. The quality of the RTH-based SPEI is the key to determining whether it can be further used. Therefore, the accuracy of RTH-based SPEI is initially evaluated and compared with TH-based SPEI when the PM-based SPEI is considered as reference.
The RTH-based SPEI is initially compared with TH-based SPEI at 88 GNSS stations over the period of 1979-2016, and the average SPEI differences between RTH-PM and TH-PM under 1-, 3-, 6-, and 12-month scales are presented in Figure 9. Those multi-month scales are selected and compared Sensors 2019, 19, 5566 9 of 17 because they correspond to the monthly, seasonal, semi-annual, and annual scales. It can be observed from Figure 9 that the SPEI differences of RTH-PM are smaller than those from TH-PM, especially for the 1-, 3-, and 6-month scales, which further indicates the good performance of the RTH-based SPEI at those month scales.
Pearson's correlation is then introduced in comparison to further analyze the relationship between SPEI calculated based on different models at different multi-month scales. Figure 10 shows the Pearson's correlations of SPEI calculated based on different models under different multi-month scales from 1 to 24. The correlation between RTH-and PM-based SPEI is high at different month scales, especially for the 12-, 23-, and 24-month scales. Such a result conforms to the conclusion derived from Figure 9. In addition, it can be observed that the correlation values between TH-PM and RTH-PM are similar at 12-, 23-, and 24-month scales. This is because the monthly climatic water balance is standardized during the process of calculating SPEI (Appendix B). In those three scales, the changing trend of accumulated climatic water balances derived from TH and RTH models are both similar to those from the PM model. However, the values of the climatic water balance derived from the RTH model are closer to that from the PM model when compared to the TH model.   Pearson's correlation is then introduced in comparison to further analyze the relationship between SPEI calculated based on different models at different multi-month scales. Figure 10 shows the Pearson's correlations of SPEI calculated based on different models under different multi-month scales from 1 to 24. The correlation between RTH-and PM-based SPEI is high at different month scales, especially for the 12-, 23-, and 24-month scales. Such a result conforms to the conclusion derived from Figure 9. In addition, it can be observed that the correlation values between TH-PM and RTH-PM are similar at 12-, 23-, and 24-month scales. This is because the monthly climatic water balance is standardized during the process of calculating SPEI (Appendix B). In those three scales, the changing trend of accumulated climatic water balances derived from TH and RTH models are both similar to those from the PM model. However, the values of the climatic water balance derived from the RTH model are closer to that from the PM model when compared to the TH model.
The comparison of SPEI calculated using different models is also performed at each meteorological station under multi-month scales. The RMS values are presented in Figure 11. It can be observed that the RMS values of RTH-based SPEI are smaller than those of TH-based SPEI at all stations. Table 1 shows the statistical result of the average RMS and MAE of TH-and RTH-based SPEI when compared with the PM-based SPEI. The result shows that good performance can be obtained for the RTH-based SPEI. Finally, the average RMS improvement rate of the RTH-based SPEI compared with the TH-based SPEI is analyzed and presented in Figure 12 when the PM-based SPEI is considered as reference. It can be concluded that the accuracy of SPEI is improved to different degrees at different stations. The statistical result reveals that the average improvement rate of the RMS derived from the RTH-based SPEI is approximately 48.0%, 49.3%, 43.2%, and 9.1% under 1-, 3-, 6-, and 12-month scales in the LP region, respectively.   The comparison of SPEI calculated using different models is also performed at each meteorological station under multi-month scales. The RMS values are presented in Figure 11. It can be observed that the RMS values of RTH-based SPEI are smaller than those of TH-based SPEI at all stations. Table 1 shows the statistical result of the average RMS and MAE of TH-and RTH-based SPEI when compared with the PM-based SPEI. The result shows that good performance can be obtained for the RTH-based SPEI. Finally, the average RMS improvement rate of the RTH-based SPEI compared with the TH-based SPEI is analyzed and presented in Figure 12 when the PM-based SPEI is considered as reference. It can be concluded that the accuracy of SPEI is improved to different degrees at different stations. The statistical result reveals that the average improvement rate of the RMS derived from the RTH-based SPEI is approximately 48.0%, 49.3%, 43.2%, and 9.1% under 1-, 3-, 6-, and 12-month scales in the LP region, respectively.

Case of Spatial Analysis of RTH-Based SPEI
The temperature and precipitation derived from the Yang Kun dataset is initially validated using the corresponding data derived from CMA at 88 meteorological stations. Figure 13 shows the scatter plots of temperature and precipitation over the period of 1979-2016 in the LP region. It can be observed that the temperature derived from the Yang Kun dataset agrees well with that from CMA, whereas the precipitation is relatively poor with that from CMA. The statistical result reveals that the RMS and bias of temperature and precipitation from the Yang Kun dataset are 3.85/−0.92 °C and 0.39/−0.34 mm, respectively. Therefore, the RTH-based SPEI can be calculated at GNSS stations using the corresponding data from the Yang Kun dataset and GNSS-derived PWV. Figure 14 presents an example of a comparison of the RTH-and TH-based SPEI at XNIN Station over the period of 1999-2014, where the PM-based SPEI cannot be calculated due to lack of corresponding data.

Case of Spatial Analysis of RTH-Based SPEI
The temperature and precipitation derived from the Yang Kun dataset is initially validated using the corresponding data derived from CMA at 88 meteorological stations. Figure 13 shows the scatter plots of temperature and precipitation over the period of 1979-2016 in the LP region. It can be observed that the temperature derived from the Yang Kun dataset agrees well with that from CMA, whereas the precipitation is relatively poor with that from CMA. The statistical result reveals that the RMS and bias of temperature and precipitation from the Yang Kun dataset are 3.85/−0.92 • C and 0.39/−0.34 mm, respectively. Therefore, the RTH-based SPEI can be calculated at GNSS stations using the corresponding data from the Yang Kun dataset and GNSS-derived PWV. Figure 14 presents an example of a comparison of the RTH-and TH-based SPEI at XNIN Station over the period of 1999-2014, where the PM-based SPEI cannot be calculated due to lack of corresponding data.

Case of Spatial Analysis of RTH-Based SPEI
The temperature and precipitation derived from the Yang Kun dataset is initially validated using the corresponding data derived from CMA at 88 meteorological stations. Figure 13 shows the scatter plots of temperature and precipitation over the period of 1979-2016 in the LP region. It can be observed that the temperature derived from the Yang Kun dataset agrees well with that from CMA, whereas the precipitation is relatively poor with that from CMA. The statistical result reveals that the RMS and bias of temperature and precipitation from the Yang Kun dataset are 3.85/−0.92 °C and 0.39/−0.34 mm, respectively. Therefore, the RTH-based SPEI can be calculated at GNSS stations using the corresponding data from the Yang Kun dataset and GNSS-derived PWV. Figure 14 presents an example of a comparison of the RTH-and TH-based SPEI at XNIN Station over the period of 1999-2014, where the PM-based SPEI cannot be calculated due to lack of corresponding data.  According to the recordings of historical drought data of China meteorological network (https://cmdp.ncc-cma.net/cn/index.htm), a large surface drought occurred in the LP region in April 2013. Therefore, this month is selected to analyze the SPEI calculated based on the RTH model. Figure 15 shows the comparison result of the SPEI calculated using different models at GNSS and meteorological stations in April 2013. Due to the lack of corresponding meteorological data at GNSS stations, the PM-based SPEI cannot be obtained at the GNSS stations, which is also the disadvantage of the PM-based SPEI. It can be observed from Figure 15 that the RTH-based SPEI is superior to TH-based SPEI under different month scales when the PM-based SPEI is considered as reference, especially in 3-, 6-, and 9-month scales. This result indicates the good performance of the ET estimated based on the RTH model in this study. According to the recordings of historical drought data of China meteorological network (https://cmdp.ncc-cma.net/cn/index.htm), a large surface drought occurred in the LP region in April 2013. Therefore, this month is selected to analyze the SPEI calculated based on the RTH model. Figure 15 shows the comparison result of the SPEI calculated using different models at GNSS and meteorological stations in April 2013. Due to the lack of corresponding meteorological data at GNSS stations, the PM-based SPEI cannot be obtained at the GNSS stations, which is also the disadvantage of the PM-based SPEI. It can be observed from Figure 15 that the RTH-based SPEI is superior to THbased SPEI under different month scales when the PM-based SPEI is considered as reference, especially in 3-, 6-, and 9-month scales. This result indicates the good performance of the ET estimated based on the RTH model in this study. According to the recordings of historical drought data of China meteorological network (https://cmdp.ncc-cma.net/cn/index.htm), a large surface drought occurred in the LP region in April 2013. Therefore, this month is selected to analyze the SPEI calculated based on the RTH model. Figure 15 shows the comparison result of the SPEI calculated using different models at GNSS and meteorological stations in April 2013. Due to the lack of corresponding meteorological data at GNSS stations, the PM-based SPEI cannot be obtained at the GNSS stations, which is also the disadvantage of the PM-based SPEI. It can be observed from Figure 15 that the RTH-based SPEI is superior to THbased SPEI under different month scales when the PM-based SPEI is considered as reference, especially in 3-, 6-, and 9-month scales. This result indicates the good performance of the ET estimated based on the RTH model in this study.

Conclusions
An RTH model is proposed in this study and the RTH model is established based on the analysis of the relationship between PWV/T and ET residual. The initial value of ET for the RTH model is calculated using the TH model, and the time series of ET residual between TH and PM models are fitted based on multiple linear regression. An experiment is performed in the LP region, and 88 meteorological stations that are evenly distributed in the LP region over the period of 1979-2016 are selected. 31 GNSS stations derived from CMONOC located in this region are also used to analyze the performance of the RTH-based SPEI.
Comparison results show that the ET estimated based on the RTH model is superior to the traditional TH model, and the statistical result reveals that the average improvement rate of RMS derived from RTH-based ET is approximately 61.6% in the LP region. The RTH-based SPEI is validated and compared with the TH-based SPEI when the PM-based SPEI is regarded as a reference. The numerical result indicates the good performance of the RTH-based SPEI, whereas the quality of the TH-based SPEI is relatively poor at the selected meteorological stations. The statistical result reveals that the average improvement rate of the RMS derived from the RTH-based SPEI is approximately 48.0%, 49.3%, 43.2%, and 9.1% under 1-, 3-, 6-, and 12-month scales in the LP region, respectively. Finally, the RTH-based SPEI is tested in a case of April 2013 using the meteorological and GNSS stations. The results obtained above indicate the improved accuracy of the RTH-based SPEI, and further, verify the superiority of the proposed RTH model in this paper.

Acknowledgments:
The author would like to thank Zhang et al. (2017) and ECMWF for providing highly precise GNSS PWV data. The CMA and Yang Kun are also acknowledged for proving the corresponding meteorological data.

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

Appendix A
ET value can be calculated on the basis of the monthly mean temperature based on the empirical model proposed by [16], which is usually called the TH model and can be expressed as where T is the monthly mean temperature in • C; I is the heat index, which is calculated as the sum of 12 monthly index values i; i is calculated using the mean monthly temperature by the first formula in Equation (A2); m is a coefficient that is dependent on I, as described by the second formula in Equation (A2); and K is a correction coefficient computed as a function of the latitude and month, which can be expressed by the third formula in Equation (A2).
where NDM is the number of days of the month; and N is the maximum number of sun hours, which is usually obtained as follows: N = (24/π) · w s (A3) where w s is the hourly angle of the sun rising, which is calculated by w s = arccos(− tan ϕ tan δ) (A4) where ϕ is the latitude in radians; and δ is the solar declination in radians, which can be obtained by where J is the average Julian day of the month. In addition to calculating the rough ET value based on the TH model, an accurate ETo value can be calculated using the PM model [17]. The PM model can be expressed as [47,48]: where T and u 2 are the 2 m average air temperature in • C and wind speed in m s −1 , respectively; R n is the crop surface net radiation in MJ m −2 d −1 ; G is the soil heat flux in MJ m −2 d −1 ; e s is the saturation vapor pressure in kPa; e a is the actual vapor pressure in kPa; ∆ is the slope of the vapor pressure curve in kPa • C −1 ; γ is a psychrometric constant in kPa • C −1 ; 900 is the coefficient for the reference crop; and 0.34 is the wind coefficient for the reference crop. The preceding eight variables can be calculated using the daily maximum and minimum temperature, relative humidity, sunshine duration, and 2 m wind following the formulas of [41].

Appendix B
SPEI can be calculated following the formula proposed by [5], where the monthly climatic water balance D i of month i is initially computed using the difference between precipitation P i and PET and expressed as follows: The calculated D i values are aggregated at different time scales. A good relationship exists between SPEI and the log-logistic distribution based on the Kolmogorov-Smirnov test under multi-month scale on a global scale [13], therefore, SPEI is calculated using the three-parameter log-logistic distribution based on the standardized D series. The probability distribution function of the log-logistic distribution for D series can be expressed as where α, β, and γ are the scale, shape, and origin parameters, respectively, which are obtained using the L-moment procedure [5].