GNSS Meteorology for Disastrous Rainfalls in 2017–2019 Summer in SW Japan: A New Approach Utilizing Atmospheric Delay Gradients

We studied disastrous heavy rainfall episodes in 2017–2019 summer in SW Japan, especially in the Kyushu region using tropospheric delay data from the Japanese dense global navigation satellite system (GNSS) network GEONET (GNSS Earth Observation Network). This region often suffers from extremely heavy rains associated with stationary fronts during summer. In this study, we first analyze behaviors of water vapor on July 6, 2018, using tropospheric parameters obtained from the database at the University of Nevada, Reno. The data set includes tropospheric delay gradient vectors (G), as well as zenith tropospheric delays (ZTD), estimated every 5 min. At first, we interpolated G to obtain those at grid points and calculated their convergence, similar to the quantity proposed by Shoji (2013) as water vapor concentration (WVC) index. We obtained zenith wet delay (ZWD) from ZTD by removing zenith hydrostatic delay. The raw ZWD values do not really reflect the wetness of the atmosphere above the GNSS station because they largely depend on the station altitudes. To study the dynamics of water vapor before heavy rains, we estimated ZWD converted to the values at sea level. In the inversion scheme, we used G at all GEONET stations and ZWD data at low-altitude (<100 m) GEONET stations as the input. Then we assumed that spatial change of ZWD is proportional to G (e.g., Gx = H ∂ZWD/∂x, where H is the water vapor scale height) and estimated sea-level ZWD at grid points all over Japan. At last, we tried to justify our working hypothesis that heavy rain occurs when both WVC and sea-level ZWD are high by analyzing hourly water vapor distributions in all the days in July 2017, July 2018, and August 2019. We found that both two values showed maxima in the studied period when the three heavy rainfall (>50 mm/h) episodes occurred, that is, July 5, 2017, July 6, 2018, and August 27, 2019. Next, we performed high time resolution analysis (every 5 min) on the days of heavy rain. The results suggest that both WVC and sea-level ZWD go up prior to the onset of the rain, and ZWD decreases rapidly once the heavy rain started. It is a future issue, however, how far these two quantities contribute to forecast heavy rains.


INTRODUCTION
Disastrous heavy rains in summer 2017 to 2019 in SW Japan caused a lot of damage to property and human lives. The Japan Meteorological Agency (JMA) officially named the extreme rainfall event in 2018 July as "The Heavy Rain Event of July 2018." Precipitation records at meteorological stations show extreme rainfalls from June 28 to July 8, 2018, especially in the northern part of the Kyushu District. Figure 1 shows the precipitation rate in the high-resolution nowcast rainfall map at 17:00 JST (08:00 UT), July 6, 2018, from JMA. This was obtained by a weather radar with a 250-m resolution every 5 min. Such meteorological radars have been operated by JMA at 20 stations throughout Japan. The heavy rain occurred over a patchy region elongated in NE-SW and overlap with the stationary front. Water vapor transported along the front from SW is thought to have caused the heavy rain.
The concept of ground-based global navigation satellite system (GNSS) meteorology was proposed initially by Bevis et al. (1992), and meteorological utilization of the Japanese GEONET has been sought (e.g., Tsuda et al., 1998). Nowadays, GNSS meteorology has become one of the essential means to observe precipitable water vapor (PWV), and PWV data from GEONET have been assimilated in the mesoscale model of JMA to improve weather forecast accuracy since 2009 (e.g., Shoji, 2015). In this study, we apply a new method of GNSS meteorology to utilize atmospheric delay gradients, reflecting azimuthal asymmetry of water vapor (MacMillan, 1995), for the 2017-2019 heavy rain cases in SW Japan. Miyazaki et al. (2003) focused on such atmospheric delay gradients and showed that the temporal and spatial variations of the gradients were compatible with the humidity fields derived from zenith wet delay (ZWD) and with the meteorological conditions in 1996 summer over the Japanese Islands (especially during a front passages). Shoji (2013) and Brenot et al. (2013) demonstrated the important role of the atmospheric delay gradients to detect smaller scale structures of the troposphere than ZWD.
Recently, Zus et al. (2019) have successfully processed the Central Europe GNSS network data to show that the interpolation of ZWD observed with a sparse network can be improved by utilizing tropospheric delay gradients. They showed significant accuracy improvement for the simulation of the numerical weather model, and for the agreement of the simulation results with real observations, relative to the cases without utilizing tropospheric delay gradients.
In this study, we propose a new method to use tropospheric delay gradients to study heavy rain phenomena in Japan using the data from the GEONET data. At first, we analyze behaviors of water vapor on July 6, 2018, using tropospheric parameters obtained from the database at the University of Nevada, Reno (UNR; Blewitt et al., 2018). The data set includes tropospheric delay gradient vectors (G), as well as zenith tropospheric delays (ZTDs), estimated every 5 min. We interpolate G to obtain those at grid points and calculated their convergence, similar to the index proposed by Shoji (2013) as water vapor concentration (WVC). Then, taking advantage of the dense GNSS network, we try to reconstruct ZWD in inland regions (especially in highaltitude regions) by spatially integrating G.
The purpose of this study is to show the implication of utilizing tropospheric delay gradients, in addition to ZWD, to improve our understanding of water vapor dynamics during heavy rains. For this purpose, we reconstruct ZWD and calculate WVC and compare them with in situ rainfall data from AMEDAS (Automated Meteorological Data Acquisition System) stations of JMA. Finally, we study the behaviors of these quantities common for the three heavy rain cases in SW Japan in summer 2017-2019.

DATA AND METHODS: CASE STUDY FOR THE 2018 HEAVY RAIN
We use the data from the dense GNSS network GEONET for the entire country operated by the Geospatial Information Authority (Tsuji and Hatanaka, 2018). It consists of more than 1,300 continuously observing stations in Japan with a typical separation distance of 15 to 30 km. Because its official solution (F3 solution) does not include tropospheric parameters in high temporal resolution (Nakagawa et al., 2009), we used tropospheric delay data from the UNR database (Blewitt et al., 2018). They estimated tropospheric parameters using the GIPSY/OASIS-II version 6.1.1 software with the Precise Point Positioning technique (Zumberge et al., 1997) using the products for satellite orbits and clocks from Jet Propulsion Laboratory. They employ the elevation cutoff angle of 7 • and estimate ZTD and the atmospheric delay gradients every 5 min (Vaclavovic and Dousa, 2015). They follow the 2010 IERS convention (Petit and Luzum, 2010), and used the Global Mapping Function (Böhm et al., 2006), and troposphere gradient model of Chen and Herring (1997).

ZWD and Tropospheric Delay Gradients
The equation below is formulated by MacMillan (1995) and explains that the slant path delay (SPD) at the elevation angle θ and the azimuth angle ϕ measured clockwise from north can be expressed as follows.
SPD(θ, φ) = m(θ) · [ZTD + cotθ(G n cosφ + G e sinφ)] + ε (1) There, m(θ) is the isotropic mapping function that shows the ratio of SPD to ZTD, and G n and G e are the north and east components, respectively, of the tropospheric delay gradient vectors G, and ε is the modeling error. ZTD is the refractivity of the atmosphere integrated in the vertical direction and is the sum of zenith hydrostatic delay (ZHD) and ZWD. In this study, we calculated surface pressure at the GNSS stations assuming 1 atm at the sea level. Then, we calculated ZHD and subtracted it from ZTD to isolate ZWD. Considering average variability of surface pressure, errors by this approximation for summer ZWD remain within a few percent.
Tropospheric delay gradients are also the sum of hydrostatic and wet contributions. Because we analyze summer data, we assumed that the latter is dominant and used the gradient as those representing the water vapor. Because of the low scale height of water vapor (∼2.5 km), ZWD highly depends on the FIGURE 1 | High-resolution map of the rainfall rate at 17:00 JST (08:00 UT) July 6, 2018, from JMA (https://www.data.jma.go.jp/fcd/yoho/meshjirei/jirei01/high resorad/index.html).

FIGURE 2 | (A)
Atmospheric delay gradients at grid points at 17 JST (08 UT) on July 6, 2018 (raw gradient data are shown in Figure 4A). Using these gradient vectors at grid points, we calculated their convergence shown in (B), similar to WVC index defined by Shoji (2013).
Frontiers in Earth Science | www.frontiersin.org station altitude; that is, small ZWD observed at highland does not always imply low humidity of the troposphere above that station. On the other hand, atmospheric delay gradients are observed by directly comparing atmospheric delays in different azimuths and hardly suffer from the station altitude (Shoji, 2013). This is the reason why we use troposphere gradients to reconstruct ZWD converted to sea level. Shoji (2013) suggested that two new quantities, WVC and WVI (water vapor inhomogeneity) indices, provide valuable information on non-uniform distribution of atmospheric water vapor in meso-β and meso-γ scales in addition to PWV representing water vapor distribution in meso-α scale. WVI is derived from the post-fit phase residuals in the processing of geodetic GNSS data analysis, and it is impossible to derive such information in the UNR database. WVC indices can be derived as the spatial derivative of ZWD. However, as explained in the previous section, low scale height of water vapor makes ZWD highly dependent on station height, and it is difficult to calculate WVC from ZWD in the mountainous Japanese Islands owing to large topographic slopes throughout the region. Shoji (2013) suggested that WVC can be calculated directly using the atmospheric delay gradient vector field, which is immune from the height problem and is readily available in the UNR database.

Water Vapor Concentration
Following Shoji (2013), we use WVC to discuss its relationship with heavy rainfalls. WVC index expresses the degree of divergence/convergence of the atmospheric delay gradient and represent the short-wavelength concentration of PWV.
∇PWV is the spatial gradient of PWV. We here use the observed tropospheric gradient G, divided by the scale height H of water vapor.
is the coefficient to convert ZWD to PWV and is a function of temperature. ∇ZWD is the spatial gradient of ZWD and can be expressed as G/H (Ruffini et al., 1999). We here assumed 2.5 km for H. In our study, first, we obtained the vector G at grid points with east-west separation of 0.20 • and north-south separation of 0.15 • all over the country. Here we calculated G at grid points as the weighted mean of G at all the GEONET stations. Larger weight was given to nearer stations. The weight was taken proportional to the Gaussian function of the interstation distance with one-sigma of 20 km.
Unlike the original definition of WVC by Shoji (2013), we did not convert the wet delay to PWV by multiplying with in (3). So, our WVC index is actually −∇ 2 ZWD and is calculated numerically as the convergence of G, that is, using G at grid values. Figure 2 shows the tropospheric delay gradient vectors at the grid points and its convergence (WVC) calculated by Eq. (4).
As the first case study, we investigate water vapor in the heavy rainfall episode at 08 UT (17 JST) July 6, 2018. There are strong southward gradients in southern Hokkaido and central Tohoku, suggesting southwestward increase of water vapor. In northern Kyushu, we can see large gradient vectors line up, suggesting high concentration of water vapor there.

Sea-Level ZWD
In our study, in addition to WVC, we also reconstruct sealevel ZWD using the observed tropospheric delay gradients and ZWD at low elevation (<100 m) stations (they are assumed to be identical to sea-level ZWD there). Although the concept of WVC is adopted from Shoji (2013), it is our original technique to estimate the sea-level ZWD using the gradient vector data as the input. Below we provide the detailed explanation of this new inversion technique.
In the inversion, we estimate x(i, j; i = 1,2,. . ., m e , j = 1,2. . .,m n ), the sea-level ZWD at the grid point with coordinates (i, j; i = 1,2. . . m e , j 1,2,. . ., m n ). The input data are G (G n , G e ) at all the GEONET stations. Suppose the situation in Figure 3, that is, the k'th GNSS station is located within a rectangle made of four corners (i, j-1; i + 1, j-1; i, j; and i + 1, j). Then, the observation equation to relate G at this GNSS station G(k) to the sea-level ZWD at grid points x(i, j) is,   To regularize the inversion, we constrain the sea-level ZWD at the grids closest to the GNSS stations with altitude less than 100 m to the observed ZWD. They indicate that the sea-level ZWD at the grid point x(i, j) is the same as the ZWD y(k) observed at the k'th GNSS station, that is, We then performed inversion also applying a continuity constraint. When we estimate the sea-level ZWD over the entire Japanese Islands, the number of parameters is ∼1,600, and the number of observations is twice as large as the number of GNSS stations, ∼2,600. The number of GNSS stations with altitudes of less than 100 m is ∼100. In Figures 4, 5, we show inputs and outputs of the inversion, respectively. In Figure 5B, we confirm that the atmospheric delay gradient vectors at GNSS stations calculated using the estimated distribution of sea-level ZWD. They are similar to the input shown in Figure 4A, suggesting that the inversion is successful. The root-mean-square error between the observed and calculated atmospheric delay gradient is ∼0.65 mm in this case.   Figure 5A shows high ZWD throughout SW Japan, but such a ZWD map still lacks spatial resolution to pinpoint heavy rain as indicated in Figure 1. At a glance, WVC in Figure 2B shows good coincidence with the heavy rain given in red and orange colors in Figure 1. However, there are regions where WVC is high but heavy rain does not occur. This suggests that both two quantities need to be high for the occurrence of heavy rains. In this section, we study long-and short-term behaviors of the two quantities in several recent heavy rain episodes in SW Japan. For the former, we see hourly changes of WVC and sea-level ZWD over 1-month periods including the heavy rain episodes. For the latter, we study the change of these quantities every 5 min over the days of heavy rain.

The 2017 and 2019 Heavy Rain Cases
In Figures 2B, 5B, we show the WVC index and sea-level ZWD at 08.00 UT (17 JST) when heavy rain occurred. This can be compared with the rain images from the JMA, as shown in Figure 1. By comparing Figures 1, 2B, we can see that the WVC index successfully pinpoints the heavy rain. This is also consistent with the detailed report of this heavy rain episode compiled by Japan Meteorology Agency [JMA] (2018).
Here we perform the same calculation for a heavy rain episode on July 5, 2017, at 07 UT (16 JST) and show the results in Figure 6. These figures also show, like in the 2018 case, that both the sealevel ZWD and WVC are high where it rains heavily as shown in Japan Meteorology Agency [JMA] (2017).
Likewise, Figure 7 shows the sea-level ZWD and WVC for a heavy rainfall episode in August 2019. WVC index pinpoints heavy rainfall, consistent with the information compiled by Japan Meteorological Agency [JMA] (2019).
From these results, we hypothesize that heavy rainfall occurs when both the ZWD and WVC are high. Next, we try to test the hypothesis by studying time series of the two quantities and rain rate.

WVC, Sea-Level ZWD, and Heavy Rain
Here we try to justify our working hypothesis that heavy rains occur when the WVC index and sea-level ZWD are both high by analyzing the distribution of water vapor every hour in the three cases, July 2017, July 2018, and August 2019. We present the results in Figure 8, which explains the scatter plot of sea-level ZWD and WVC, along with hourly rainfall data based on observations at the AMEDAS station of JMA.
We selected the AMEDAS station showing the largest rain for the three episodes, that is, the Asakura station for 2017, the Yanase station for 2018, and the Shiroishi station for 2019. These Frontiers in Earth Science | www.frontiersin.org three stations are all located in northern Kyushu and recorded heavy rains exceeding 50 mm/h. Then we picked up the grid point closest to these AMEDAS and compare the three quantities over a month, that is, July 2017, July 2018, and August 2019. We can see that ZWD values are high up to 400 mm, and WVC index goes up to 40 × 10 3 mm/km 2 when heavy rain exceeding 50 mm/h occurred on July 5, 2017, July 6, 2018, and August 27, 2019 (Figure 8).
Next, we look for data from more AMEDAS stations in Japan on the same days showing hourly rain rates exceeding 50 mm/h and found 12 stations for the 2018 July and 5 stations for the 2019 August, the months including heavy rain episodes. The 2017 heavy rain was quite local (limited to northern Kyushu); we could not find enough number of AMEDAS stations showing rains exceeding 50 mm/h. From Figure 9A, we find that the probability of the heavy rain (>50 mm/h) was 14% for the range of ZWD (400-450 mm) and WVC index (35-50 × 10 3 mm/km 2 ), for July 2018. Likewise, for August 2019, the heavy rain occurred 50% for the same range of ZWD and WVC index. If we count rains exceeding 20 mm/h, then the percentages go up to 71% for July 2018 and 78% for August 2019. These results indicate that both the WVC index and sea-level ZWD are high when heavy rains occur. The results also suggest that heavy rains may not occur even when these two quantities are high. Next, we will show time series of these quantities.

Time-Series Analysis
In Figure 10, we plot the quantities in Figure 8 as the functions of time. The two quantities WVC and sea-level ZWD are given at the top, whereas the rain rate is given at the bottom. The bottom panels also include the time series labeled as "deviation." This quantity indicates the sum of the deviations of WVC and sea-level ZWD from their averages normalized by their standard deviations. For example, if both quantities deviate from the means by 2σ, the deviation is 4. Figure 10 clearly shows that the two quantities show large deviations from the average values whenever heavy rains occur. We also see sometimes that high deviation does not coincide with a heavy rain, for example, July 20, 2017, and July 23, 2018. Regarding the time sequence, the high WVC/ZWD times seem to "coincide" with heavy rains rather than "precede" them. Hence, the usefulness of monitoring these quantities for weather forecast is not clear from this figure.
Lastly, we analyzed high time resolution (every 5 min) behaviors of WVC, ZWD, and rain rate, for the three heavy rain days in July 5, 2017, July 6, 2018, and August 27, 2019, selecting two AMEDAS stations from each episode (Figure 11). We estimated the sea-level ZWD and WVC every 5 min for these 3 days. The results show that both WVC indices and sea-level ZWD show large values at the start of the heavy rains. However, for many cases (e.g., Figures 11A,E,F), heavy rains already started before the two quantities reach their peaks, suggesting limited applicability of these quantities for weather forecast. ZWD often showed rapid declines after the start of the heavy rainfalls possibly because of the transformation of water vapor to liquid water.

CONCLUSION
The atmospheric delay gradients are little influenced by station altitudes and provide useful information of spatial gradient of relative humidity in lower atmosphere. By using atmospheric delay gradients, we could estimate distributions of sea-level ZWD, reflecting relative humidity of the air column above the station. This is possible only in regions with dense GNSS networks like Japan, because we relate tropospheric delay gradient to horizontal spatial gradient of ZWD. We also showed the importance of the WVC index reflecting small-scale enhancement of water vapor. As seen in the scatter diagrams in Figures 8, 9, sea-level ZWD and WVC are not correlated; that is, high sea-level ZWD does not always mean high WVC and vice versa. We showed that heavy rain occurs only when both quantities show high values.
The results suggest that monitoring these quantities is useful for the nowcast of heavy rains. However, their potential to forecast heavy rains is yet to be studied. As the next step, we will need to explore the way to use these two quantities for weather forecast, for example, by putting them to numerical weather models by some means. Computation time of sea-level ZWD and WVC is short, and we can convert ZTD values estimated by GNSS data analysis in near real time to sea-level ZWD and WVC.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author. Tropospheric parameters used in this study are available from geodesy.unr.edu.