GPS Meteorology : Error in the estimation of precipitable water by ground based GPS system in some meso-scale thunderstorms-A case study

Remote sensing by ground based GPS receivers provide continuous and accurate measurement of precipitable water (PW) of an order of 1.5 mm comparable to radiosondes and water vapour radiometers. In the present work we have examined the amount of PW variation in three thunderstorms accompanied with rain shower that occurred over the GPS station. In all the three thunderstorms event heavy rain was reported. However on comparison of observed rainfall with GPS estimated precipitable water (hourly) in real time, it is observed that among the three, in one event the amount of precipitable water (PW) is much less (~20mm) for the same amount of rainfall. After analysing and taken into account various source of error, we suspect that in a mesoscale thunderstorms or squall lines associated with heavy rainfall, discrepancies arise because the wet mapping functions that used to map the wet delay at any angle to the zenith do not represent the localized atmospheric condition particularly for narrow towering thunder clouds and non-availability of GPS satellites in the zenith direction. On the other hand for the larger thunder cells the atmosphere is very nearly azimuthally symmetric with respect to GPS receiver, the error due to the wet mapping function is minimal.


Introduction
Water vapour is an important atmospheric gas. The concentration of water vapour in the atmosphere is highly variable both spatially and temporally. Active weather is strongly correlated to the water vapour distribution in the atmosphere. The conventional method of measurements of water vapour does not normally have a resolution high enough to resolve these variations. Accurate and timely measurement is very important when making weather forecasts and now casting.
The prevention of water disasters is one of the most important issues in urbanised areas. Cumulonimbus when accompanied by a large hailstones, severe squalls or tornadoes are the most intense or meteorological phenomena and yet they are one of the most difficult to forecast. The number of convective strong rainfalls has been increased due to urbanisation and partially due to urban heat island (Kanda et al., 2001). Compared to synoptic scale disturbances, the spatial and temporal behaviour of convective cells in tropical region is more or less chaotic and thus it is still difficult to predict and detect it even by using the most advanced numerical techniques. Although Doppler-radars are powerful to detect raindrops, they cannot monitor moisture and have blind regions near the ground where the initiative convective cells are generated and "cone of silence" in which it does not detect any cloud development.
In NWP modelling, weather is an initial value problem. This means that, to predict changes in the weather it is necessary to have an accurate picture of the atmospheric conditions at the beginning of the forecast period. Moisture is one of the important variables in the NWP model. The serious deficiency in the current observing system is the absence of moisture profiles over the oceans and over land with sufficient horizontal, vertical and temporal resolution. This is especially important and crucial for short-period precipitation and forecasting.
The technology of ground-based GPS and GPS Radio Occultation (GPS RO) are valuable observations for operational meteorology and has a great advantage of monitoring precipitable water vapour (PWV) even before cloud formation free of clutters due to surface obstacles and thus it can be largely expected to detect the initiation of convective cells. The qualities of GPS RO temperature/ humidity observations are generally very accurate with good global coverage; these observations are therefore assimilated operationally in global NWP models. Because of their random distribution over the globe, the footprint of about 300 km and the irregular observation times these data are less suitable for now-casting. For short-range weather forecasting this data could be valuable when there are multiple GPS RO receivers operational. On the other hand GPS Zenith Total delay (ZTD) and PW can be estimated accurately and continuously with a temporal resolution of 10 minutes or less for a surface network of GPS receivers with a horizontal resolution between 50 and 100 km suitable for both nowcasting and NWP applications (Poli et al., 2007).
In recent decades ground based GPS PW measurement attracted the attention of meteorologists and network of continuously operating GPS being operational in many countries for various applications. In our present work we have derived GPS PW in three thunderstorm events that occurred during Indian summer monsoon and studied the utility and various source of errors involving in the measurements.

GPS sites and data set
India Meteorological Department (IMD) established 5 GPS site at New Delhi, Chennai, Mumbai, Guwahati and Kolkata and starts monitoring PW from April 2007 onwards. All GPS sites use Leica AX 1202 antenna connected to Leica GRX 1200 dual frequency receiver via a 30 m low loss cable. The antenna mounts were on concrete pillars fixed on the roof top of the building (approximately 16 to 20 meters height from the ground). A meteorological sensor (Met3A) is also co-located near the antenna mount for temperature, pressure and humidity measurements.
GPS data were processed using GAMIT-10.31 software. The Precise Point positioning (PPP) strategy was used for the data processing. This strategy estimates all station related parameters relative to the precise JPL satellite orbits and satellite clocks that have high accuracy. In order to reduce the multipath effects a 10° elevation cut-off angel was fixed. The New Mapping functions (NMFs) (Niell et al., 1996) were used to map the signal path delay from zenith direction to other elevation angels.
The GPS ZTD estimates are obtained using the GAMIT software at a 30/60 minute interval with a processing window (sliding window) of 9 hours (Foster et al., 2005). Zenith Hydrostatic Delay (ZHD) estimates at each GPS receiver were calculated using surface pressure data measured from the surface Meteorological sensor assuming the hydrostatic equilibrium, is given by (Davis et al., 1985), (1) with ZHD in mm and P s in hPa. The term accounts for the variation in gravitational acceleration with latitude φ and the height H of the surface above the ellipsoid (in kilometres). The ZWD is obtained by subtracting the ZHD from ZTD, the PW (mm) estimates were then derived by scaling the ZWD with the multiplication factor Π given by: where, ρ is the density of liquid water, R v is the specific gas constant for water vapour, k 3 and k 2 ' are physical constant and given by [Bevis et al., 1994]. The water-vapour weighted mean temperature of the atmosphere is defined and approximated as (Davis et al., 1985); can be calculated from vertical profiles of water vapour pressure P v and temperature T, usually derived from NWP analysis fields (Wang et al., 2005). Assuming a linear relation with surface temperature it is also possible to approximate T m from station surface temperature observations T s (Bevis et al., 1992).
The PW thus estimated was validated with nearest Radio sonde data (0000 and 1200 UTC) observed between the year 2007 and 2009 (Puviarasan et al., 2010). The result shows a standard deviation of 6.23 mm and RMSE of 7.09 mm which is mainly attributed to error associated with Radio sonde humidity sensor. It is worth to mention here that after the induction of Graw GPS sonde in IMD, the error associated with humidity sensor reduced considerably and both the observations well agrees with each other. This result is yet to be published.

Event-1
On 1 st August, 2007, during Indian summer monsoon, a thunderstorm with rainfall activity started around 2000 UTC over Delhi region and that continued till 0300 UTC of next day. Total rainfall of 182 mm with maximum rainfall intensity of 60 mm/hr was observed at around 2300 UTC. Since 0000 UTC of 1 st August before the development of thunder clouds, GPS PW steadily increased from 56.82 mm at 0000 UTC to 69 mm at 1500 UTC [ Fig. 1(a)] and decreased thereafter, however the average PW is maintained around 65 mm throughout the day. The PW estimated from RSRW at 0000 UTC and 1200 UTC was 55.14 and 65.26 mm respectively and is well agreed with GPS PW observation.

Event-2
During southwest monsoon, on 27 th July, 2009 over Delhi region a thundershower accompanied with hailstone between 1415-1515 UTC with Rainfall intensity of 80 mm/hr and total rainfall of 137.1 mm was recorded.
[ Fig. 1(b)] shows the time series of PW along with rainfall recorded at the GPS site. In this event also the GPS PW was increased steadily and reached a maximum of 70 mm before the rainfall activity starts. At 0000 UTC of RSRW observation, the PW was around 70 mm and at the same time it was 62 mm in GPS observation. The difference of 8 mm between the two independent observations is attributed to error introduced by humidity sensor that some time become saturated while passing through the clouds. AT 1200 UTC of RSRW observation PW was 69 mm while GPS was 63 mm (difference of 6 mm). Through the RSRW measurement showed a higher value of PW at 0000 UTC, no thunderstorm or rainfall activity near the RSRW or GPS site. This was confirmed from 0000 UTC Kalpana satellite image [ Fig. 1(c)] fin which no significant cloud present over Delhi region. This confirms the validity and accuracy of GPS measurements. At around 0900 UTC the GPS PW increased to 65.5 mm and followed by 35 mm of rainfall it reduced to 58 mm at 1100 UTC. Kalpana image of 1100 UTC [ Fig. 1(d)] shows presence of cloud over Delhi region. AT 1200 UTC the RSRW PW was 69 mm whereas GPS PW was only 63.5 mm and thereafter increases steadily and reached a value of 70 mm at 1500 UTC followed by thunderstorms accompanied by hailstones and heavy rainfall. In Figs. 2(a&b) MODIS image at 1645 UTC showing thick convective cloud spread over Delhi and adjoining region.

Event-3
On 21 st August, 2009 between 1045-1145 UTC a thunderstorm followed by hailstone and intense shower of 78 mm/hr was occurred over Delhi GPS site. This event also took place during the summer monsoon at Delhi GPS site. .5 mm at 1030 UTC just before the occurrence of thunderstorm. This event somewhat differ from the first two cases in which GPS PW increased to 70 mm before storm activity. In all the three events the rainfall intensity is nearly equal; however the measured GPS PW in the third event is relatively much smaller. In fact GPS observation does not show any appreciable increase in PW during or before rainfall activity starts. It is worth to note that PW of 45 mm in break monsoon season is sufficient for thunderstorm with heavy rain if there is any lifting mechanism. Nowadays ground based GPS water vapour monitoring is one of the tool available for nowcasting and input to NWP applications, so it is important to find the reasons for lower value of PW by GPS in some mesoscale thunder storms event

Error due T m
Water vapour weighted vertically averaged mean temperature of the atmosphere T m is an important parameter of the relationship between total precipitable water and the zenith wet delay because, the accuracy of GPS estimates of precipitable water is directly related to the accuracy of T m (Ross et al., 1997). The water-vapour weighted mean temperature of the atmosphere is defined and approximated as (Davis et al., 1985); ∆ T m is function of vertical profile of atmospheric temperature and humidity.
Using equations (4) and (7), the relative error of PW due errors in T m can be given by Wang et al. (2005).
k k is small (~5.9 × 10 -5 ). Hence the relative error of PW approximately equals to that of T m . On the basis of equation (7), for T m values 240 K to 300 K, the (8) 1% and 2% accuracies in PW require errors in T m less than 2.74 K and 5.48 K on average respectively (Wang et al., 2005). T m can be calculated from soundings which are usually unavailable at high temporal and spatial resolution. A more powerful approach would be to use operational meteorological models to predict the actual value of T m . In the event that operational model is not available, T m solely on the basis of observed surface temperature. But this approach suffers from systematic T m overestimation at mid and high latitudes (up to 5 K) and underestimation at low latitudes (up to 6 K) (Wang et al., 2005). Nevertheless, GPS PW is rather robust against T m uncertainties. In Table 1, some of the relations available in the literature to derive mean atmospheric temperature (Jade et al., 2005) were shown. In Table 2, T m derived from respective regression equations from Table1 and from RSRW for the three events is shown. It is understood that the estimated T m for the three events using the regression equations is within 5 K with the observed T m from RSRW except Schueler et al. (2001, Eqn. 11).
This T m uncertainty of 5 K corresponds to 1.7 -2.0% in PW only (Hagemann et al., 2003). Thus for the event 3, the maximum error due to T m would be 2% and may lead to an error in PW of 0.93 mm only. Nevertheless a dip in GPS PW at 1100 UTC [Figs. 1(b&g)] may be attributed to sudden decrease of surface temperature [Figs. 1(e&f)] due to heavy rain that causes significant error in determining mean atmospheric temperature based on T s .

Error in estimation of ZHD
In GPS meteorology it is assumed that the atmosphere is in hydrostatic equilibrium. It does not include the effect due to non-equilibrium conditions. It is in fact, difficult to assess this effect without actually integrating vertical profiles of vertical wind acceleration (Davis et al., 1985). Fleagle and Businger (1980) stated that only under extreme weather conditions (thunderstorm or heavy turbulence) to these vertical accelerations reach 1% of gravity, corresponding to an error in ZHD of about 20 mm/1000 hPa or approximately 3 mm in PW. Fig. 2(e) shows the wind speed and wind direction (time is in IST) along with pressure variation at GPS site. In most circumstances atmospheric pressure is closely approximated by the hydrostatic pressure caused by the weight of air above the measurement point. However during thunderstorms the falling rain creates downdrafts. The downdraft will push down out of the thunderstorm and when reach the surface they spread out and turn into the destructive straight-horizontal winds. The sudden

Model
Mean temperature Bevis et al. (1992Bevis et al. ( , 1994 Eqn. (7 Fig. 1(f)] from 974 hPa to 980 hPa between 1030 to 1130 UTC is not due to atmosphere above the measurement point but due to gusty winds. This dynamical pressure produce a temporary increase of surface pressure thus ZHD, results a dip in PW measurements of the order 0.35 mm/hPa. So in our case the surface pressure was increased 3 hPa on an average between 1030 to 1130 UTC which corresponds to 1 mm dip in PW. Thus the sudden dip in PW at 1100 UTC in Fig. 1(e) is not only by fall in surface temperature but also by dynamical pressure.

Error due to mapping function
When estimating ZTD from GPS signals it is assumed that the atmosphere is symmetric in the horizontal plane and these delays in the signal from GPS satellite to receiver at arbitrary elevation angle θ, can be mapped onto delays in zenith direction through a known mapping function. Thus total delay for a path with an elevation angle θ can be computed from the zenith hydrostatic and zenith wet delays via: where, ) ( M h  is the hydrostatic mapping function and ) M w  is the wet mapping function. Due to different behaviour of dry and wet components of the atmosphere,  the mapping functions are different for hydrostatic and the wet parts. These mapping functions were usually developed based on the assumption that the atmosphere is azimuthally isotropic. Many mapping functions have been developed based on the use of different input parameters, e.g., meteorological data or site coordinates. For example Chao's (1972) mapping functions use one each for h M and w M , make no reference to meteorological conditions, whereas the hydrostatic mapping function derived by Davis et al. (1985) incorporates surface temperature, pressure and relative humidity, the height of the tropopause and the tropospheric temperature lapse rate. Similarly in new Niell's mapping function (NMFs) (Niell, 1996), the delays are modelled separately for the hydrostatic and wet components. Nearly all of the mapping functions that have been suggested in the literature assume no azimuthal variation in path delay (Bevis et al., 1992). The assumption of azimuthal symmetry at 15° elevation typically produce a rms variations of 7 mm, but this effect may go 5 times when the local troposphere has large

Figs. 3(a-c).
Three stages of thunder storms.http://en.wikipedia.org/wiki/Thunderstorm) lateral temperature, pressure or humidity gradients usually accompany strong frontal weather systems (Bevis et al., 1992). Fig. 2(c) shows the 21 st August, 0500 UTC MODIS (Aqua) 11 µm image. Fig. 2(d) shows the variation of brightness temperature (BT) of 11 µm IR channel in MODIS from Aravalli region (Western India) to New Delhi. Along this line it is observed that the BT is decreases and at some point there is a sudden dip in BT of the order of 8° K -A frontal surface separating a shallow, warm moist layer from a deeper, hot and dry layer. Some asymmetry also found in temperature and humidity gradient in east-west direction on 21 st August, 2009. It is worth to note that the presence of dry line is favourable conditions for the formation of thunderstorm.

Error due to satellite positions
In above, we have discussed all possible error sources and if this has been taken into consideration, the observed PW is still much less than the expected one. If one accepts that GPS PW can be one of a tool for a real time forecaster then the questions arises whether the initial 46 mm of PW would give 76 mm of rainfall? If it is so, then most of the monsoon days the PW is above 45 mm and many occasions the PW reaches 60 to 62 mm without any rainfall. In the first two events the PW reached its peak value of 70 and 69 mm before the intense rainfall. However, in third events it is not so. This means that the GPS measurements still does not representing the actual atmosphere completely. If so under what situations?
Now consider the situation as in Figs. 3(a-c). Assume that the GPS antenna is at the center of towering cumulus cloud [ Fig. 3(a)], so that the atmosphere is symmetric with respect of GPS receiver on the ground. Now for a 6 km towering cloud, to capture it by GPS receiver, there should be at least one GPS satellite in the Zenith Direction. Or if the towering cumulus extends ~ 3 to 4 km radius in the horizontal and 6 km in the vertical then, any GPS satellites between 90 to 60 degree elevations can represent the cloud fully. If the cloud extends up to 24 km (12 km radius) in the horizontal direction then satellites up to 30 degree elevation can track the clouds. Below 30 degree elevations satellite signals pass through only part of the thunderstorm and do not depict the atmosphere in the Zenith properly. In the same way during mature stage [ Fig. 3(b)], the thunder cloud reach up to height of 18 km, to perceive this by GPS satellites at 30 degree elevation, the cloud should extend up to 44 to 45 km or radius of 22 km in the horizontal direction. For 60 degree elevation, the thunder cloud should extend at least 30 km (15 km radius) in the horizontal direction. However in tropics quite often small scale thunder clouds or squalls do not have large horizontal coverage. During dissipating stage also the base of the cloud [ Fig. 3(c)] is at a height of 6 km, most of the low elevation GPS signals reach the antenna without seeing the cloud. Thus non availability of GPS satellite in near zenith direction always underestimates PW over GPS station. Figs. 4(a&b) Fig. 4(b)]. One may also notice in the sky plots [ (Fig. 6(a)], below 60 degree elevation towards north of the plot center between 315-040 degree azimuth, no satellites covering this region, thus GPS signals does not detect presence of clods in these regions. This is due to inclination of GPS orbital plane with respect to equatorial plane of earth. In radar meteorology there is a "cone of silence" directly above the radar, in which radar does not detect any object that comes under cone region. Similarly GPS measurements do have blind region that does not detect presence of any thunder clouds. In first two events this zone does not have any effect in PW measurements because the clouds were extended both horizontally and vertically over large region, thus all GPS signals pass through clouds before reaching the antenna. On the other hand for event 3, the thunder clouds develops and enter the GPS site from north north-westerly direction and thus most of the satellite signals reach the GPS antenna without seeing the cloud [ Fig. 4(b)]. This is evident from Kalpana satellite image (Fig. 5)  the uncovered zone is relatively less in Chennai which is at lower latitude (~13°) than at Delhi which is at higher latitude (~28°). This error however can be minimised by incorporating other constellation such as GLONASS, Galileo and COMPASS. Fig. 6(c) shows that, after inclusion of GLONASS with GPS observation, the blind region is reduced considerably. This is because the GLONASS orbital plane is inclined 65 degree whereas GPS orbit is inclined at 55 degree with respect to equatorial plane of the earth.

Conclusions
Knowledge about the amount of atmospheric water vapour is important in many areas like nowcasting application, hydrological cycle, NWP and climate modelling etc. Several techniques have been developed and being in use to obtain precipitable water in the atmosphere. However each technology has some advantages and disadvantages and some limitations. Remote sensing of amount water vapour based on ground based GPS systems is an independent method, works under all weather conditions, gives continuous measurement of water vapour at every 10 minutes or less, more accurate, stable and cost effective and requires less maintenance.
Though GPS can be useful to monitor PW changes, but at times does not capture meso scale thunder clouds or narrow vertical cloud. The wet mapping functions that used to map the wet delay at any angle to the zenith do not represent the localized atmospheric condition particularly for narrow towering thunder clouds. Presence of blind region and non availability of GPS satellite in the vicinity of zenith direction are a major source of error. For the first two events this has no effect, since the thunder clouds have large vertical and horizontal coverage so almost all GPS signals seen the cloud before reaching the GPS antenna. However for the third event GPS does not able to pick up the cloud and measured PW (46.5 mm) was much lower than the expected PW (70 mm). This is mainly due to the fact that very narrow thunder clouds were developed to the north of GPS stations where GPS observation has some blind region below 60 degree elevation and approached the GPS site from north/north westerly direction and large asymmetry components. The blind region is higher at higher latitude and low at lower latitude due to inclination of GPS orbit with respect to equatorial plane. However the blind region can be reduced by incorporating other constellation such as GLONASS, Galileo and compass. Thus inclusion of other constellation with existing GPS can give better estimation of PW. The dip in PW at 1100 UTC for event 2 is attributed to sudden decrease of surface temperature while for event 3 is due to both surface temperature and dynamical pressure. Although the GPS not the 45 mm of PW is sufficient for thunderstorm with heavy rain particularly break monsoon condition Instead of estimating PW in the Zenith direction, estimation of Slant Wet Delay (SWD) along the signal direction and computing PW in each quadrant or octant (by dividing the atmosphere) would give better understanding of PW distribution and useful for nowcasting. Also fast developing thunderclouds PW estimation should be at every 10 to 15 minutes instead of 30/60 minute since GPS signals are less sensitive to water droplets.