Atmospheric Water Vapour Determination Using GPS Signals for Numeric Weather Prediction in Tanzania

Numeric weather prediction (NWP) is modern and best technique for accurate weather forecasting which impacts a lot of human activities and survival. NWP models use different input parameters including atmospheric water vapour in which its measurement is still faced with numerous challenges locally and at a global scale [1]. Atmospheric water vapour is an important greenhouse gas that plays a vital role in maintaining Earth’s energy balance and is also a vital meteorological parameter in NWP [2]. The amount and distribution of water vapour in space (horizontal and vertical) is a major parameter in the development of NWP models and its importance cannot be underestimated. The measurement of such parameters have performed different with different techniques including convection methods like radiosondes (RAOBs), microwave radiometers (MWRs) and sun photometers. However, these existing systems are inadequate in global scale and hence there is clear need of developing a superior system that is high accurate with high temporal and spatial resolution. That system should be free from meteorological conditions, low maintenance and most of all relatively cheaper [1]. GNSS meteorology technique in retrieving atmospheric water vapour was first suggested by [3] and later clarified more by Solheim and Ware, The concept is based on tracking propagation delays due to neutral atmosphere as the magnitude of signal delay is treated to be proportional to the amount of water vapour in atmosphere. Given the straightforward GNSS ZTD computations as well as the low initial cost, this technique suppresses most of challenges faced by convection techniques. Tanzania Meteorological Agency (TMA) has 26 synoptic stations with only four Radiosondes across the country. Although some regions have varying and challenging weather, these are the only system used to measure atmospheric water vapour in all Tanzanian regions. This deficit has caused the TMA to apply most products from global numerical weather models for their predictions. Some of the global models are computed by different programs around the world example the European Centre for Medium-Range Weather Forecasts (ECMWF) and the Global Continuous Observing System of the World Weather Watch (GCOS) Upper-Air Network –GUAN (Figure 1). These global predictions are very precise in places where the station are available. Given a sparse GUAN instrumentation in Tanzania, the TMA to lagging behind in providing accurate predictions. Since Ground-based GNSS CORS stations form part of the GUAN, this paper utilizes GNSS meteorology concepts to investigate the GNSS CORS at Dodoma. This CORS was used in determining atmospheric water vapour (precipitable water) which is the key requirement in weather prediction. The reason for choosing GNSS CORS at Dodoma is its proximity to one of the GUAN Station.  The correlation results between the GUAN prediction and the computed CORS precipitable water vapor will be integrated to other CORS sites. This integration will thereby improve infrastructures for NWP in Tanzania and help TMA in weather forecasting. 

Figure 1: GUAN Ground based GPS stations as of February 2021 as recorded by ECMWF. Stations that contain both daytime and nighttime data are denoted by green crosses, whereas daytime-only stations are denoted by red diamonds, and nighttime-only stations are denoted by blue triangles.

The error due to troposphere on GNSS signals is non dispersive and can be split into two components, hydrostatic and wet delays. Although satellites are located at different angles from the horizon, these delays are mapped onto zenith direction to get zenith hydrostatic delays (ZHD) and zenith wet delays (ZWD) by using mapping functions. The sum of the mapped delays is called Zenith Total Delay (ZTD). These tropospheric effects can introduce errors from 2.5 m to as high as 25 m depending on the elevation of the satellite [4]. The ZTD and ZHD are determined directly from GNSS processing using GAMIT/GLOBK and gLAB software. The GAMIT/GLOBK software uses in-built algorithms to compute the zenith delay parameters at time interval set in the processing using a piecewise-linear (PWL) function. The correction for the hydrostatic delay is applied to the resulting zenith delay estimates and later these residual wet delays are then converted into precipitable water (PW). The resulting PW estimates become the input in computing the required ZWD at the set time intervals [5]. The gLAB on the other hand uses Precise Point Positioning (PPP) technique to determine ZTD using different algorithms that compute the ZHD.  The gLAB first approach uses algorithms and methodologies explained example by [3]. In this method, the ZHD is computed using the observed surface pressure above the certain height (H).

 The second approach computes ZHD from non-weather observation instead it uses station-specific parameters [6]. These parameters may include the 3D coordinates, zenith tropospheric delays, phase bias and clock parameters. In this approach the exponential nature of ZHD is described by its elevation in the relation elaborated in equation 3

 Where a = 2.3m, b = 0.116, and H is station height above sea level.

The combination of the three equations (1, 2, 3) are implemented in the gLAB software to compute the ZHD. The ZWD is computed from the directed determined ZTD and the computed ZHD using equation 1. In the case of GAMIT/GLOBK software the estimated atmospheric values (ZWD and PW) for each file at given epoch are computed as a result of processing using a piecewise-linear (PWL) function following equations defined in [3] 

The parameters used in this computations are the atmospheric water vapour in terms of Integrated water vapour (IWV) in kg/m^2, which is the quantity of the atmospheric water vapour over a specific location and PW used to express the height of an equivalent column of liquid water in units of length. Where p_v is the partial density of water vapour in kg/m³; the height h in meters and R_w is the specific gas constant for water vapour in J/(kgK). PW relates to IWV by dividing with the density of liquid water, ρw. PW= IWV/ρw Again IWV is related to the ZWD using a dimensionless quantity as conversion factor, Π: The conversion factor is given by Bevis et al, (1994) as

With  T_m = 0.72 T_s + 70.2 hereby T_s is absolute surface temperature, p_w is density of water, R_w is the specific gas constant for water vapour in J/(kgK), k'₂  and k₃ are constants based on laboratories estimates defined by [3]. The computed  ranged from 0.160225 to 0.163567 and the computed mean is 0.161802. The dimensionless constant is a function of season and location and has been computed in previous studies and found to be be approximately 0.15 [3]. In this study the dimensionless constant falls within the 20% recommended range. This show a good agreement between the computed value from this study to the previous studies as it was also identified in the study by [7,4].

The main input in this study are one year long daily observation and navigation files for Dodoma CORS (DODM) site downloaded from UNAVCO achieve. The files are stored in hatanaka compressed and required to converted into RINEX format in order to be read by gLAB software.  Other data include the meteorological files for one year long of DODOM, IGS final orbit, absolute antenna phase center corrections for both receivers and satellites and the total amount of water vapour present in a vertical atmospheric column (IWV) downloaded from the ECMWF.The IWV was used for validation and statistical analysis of PW from the DODM CORS.  Meteorological files were the inputs in computation of ZWD, PW and IWV and the IGS final orbit and absolute antenna corrections were used for ZTD. The Meteorological data were spliced into daily file (24 hrs) in order to correspond to observational and navigation files. In the case of IWV from ECMWF, the netCDF formats downloaded convert mean daily value at each 0.125⁰ grid which were later interpolated to get most probable value of IWV at experiment site DODM. The one-year data used started from day 168 of 2017 to 168 2018 with exception of day 291-304 of 2017 and 111-164 of 2018, in which DODM did not have data. 

• Data Processing 

In the determination of PW two software with complete different processing strategies are implemented using GAMIT/GLOBK and gLAB software. The gLAB is based on PPP in which Precise satellite orbit and clock corrections produced by the IGS are used in computing solutions while GAMIT/GLOBK uses double differencing in which clock error and integer ambiguity are eliminated when computing the final solution. The following sections explain the processing strategies for each software. 

• gLAB Software 

The input data were decimated at every 30 sec with elevation mask of 5 degrees and cycle slip detection at every 40 sec data gap. These data were modeled using Klobuchar model for ionospheric correction and corrections for relativistic, antenna phase centres and ocean tidal loading effects [8]. For the case of tropospheric correction which is of more interest, this study used UNIB-3 Nominal model which computes dry (T_dry) and wet (T_wet) delays from the receiver’s height and estimates five meteorological parameters: pressure, temperature, water vapour pressure, temperature lapse rate and water vapour lapse rate according to procedures described in RTCA-MOPS. This procedure used the more refined Niell mapping model described in [9] which considers different obliquity factors for the wet and dry components and do not equire any surface meteorological measurements.  The results is the determined ZTD at each 30sec which was later averaged into mean daily ZTD using their respective standard deviation. PW was computed using the resulting ZTD following the methodologies described in Equations 1 to 5. The results are PW_gLABwd for weather dependent and PW_gLABwf for weather free. 

• GAMIT/GLOBK Software 

This study analyzes the the yearlong GPS data using the GAMIT-GLOBK software [5]. Usually a full process follows the two steps strategies which is GAMIT and later combination through GLOBK. These steps are partly described by e.g. [10-14] However, in this process we used on one step as we only require the estimated atmospheric values. In this process this study uses for each day the doubly differenced GPS phase observations to estimate daily station coordinates, satellite state vectors, 7 tropospheric delay parameters at each station per day, 2 horizontal tropospheric gradients per day, and phase ambiguities applying IGS final orbits and Earth Orientation Parameters (EOP). This study in addition, applied the absolute phase center correction using the IGS tables [15], the current ocean tide model [16], solid Earth and polar tide correction to obtain the loose daily solution vector and its variance covariance matrix for station and orbital elements as quasi-observations. Through this process the estimates are atmospheric values for the site DODM at each hour and later average it for each day. The GAMIT incorporated operations allows for extraction of zenith delay estimates (ZTD), apply for hydrostatic delay (ZHD) and convert the residual wet delays to PW. The source of the applied hydrostatic corrections is the measurement of station pressure recorded in the Receiver Independent Exchange Format - RINEX met file [5] . The result from this process are the mean daily values for ZTD and PWD.

The computed PW values for each day were analyzed on daily basis following the results from gLAB and GAMIT/GLOBK. The analysis based on their correlation and their descriptive statistics for each respective dry and wet season. According to weather records in Dodoma the wet season is considered from November to May and dry from End of May to stand of November each year. First the PWs from GAMIT (PW_GAMIT) was plotted against the PWs from gLAB - PW_gLABwd and PW_gLABwf­ (Figure 2). The average PW difference ranges between 2.5mm to 0.5mm in dry season and between 3.5mm to 0.5mm for wet season. Given the ranges, better agreement is observed during dry season and much weaker agreements in wet season, however time series pattern looks similar in both cases.  The computed difference between the two values is less than 3mm showing that there is a great agreement between the two independent GNSS computed PW values.

Figure 2: Comparison of GNSS computed PWs for site DODM. Yellow is gLAB computed PW weather dependent, gray is GAMIT computed PW and red is gLAB computed PW weather free.

To compare the agreement between GPS derived PWs more precisely, we decided to quantitatively investigate the consistence of the two GPS solution from the original data by determine a measure of the strength and direction of the linear relationship between the two time-series. This analysis is used to quantify the degree by which two variables are related and that the evaluation extends to tell how much one variable cause changes to other variable and provide a coefficient that show the degree of this relationship. We use the linear correlation coefficient computed on the bases of least squares fitting of the two data sets in each component e.g. [17]. We use the so called Persons correlation coefficient applied to a sample of data (not the entire population), since our comparison is based on only one year of GPS observations. The correlation is given by:

Where X, Y and X, Y are GNSS derived PW from gLAB and GAMIT and their corresponding mean respectively for the entire time-series of n data points. For the perfect correlation the value of is is such that r= +1 for perfect positive correlation and r = -1, for perfect negative correlation (Figure 3A,3B,3C).

Figure 3A: The Correlation between gLAB computed PW weather free against the GAMIT computed PW. The computed Pearson correlation coefficient is 0.9867

Figure 3B: The Correlation between gLAB computed PW weather dependent against the gLAB computed PW weather free. The computed Pearson correlation coefficient is 0.9967

Figure 3C: The Correlation between gLAB computed PW weather dependent the GAMIT computed PW. The computed Pearson correlation coefficient is 0.9922 

As it was expected better correlation is observed between two PW solutions from gLAB, however the weather dependent PW solution seems to correlate better with GAMIT computed PW than the weather free solution.  According to the results all solution correlates better

To compare agreement between the two datasets, the GNSS computed PW daily values for Dodoma GPS Site (DODM) time is interpolated from global atmospheric reanalysis center (ERA-Interim).  This study first resample the GPS computed PW from both gLAB and GAMIT estimates with the global interpolated PW values for DODM when both GPS and Global data are available. The extracted values are the compared graphically. The comparison show a similar trend in a sense that both had relatively smaller values during dry season and larger values during wet season. In comparison with PW_ERA-Interim, the _­PW_glaBwd have relatively smaller values in both seasonscwith more vivid departure during wet season (Figure 4). Contrary, weather-free model values (PW_gLaBwf) had larger values than PW_ERA-Interim in both seasons (Figure 5). The values from PW_GAMIT had mainly smaller values than PW_ERA-Interim during dry season contrary to wet season in which PW_GAMIT have oversampled PW_ERA-Interim in majority (Figure 6).  The mean seasonal difference are shown which supports the trend shown on graphs (Figure 7a,7b), whereby during both seasons, the differences were minimum when results from GAMIT (PW_GAMIT) were used.

Figure 4: Comparison of global interpolated PW and gLAB computed PW weather dependent for site DODM .  Blue is gLAB computed PW weather dependent and Green is global interpolated Pw from ERA-Interim.

Figure 5:  Trend and differences between global interpolated PW and gLAB computed PW weather dependent for site DODM with respect to time

Figure 6: Comparison of global interpolated PW and gLAB computed PW weather free for site DODM .  Orange is gLAB computed PW weather dependent and Green is global interpolated Pw from ERA-Interim

Figure 7a:  Trend and differences between global interpolated PW and gLAB computed PW weather free for site DODM with respect to time

Figure 7b:  Trend and differences between global interpolated PW and GAMIT computed PW for site DODM with respect to time

This study further increases the analysis by quantitatively investigate the consistency of the two solution from the original data by determine a measure of the strength and direction of the linear relationship. We use the same analysis of least squares fitting of the two data sets as explain by [17].  The slope = 1 or -1 is an ideal case were the two datasets are perfectly correlated. The correlation was analyzed between global interpolated PW for DODM site against the GNSS computed PW from gLAB and GAMIT software (Figure 8A,8B,8C).

Figure 8A: The Correlation between gLAB computed PW weather dependent with the global interpolated PW for DODM site. The computed Pearson correlation coefficient is 0.9688

Figure 8B: The Correlation between gLAB computed PW weather-free with the global interpolated PW for DODM site. The computed Pearson correlation coefficient is 0.9588

Figure 8C: The Correlation between GAMIT computed PW with the global interpolated PW for DODM site. The computed Pearson correlation coefficient is 0.9672

The relationship is such that, if the GNSS computed PW and global interpolated PW values for DODM are perfectly correlated the values should be aligned with the ideal line. Our result shows a strong positive correlation in all cases, however, a better agreement was observed between gLAB computed weather-dependent PW, followed by the GAMIT computed and finally the gLAB computed weather-free (Figure 8A, 8B, 8C). The better agreement of GAMIT and gLAB weather-dependent PW to the globally interpolated ERA-Interim may be attributed by its dependence to weather signals recorded in the RINEX file as compared to weather-free computations. The obtained correlation coefficient values (r > 0.96 values) for this study are larger than correlation values obtained on similar study for GNSS site in Africa [2,18]. The computed PW estimate for KNUST site located in Ghana compared with the global interpolated PW values from ERA-interim show a correlation r < 0.8 [2]. Similar study was conducted to compare nine (9) International GNSS service station (IGS) with global re-analysis, radiosondes and AERONET data spanning 3 years period [18]. An averaged correlation coefficient values of 0.81 and 0.67 were obtained for ERA-40 and NCEP respectively which is way smaller than the one obtained in this study.

This study has outlined how Precipitable water (PW) can be determined from GNSS observed raw data in either Receiver Dependent or Receiver Independent Exchange Format (RINEX). The study was conducted on a long serving CORS station at Dodoma (DODM) using two software’s. gLAB and GAMIT. Based on a yearlong data from DODM site, different result were obtained using different approach obtained in gLAB as well as GAMIT software. The results from the software used were compared as well as the results from the global interpolated values of PW from ERA-Interim. The analysis of GNSS computed PW values shows a good agreement with PW retrieved from Global reanalysis data of ERA-Interim with strong positive correlation and a small difference in mean daily values for both seasons.  A good agreement observed between GNSS computed PW values shows a correlation <0.98 while between GNSS and global interpolated values of above 0.96 which is close to 1 (Perfect correlation). Although GAMIT results show better correlation to the global interpolated PW values does not mean the gLAB PW values exhibits less comparisons as the differences are very small.  Given these correlations, this study provides great indication of how GNSS data can be used to retrieve the key meteorological values. When computed at different locations, the meteorological data can be used to improve weather forecasting in the sub-region with minimum cost and high spatial and temporal resolution.  Given these results, better weather forecasting across Tanzania can be improved if GNSS data are incorporated in the weather forecasting. Since the GNSS metrology is increasing between geodetic community there is a need for a collaborative effort between Tanzania Surveying and Mapping Division (SMD) and TMA. The first basic study can be done using the available CORS stations across the country and later increase densify GPS base stations so as to fully exploit meteorological advantages from GPS technology in Tanzania.

This work has been carried out at the Department of Geospatial Sciences and Technology (SGST) Ardhi University Dar Es Salaam as part of undergraduate dissertation. We wish to thank Ardhi University for providing facilities to process and analyze GNSS data. We also wish to acknowledge the use of the freely available GNSS data for DODM from UNAVCO as well as the Global reanalysis model of ERA-Interim. Most of the figures were plotted using Microsoft excel and ArcGIS Pro software.

1. Wolfe DE, Gutman SI (2000) Development of the NOAA/ERL Ground-Based GPS Water Vapor Demonstration Network: Design and initial results. Journal of Atmospheric and Oceanic Technology, 17: 426-440.
2. Acheampong AA, Fosu C, Amekudzi LK, Kaas E (2017) Precipitable Water Comparisons Over Ghana using PPP Techniques and Reanalysis Data. South African Journal of Geomatics 6.
3. Bevis M, Businger S, Chiswell S, Herring TA, Anthes RA, et al. (1994) GPS meteorology: mapping wet delays onto precipitable water. J Appl Meteor 33: 379-386.
4. Seeber G (2003) Satellite geodesy, 2nd edition, Walter de Gruyter Berlin New York.
5. Herring TA, King RW, McClusky SC (2015) Documentation for the GAMIT GPS software analysis, release 10.71, Mass. Inst. of Technol. Cambridge.
6. Zumberge JF, Heflin MB, Jefferson DC, Watkins MM, Webb FH (1997) Precise Point Positioning for the efficient and robust analysis of GPS data from large networks, Geophys. Research 102: 5005-5017.
7. Liou K, Meng CI, Newell PT, Lui ATY, Reeves GD, et al. (2001)Particle injections with auroral expansions. Journal of Geophysical Research 106: 5873-5881.
8. Subriana JS, Zornoza JMJ, Hernandez-Pajares M (2013) GNSS Data Processing, Vol 1: Fundamentals and Algorithms. ESA Communications, Noordwijk, Netherlands.
9. Niell AE (1996) Global mapping functions for the atmospheric delay at radio wavelengths. J. Geophys. Research 101: 3227-3246.
10. Dong D, Dickey JO, Chao Y, Cheng M (1997) Geocenter variations caused by atmosphere ocean and surface ground water, Geophys. Res. Letters 24: 1867-1870.
11. McClusky S, Ergintav S, Hamburger M, Kastens KA, Kotzev KA, et al. (2000)GPS constraints on plate motions and deformations in eastern Mediterranean and Caucasus,  geophysical Research 105: 5695-5719.
12. Nocquet JM, Willis P, Garcia S (2006) Plate kinematics of Nubia-Somalia using a combined DORIS and GPS solution, J. Geod. 80: 591-607.
13. Stamps DS, Calais E, Saria E, Hartnady C, Nocquet JM, et al. (2008) A kinematic model for the East African Rift, Geophys. Res. Letters, 35: 05304.
14. Saria E, Calais E, Willis AP, Farah H (2013) A new velocity field for Africa from combined GPS and DORIS space geodetic solutions: Contribution to the definition of the African reference frame (AFREF), J. Geophys. Res. Solid Earth 118: 1677-1697.
15. Saria E, Calais E, Willis AP, Farah H (2013) A new velocity field for Africa from combined GPS and DORIS space geodetic solutions: Contribution to the definition of the African reference frame (AFREF), J. Geophys. Res. Solid Earth 118: 1677-1697.
16. Lyard F, Lefevre F, Letellier T, Francis O (2006) Modelling the global ocean tides: Modern insights from FES2004, Ocean Dyn 56: 394-415.
17. Cohen J, Cohen P, West SG, Aiken LS (2002) Applied multiple regression / correlation analysis for the behavioral sciences (3rd ed.). Psychology Press.
18. Bock O, Bouin MN, Doerflinger E, Collard P, Masson F, et al. (2008) The West African Monsoon observed with ground-based GPS receivers during AMMA, J. Geophys. Res 113.