Determining precipitable water in the atmosphere of Iran based on GPS zenith tropospheric delays

Precipitable water (PW) is considered as one of the most important weather parameters in meteorology. Moreover, moisture affects the propagation of the Global Positioning System’s (GPS) signals. Using four different models, the current paper tries to identify the best relationship between the atmospheric error known as zenith wet delay (ZWD) and PW. For that matter, based on 54,330 radiosonde profiles from 11 stations, two different models i.e. linear and quadratic have been derived for Iran. For analyzing the accuracy of these models, ZWDs of three permanent GPS stations located in the cities of Tehran, Ahvaz and Tabriz have been used. Applying the aforementioned models as well as those already developed for Europe and the U.S., PWs are derived at the position of these stations in Iran. Further, in this research, root mean square error (RMSE) and bias are the measures for selecting the optimal model. Here, the bias and the RMSE (between GPS and radiosonde derived PWs) for the proposed linear model for Iran is 1.44 mm and 4.42 mm, and for quadratic model 2.18 mm and 4.74 mm respectively while, the bias and the RMSE for Bevis’ linear model is 2.63 mm and 4.98 mm and for Emardson and Derk’s quadratic models are 2.80 mm and 5.08 mm respectively. As such, it is observed that the bias of the proposed linear model for Iran is 1.19 mm and 1.36 mm less than the Bevis’ and Emardson and Derk’s models. In addition, the RMSE of the proposed linear model is 0.56 and 0.66 mm less than the RMSE of the later ones. This emphasizes that the estimation of the model coefficients must be based on regional meteorological


Introduction
Water vapor plays an important role in atmospheric processes.This parameters is one of the most variable components of the Earth's atmosphere.Its non-uniform distribution, which is due to the atmospheric phenomena above the Earth surface, depends both on space and time.Due to the limited spatial and temporal coverage of observations, estimates of water vapor are always shown with some uncertainties [Emardson et al. 1998, Jarlemark et al. 1998].
Being a key element in various atmospheric phe-nomena such as precipitation, the amount of water vapor is estimated by applying integrated water vapor (IWV) and precipitable water (PW).In other words, IWV signifies the amount of water vapor at a specific direction and the PW is the height of the equivalent column of the liquid water represented by , where t w is the liquid water density [Bevis et al. 1992].
The PW is used as one of the inputs for the numerical weather prediction models [Vedel and Huang 2003] as well as the pattern analysis and forecast of rain which is considered as the most important and fundamental hydrological cycle [Seco et al. 2012].In addition, the PW is also a useful tool for determining the trends in climate change [Bevis et al. 1992].The existing restriction in the water vapor analysis is a major source of error in shortterm (0-24 hours) forecasts of precipitation [Bevis et al. 1992].A linear power structure has been reported for the time variation of the PW [Hogg et al. 1981, Jarlemark et al. 1998].The spatial structure of the PW analyzed by Emardson et al. [1998] demonstrated the significance of local models in estimating the precipitable water.
The ground-based GPS stations have improved both spatial and temporal resolutions of measurements required for estimating the water vapor [Rocken et al. 1991].For this purpose, the propagation delay experienced by GPS signals in their transition through troposphere is used.In other words, the delay in this part of the atmosphere is partially controlled by the existing gaseous water.Askne and Nordius [1987] highlighted the required linear relation between the wet part of tropospheric delay and the PW.Over a decade later, Emardson et al. [1998] revised this ratio for the Scandinavian region.Following years, the same technique was also applied for estimating the PW in near real time [Karabatic et al. 2011].Due to their reasonable accuracy, the PWs PW IWV w t
were computed using GPS in rain pattern analysis and forecast models [Seco et al. 2012].
Taking into account the ZWDs, the current research primarily intends to assess the necessity and generate a local (conversion) model for the accurate estimation of the PW in Iran.Since radiosondes are expendable instruments, their high cost often restrict the number of launches (normally once to twice a day at 00,00 and 12,00 in UTC) at each station.It must be remembered that Iran possesses more than 107 permanent GPS stations hence; the possibility of computing the PW from ZWDs provides a unique opportunity for improving both spatial and temporal resolutions required for estimating the PW in this country.The paper analyzes the validity of the existing models for estimating the PW parameter through ZWDs in Iran.
The next section discusses the conversion models which have already been published in other parts of the world.The paper highlights the development process of these models as well as the PW computation through ZWDs.The third section introduces the most appropriate model to the climatic conditions of Iran.For this purpose, the PWs estimated through different conversion models are compared with the precipitable water directly extracted from the radiosonde data.The most appropriate model is then selected using the RMSE and the bias of GPS and radiosonde derived PWs.

Estimating precipitable water from zenith wet delays
The ZWD of the GPS tropospheric signals is due existing gaseous waters in the atmosphere.Taking Π as the required model for converting the ZWD to the PW, one can write: PW = ZWD/Π.To estimate ZWD, a slant tropospheric delay (STD) is computed along the propagation path of the GPS signals.This parameter is the signal delay between a satellite transmitting the signal and a ground station.When the STD is mapped from the line of sight to the zenith direction, it is called zenith tropospheric delay (ZTD).The ZTD is comprised of hydrostatic and wet components.The ZHD is modeled using temperature and pressure at the observing site.The ZWD is computed as the difference between the observed ZTD and the ZHD [Jin and Luo 2009]: (1) Various models are available for computing the ZHD [Hopfield 1969, Saastamoinen 1972, Goad and Goodman 1974].For example, according to Saastamoinen [1972]: (2) where P 0 , H 0 and { are the surface pressure in hPa, orthometric height of the station in meters and the latitude of the observing site respectively.
When no pressure and temperature sensors are available at the GPS station, the data can be extrapolated from the nearest meteorological site.As the vertical temperature gradient is used to transfer the temperature, Equation (3) can be used to transfer the pressure [Karabatic et al. 2011]: where P syn , T syn , h syn , are the pressure, temperature and height at meteorological station respectively.Parameter c is the temperature gradient and R is the gas constant.
Since, the zenith tropospheric delay is normally estimated through analyzing the GPS measurements, a time series of ZWD can be computed at each permanent GPS station.Using the ZWD estimates as well as the conversion model, the PW amount can be calculated independently from the radiosonde data profiles.
As the PW depends on the temporal and spatial variations of the atmosphere, the conversion model Π is also not a constant parameter.Moreover, the best estimate for this parameter is obtained when the climate condition is taken into account [Emardson and Derks 2000].
To compute PW from ZWD, there are various existing models which are classified into two groups of linear and quadratic here.The following subsections introduce these models in further detail.

Linear model
The ratio between the ZWD and the PW parameters [Bevis et al. 1992] is as followed: Here, R v is the specific gas constant of water vapor, that is 461.45JKg −1 K −1 , and the atmospheric refractivity constants k 3 and are approximately 3.7×10 5 K 2 mbar −1 and 17 Kmbar −1 respectively.T m is the weighted average of the atmospheric temperature.The value of T m is calculated using Equation ( 5).Therefore, computation of this parameter at each point requires water vapor pressure and temperature along a vertical profile: In the equation above, T is the temperature in Kelvin, z is the vertical coordinate and e w is the water vapor pressure.Further, this parameter is obtained using the following equation [Holton 2004]: Where L v is the latent heat of vaporization and is equal to 2.5×10 6 JKg −1 and T d is the dew point temperature.To calculate the distance between any two consecutive recorded levels, after taking logarithm of the equation proposed in [Nafisi et al. 2012] can be used: Here P i is the pressure on level i in hPa, R d is the gas constant for dry air and is equal to 287.05 JKg −1 K −1 , local gravity g is a function of latitude and height of the station and is the mean virtual temperature (T v ) between the two consecutive levels.The virtual temperature is defined as the one that a hypothetical system of dry air would have, in relation to the actual condition of (moist) air at the same density and pressure.This parameter is expressed as [Wallace and Hobbs 2006]: where M v and M d are the molecular mass of water vapor and dry air respectively.Their ratio is defined as follows [Andrews 2010]: T m can be estimated as a linear function of the temperature at 2 meters high above the Earth surface (T s ) [Bevis et al. 1992]: Analyzing more than 8000 vertical profiles of radiosonde stations in order to compute the unknown coefficients a i , i = 0,1 in the United States, Bevis et al. proposed the following functional relation between T m and T s :

Quadratic model
According to Bevis et al. [1992], the conversion model Π and the temperature behave similar to the parameter T m with respect to time.Based on the same, Emardson and Derks [2000] later proposed the following functional relation between the conversion model Π and the difference of surface and mean surface temperatures: Where T Δ is the surface temperature minus the mean surface temperature The mean surface temperature is computed by averaging the temperature obtained from radiosonde profiles.Further, the ZWD and PW estimates are required for computing the unknown coefficients a i ,i = 0,1,2.Based on that as well as using 120,000 radiosonde profiles in Europe, Emardson and Derks [2000] proposed the following functional relation [Emardson and Derks 2000]: To compute the zenith wet delay from the radiosonde profiles, one has [Rocken et al. 1995]: Here, this integral is computed from the surface (position of the station) to the upper troposphere where T and e w are derived from radiosonde profiles.Practically, since these parameters are not available in continuous form, the discrete form of the integral is used: In Equation (15), n is the number of levels in a radiosonde profile.Rocken et al. [1995], computed the PW using the following equation: Where t w is the density of liquid water.This integral is also computed from the surface (position of station) to the upper most tropospheric layer in the zenith direction.And for a similar aforementioned reason, the discrete form of this equation was used [Hagemann et al. 2003]: The PW, for each profile, can be calculated using the temperature and water vapor pressure derived from the radiosonde profiles.Then, the PW and ZWD for every profile can be used to calculate the conversion model Π.Thereafter, unknown coefficients as in Equation (12) are computed using the corresponding surface and the mean surface temperatures derived from the radiosonde profiles.

Results and discussion
The proposed linear and quadratic models briefly introduced in Sections 2.1 and 2.2 are computed first time for Iran.For this purpose, the current study could use profiles of 11 radiosonde stations.The acquired data covers the time interval of 1996 to 2012.These profiles are obtained from radiosonde lunches which were made once or twice a day at 00,00 and/or 12,00 in UTC.The corresponding profiles are available for download in the text format through http://www.esrl.noaa.gov/raobs/.
Existing gaps in the radiosonde data and the small number of lunching stations compared to the much larger density of the GPS as well as the continuous data provided by these equipments emphasize on the necessity of developing a suitable conversion model for Iran. Figure 1 illustrates the position of the radiosonde and the GPS stations in this country.

Derived models
In the current research, following Section 2, the coefficients of Equations ( 10) and ( 12) were locally computed using the data introduced in the previous section.This section briefly highlights the obtained results.
3.1.1.Empirical linear model P in the Equation ( 4) is only a function of the parameter T m .In order to compute this parameter in the proposed area of research, model of degree one, similar to the Equation ( 10), is used.Here, 54,330 radiosonde profiles have been taken into account.Further, the unknown coefficients of the proposed linear model are computed from the inputs given by the above radiosonde profiles.In other words, the inputs are the parameters T m and T s which are computed and derived respectively from the radiosonde data.Here, T m is computed using Equation ( 5).The obtained polynomial model is given in Equation ( 18).The estimated standard deviations for the interception and slope parameters of this model are ±0.1 degree of Kelvin and ±0.003 respectively.

Empirical quadratic model
In order to find a conversion model of the second type (see Equation 12), the ZWD and PW parameters are computed using the equations given in Section 2.2.For this purpose, the radiosonde data introduced above are used.The same number of radiosonde profiles is also used for estimating the aforementioned polynomial model (see Equation 18).The estimated model coefficients are given in Equation ( 19) where the root mean square error stands at 0.09.Moreover, the estimated standard deviations for the parameters of this model are ±0.001,±0.0001 and ±0.00001 for zero, first and second order coefficients, respectively.

Computing precipitable water from zenith wet delays
In the current study, for analyzing the accuracy of the developed models, three GPS stations located in climatologically different areas including Tehran, Ahvaz and Tabriz were selected.These three cities, respectively, have semi-dry, wet and semi-wet climatic conditions.
A comparison of PWs estimated from ZWDs with those directly derived from radiosonde profiles helps analyze the accuracy of conversion models.In fact, PWs estimated directly from the radiosonde data are used as benchmark for comparing the existing and the proposed models in this research.RMSE and bias are used to select the optimum result.
Thereafter, the PW time series obtained from the GPS zenith wet delays was analyzed using the least squares harmonic estimation technique.That analysis provided an insight into the frequency content of this time series and therefore could be considered as an independent approach for evaluating the proposed models.
To analyze the obtained results, two distinct time (18) (19) intervals were taken into account with two signified low and high rainfall, respectively.In fact, the Meteorological Organization of Iran reported the minimal precipitation at Tehran, Ahvaz and Tabriz stations in January, July and September of 2010, 2012 and 2011, respectively.The maximum precipitations were reported for these stations in May, February and December of the same years, respectively.As such, the zenith tro-pospheric delay has been computed for every period using the GPS measurements of the above stations.It is worth mentioning that the distance between the GPS stations of this study and the radiosonde stations located in Tehran, Ahvaz and Tabriz is 2.3, 36.4 and 11.5 kilometers.The height differences the GPS and the corresponding radiosonde stations for Tehran, Ahvaz and Tabriz are approximately 2.5, 1.0 and 35 meters re-PW DETERMINATION FOR IRAN  spectively.Charts 2 and 3 compare the RMS and the bias of the models for the GPS stations.
The analysis of the obtained results shows that the models found for Iran perform better than those given by Equations ( 4) and ( 13).Moreover, the RMSE and the bias for the proposed linear model is less than the quadratic one both for dry and wet months on all of the above GPS stations.Figures 4 and 5 illustrate the degree of misfit between the proposed linear model and the PW parameters which are computed directly from the high resolution radiosonde profiles as well as the time series of the estimated precipitable water at the Tehran station.These figures show the high and the low precipitations for May 2010 and January 2010, respectively.Dots are the precipitable water computed from the radiosonde records.
The frequency content of the PW time series derived from the GPS zenith wet delays can be considered as an independent approach for the evaluation of the linear model of this research.The ZWDs of Tehran, again a typical station, are used for computing the PW time series elements.This time series starts at 2010 and ends at 2012 (see Figure 6) The frequency content of this time series can be analyzed using one of the existing frequency analysis techniques.Here, the least squares harmonic estimation has been used [Amiri-Simkooei 2007] because, this method has been efficiently used for analyzing the frequency content of various time series i.e.GPS observations [Amiri-Simkooei 2007], time series of Total Electron Content [Amiri-Simkooei and Asgari 2012], tidal data [Mousavian and Mashhadi-Hossainali 2012] and coordinate time series for detecting slow slip events [Mousavian and Mashhadi-Hossainali 2014].
The application of the least squares harmonic estimation suggests annual and seasonal frequencies in the above PW time series.In other words, the annual and seasonal frequencies are the first that approved by the hypothesis test of this method (see Amiri-Simkooei [2007] for further details).Moreover, since the least squares harmonic estimation is sensitive to the length of the input data [Mousavian and Mashhadi-Hossainali 2012], a detailed analysis of the frequency content requires a PW time series of much larger length.Further, the existence of annual and seasonal frequencies in the PW time series indicates that the adopted linear model does not distort the expected behavior of the PW parameter.As a matter of fact, this can be considered as independent evidence on the applicability of the proposed linear model in the field of meteorology and climatology in Iran.

Conclusion
As discussed, the PW is a useful tool for determining the climate change.It is also used for producing the rain, as the most important and fundamental   hydrological cycle, prediction models.However, restriction in spatial and temporal resolution of the PW estimates seems to be the major source of error which can significantly and reliably be sorted out through the GPS-derived PW parameters.
The RMSE and bias of the proposed models are slighter than the existing ones.As the obtained results is valid for both dry and wet climatic conditions in Iran, the application of these models seems proper for computing the PW from the ZWDs in this country.In addition, the empirical linear model investigated in this research show that the RMSE and bias of the proposed model is less than the quadratic one.As such, the application of this model is also efficiently superior for computing the PW from ZWDs in Iran.
The accuracy analysis of the PW parameters obtained from the GPS zenith wet delays and the linear model of this study, as well as their frequency analysis approves the reliability of the GPS-derived precipitable water.Therefore, the linear model which has been determined here can be used for estimating the PW parameters together with their direct estimates from radiosonde profiles.This improved both spatial and temporal resolutions of this parameter in the country.

Figure 1 .
Figure 1.Position of the radiosonde stations (solid) and the permanent GPS stations in Iran (circles).The position of the IGS stations TEHN, AHVZ and TABZ is shown by a solid circles and names here.

Figure 2
Figure2(top to bottom).The RMS and bias of the proposed as well as the existing models in the rainy months: i) linear model proposed inBevis et al. [1992], ii) linear model proposed in this research, iii) quadratic model proposed inEmardson and Derks [2000] and iv) quadratic model proposed in this research.

Figure 3
Figure3(top to bottom).The RMS and bias of the proposed as well as the existing models in the dry months: i) linear model proposed inBevis et al. [1992], ii) linear model proposed in this research, iii) quadratic model proposed inEmardson and Derks [2000] and iv) quadratic model proposed in this research.

Figure 4 .
Figure 4. Time series of the estimated pecipitable water at station Tehran in January 2010.

Figure 5 .
Figure 5.Time series of the estimated precipitable water at station Tehran in May 2010.

Figure 6 .
Figure 6.Time series of PW based on the GPS ZWD and Π calculated using linear model of Iran during 3 years.