Generation of real‐time mode high‐resolution water vapor fields from GPS observations

Pointwise GPS measurements of tropospheric zenith total delay can be interpolated to provide high‐resolution water vapor maps which may be used for correcting synthetic aperture radar images, for numeral weather prediction, and for correcting Network Real‐time Kinematic GPS observations. Several previous studies have addressed the importance of the elevation dependency of water vapor, but it is often a challenge to separate elevation‐dependent tropospheric delays from turbulent components. In this paper, we present an iterative tropospheric decomposition interpolation model that decouples the elevation and turbulent tropospheric delay components. For a 150 km × 150 km California study region, we estimate real‐time mode zenith total delays at 41 GPS stations over 1 year by using the precise point positioning technique and demonstrate that the decoupled interpolation model generates improved high‐resolution tropospheric delay maps compared with previous tropospheric turbulence‐ and elevation‐dependent models. Cross validation of the GPS zenith total delays yields an RMS error of 4.6 mm with the decoupled interpolation model, compared with 8.4 mm with the previous model. On converting the GPS zenith wet delays to precipitable water vapor and interpolating to 1 km grid cells across the region, validations with the Moderate Resolution Imaging Spectroradiometer near‐IR water vapor product show 1.7 mm RMS differences by using the decoupled model, compared with 2.0 mm for the previous interpolation model. Such results are obtained without differencing the tropospheric delays or water vapor estimates in time or space, while the errors are similar over flat and mountainous terrains, as well as for both inland and coastal areas.


Introduction
Time-varying two-dimensional precipitable water vapor (PWV) fields (maps) are used in meteorological nowcasting, including the identification of events dominated by horizontal advection [Benevides et al., 2015], for assessing moisture transport in the lower troposphere [e.g., Mengistu Tsidu et al., 2015], for relating humidity fields to precipitation events [e.g., Boniface et al., 2009], for assessing the severity of tropical cyclones [e.g., Shoji et al., 2011], and for assessing the impact of new assimilated observations for forecasting precipitation [e.g., Yan et al., 2009]. Such maps are also essential for correcting synthetic aperture radar (SAR) images for atmospheric effects to enable small (and long wavelength) geophysical signals to be measured, including interseismic strain accumulation and postseismic motion, observations of which not only give insight into the mechanics of a fault but also play key roles in estimating the likelihood of future earthquakes [Wright et al., 2004;Gourmelen and Amelung, 2005;Fialko, 2006;Walters et al., 2013;Wen et al., 2012;Li et al., 2009a]. Dense PWV fields also enable Global Navigation Satellite Systems (GNSS) Network Real-time Kinematic (RTK) observations to be corrected for signal delays due to water vapor on propagating from space through the Earth's neutral atmosphere ("troposphere") to a ground-based receiver. Such corrections are essential for centimeter-level positioning, particularly heights, and enable (subject to sufficient GNSS base station coverage) Network RTK to be used for geophysical and engineering applications that have normally only used local base station RTK, such as river channel mapping [e.g., Brasington et al., 2003], glacier flow and debris mapping [e.g., De Paoli and Flowers, 2009], coastal erosion [e.g., Lee et al., 2013], crustal deformation [e.g., Genrich and Bock, 2006] and structural [e.g., Im et al., 2013] monitoring, precision farming [e.g., Pérez-Ruiz et al., 2011], embankment instability, and landslide monitoring [e.g., Gili et al., 2000;Zabuski et al., 2015]. oceanic areas with Sun glint, and/or above clouds over both land and ocean [Gao and Kaufman, 2003]. MODIS swaths do not coincide in time with images from SAR satellites such as ERS-1/2 (nominally 30 min of time difference), Envisat (nominally 30 min of time difference), or Sentinel-1 (nominally 4.5 h of time difference); therefore, for SAR atmospheric corrections, temporal interpolation errors will always prevail. The Envisat SAR satellite included the Medium-Resolution Imaging Spectrometer (MERIS), a passive push-broom imaging instrument measuring the solar radiation reflected from the Earth's surface and clouds in the visible and near-IR spectral range during the daytime. MERIS provided near-IR water vapor products above land or ocean surfaces under cloud-free conditions  or above the highest cloud level under cloudy conditions [Albert et al., 2001] at two nominal spatial resolutions: 0.3 km for full-resolution mode and 1.2 km for reduced-resolution mode. Both MODIS and MERIS near-IR PWV measurements have been validated against radiosonde and GPS-based products at discrete points with 1 mm and 1.2 mm respective RMS agreements [Li et al., 2009b;Li et al., 2006b;Li et al., 2005;Li et al., 2006c]. However, neither the MODIS nor MERIS instruments provide reliable PWV values above land or ocean surfaces when clouds are present, which restricts their use. Similarly, the Atmospheric Infrared Sounder (AIRS) instrument, also on board the Aqua satellite, not only provides global coverage amounting to 324,000 PWV profiles with a vertical resolution of 2 km and a spatial resolution at nadir of 45 km but also only operates in clear skies and in partly cloudy conditions [Aumann et al., 2003]. Comparisons between AIRS products and GPS PWV showed RMS differences of about 3-4 mm and 15% differences in the 2 km layers between AIRS and radiosondes [Divakarla et al., 2006;Raja et al., 2008]. PWV fields may also be obtained from global atmospheric models (both reanalysis and rapid update cycle predictions) but suffer from their coarse spacing, low temporal resolution (e.g., 3 or 6 h for the European Centre for Medium-Range Weather Forecasts (ECMWF) model), and failures during periods of atmospheric turbulence [Jolivet et al., 2014].
GPS provides a means of obtaining PWV with continuous temporal resolution without any cloud or weather dependence, albeit at discrete points where the GPS receivers are located. First, the time-varying tropospheric zenith total delay (ZTD) is estimated per receiver by using the GPS data alone, and then the zenith hydrostatic "dry" delay (which may be accurately modeled by using surface pressure data) is subtracted to leave the zenith wet delay (ZWD), which is spatially and temporally much more variable. The ZWD may then be readily converted to PWV by using estimates of the mean temperature of the atmosphere, based on surface temperature measurements [Bevis et al., 1992]. The pointwise GPS PWV measurements agree to those from radiosondes and microwave radiometers with standard deviations of about 1-2 mm [e.g., Ohtani and Naito, 2000;Liou et al., 2001; and may then be interpolated to provide PWV fields, which has been attempted by using different models. For example, Jarlemark and Emardson [1998] evaluated a gradient model, a linear regression in time model that ignored observational directions, and a turbulence model that yielded at least 10% improved RMS error than the other two models. Williams et al. [1998] used a structure function to model the water vapor variation in space, but with respect to a reference value, while Janssen et al. [2004] not only found that Inverse Distance Weighted (IDW) and Ordinary Kriging interpolation models perform better than spline interpolation but also only considered double-differenced ZTDs as they were considering interferometric synthetic aperture radar (InSAR) atmospheric corrections only. A deficiency of all these models is that they did not consider the terrain elevation dependence of water vapor, and hence, the interpolated values may contain large errors in regions with highly varying topography [e.g., Walters et al., 2013]. To deal with the water vapor stratification, Emardson and Johansson [1998] incorporated a height scaling function with a best linear unbiased estimator and suggested an interpolated ZWD accuracy of about 1 cm, but only one station in Sweden was considered and the height variation across the network was only about 200 m. Li et al. [2006a] proposed a GPS topography-dependent turbulence model but reported that interpolation models should be applied to ZTD/ZWD values differenced in time rather than the absolute ZTD/ZWD values themselves, as this can reduce the influence of topographic effects on the ZTD/ZWD variations. For interpolating undifferenced GPS ZWD point values, Onn and Zebker [2006] utilized a frozen flow hypothesis to model the water vapor variation in time. Then Xu et al. [2011] showed that incorporating this height scaling function approach with an interpolator model based on the estimator of simple Kriging with varying local means (we will refer to this model as SKlm + Onn) improved the ZWD interpolation RMS accuracy by 29% compared with using the Berrada Baby et al. [1988] semiempirical height scaling function. In a different approach to account for variations with topography, Bekaert et al. [2015] employed an InSAR phase observation-based power law correction model which used a fixed reference at the relative top of the Journal of Geophysical Research: Atmospheres 10.1002/2016JD025753 troposphere and described how the phase delay varies with altitude. To separate deformation and tropospheric signals, a frequency band insensitive to deformation is required. Benevides et al. [2016] also attempted to constrain GPS PWV with InSAR-derived PWV maps containing the topography signal. However, these models did not take into account that the InSAR measurements themselves have uncertainties of up to several centimeters and are susceptible to not detecting geophysical signals such as volcano inflation/deflation and interseismic slip rate [Williams et al., 1998]. Hence, we consider the SKlm + Onn model to represent the current state of the art for the generation of PWV maps.
Several previous studies have noticed the coupling effect of the tropospheric turbulence and terrain elevation dependency [e.g., Treuhaft and Lanyi, 1987;Li et al., 2006a;Xu et al., 2011;Benevides et al., 2016]. However, in the presence of strong atmospheric turbulence, previous models are inadequate for correcting SAR images to be used for subcentimeter-level ground motion monitoring [e.g., Walters et al., 2013;Fattahi and Amelung, 2015] or for the highest Network RTK positional precisions when such variations are not eliminated by data differencing. The aim of this study is therefore to improve the accuracy of GPS interpolated tropospheric water vapor maps by accounting for the coupling effect of both the terrain elevation dependency and tropospheric turbulence, demonstrate this over varying terrain, and compare with the SKlm + Onn model. We generate temporally continuous high spatial resolution maps by using undifferenced (neither spatially nor temporally) ZTD values estimated pointwise at GPS stations by using the precise point positioning processing (PPP) technique in real-time mode. Hence, they are directly applicable for correcting SAR images, enabling rapid response to ground motion from large earthquakes and volcanoes, for removing tropospheric effects in Network RTK GPS processing, and for meteorological nowcasting via potential assimilation in numerical weather models, noting that ZTD is the form of GPS-based water vapor observation used operationally by, for example, the Met Office [Bennitt and Jupp, 2012] and Météo-France [Poli et al., 2007]. We validate our interpolated absolute ZTD fields first by cross validation with the GPS pointwise ZTDs, and then by converting to PWV and comparing with the MODIS near-IR water vapor product, which we choose due to its high spatial resolution and accuracy . The MERIS near-IR water vapor product has similar quality with an even higher spatial resolution [Li et al., 2006c], but it is only available during the period from May 2002 and April 2012 before the International GNSS Real-time Service began (April 2013). We use data from all of 2015 (hence covering all four seasons) for a study region of California of about 150 km × 150 km, i.e., 37°06 0 to 38°30 0 N and 122°40 0 to 121°00 0 W, incorporating 41 Plate Boundary Observatory (PBO) GPS  Figure 1. This region was selected as there is a dense spacing of continuous GPS stations whose data have been extensively analyzed previously (including postprocessed ZTDs readily available for quality control), it is comparable to the size of a SAR image swath (e.g., 100 km for Envisat stripmap mode and 250 km for Sentinel-1 Interferometric Wide Swath mode), there are large variations in topography, and it experiences relatively cloud-free conditions as needed for validation with MODIS.

Iterative Tropospheric Decomposition Model
Tropospheric delays, especially the part due to atmospheric water vapor, vary both vertically and laterally over short distances and are often considered as the sum of (i) a stratified component highly correlated with topography which therefore delineates the vertical tropospheric profile and (ii) a turbulent component resulting from disturbance processes (e.g., severe weather) in the troposphere which trigger uncertain patterns in space and time [Williams et al., 1998;Cho et al., 2003;Li et al., 2006a]. We present an iterative tropospheric decomposition (ITD) model to effectively separate the turbulent-and elevation-dependent ZTD components. The ITD model decouples the ZTD into stratified and turbulent delays, which enable the more accurate interpolation of dense ZTD fields from pointwise values from a set of GPS reference stations across a region. It is defined as where, for the ZTD at location k, T represents the turbulent component and x k is the station coordinate vector in the local topocentric coordinate system; the stratified component is represented with an exponential function with coefficient β in which L 0 is, for the region considered, the stratified component delay at sea level and Þis the scaled height; ε k represents the remaining unmodeled residual errors, including unmodeled stratified and turbulent signals. The turbulent component usually consists of mediumto-long wavelength signals that can be interpolated by an IDW method. If n GPS stations are used in the region considered, then the IDW model reads as where w ui denotes the interpolation coefficient; u and i are indices for the user and reference stations, respectively; and d ui represents the horizontal distance from the user to reference station. Reference stations at distances more than 100 km from the user station show limited correlation [Emardson and Johansson, 1998], so they are not used.
In order to decompose the ZTD into stratified and turbulent components, which can account for substantial amounts of the ZTD but behave very differently, an iterative separation procedure was used: 1. The ZTDs from all GPS stations within the user to reference station decorrelated range limit are used to estimate initial values for the exponential coefficients β and L 0 , assuming that the turbulent component values in equation (1) are zero. 2. The residuals ε k , which are the summation of the unmodeled errors and the turbulent component, are computed by subtracting per station the stratified delay (as modeled by the estimated exponential coefficients) from the ZTD. 3. The turbulent component, T in equation (1), is computed per reference station from the residuals ε k by using the IDW function w ui given in equation (2): Steps 2-4 are repeated until the exponential coefficients β and L 0 converge. The final outputs are the exponential coefficients for the decorrelated range limit considered, plus the turbulent delay component and residuals per reference station. 6. The ZTD at the location of interest is obtained by interpolation of the reference station turbulent component and residuals and added to the stratified delay computed by using the final values of the exponential coefficients.
It should be noted that the assumption of the ITD model is that the turbulent component obeys the IDW interpolation law and the stratified component obeys the exponential law and, importantly, that these two components are not tightly coupled together. We later show (section 6 and Figure 9) that the two components are indeed separable and the convergence state can be reached rapidly. Although here we present the ITD model by ZTDs in zero difference mode, it is also suitable for interpolating differenced ZTD or PWV/ZWD. Hence, the input to the ITD model may be, as well as GPS ZTDs, other water vapor products such as those from the ECMWF numerical weather model, MODIS, or their combination.

Real-Time Mode GPS Precise Point Positioning Method
The pointwise ZTD values for the 41 study region GPS stations used for the generation of the interpolated maps were estimated by using real-time mode GPS PPP processing. We used a PPP software which is a highly self-modified version of RTKLIB (www.rtklib.com), employing an extended Kalman filter [Gelb, 1974] to estimate in the state vector the constant ambiguities and time-varying receiver coordinates, receiver and satellite clocks (considered as white noise), while the ZWD was estimated as a random walk parameter as a correction to an a priori ZTD from the UNB3 global empirical model [Leandro et al., 2006], employing the Global Mapping Function [Boehm et al., 2006a], and east-west and north-south tropospheric gradients were estimated. We used the ionospherically free pseudorange and carrier phase observables and applied absolute International GNSS Service (IGS) satellite and receiver antenna phase center offset corrections. We also applied corrections for antenna phase wind up [Beyerle, 2009], relativistic effects [Kouba, 2009], pseudorange differential calibration delays [Kouba, 2009], Earth tide [McCarthy, 1996], and ocean tide loading effects by using FES2004 coefficients obtained from http://holt.oso.chalmers.se/loading. Uncalibrated phase and pseudorange hardware delays were assumed to be absorbed by the (float) ambiguity parameters and estimated receiver clocks, respectively.
PPP relies on highly accurate satellite orbits and clocks [Zumberge et al., 1997], which are usually held fixed in postprocessed PPP. For our real-time mode processing, we used fixed satellite orbits from the International GNSS Service (IGS) Real-time Service, which were generated by decoding the IGC01 solution streams (products.igs-ip.net), every 15 s to match the GPS observation sampling rate of the PBO stations used. However, the satellite clocks have unpredictable behavior which makes their real-time prediction challenging , so we did not fix these to the real-time product values but estimated corrections to them by using satellite clock parameters with the Gundlich and Koch [2002] robust estimation method. Additional constraints were introduced to overcome the rank deficiency of the normal equations, namely where dt k,RTS is the initial value of the satellite clock given by the real-time product and acts as a pseudoobservation; res(dt k ) and dt k are the satellite clock residual and value, respectively. The satellite clock parameters were estimated as white noise parameters with a sigma of 0.001 ns, and the error messages contained in the real-time satellite clock product were used to determine the weights of the pseudo-observations in equation (4). An iterative process was used to identify some clock outliers which were hence ignored or assigned less weight in subsequent iterations [Gundlich and Koch, 2002].
We processed the data from 1 January to 31 December 2015 per GPS station in daily, discrete 24 h batches in real-time mode, with an elevation angle cutoff of 10°. The tropospheric delay was estimated every 5 min by using a process noise of 5.0e À 8 km/ ffiffi s p . To enable the fastest PPP solution convergence and separation of the ambiguities from the other estimated parameters, which is particularly problematic when using real-time satellite orbits and clocks [e.g., Yao et al., 2014], we applied loose constraints to a priori receiver coordinate values obtained from the PBO GPS Station Position Time Series (http://pbo.unavco.org/data/ gps). Nearly 70% of daily solutions converged within 30 min (convergence time here means from the beginning epoch to an epoch whose horizontal component bias is less than 10 cm and height component bias is less than 15 cm, and the overall standard deviation of its next 20 consecutive epochs also satisfies this requirement), with 90% of daily solutions converging within 50 min. The results presented hereafter are based only on the ZTD values after convergence was attained, according to these criteria.

Validation of Real-Time Mode GPS Pointwise ZTDs
To validate our real-time mode PPP (RTPPP) GPS ZTD estimates per station, we compared with postprocessed 5 min "truth" values computed by the Geodesy Laboratory at Central Washington University by using the NASA Jet Propulsion Laboratory/Caltech GPS Inferred Positioning System (GIPSY) software version 6.2 in PPP mode with fixed IGS final orbits and clocks. The values were obtained from ftp://data-out.unavco.org/ pub/products/troposphere/, and in their processing, the VMF1 gridded tropospheric mapping function was used together with ECMWF gridded a priori ZHDs and ZWDs [Boehm et al., 2006b], while estimating the ZWD and tropospheric gradients (east-west and north-south) with process noise values of 5.0e À 8 km/ ffiffi s p and 5.0e À 9 km/ ffiffi s p , respectively. We computed the differences between our RTPPP and GIPSY ZTDs at the common 5 min epochs, excluded all the outliers greater than three times the standard deviation, then for each station computed per day the mean difference and the root-mean-square (RMS) of the differences to assess the quality of the real-time mode processing. These are shown for a sample day (2 September 2015) in Figure 1, chosen as it is indicative of the median differences for all days of 2015. The mean of the per station differences across the whole network is 1.9 mm for the sample day, indicating that no significant systematic error exists between the RTPPP and GIPSY ZTDs, including stations in mountainous areas. The mean RMS is 10.1 mm, and more than 80% of the stations have an RMS value smaller than 12 mm, which is deemed sufficient quality for assimilation into real-time meteorological models [Shoji et al., 2011].
Apart from considering the spatial distribution of the differences, it is also important to assess the RTPPP performance over time. We therefore in Figure 2 show the 5 min RTPPP ZTDs plotted against the common epoch 5 min GIPSY ZTDs, from all 41 stations for the entire year, and plot the differences as a histogram. A linear regression fit gave GIPSY ZTD = 0.989(±0.002) × RTPPP ZTD + 0.024(±0.003), and the correlation coefficient between them was 0.99, demonstrating high consistency between the RTPPP and GIPSY ZTDs not just spatially but also temporally. About 82% of solutions show differences smaller than 15 mm with 73% below 12 mm. The RMS difference is 12.5 mm, commensurate with the spatial RMS difference and further indicating

Cross Validation of Interpolated ZTD at GPS Stations
Cross validation was used to evaluate the performance of the ITD model for interpolating ZTDs and compared with the SKlm + Onn model. While Xu et al. [2011] applied the SKlm + Onn model to ZWD, as ZWD also dominates the spatiotemporal variations of ZTD (with the ZHD being readily determined with surface pressure), we may also apply it to ZTD interpolation. Since the dry and wet components are both important for some applications such as InSAR atmospheric corrections [Elliott et al., 2008;Jolivet et al., 2014], and the separation of ZWD from GPS-derived ZTD will introduce additional uncertainties, it is better to use total delays rather than wet delays only. To provide an indication of the ZTD quality used for single-epoch SAR corrections, and to be commensurate with subsequent MODIS validations (section 5), we adopted the approach of Xu et al. [2011] and used one ZTD value per day (that at 14:00 Pacific standard time, i.e., local time) per station for the whole of 2015. The RTPPP ZTDs were interpolated for each GPS station in turn, with the 40 other GPS station ZTD values providing the input. Hence, for ITD, per epoch, L 0 and β of equation (1) were estimated for the network and the turbulent component estimated per station. Validation was carried out by comparing the interpolated ZTDs with the RTPPP ZTD estimates themselves at 14:00 local time and computing the RMS difference and mean absolute error (MAE) for all stations for each day. This was repeated by using interpolation with the SKlm + Onn model.
The time series of the cross-validated daily MAE and RMS differences from the ITD and SKlm + Onn ZTD interpolation models are shown in Figure 3. It is clear that the ITD model leads to both lower annual mean MAE and RMS difference values than SKlm + Onn; i.e., the MAE and RMS reduce from 6.2 mm to 3.2 mm and 7.4 mm to 4.1 mm, respectively. It can also be seen from Figure 3 that the improvement of ITD is greater than that of SKlm + Onn in colder seasons (e.g., between days 0 to 100 and 280 to 365), when the medium-to-long wavelength turbulent signals are greater than the short wavelength ones and can be effectively modeled by ITD. This, in turn, enables stratified delays to be separated from the turbulent component in the ITD iterative process and better determined by using regression analysis. The performance of the two models is more similar in the summer (i.e., from around day of year 150 to 220), indicating that short wavelength water vapor effects are significant and variable and cannot be fully mitigated by either model. Figure 4 for both the ITD and SKlm + Onn interpolation models, for all 41 stations for all of 2015. As for the time series, substantial reductions in the scatter are observed for ITD compared with those for SKlm + Onn; i.e., the RMS difference decreases Journal of Geophysical Research: Atmospheres 10.1002/2016JD025753 from 8.4 mm to 4.6 mm. Small correlation and linear fit improvements are also obtained with ITD compared with SKlm + Onn: the correlation coefficient increases from 0.98 to 0.99, the linear fit improves from a slope of 1.025 to 1.012, and the intercept from 0.026 to 0.004 m, as detailed in Table 1. Furthermore, the proportion of differences under (magnitude) 10 mm increases from 61% for SKlm + Onn to 89% for ITD and increases from 32% to 69% for the proportion under 5 mm magnitude.

Cross comparisons of the daily interpolated ZTD values are shown in
When considered for all stations for the entire year, the ITD interpolation model has been shown to substantially improve on the SKlm + Onn model. However, in the cross validation, some large differences occurred (more than 2 cm magnitude for both models), which suggests that the interpolation result is influenced by the GPS reference station distribution. To investigate this, the annual MAE and RMS differences per station are plotted in Figure 5, for both the ITD and SKlm + Onn models. It can be seen that the smaller MAE and RMS differences occur where the station coverage is denser, but the SKlm + Onn MAE and RMS values show substantial degradation compared with ITD in the northwest of the region where there are fewer stations, e.g.,~12 mm MAE for SKlm + Onn compared with~5 mm for ITD. Meanwhile, the largest RMS value for any station is only~8 mm for ITD, improved from~12 mm for SKlm + Onn. Such station distribution degradations are consistent with the interpolation law and can be mitigated as far as possible by using stations with as uniform a distribution as possible. In terms of terrain effects on the MAE and RMS, for the ITD model, stations in the mountainous areas show approximately comparable precision with those at lower altitudes, whereas with SKlm + Onn, larger MAE and RMS values arise and the same applies in coastal areas. We attribute this improvement to the height scaling model working under the assumption that the turbulent delays are insignificant and have short wavelength compared with stratified delays; therefore, the height scaling is easily biased when strong tropospheric turbulence with medium-to-long wavelength signals occur.

Validation of Interpolated PWV Maps With MODIS
To provide further validation of the ITD interpolation model and its improvement over the SKlm + Onn model, including a detailed spatial resolution assessment, and to provide an accuracy assessment with an independent data set, the RTPPP GPS pointwise ZTD values were converted to PWV, interpolated to 1 km pixels across the entire study region, and compared with the MODIS near-IR PWV product. Pressure and temperature data at 5 min temporal resolution from co-located meteorological sensors were available at 4 of the 41 GPS stations and obtained from unavo.org. These were supplemented by 10 meteorological stations which were located up to 10 km outside the study region. The meteorological data were first interpolated to all 41 stations by using the  differential models, and according to their cross tests, the resulting pressure and temperature errors were less than 1 hPa and 2 K, respectively. The interpolated pressure measurements at each GPS station were used to directly compute ZHD by using the Saastamoinen [1972] model and subtracted from the RTPPP ZTD estimates, with the resulting ZWD pointwise values converted to PWV by using the Bevis et al. [1992] model with the interpolated surface temperature measurements. To enable the comparisons, 1 year MODIS Level 2 data from the Terra satellite were obtained across the study region, providing one PWV map at the Terra orbit track time of each day (during daytime, about 10:30 local time). The Level 2 data were generated at the 1 km spatial resolution of the MODIS instrument by using the near-IR algorithm [Gao and Kaufman, 2003]. About 30% of days had severe cloud conditions, so we excluded them as only a few grids can be obtained. Areas with cloud conditions or above water were also masked, and only the cloud-free land areas were used in the comparison, with the percentage of data points available and used in the comparisons shown in Figure 6e, with the average cloud-free pixels over the entire year being 35%.
As the first step, cross validation was carried out at each GPS reference station, using one GPS PWV per day per station, taken at the MODIS acquisition time. We used MODIS PWV as "truth" values to validate the interpolated PWV at each GPS station (cross test), and the MODIS PWV values were averaged over boxes of 3 × 3 pixels centered on the GPS station's location if the centered pixel was missing. While in such cases the spatial resolution of the MODIS PWV maps was effectively reduced from 1 km × 1 km to 3 km × 3 km, this resolution is still much greater than the~15 km GPS station spacing across the study region, providing a compromise between the resolution and maximizing the number of MODIS pixels for comparison. PWV cross comparisons for all daily values available for each GPS station for the whole of 2015 are shown in Figures 6a and 6b, with the RMS of the differences being 1.48 mm for the ITD model and 1.73 mm for SKlm + Onn, as also listed in Table 1. The ITD model also results in a better linear regression fit, with a slope of 0.97 and an intercept of 0.33 mm compared with respective values of 0.95 and 0.63 mm with the SKlm + Onn model, and the correlation coefficient increased from 0.97 to 0.98. The SKlm + Onn model also produced a greater mean difference than ITD during times of lower PWV content (i.e., 0-10 mm), with the more effective ITD modeling in the times with less water vapor (i.e., from midautumn to midspring) being consistent with the GPS cross-validation results shown in Figure 3, i.e., the improved separation of medium-to-long wavelength turbulence signals and improved performance across mountainous areas.
The interpolated RTPPP PWV values and resulting maps were then compared spatially with the MODIS PWV maps. The results are shown (RTPPP GPS minus MODIS) in Figure 7 for MODIS images acquired on 3 September and 19 November 2015, chosen as they are sample days which are virtually free from cloud conditions across the whole study region. The height for each grid cell was resampled by the 30 m Shuttle Radar Topography Mission DEM. Some large differences are visually apparent, mostly across the areas with frequent cloud masks and near San Francisco Bay, but MODIS PWVs above water areas also involve a different retrieval algorithm compared to those above the land, resulting in differences and discontinuities at the land edge. Furthermore, any values over water areas have been removed since PWVs above water (bay, lake, or ocean) share different characteristics from PWVs over land areas which cannot be well-described by the interpolation model [Sobrino et al., 2003]. On average, the mountainous areas give more negative differences than flat terrain, showing that MODIS tends to overestimate PWV compared with GPS with increasing altitude (i.e., small PWV values), as previously found by . It can also be seen in Figure 7 that the edge areas with fewer GPS stations produce larger differences than central areas, confirming as discussed in section 4 that improved GPS station coverage will improve the quality of interpolated PWV maps. The ITD model produces smoother difference maps than SKlm + Onn and has a lower percentage of large differences. ITD also performs much better than SKlm + Onn in coastal areas where the PWV is more changeable and gives more complicated turbulent components.
To further spatially assess the improvement of ITD RTPPP PWV values with MODIS over those from the SKlm + Onn model, differences between the RTPPP interpolated PWV and all MODIS pixels (cloud free) were computed across the study region bounded by the GPS stations for the entire year 2015 and for all weather conditions. The differences are shown in Figures 6c and 6d, for each of the ITD and SKlm + Onn models. In general, most PWV pairs are located along the 1:1 line, implying a strong correlation between the RTPPPbased PWV and MODIS PWV maps. The SKlm + Onn PWV plot exhibits greater differences in comparison with the ITD PWV plot. The scatters are distributed rather unsymmetrically, especially when the PWV amounts fell between 5-15 mm (typically for mountainous areas) and 30-35 mm due to the lower scale factor (0.912 versus 0.934).
To illustrate the finer spatial detail of PWV and the performance of interpolated RTPPP PWV, in Figure 8 we plot ITD-based PWVs over mountainous areas, shown as 37°09 0 to 38°30 0 in latitude and À122°12 0 to À121°00 0 in longitude ( Figure 8a) and 37°09 0 to 37°50 0 in latitude and À122°30 0 to À121°00 0 in longitude (Figure 8b difference between MODIS and ITD PWV for the eight profiles considered is 1.51 mm, and the RMS differences for mountain (gray polygon in Figure 8) and flat areas are 1.57 mm and 1.47 mm, respectively. These agreements demonstrate that the ITD model is capable of retrieving detailed water vapor distributions over a wide region, thereby showing its potential application for monitoring local extreme weather events.

Discussion on Tropospheric Turbulence
The principal aim of the ITD model is to separate the elevation-dependent ZTD/PWV component from the turbulent component, which is the most variable and uncertain part, and can easily bias the vertical ZTD scaling, making the separation of the two components challenging. Due to the constraints of the density of GPS stations, only medium-to-long wavelength turbulent signals are expected to be successfully modeled by using ITD. To illustrate the size and variation of the turbulent component, and the importance of iterating the solution until convergence arises, a sample with three GPS stations was considered: P177, P230, and S300 (Figure 1), which are in different parts of the study region, are at different elevations, and are at varying distances from the nearest other GPS reference stations. P177 is near the ocean, while S300 and P230 are in mountainous areas with elevations of~500 m and~700 m, respectively. Three epochs from different seasons (spring, summer, and autumn) were considered, and the variation of the turbulent component and its convergence with the number of iterations is illustrated in Figure 9.
It is clear from Figure 9 that the turbulent ZTD component can reach several centimeters and can be efficiently separated from the elevation-dependent component. Although the first iteration enables the majority of the turbulent component to be determined, the subsequent iterations are needed for robust estimation. Figure 9 (fourth column) further indicates the performance improvement with increasing number of iterations, with the RMS difference (computed through the RTPPP ZTD cross validation from all 41 stations at the corresponding epoch of each row) becoming smaller and tending toward convergence. Around six iterations are typically needed for convergence, after which submillimeter RMS changes arise, with improvements of typically less than 1 mm. As the most important feature of the ITD model, the convergence tendency in Figure 9 reveals that the turbulence effect can be reduced by separating the stratified and turbulent components through iteration. The decoupling procedure acts very similarly to robust estimation, in which the ZTDs from stations exhibiting strong turbulence will contribute less in height scaling but account for more in the turbulent delay interpolation. In this way, the systematic patterns of turbulence resulting from local weather conditions or topography [Betts, 2001;Cho et al., 2003] can be modeled better. The iteration also allows for the detection of ZTD outliers, which is a not uncommon occurrence in real-time PPP due to the unpredictable behavior of the satellite clocks. Consequently, the ITD model enables both fitting of the tropospheric vertical profiles and modeling of the turbulence processes. As the fourth column in Figure 9 suggests, successfully accounting for this results in the overall RMS difference of the cross validation test reducing and converging.

Conclusions and Outlook
In this study, an iterative tropospheric decomposition model has been developed for the generation of highresolution regional tropospheric PWV fields (maps) from the interpolation of pointwise GPS ZTDs, without any data differencing. For a California study region of around 150 km × 150 km, the approach of decoupling the terrain elevation dependency and the tropospheric turbulence contribution to ZTD in an iterative procedure (typically 4-6 iterations were required) led to improved accuracy interpolated tropospheric water vapor maps over those based on previous studies, such as the tropospheric turbulence and elevation-dependent model SKlm + Onn of Xu et al. [2011]. To not only be applicable to postprocessed SAR atmospheric corrections but also facilitate SAR for rapid response to monitoring earthquakes and volcanoes, Network RTK positioning, and meteorological forecasting, we used real-time mode PPP GPS ZTD values estimated every 5 min (which were validated with postprocessed GIPSY ZTDs with an overall RMS difference of 12.5 mm for all 41 stations for all of 2015) to generate the tropospheric maps. Cross validation of the GPS ZTD values resulted in 4.6 mm RMS differences by using the ITD model compared with 8.4 mm by using the SKlm + Onn model, using one value per station per day (14:00 local time) for all of 2015. Whereas the SKlm + Onn interpolation model has degraded performance over mountainous areas, the cross-validation ITD model RMS and mean differences are similar for both mountainous and flatter terrain and also similar for both coastal and inland areas. The cross-validation improvements using ITD are smallest in the summer months. Spatially, we generated PWV values for 1 km pixels for all land-covered parts of the region and compared with daily MODIS PWV near-IR product values, with the RMS difference for the year being improved from 1.96 mm by using the SKlm + Onn model to 1.71 mm by using ITD. Furthermore, the spatial PWV gradients using the ITD model and MODIS across a variety of topography were nearly identical to each other. The overall RMS difference between MODIS and ITD PWV profiles is 1.51 mm, and the RMS differences for mountain and flat areas are 1.57 mm and 1.47 mm, respectively. Hence, the ITD PWV fields are also able to reveal detailed water vapor information over varying terrains.
To provide an indication of the potential of the real-time mode ITD model interpolated tropospheric maps for meteorological and geodetic applications, including revealing detailed information of local weather processes, we show in Figure 10 the detailed 2-hourly PWV information during a rainfall process over the study region on 2 November 2015 (10 A.M. to 8 P.M. UTC). The arrows represent the PWV increasing (upward) or decreasing (downward) during each preceding 2 h. One important fact is that the PWV over mountainous areas decreased during the whole rainfall process, but other areas experienced increasing and decreasing PWV before and after the rainfall, respectively. We do not explain the patterns further as the focus of this paper is on showing the quality of RTPPP ZTD and PWV maps.
The generated spatially dense PWV fields with continuous, high (5 min) temporal resolution are not only suitable for correcting atmospheric effects in SAR images at the instant of acquisition but will also ensure the identification of water vapor variation from ground motion between image acquisition times [Foster et al., 2006]. What is more, the high performance of the dense PWV maps using the ITD model is especially useful for mitigating the effects of water vapor for SAR measurements in mountainous areas, which usually suffer from vertical stratification and turbulent mixing due to the orography [Wadge et al., 2002].
In the near future, more research is planned to investigate the capacity of high-resolution real-time tropospheric products using the ITD interpolation model. This includes the testing of the ITD model on ECMWF ZTD, MODIS, and radiosonde PWV in addition to the GPS interpolations considered here. The abundant data sources will contribute to generate global coverage and all weather condition dense water vapor fields and to further investigate detailed tropospheric turbulence characteristics and processes.