Applications of GNSS Slant Path Delay Data on Meteorology at Storm Scales Applications of GNSS Slant Path Delay Data on Meteorology at Storm Scales

This chapter focuses on applications of Global Navigation Satellite Systems (GNSS) slant path delay data (SPD) to obtain signals from thunderstorms or rainbands. Current operational numerical weather prediction systems (NWPs) use water vapor distributions derived by GNSS technology as vital information for predicting convective rainfall. Mostly, zenith total delay or integrated water vapor data are used at horizontal scales of several tens of kilo- meters for this purpose. Beyond such operational use, SPD can be used to obtain information on storms (cumulonimbus) at horizontal scales of less than 10 km. For instance, found that SPD represents very small-scale phenomena of less than 10 km and can be used to estimate water vapor distribution around a thunderstorm with a strong tornado, and succeeded in improving the forecast skill of a rainband at 10 km scale. This chapter reviews SPD, which is invaluable for predicting thunderstorms and/or rainbands.


Introduction
The Earth's atmosphere, including the ionosphere, affects navigation signals transmitted by Global Navigation Satellite Systems (GNSS), which causes positioning errors. It is able to eliminate the ionospheric effect using a pair of GNSS carrier waves, and then GNSS analysis estimates the signal delay by the atmosphere as an unknown parameter [1][2][3][4]. The atmospheric delay is obtained by integrating the refractivity along the ray path, which is calculated with the variables of temperature, pressure, and water vapor pressure in the atmosphere. In the 1980s, several studies showed the feasibility of sensing the atmosphere using GNSS signal delay. For instance, Askne and Nordius [5] found a relationship between signal delay in the zenith direction (zenith total delay; ZTD) and the precipitable water vapor (PWV; vertically integrated water vapor) in a formulation. Studies such as these led to the establishment of a new interdisciplinary research field (Bevis et al. [6], Businger et al. [7]). The Global Positioning System (GPS) maintained by the United States is the world's first operational GNSS, and has contributed significantly to the progress of the meteorological applications of GNSS. One of the most important contributions of GPS/GNSS meteorology is an innovative application of GNSS observations as a water vapor sensor.
Water vapor in the atmosphere plays an important role in convective rainfall systems. Particularly in the summer in Japan, convective clouds form frequently during the evening and develop into thunderstorms. Weather radar has been used in several observational studies to investigate the evolution of such convection. However, since weather radar only detects the distribution and movement of raindrops, it is still difficult to predict the precise position and time of the initiation of heavy precipitation. GNSS tropospheric delays are sensitive to water vapor, and can be used with radar to monitor the early stages of the convergence of water vapor before severe precipitation. Since the mid-1990s, procedures for retrieving precipitable water vapor (PWV) using the GNSS have advanced significantly.
In Japan, the Geospatial Information Authority of Japan (GSI) operates a nationwide permanent ground-based GNSS observation network called the Earth Observation Network (GEONET), which covers the entire Japanese archipelago with an average spacing of 17 km. It is regarded as one of the densest GNSS networks in the world. On July 13, 2012, the GSI began to provide data from the Russian GLONASS (Globalnaya Navigatsionnaya Sputnikovaya Sistema) and Japan's QZSS (Quasi Zenith Satellite System) along with GPS data. Also, from March 2016, the GSI started providing GALILEO data.
Several studies have been shown GNSS PWV as a useful tool for monitoring heavy rainfall. Kanda et al. [8] found that periods of maximum PWV tended to be a precursor of the initiation of heavy rain, and that larger hourly increases in PWV led to higher frequency of precipitation.
Niimura et al. [9] statistically investigated the relationship between precipitation, PWV, and temperature, and found the correlation between the observed PWV and surface temperature. They also found a threshold of PWV depending on the temperature; for instance, when PWV exceeds the threshold, frequency of precipitation rapidly increases. After their statistical study, it was concluded that this threshold is connected with the humidity averaged vertically above the observation point, and that the saturation level of the entire atmosphere can be estimated with PWV. Inoue and Inoue [10] also statistically studied the two-dimensional PWV field derived by GNSS in connection with thunderstorms in summer. They found that cloud-to-ground (CG) strokes were observed 1 hour after large amounts of PWV and its 30-min increments.
These studies have illustrated that observing PWV and its variation are useful for monitoring heavy rainfall. However, monitoring PWV alone is not always enough to capture information on the precursor of a severe storm. For instance, Niimura et al. [9] illustrated a good agreement between the level of saturation of the entire atmosphere given by PWV and stratiform rain, but the relationship was unclear in localized heavy rain cases. Inoue and Inoue [10] statistically showed that the maximum CG stroke only in 40% of the observed thunderstorms was observed 15-30 min after the maximum PWV. Therefore, it is necessary to find other method for the remaining 60%.
Since the development of the GPS/GNSS meteorology, the applications of slant path delays (SPDs; delays along paths from GNSS satellites to GNSS receivers) have been studied for water vapor variation in local and its role in the development of hazardous thunderstorms. One such approach is the tomographic method, which estimates structures of water vapor in 3D. For this purpose, several GPS/GNSS observational campaigns with dense networks have been conducted, and obtained successful results [11][12][13]. However, it is difficult for operational GNSS networks to retrieve 3D structures of water vapor using only GNSS SPDs, due to the low density of the network, small number of SPDs, and slow movement of GNSS satellites.
Another application is assimilation of GNSS-derived water vapor information for improving the initial fields of numerical weather prediction (NWP) models. The impact of this information on mesoscale NWP systems has been thoroughly investigated (Nakamura et al. [14], Koizumi and Sato [15], Seko et al. [16]). These studies mainly used PWV data retrieved by GNSS analysis using IGS's precise final ephemerides. Shoji [17] developed a near real time (NRT) PWV analysis system. This system provides the PWV data by the time when the Japan Meteorological Agency (JMA) assimilates the data in the operational mesoscale data assimilation (DA) system. Using PWV by the NRT system, Shoji et al. [18] succeeded to predict heavy rainfalls, which JMA's mesoscale NWP model failed operationally. Further, an impact assessment of the local-scale atmospheric phenomena on GNSS positioning was also carried out using mesoscale NWP outputs (Seko et al. [13], Shimada et al. [19], Ishikawa et al. [20]).
Other studies have investigated the assimilation of SPD data. The advantage of SPD for PWV and ZTD is that it includes information about several atmospheric parameters (pressure, temperature, and humidity) in several directions from each receiver. Therefore, the assimilation of SPD data is expected to improve the water vapor field with the thermodynamics field of the model above and around the observation points.
Ha et al. [21] assimilated virtual GNSS slant water vapor observations, which are accumulated water vapor along slant paths (SW data), as an Observation System Simulation Experiment (OSSE). The SW data is advantageous to PWV observations in its various directions from each receiver. They succeeded to simulate a realistic squall line after assimilation of the SW data with a fifth-generation mesoscale model (MM5)-4D-Var system (Zou and Kuo, [22]) at a 27-km grid spacing. Järvinen et al. [23] assimilated SPD data with a High Resolution Local Area Modeling-3D-Var system (HIRLAM; Gustafsson et al. [24]) at a 9-km grid spacing, and showed that the analysis increments of the SPD assimilation were larger that of the ZTD assimilation, and that the horizontal distribution of the SPD analysis increments was different from that of ZTD assimilations. Bauer et al. [25] tried to assimilate SPD data with the MM5-4D-Var system and succeeded to improve quantitative precipitation forecasting skills. Though they demonstrated the superior QPF scores by assimilating SPD data over a month, their horizontal grid spacing was 18 km; this is not enough to represent thunderstorms directly in the model. In the past few years, the German Weather Service (DWD) has developed an assimilation method for SPD data with their operational data assimilation system (Kilometer-Scale Ensemble Kalman Filter; KENDA) and obtained a 20% improvement on prediction of a convective rainfall system (Potthast 2017; personal communication).
The SPD data represent both vertical and horizontal atmospheric conditions, whereas the ZTD and PWV observations contain only vertical information. Thus, it is advantageous to assimilate the SPD observations at storm scales. For instance, a GNSS signal with a 30° elevation angle at a GNSS receiver travels from the top of the troposphere to the receiver with a horizontal distance of 17 km. Thus, the SPD data cover only two model grid cells, when an assimilation system with 20-km grid spacing is applied. As a result, it is expected that the assimilation effect of SPD data would be similar to the ZTD assimilation. Therefore, it is important for SPD assimilations to use assimilation systems with high grid spacings, hopefully less than 5 km.
In this chapter, the use of SPD data for the monitoring of hazardous convection is described in Section 2, and the first assimilation of SPD data with a DA system with 2-km grid spacing (storm scale) is given in Section 3. The chapter is summarized in Section 4.

Using SPD to monitor hazardous convection
GNSS-derived PWV is, so to speak, a representative value of PWV within an inverted coneshaped space above each GNSS antenna. GNSS-derived PWV is inherently unsuitable for monitoring water vapor variation under severe storms where several-kilometer-scale phenomena prevail. Is SPD suitable in this case? In this section, we first divide SPD into three components based on an SPD model in GNSS analysis, and discuss the horizontal scale of each component. Next, we propose new indices that represent finer water vapor distribution than the GNSS PWV. Finally, we introduce a procedure to analyze a several-kilometer-scale PWV distribution around each GNSS antenna using SPD.
Following MacMillan [26], the SPD between a GNSS satellite and a receiver at an elevation angle and a direction angle measured clockwise from north can be written in the following form: where ZTD, m ( θ ) , G n (G e ) and ε are the total delay in the zenith direction, the isotropic mapping function that describes the ratio of SPD to ZTD, the delay gradient parameters in the north (east) directions, and the postfit phase residual, respectively. ZTD is vertically integrated refractivity (N) of the atmosphere in the zenith direction. Refractivity is expressed by temperature (T), partial dry air pressure (P d ), and partial water vapor pressure (P w ): where n is refractive index, and K 1 , K 2 , and K 3 are constants that have been determined theoretically or by fitting observed atmospheric data. Several studies have evaluated GPS-derived gradient parameters by comparison with those observed by water vapor radiometers (WVR), and found good agreement (e.g., [27][28][29]).
The gradient parameters represent the horizontal first-order gradients of water vapor. The postfit residuals are expected to contain information on higher order water vapor inhomogeneity (HI). However, other errors that do not originate from the atmosphere are also included (e.g., antenna phase center variation (PCV), signal scattering, multipath, errors in satellite orbit, and errors in clocks of both satellite and receiver). Therefore, to estimate HI, it is necessary to remove all errors not related to atmospheric inhomogeneity. Shoji et al. [1] performed a procedure to eliminate multi-path and satellite clock error-induced residuals to reconstruct HI components. The correlation coefficient of each component retrieved at a different GNSS station and sorted by distance demonstrated that the horizontal scale of the ZTD can be considered as 644 ± 120 km, the gradient parameter G n , G e as 62 ± 23 km, and the HI as 2-3 km. This result suggests that ZTD, G, and HI relate to atmospheric motion of the mesoα , mesoβ , and mesoγ scales, respectively (Figure 1).
Sato et al. [30] compared zenith-scaled SPD using a mapping function and that retrieved from radiosonde observations. They found that the zenith-scaled SPD, in which the path is closest to a radiosonde path, exhibited better agreement than zenith total delay (ZTD) retrieved by standard GNSS analysis (i.e., a representative value of the inverted-cone-shaped space above the GNSS antenna).
Although it requires some careful effort to retrieve, GNSS SPD possesses information on local-scale atmospheric activity. Shoji [31] proposed procedures for retrieving two indices indicating the degree of inhomogeneity of water vapor using the GNSS SPDs. One index (WVC) describes the spatial concentration of water vapor (Eq. 3), whereas the other (WVI) indicates higher order water vapor inhomogeneity (Eq. 4). The horizontal scales of the two indices are considered to be approximately 60 km and 2-3 km, respectively.
where ∇ PWV G is the horizontal gradient of PWV estimated from the atmospheric gradient parameter (G).  where HI PWV is the inhomogeneity component of SWV normalized as a vertical value using a mapping function. To extract HI PWV from the postfit phase residual ( ε ), it is essential to eliminate the effects of GNSS antenna phase center variation (PCV), multipath effects, and errors in satellite orbits and clocks.
The statistical examination between these indices and precipitation (Figure 2) illustrate that the inhomogeneity indices show stronger correlation with rainfall amount than with PWV. It seems that PWV related to precipitation of less than 10 mm h −1 , but was not connected to precipitation greater than 10 mm h −1 . It is also true for both present and imminent precipitation. Furthermore, the spatiotemporal variations of the indices were also examined on a particular thunderstorm on August 11, 2011. Both the WVC and WVI indices increased ahead of the convective initiation (Figure 3). The WVI index is based on the standard deviation of the SPDs measured at a GNSS station, so the directional information for each SPD is neglected.
Shoji et al. [2] proposed a new method for estimating PWV distribution around each groundbased GNSS station on a scale of several kilometers (Figure 4).    In this approach, the following three assumptions were set.

1.
The horizontal gradient of dry refractivity (N dry ) is small enough to be negligible.

2.
The difference between ZTD EST (retrieved value through standard GNSS analysis) and ZTD SPD (normalized SPD into the zenith direction using a mapping function) is caused by the several-kilometer-scale horizontal gradient of water vapor refractivity (N wet ) alone.

3.
The horizontal N wet gradient (g wet ) decreases exponentially with height.
where g wet (z) is the horizontal water-vapor gradient at altitude Z, and H is the scale height of g wet .
Under these assumptions, the horizontal PWV gradient can be expressed as the following equation.
where Π is the proportionality coefficient. The horizontal gradient of PWV is expressed as a function of the ZTD difference, elevation angle, and scale height.
For practical monitoring of severe convection, we must carefully assess the accuracy of PWV distribution using the new method. Shoji et al. [3] quantitatively evaluated the method of Shoji et al. [2] using the simulation results of a high-resolution NWP model performed by Mashiko [32] for a tornadic supercell case, which generated an F3 tornado.  Figure 5, which illustrate distances at where the RMS difference was the smallest, show the distance increases as the elevation angle decreases. The RMS was the smallest within 1 km range from a virtual GNSS station for PWV SPD with a 77.6° elevation angle. In case of PWV SPD with a 17.6° elevation angle, RMS was the smallest of 1.5 mm at a 4.5 km distance. The distance and elevation dependency of minimum RMSE are illustrated by the red dots ( Figure 5). Overall, the conventional procedure causes about RMS of 0.5 mm at the GNSS site, and the error by simple extrapolation increases with distance, reaching 1.5 mm at 1 km. The distance dependency of PWV errors in PWV SPD differs in each elevation angle. From this result, we can estimate PWV with RMSE of less than 2 mm within 6 km from a GNSS station by using SPD with an elevation angle over 15°. Essentially, with PWV SPD , it is able to estimate the PWV distribution around each GNSS station with better than half the RMSE of that obtained by the conventional GNSS PWV retrieval method (PWV EST ).  The new method demonstrated the capability to capture a strong PWV gradient associated with the parent storm of the F3 tornado that struck Tsukuba City in Ibaraki Prefecture, Japan, on May 6, 2012, with a numerical model simulation at a super-high resolution of 50 m (NHM50m; Figure 6). An area of large PWV contrast centered on strong precipitation (Figure 6a) implies a strong upward wind in front of and a strong downdraft behind the parent storm. The PWV contrast toward the tornado is also seen in Figure 6b, whereas, no such PWV gradient is seen at all in Figure 6c.

Data assimilation system
Data assimilation (DA) techniques optimize differences between first guess fields provided by numerical models and observations. One advanced DA theory is a variational method (Sasaki [33,34]), and the NHM-4DVAR is a variational DA system for predicting thunderstorms (Kawabata et al. [35,36]). The NHM-4DVAR is based on the Japan Meteorological Agency Nonhydrostatic Model (JMANHM; Saito et al. [37]), which includes three-ice bulk cloud microphysics, as the forward model. In the first version of NHM-4DVAR, the adjoint model considered only dry dynamics and advection of water vapor. Kawabata et al. [36] implemented an additional warm rain process for assimilating radar reflectivity. The horizontal resolution of NHM-4DVAR is 2 km. This 2-km spacing is known as "storm scale" because in this range it is possible to explicitly represent cumulonimbus and thunderstorms by numerical models. The control variables are the three wind components, potential temperature, surface pressure, nonhydrostatic pressure, total water (water vapor and cloud water), the relative mixing ratio of rainwater, and the pseudo-relative humidity (only for lateral boundary conditions).
This method optimizes the difference between observations and numerical models using their error statistics. To measure the difference, we need observation operators to convert the numerical model fields to observation spaces. For SPD DA, Eq. (2) is used as the operator and provides delays of GNSS radio waves from GNSS antennas to the model top boundary as Here, ∆L denotes the atmospheric delay (m), n denotes the refractive index along the path in a grid cell, and ds is the path length (m) in each model grid cell. Because the top boundary of numerical models is limited to a certain height level (usually 20-40 km), we need an assumption to add a delay from the top boundary to a GNSS satellite.
The authors assumed that the delay decreased to 1/e every 10 km from the top boundary to 200 km height (Bean and Thayer [38]) and that above 200 km the amount of delay was zero. In addition to this assumption, a straight-line assumption was adopted; the bending effect was thus eliminated. To account for these assumptions, a relatively large observational error was implemented. To calculate SPD in the model, first, the linear path is determined from the receiver to the GNSS satellite. Then, the middle point of the path within the model grid cells are set down. Delays are interpolated and averaged from eight points surrounding the middle point, with each weight set according to the distance. Finally, the slant path delays are calculated by integrating each delay in model grid cells from the receiver to the top boundary of the model.
The observations over elevation angles of 5° or more were assimilated in their study. The operator works on the World Geodetic System 1984 (WGS84; National Imagery and Mapping Agency 1997) [39].

Single set of SPD observations at a single site
First, three experiments were performed using NHM-4DVAR with a 2-km horizontal grid spacing in which three SPD observations (SO_SPD), one ZTD observation (SO_ZTD), or one PWV observation (SO_PWV) from a single observation site were assimilated. Note that the ZTD and PWV observations were derived from the SPD observations originally. These experiments were conducted to confirm the effects of SPD assimilation on a single analysis and to examine the differences between SPD, ZTD, and PWV assimilations. The assimilation window was set at 10 min, and the observations were assimilated every 5 min (at 0, 5, and 10 min in the assimilation window). The observational data set was chosen by considering the horizontal distributions, elevation angles, and the first guess field from an experimental data set.
The propagation paths of GNSS signals from three satellites to a receiver in the model atmosphere in both the horizontal (Figure 7a) and the vertical plane (projected from the south; Figure 7b) are illustrated. Path I with the smallest elevation angle is also the longest, while path II at a near 90° elevation angle is the shortest. The large amount of delays is illustrated mostly in low altitudes area by warm colors.
The distributions of the increments (analysis minus first guess) of PWV in SO_SPD (Figure 8a) are different from that in SO_ZTD (Figure 8b) at the end of the assimilation window. In SO_ZTD, the increment distribution (Figure 8b) is axisymmetric and elliptical mostly (i.e., isotropic). Although 4D-Var captures analysis increments shaped flow-dependent (anisotropic) in general, it is not  long enough for the 10-min assimilation window to capture such flow dependency. As a result, the shape of the analysis increment seen in SO_ZTD is isotropic mostly. In contrast, the inhomogeneous distribution of the analysis increment in SO_SPD (Figure 8a) is given by the various distribution of the slant paths. Moreover, the maximum value of the increments is much larger in SO_SPD (10 mm; Figure 8a) than that in SO_ZTD (6 mm, Figure 8b). The increment distributions in SO_PWV and SO_ZTD were close to each other (not shown).
Seeing vertical cross sections of the mixing ratio of water vapor (Qv) along path III (Figure 9), the distributions of analysis increment in SO_ZTD and SO_PWV are close, and the magnitudes of the increment are quite similar: The distributions of Qv reach at high altitudes from 1 to 5 km in vertical and extend to 5-8 km away in horizontal. The increment along path III in SO_SPD is distributed along a distance of 15 km in horizontal and over 8 km in height; furthermore, the magnitudes of the increment are much larger than others, at low altitude especially (around 3 km), because all of the slant paths are within a narrow area of the lower troposphere above the observation site. The total weights in the cost function show good agreement among the assimilations of SO_SPD, SO_ZTD, and SO_PWV, but their effects are seen at different places.

Actual observational data assimilation
Next, an actual observation assimilation and forecast experiments were conducted using NHM-4DVAR with the assimilation window of 30 min. Three-hour forecasts were performed after the 4D-Var analysis with multiple observations. The assimilation was started at 1100 JST on August 19, and the forecast was performed from 1100 to 1400 JST. The SPD, ZTD, or PWV data were assimilated every 10 min (hereafter, A_SPD, A_ZTD, and A_PWV, respectively).
Only the observations listed above, and no other kinds, were assimilated. No cumulus parameterization was applied in the experiments. The experimental domain and the distribution of GNSS receivers are illustrated in Figure 10. Unlike the PWV and ZTD data, it can be said that SPDs in the model atmosphere (Figure 10a) contain a great deal of horizontal information.  The horizontal distributions of rainfall in A_ZTD and A_PWV (Figure 11c and d) were close to each other. In A_SPD (Figure 11b), the initiation of intense rainfall enhanced over southern Okinawa at 1200 JST, and the maximum intensity reached 47 mm h −1 (not shown). By 1300 JST, a rainband had formed over the island, and the distribution and intensity resembled the radar observation (Figure 11a). Therefore, it can be concluded that the assimilation of SPD data showed improvement of the rainfall forecast with respect to both timing and intensity compared with the assimilation of PWV and ZTD data. This improvement was obtained through modifications in the atmospheric profile in the simulation after the SPD assimilation (not shown).

Summary
GNSS analysis is performed with a number of models for solid earth tide, ocean tide, earth rotation, and atmospheric delay. However, conventional GNSS analysis inevitably has errors caused by local-scale variations in the atmosphere, even if estimates of the gradient parameters are applied. These errors affect positioning as well as atmospheric delays. Spatial correlation in screened post-fit phase residuals (i.e., the higher order inhomogeneity part) must be caused by local-scale variations in the atmosphere, which are not resolved by ordinary ZTD and gradient parameter analysis. Thus, GNSS postfit residuals can be used to detect the inhomogeneous distribution of water vapor on a local scale of less than 10 km caused by weather disturbances like convective thunderstorms. In this chapter, we introduced two new indices and one experimental method to retrieve a several-kilometer-scale PWV distribution. These approaches may be a preliminary way to serve operational monitoring of cumulus convection.
Currently, multiple GNSSs (e.g., the Russian GLONASS, EU's GALILEO, China's BeiDou, and Japan's Quasi-Zenith Satellite System (QZSS)) can serve to more precisely resolve localscale water vapor variation. In addition, the advancement of real-time GNSS analysis has been progressing rapidly. Real-time orbit and clock corrections have offered officially by the International GNSS Service (IGS) since April 2013. Moreover, a Multi-GNSS orbit and clock estimator called MADOCA (Multi-GNSS Advanced Demonstration tool for Orbit and Clock Analysis) has been developed by the Japan Aerospace Exploration Agency (JAXA) and they started to provide real-time ephemerides via the Internet (https://ssl.tksc.jaxa.jp/madoca/public/public_index_en.html) in September 2015. These rapid advancements in GNSS technology should speed the reality of new usage of GNSS-derived SPD.
Another application of SPD data at a local scale is data assimilation (DA). As shown in this chapter, a promising avenue for DA applications is to improve the initial conditions of a numerical weather model in high-resolution simulations. NHM-4DVAR with a 2-km horizontal grid spacing was used for assimilating GNSS slant path delay data. Compared with simulations after assimilating PWV and ZTD data, the assimilation of SPD showed inhomogeneous increments. Moreover, the analysis increments in the assimilation of SPD at low altitudes were larger than in the others. Conducting actual observation assimilations, the SPD assimilation clearly improved the forecast of both the timing and intensity of the rainband. This improved forecast was given through the decreased atmospheric stability over Okinawa Island.
The SPD data include information on temperature, dry atmospheric pressure, and water vapor pressure. In addition, the data contain both horizontal and vertical information on those atmospheric parameters. Because atmospheric inhomogeneity is greatest in the lower troposphere, assimilating SPD data is a promising way to improve forecast of a rainband. A GNSS signal along a path with a 15° elevation angle propagates about 11 km in horizontal and travels 3 km in vertical. Therefore, it is concluded that data assimilation systems at storm scales like used in this study are necessary to assimilate SPD data for taking advantages.