Utilization of GPS Tropospheric Delays for Climate Research

The tropospheric delay is one of the main error sources in Global Positioning Systems (GPS) and its impact plays a crucial role in near real-time weather forecasting. Accessibility and accurate estimation of this parameter are essential for weather and climate research. Advances in GPS application has allowed the measurements of zenith tropospheric delay (ZTD) in all weather conditions and on a global scale with fine temporal and spatial resolution. In addition to the rapid advancement of GPS technology and informatics and the development of research in the field of Earth and Planetary Sciences, the GPS data has been available free of charge. Now only required sophisticated processing techniques but user friendly. On the other hand, the ZTD parameter obtained from the models or measurements needs to be converted into precipitable water vapor (PWV) to make it more useful as a component of weather forecasting and analysis atmospheric hazards such as tropical storms, flash floods, landslide, pollution, and earthquake as well as for climate change studies. This paper addresses the determination of ZTD as a signal error or delay source during the propagation from the satellite to a receiver on the ground and is a key driving force behind the atmospheric events. Some results in terms of ZTD and PWV will be highlighted in this paper.


Introduction
One of the most atmospheric parameters that play a key role in the Earth's atmosphere and a highly variable atmospheric constituent in determining the accuracy of weather forecasts as well as the positioning accuracy in altitude determination is the total tropospheric delays. Tropospheric delay value represents the total amount of signal delay that would result from dry gases (~90%) presents at the troposphere (N 2 , O 2 , Ar, etc.) and the water vapor condensed in form of clouds (~10%) at a particular time and over a given location from the surface of Earth to the top of the atmosphere. On the other hand, the tropospheric delay amount depends on the signal path through the neutral atmosphere. The current state analysis revealed that the tropospheric delay is retrieved from the Global Navigation Satellite System (GNSS) receivers, which known as the total troposphere zenith path delay (ZPD), or often referred to as zenith tropospheric delay (ZTD) [1]. In this method, GNSS (e.g., Global Positioning System (GPS) for simplicity) satellite sends electromagnetic signals through the atmosphere to a receiver on the ground at a fixed location. Total delay in the GPS signals is measured, and the ZTD is obtained from a summation of dry and wet components [2].
Since the ZTD data is strongly significant in meteorological studies such for improving numerical weather prediction (NWP) models, to benefit, their value needs to be converted into the amount of water vapor in the atmosphere, which so-called the precipitable water vapor (PWV) [3]. PWV is one indicator of the critical components of our atmosphere, where their distribution and content are critical variables that can explain the evolution of various physical processes in the atmosphere. With the availability of PWV data at fine temporal and spatial resolution, unpredictable phenomena such as extreme weather and 2

1234567890
The climate change at a specific location can be mitigated its impacts before they happen, and this, in turn, will maintain the resilience of the system either in local, regional, and global scales. Hence, this paper aimed at determination of ZTD for climate studies. Conversion of ZTD into PWV will be discussed. Quantification of both parameters and interpretation of their physical characters are crucial in understanding on how the global climate system such as global warming changing over time. Based on that basis, some results from Antarctica research will be presented.

Determination of Tropospheric Delays
When the radio signals traveling through the atmosphere to a receiver, it is retarded and affected significantly by the presence of free electrons in the ionosphere and the refractivity of gases, hydrometeors, and other particulates in the neutral atmosphere. During the propagation, the time at which the signal received is given by clock in a receiver, and the difference between transmits time and receiver time is so-called pseudorange. Figure 1 shows the GPS signals propagates through the atmosphere. The difference between the actual path of the carrier (S) and the straight-line (rectilinear) path in a free-space (G) will obtaine dt, where dS and dG are an infinitesimal part of the path length, v and c are speed of the radio signals in medium and in free-space respectively. From the Fermat's principle, the transmission will follow the path between two points (P 1 and P 2 are the initial and end point of the S, respectively). Let D be the delay of the radio wave, with θ is the elevation angle of the direction from which the wave arrives, S be the true path along which the radio wave propagates, and G the shortest geometric path along which the signal would transverse, assuming as n = 1. Then In the first term Eq. (2), the integral is performed along the line increment dS of straight ray path (excess path delay) from a receiver to a GPS satellite. The second term indicates the geometry delay due to ray path bending. As the electromagnetic wave propagates in the atmosphere, it is continuously refracted due to the varying index of refraction of the air between the top of the atmosphere and the ground. Basically, it has two effects on a ray path: bending and retarding, which both produce an excess path length with respect to propagation in a free-space. Excess path length due to signal retarding in the troposphere is expressed as [5] [ ] where ∆D is the difference delay between S and G. However in the neutral atmosphere several simplifying assumptions can be made to simplify Eq. (3). First, the atmosphere is assumed to be spherically symmetric, where the Earth is a sphere and the properties of the atmospheric vary only with geometric radius. In this way the atmosphere can be considered layered with a refractive index characterizing each layer. Second, the atmosphere is usually assumed to be azimuthally symmetric; that is with no variation in azimuth in each layer. In this way the electromagnetic ray is confined to a plane defined by the start and end points of the ray and the geocentric. These assumptions allow us to represent the refractive index profile as a function of geocentric radial distance only, n(h). With the application of the refractivity expression (3), Eq. (2) for S can then written as where the refractive index is integrated along the path between points h 0 and h 1 which are the geocentric distance of the user's antenna and the geocentric distance of the 'top' of neutral atmosphere, respectively.
Angle β z is the actual (refracted) zenith angle of the ray path at distance h (see Fig. 1). The path delay is caused by variation of n from unity, hence: This gives in the first term the excess delay equivalent path length and in the second term the geometric length of the curved path. To obtain the total tropospheric delay, we can subtract the delay in free-space distance to get the following integral equation [6]: To summarize, the first integral accounts for the difference between the electromagnetic and geometric distances along the ray path and the bracketed integrals account for the curvature of the ray path, i.e. the difference between the refracted and rectilinear geometric distances. Angle α z is the true (unrefracted) satellite zenith angle and hence constant along the unrefracted path. The total tropospheric delay, can be simplified as where δ is the tropospheric correction (recently, known as gradient tropospheric delay, if δ = 0, the asymmetric components are neglected). Propagation delays at arbitrary elevation angles are determined from the zenith delays and so-called "mapping functions". Tropospheric delays increase with decreasing satellite elevation angle. This is accounted for multiplying the zenith delays by a correction factor, m. In general, the total tropospheric delay from (8) following Davis et al. [5] can be rewritten as where ZTD is the zenith tropospheric delay or equal to the total tropospheric delay (  Figure 2 shows the flow diagram for the determination of PWV from GPS signals and surface meteorological measurements [7]. The PWV is extracted from ZWD. ZTD is a main key in determining the PWV, which is composed of ZHD and ZWD.

Determination of Precipitable Water Vapor
where for the hydrostatic component subscript j is replaced by h and for the wet component subscript j is replaced by w respectively. Trop j N 0 , is the total refractivity at the surface of the Earth and is given as: where P is the measured total atmospheric pressure in mbar. The refraction constants where T and T K are the surface air temperatures and both are in Celsius and Kelvin respectively. The P w (in mbar) is partial pressures of water vapor and is obtained from relative humidity (H) with using the formula as recommended by WMO Technical Note No. 8.

1234567890
The The ZHD is calculated using Saastamoinen (SAAS) model [10]. This model is most precise that uses the surface pressure measurements and a correction factor to correct the local gravitational acceleration at the center mass of the atmospheric column. The ZHD from SAAS can be expressed as follows and is close to unity, where λ is the station latitude (in degree) and h is the height of the site above the ellipsoid (in km) measured by GPS. In this work, the hydrostatic Niell mapping function, m h (θ) is employed to reduce the dependence of the zenith delay (ZTD) to the satellite elevation angle. This mapping function has been normalized to yield a value of unity at the zenith [11]: where h is the height of the site above geoid in meters. Here, ( ) represents the three-term continued fraction in the Marini mapping function, In the eq. (16), the coefficients a ht = 2.53 x 10 -5 , b ht = 5.49 x 10 -3 and c ht = 1.14 x 10 -3 were determined by a least-squares to fit the height corrections at nine elevation angles. In modern mapping function, Boehm et al. [12] has updated the "b" and "c" coefficients of the Marini continued fraction form, as shown in (16). However, the hydrostatic "a" coefficients are still valid for zero heights. This mapping function is later called the hydrostatic Vienna Mapping Function (VMF1). To increase the precision, it is using the European Centre for Medium-Range Weather Forecasts (ECMWF) instead of the Numerical Weather Models (NWMs). More than that, the coefficient c shown in (15) is now modeled to remove the systematic errors as (17) where DoY is the day of the year and 28 January has been adopted as the reference epoch [11], ϕ is latitude (deg), and ψ specifies the northern (0) and southern hemisphere (π). On the other hand, the ψ value was selected π/2, c 0 = 0.062, c 10 = 0.0015, and c 11 = 0.006. Detail parameters for c 0 , c 10 , c 11 and ψ needed for computing the coefficient c of the hydrostatic mapping function can be found in [12]).
Since the ZWD cannot be modeled by using meteorological parameters, it is calculated from subtracting ZHD from ZTD. PWV now can be calculated as given by Bevis et al. [8] ( ) where lw is the liquid water density in 1000 kg m -3 , v is the universal dry gas constant (287.054 J Kg -1 K -1 ), 2 ′ (22.1 K mbar -1 ) and 3 (3.701 × 10 5 K 2 mbar -1 )) are the refraction constants. refraction constants and temperature T m are estimated by Bevis et al. [8]. T m for this work is calculated from the surface temperature (T) measured at the particular site. Figure 3 shows the example of ZTD results for Davis station (DAV1: 68.58°S, 77.95°E, and ellipsoid of 44.4 m) in Antarctica at one hour resolution. The results are compared to the Crustal Dynamics Data Information System (CDDIS) NASA. The ZTD in this work was computed using a Tropospheric Water Vapor Program (TroWav) developed using Matlab source codes [4,7,13]. The program is capable of computing ZHD, ZWD, ZTD, and PWV. As shown in the figure, a small drift between TroWav ZTD and NASA ZTD exhibited for two peaks with one valley, where TroWav ZTD reached its highest value on the DoY = 363. This drift for three days is about 1.02 mm higher from their mean difference. The GPS ZTD at DAV1 exhibit almost similar fluctuation compared to the ZTD from NASA and shown a strong correlation (R 2 = 0.70 at the 99% confidence level). The mean difference values between TroWav ZTD and NASA ZTD is 2.04 mm. On the other hand, the ZTD average from TroWav at DAV1 for the period of 29-31 December 2015 is 2.246 ± 0.009 m while ZTD obtained from NASA is 2.262 ± 0.011 m. Another work to compare the ZTD results between TroWav and GeoForschungsZentrum (GFZ) and Center for Orbit Determination in Europe (CODE) products has reported by Suparta and Kemal [1]. One can be found from this observation is that the variation of both ZTDs in this region is strongly following the pattern of surface pressure.

PWV Results
After the ZTD computed by TroWav program and the accuracy obtained are reasonable, then the PWV value can be retrieved. Figure 4 shows the PWV variability in monthly average for four selected stations in Antarctica. In this paper the data is collected and processed for Davis