GNSS-based PWV Application for Short Term Rainfall Prediction in Mountainous Region

Rainfall in the mountainous region has unique characteristics related to time-varying and spatial distribution. In Mt. Merapi region, located at the border of Special Region of Yogyakarta Province and Central Java Province, Indonesia, rainfalls are typically classified as the deep-convective type, which occurs in a short period with high intensity. Therefore, short term rainfall predictions in a proper way remain challenging tasks. The use of remote monitoring instruments such as the GNSS (Global Navigation Satellite System) is believed to provide a better measurement accuracy through the identification of water vapor variation in the process of deep convection weather. GNSS observes the geodetic position of the GNSS antenna or receiver while it broadcasts microwave signals continuously through the atmosphere to the ground-based receivers. As they travel through the atmosphere, the microwave signals are mostly influenced by ionospheric and neutral atmospheric effects, which cause some delays. By using a sufficiently dense network of GNSS receivers, the impact of the neutral atmosphere delay can be estimated as a by-product of the geodetic processing. These delays can be regarded as an integrated water vapor along the path, namely Precipitable Water Vapor (PWV), which can be indirectly measured by Zenith Total Delay (ZTD). By studying the relationship between time-varying PWV and rainfall, it can be found that the PWV level increases sharply before raining. Through the deployment of GNSS receivers, the spatial feature of rainfall characteristics is also depicted. The initial results showed that the increase of PWV is strongly correlated to rainfall occurrence based on the rain gauge measurement around Mt. Merapi region. The results show that the correct forecasted rate is about 47%-62% with the PWV increment time is three hours.


Introduction
Water vapor plays an important role in atmospheric processes such as rain, adverse weather, and global climatological condition [1]. An emerging method of Global Navigation Satellite System (GNSS) for meteorology triggers a new powerful approach to the remote sensing of atmospheric water vapor. The GNSS consists of satellites constellation, which transmits radio signal and the mounted GPS antenna, which receives the signal from the ground. Along the propagation path, the transmitted signal is mostly interrupted by the presence of atmospheric water vapor, causing a delay. The water vapor contained in the atmosphere is expressed as the height of an equivalent column of liquid water, namely Precipitable Water Vapor (PWV), which can be indirectly measured by the integration of a Zenith Total Delay (ZTD). The ZTD has a contribution from pressure, and thus, it has strong constrained with topography, while PWV has much less effect of altitude [2]. Therefore, PWV is more reliable for rainfall prediction in the mountainous region.
Researches related to PWV and rainfall occurrence relationships are mainly studied on the Mediterranean or humid subtropical region with a uniform topography [3][4][5]. However, the region, which has a tropical wet climate with complex terrain condition, has not been thoroughly studied. Rainfalls in the mountainous region are mostly influenced by the confusion of the hydrodynamic environments such as prevailing wind direction, air temperature collision, and many others. This circumstance features the characteristics of rainfall in the mountainous region, which mainly categorized as the deep-convective type where rainfalls occur in a short period with high intensity. The Mt. Merapi region, located at the border of Special Region of Yogyakarta Province and Central Java Province, Indonesia, has a typical tropical climate characteristic. Its relative humidity rate is about 70% to 90% with extreme variations in rainfall due to west monsoons of more than 2,000 mm/year [6]. With this condition, rainfall in Mt. Merapi region may hold a beneficial potency and a hazardous risk at the same time. Therefore, this study aims to investigate the characteristics of rainfall in a mountainous region with particularly having tropical climate environments through the application of GNSS-derived PWV. This study is intended to find the appropriate threshold of PWV to rainfall occurrence or no-occurrence so that earlier detection of rainfall can be well predicted.

Instrument specifications
The GNSS receivers in Mt. Merapi region has been set up since November 2017 with the primary purpose to measure the deformation of landform around mountain area during the volcano activities. There are four instruments of GNSS receivers with Leica GR10 type, divided into three instruments installed on the ledge area and one instrument installed on the ridge area of the mountain area. Besides, there is nine Automatic Rainfall Recorder (ARR) stations, which are sparsely installed on this study site.   figure 1, the topography condition of Mt. Merapi region with the distribution of rain measuring instruments installed is well depicted. By considering the coincide position between the GNSS receivers and the ARR stations and considering the data available from both instruments, there are two GNSS receivers and three ARR stations remaining that will be used for analyses. The GNSS receivers selected are Deles and BPTK stations and the ARR stations selected are BE-D4, GE-Kaliadem, and BO-CO1 stations. Here, Deles station position corresponds to BE-D4 and GE stations, while BPTK station position corresponds to BO-CO1 station. The ARR stations are selected by taking into account the influence of distance factors from the GNSS receivers. The altitude position of Deles station and BPTK station is located at +1,400 m and +112 m, respectively. The interval data of the GNSS product is 2 seconds, while the interval data of the ARR product is 5 minutes. For further analysis, both data products are converted into an hourly data interval.

Meso-Scale Model
The mesoscale numerical prediction system has been operated since 2001 to provide information for disaster prevention and aviation safety. The MSM product [7] includes atmospheric components such as temperature, surface pressure, relative humidity, and oceanic components. The horizontal resolution of the MSM is 10 km by covering 40 vertical layers. The time resolution produces 18-hour forecasts every 6 hours at 00, 06, 12, and 18 UTC. The MSM generated by Japan Meteorological Agency (JMA) is used to generate the atmospheric surface pressure, the surface temperature, the wind speed, and the wind direction. The atmospheric condition retrieved from MSM that covers the study site is shown in figure 2 and figure 3. Figure 2 and figure 3 shows surface pressure and relative humidity condition on December 20, 2017, respectively. These data are then used as an input to calculate the zenith hydrostatic delay and zenith wet delay conditions.

PWV retrieval
The PWV retrieval based on GNSS is mainly following the procedure (see figure 4): -obtaining the GNSS raw data from the GNSS receiver. The raw data products include Glonass Nav data (g.file), Navigation data (n.file), and Observation data (o.file), processing the acquired data using RTNet Software, -extracting the ZTD values, and -calculating the PWV.
The RTNet is designed for processing of GPS networks in the real-time application by using a squareroot filter processing. Real-time estimation of satellite clocks for regional or global networks is possible without the requirement of a stable reference clock. Therefore, the RTNet can estimate a large number of parameter types such as the characteristics of the ionospheric and tropospheric biases [8]. In order to measure the PWV contents, some factors that cause signal transmission delay are distinguished into the hydrostatic delay, which influenced by ionospheric factors and the wet delay, which influenced by atmospheric water vapor factors. The relationship among the parameters can be expressed as: Where ZTD is zenith total delay, H is zenith hydrostatic delay, and is zenith wet delay. All parameters are in the meter unit.
Where P is the station total pressure, is station latitude, and H is station height. Then integrated water vapor along the path, namely Precipitable Water Vapor (PWV) can be calculated by:

Evaluation Methods
The GPS data were processed using the RTNet software in an hourly interval. The ZHD was modeled throughout a global grid containing surface pressure and temperature data with 6-h output intervals based on the MSM analysis. The analysis was taken based on GNSS data retrieved from November 4, 2017, up to April 30, 2018. During this period, the number of rainfall events recorded at BPTK and Deles station was 32 and 53, respectively, which vary from moderate to very heavy rainfall (2 -60 mm/h) [9].
In general, the PWV's trend has a strong correlation with the rainfall occurrence of which the PWV increases before rainfall and then decreases when rainfall begins to weaken. In this section, three factors -PWV variation, the rate of change of PWV, and PWV-increment time are selected as rainfall forecasting factors. The PWV variation is described as a difference between the maximum PWV before rainfall and the minimum PWV before the adjacent maximum PWV. The PWV rate of change is described as a derivation of the PWV variation by increment time of the minimum to maximum adjacent PWV correspondingly. For the appropriate time forecast, the duration of the PWV increment time is used.
First, the rainfall occurrence and no-occurrence are classified (see table 1). In this table, "a" stands for "true positive" statement, "b" stands for "false positive" statement, "c" stands for "false negative" statement, and "d" stands for "true negative" statement. Second, the forecast indicators (i.e., correct rate, false alarm rate, and KSS value) are calculated. The correct rate has a range of 0 to 100% with 100% representing a perfect forecast. The equation is described as follows.
The false alarm rate is simply the fraction of observed non-events that are false alarms. The best score for the false alarm rate is close to 0, which means have as few false alarms as possible. The equation is described as follows.
To indicate whether the forecast occurrence or no-occurrence, the decision is made based on the Hanssen-Kuipers Score (KSS). This score measures the ability of the forecast to distinguish between occurrence and non-occurrence events (10). The best score is 1, which is obtained when the correct rate is 100%, and the false alarm rate is 0%. The equation is described as follows. Finally, the event forecast threshold is determined based on the forecast indicators score. The threshold of each PWV parameters is selected when the KSS value is closest to 1.

PWV and rainfall characteristics
Study of the relationship between PWV characteristics and rainfall occurrence is conducted at two GNSS receiver stations (i.e., BPTK and Deles station), which of each represents flatter and mountainous area. The analysis is conducted during the rainy season in 2017-2018. During this period, the rainfall events vary from medium to very heavy rainfall with the rain rates ranging from 2 to 60 mm/h, and the rainfall duration recorded is about 4 hours on average. The total events that are counted, including rain and norain occurrence are 43 events at BPTK station and 76 events at Deles station. The results of the PWV value and the corresponding rain rates are shown in figure 5 and figure 6. The average value of PWV at BPTK station is 61 mm and at Deles station is 35 mm. The altitude factor of each station installed mainly causes the difference value of PWV estimated at both stations. The PWV value represents the accumulated water vapor at various height. Therefore, the higher altitude tends to have lower PWV value and vice versa.
The results are shown in figure 5 and 6 reveal that the PWV will rapidly increase and aggregate in a few hours before rainfalls. However, the variation of PWV increase may give a different response to rain rate, which makes difficult to determine a direct correlation between the PWV value and the rain rate. As an example, at BPTK station, rainfall event on January 20, 2018 has a maximum rain rate 66 mm/h and a corresponding PWV value is 62 mm. With the same PWV value, the rainfall event on January 31, 2018 has a maximum rain rate 21 mm/h. The similar case is found at Deles station. Rainfall event on November 8, 2017 has a maximum rain rate value 47 mm/h, and a corresponding PWV value is 30 mm. With the same PWV value, the rainfall event on November 9, 2017 has a maximum rain rate 23 mm/h. At this location, it is found that in some cases, the PWV value tends to decrease before rainfalls occurred. It reveals that the PWV characteristics may give a different response to rainfall pattern depending on the position of station installed. From both GNSS receiver stations, we can distinguish the effect of mountainous and flatter area condition in terms of rainfall forecast  The characteristics of PWV-variation (∆PWV) before and after rainfall occurred can be seen from figure 7. The trend of ∆PWV at BPTK station is rapidly increased from time 12:13 to 14:13 before rain rate reaches the maximum value and it gradually increases as long as rainfall subsides, then it continuously decreases when the rainfall tends to stop. On the other hand, the trend of ∆PWV at Deles station has a gentler slope, which gives different rainfall pattern and intensity. The trend of ∆PWV is rapidly increased from time 15:13 to 17:13 before rain rate reaches the maximum value and it gradually increases as long as rainfall subsides, then it continuously decreases when the rainfall tends to stop.

PWV threshold determination
The threshold of PWV is determined based on the parameter of PWV-variation, PWV-rate of change, and PWV-increment time. Here, the interval of PWV-variation and PWV-rate of change is set to be 0.5 mm and 0.1 mm/h, respectively. While the PWV-increment time is calculated based on the selected threshold of PWV-variation and PWV-rate of change. Each interval is counted for all possible rainfall events, as stated in table 1. For example, for PWV-variation value 0.5 mm, every PWV-variation that has value ≥ 0.5 mm will be counted either as true positive, false positive, false negative, or as true negative. For PWV-variation value 1 mm, every PWV-variation that has value ≥ 1 mm will be counted either as true positive, false positive, false negative, or as true negative, and so forth. Yield, we can determine the optimum threshold for rainfall forecast. In the case of PWV-rate of change, the classification of rainfall occurrence and no-occurrence is conducted through similar steps. Once classifications have been made, the quantification of possible forecast indicators (i.e., correct rate and false alarm rate) can be calculated based on equation (5) and (6). Afterward, the threshold selection is determined based on equation (7). In this study, unclassified rainfall events from moderate to very heavy rainfalls are included. The result of the correct rate and false alarm rate can be seen in figure 9 and figure 10. From the figures, we can see the different response of correct rate and false alarm rates. The correct rate and false alarm rate difference at BPTK station is more obvious than at Deles station. Moreover, with a prevailing false alarm rate occurs frequently, it may cause bias in determining intersection area between the correct rate and the false alarm rate.
In figure 9(a), the maximum KSS value obtained is 0.39, which corresponds to the value of PWVvariation 0.5 mm. The corresponding correct rate and the false alarm rate is 75% and 36%, respectively. In figure 9(b), the maximum KSS value is 0.27, which corresponds to the value of PWV-rate of change 0.1 mm/h. The corresponding correct rate and the false alarm rate is 81% and 54%, respectively. The average increment time of PWV from an initial condition to the maximum is about 5 hours. In figure 10(a), the maximum KSS value obtained is 0.12, which corresponds to the value of PWV-variation 1 mm. The corresponding correct rate and the false alarm rate is 47% and 34%, respectively. In figure 10(b), the maximum KSS value is 0.36, which corresponds to the value of PWV-rate of change 0.3 mm/h. The corresponding correct rate and the false alarm rate is 62% and 26%, respectively. The average increment time of PWV from an initial condition to the maximum is about 3 hours.
( a )  Figure 9. The relationship between the correct rate, the false alarm rate, and the KSS value at BPTK station in case of (a) PWV-variation, and (b) PWV-rate of change.

KSS
(b) Figure 10. The relationship between the correct rate, the false alarm rate, and the KSS value at Deles station in case of (a) PWV-variation, and (b) PWV-rate of change.  Based on the results, the threshold is determined for the following conditions. In case of the mountain area (Deles station), when the PWV increments have reached to 3 hours with the PWV-variation has reached to 1 mm, it is counted as rainfall will occur. However, if the PWV increments have not reached to 3 hours even the PWV-variation has reached to 1 mm or larger, it is counted as rainfall will not occur. In case of the flatter area (BPTK station), when the PWV increments have reached to 5 hours with the PWV-variation has reached to 0.5 mm, it is counted as rainfall will occur. However, if the PWV increments have not reached to 5 hours even the PWV-variation has reached to 0.5 mm or larger, it is counted as rainfall will not occur.

Conclusion
The parameters of PWV-variation, PWV-rate of change, and PWV-increment time show high sensitivity to rainfalls. Therefore, the threshold of GNSS-derived PWV for short-term rainfall prediction is well defined. The use of correct rate factor, false alarm rate, and Hanssen-Kuipers Score (KSS) as forecast indicators can distinguish the possibility of occurrence and no-occurrence rainfall events and can determine the appropriate thresholds, especially in a mountainous region. The results have revealed that the PWV-rate of change in the mountain area rises three times quicker than in the flatter area. Therefore, the time to forecast is much shorter. However, the result of the correct rate of 47%-62% is considered good. While, for flatter area, the result of the correct rate of 75%-81% remains to compensate for other studies (4,5). The limitation of this study is that the thresholds selected only give us information about the occurrence or no-occurrence rainfall events without knowing which class of rainfalls will occur. It requires more data trained to provide better short-term rainfall prediction based on rainfall class.