High Spatial Resolution mapping of Precipitable Water Vapor using SAR interferograms, GPS observations and ERA-Interim reanalysis

A high spatial resolution of the Precipitable Water Vapor (PW V) in the atmosphere is a key requirement for the 10 short-scale weather forecasting and climate research. The aim of this work is to derive temporally -differenced maps of the spatial distribution of PW V by analyzing the atmospheric delay “noise” in Interferomet ric Synthetic Aperture Radar (InSAR). A t ime series maps of differential PW V were obtained by processing a set of ENVISAT ASAR images cover the area of Southern Californ ia, USA from 06 October 2007 to 29 November 2008. To get a more accurate PW V, the component of dry delay is calculated by using ERA-Interim reanalysis products. In addition, the ERA-Interim was used to compute the 15 conversion factors required to convert the zen ith wet delay to water vapor. The InSAR-derived differential PW V maps were calibrated by means of the GPS PW V measurements over the study area. We validated our results against the measurements of PW V derived from the MEdium Resolution Imaging Spectrometer (MERIS) which is located together with ASAR sensor onboard the ENVISAT satellite. Our comparative results show strong spatial correlations between the two data sets. The difference maps have Gaussian distributions with mean values close to zero and standard deviations below 2 mm. The 20 advantages of the InSAR technique is that it provides water vapor distribution with a spatial resolution as fine as 20 m and an accuracy of ~2 mm. Such a high spatial resolution maps of PW V could lead to much greater accuracy in meteorological understanding and quantitative precipitation forecasts.

to compute the conversion factors required to convert the zenith wet delay to water vapor. The 23 InSAR-derived differential PWV maps were calibrated by means of the GPS PWV 24 measurements over the study area. We validated our results against the measurements of PWV 25 derived from the M Edium Resolution Imaging Spectrometer (M ERIS) which is located together 26 with ASAR sensor onboard the ENVISAT satellite. Our comparative results show strong spatial 27 correlations between the two data sets. The difference maps have Gaussian distributions with 28 mean values close to zero and standard deviations below 2 mm. The advantages of the InSAR 1 technique is that it provides water vapor distribution with a spatial resolution as fine as 20 m 2 and an accuracy of ~2 mm. Such a high spatial resolution maps of PWV could lead to much 3 greater accuracy in meteorological understanding and quantitative precipitation forecasts. With 4 the launch of Sentinel-1A satellites, every few days (12 days) a new SAR images can be 5 acquired with a wide swath up to 250 km, enabling this a far unique operational service for 6 InSAR-based water vapor maps with unprecedented spatial and temporal resolution. Although the water vapor contributes only about 10% of total atmospheric delay, this source of 13 error is not easily eliminated due to its high spatial and temporal variability. Our aim in this 14 paper is to investigate the tropospheric delay "noise" of InSAR as a meteorological signal to 15 measure the water vapor content in the atmosphere. We will present a new approach for accurate 16 water vapor constructing with a high spatial resolution by combing InSAR observations, GPS 17 data, and a Global Atmospheric M odel (ERA-Interim), and evaluate the performance. short-term (0-24 hour) precipitation prediction. The advantage of satellite-based InSAR, a 29 relative new tool for measuring water vapor content, is that it could provide maps of water 30 vapor with a spatial resolution as fine as 10-20 m over a swath of ground about 100 km wide. 31 The Sentinel-1A satelllite, a new SAR mission launched in April 2014, providing a wide swath 1 up to 250 km that gives this technique a more intersting perspective for meteorology studies. 2 3 In this paper, we used the InSAR data in combination with GPS measurements and ERA-4 Interim reanalysis products to prescisely estimate the water vapor content in the atmosphere. 5 The main concept of InSAR for constructing water vapor maps is that the tropospheric phase 6 delay is considered as our interested signal to be extracted and the other phase components are 7 treated as noise to be removed. The tropospheric phase delay mainly consists of two 8 components: hydrostatic delay and wet delay. The hydrostatic delay varies with local 9 temperature and atmospheric pressure which is smooth in time and space, while the wet delay 10 varies with water vapor partial pressure which is more spatially and temporally varying. Within 11 a typical interferogram area of 100×100 km, the pressure usually varies less than 1hPa, while 12 a significant changes of the water vapor partial pressure are common. Consequently, the wet 13 delay in the interferogram is much greater than the hydrostatic delay. Therefore, most studies 14 have focused on estimating the wet delay and neglected the hydrostatic delay. However, recent 15 studies also show that hydrostatic delay varies significantly at low elevation and cannot be 16 neglected (Doin et al., 2009;Jolivet et al., 2014). Thus, to obtain accurate PWV maps, 17 hydrostatic delay in InSAR must be taken into account. In this work, we compute the 18 component of hydrostatic delay by using ERA-Interim reanalysis products. Using the water 19 vapor conversion factor, the InSAR-derived zenith wet delay is then mapped onto Precipitable 20 Water Vapor (PWV), a quantity representing the water vapor content in the atmosphere. In this 21 study, the outputs of temperature and specific humidity from ERA-Interim model are used to 22 estimate this water vapor conversion factor. It should be noted that water vapor maps from 23 InSAR are derived from the difference between the water vapor variations during two Synthetic 24 Aperture Radar (SAR) images with a temporal separation of one or more days, we call ∆PWV 25 hereafter. The temporal interval relies on the space-borne InSAR mission: 1 day (tandem ERS-26 1/2), 11 days (TerraSAR-X, Cosmo-SkyM ed), 12 days (Sentinel-1), 35 days (ENVISAT-ASAR, 27 RADARSAT) and 46 days (ALOS-PALSAR). The main problem of this operation is that the 28 ∆PWV maps is a relative measurement and requires absolute observations to calibrate each 29 ∆PWV map. The calibration procedure was implemented by using absolute measurements of 30 PWV from a few GPS stations in our study area. After that, the calibrated ∆PWV maps were 31 evaluated by comparing with the ∆PWV from the collocated GPS stations. Finally, we made 32 a comparative analysis of ∆PWV maps from InSAR and M ERIS pixel by pixel, and by 1 inspecting the spatial properties. 2 2 Study area and data sets 3 We carried out the study using data sets collected in the Los Angeles basin located in Southern 4 California, USA. This study area neighbors the Pacific Ocean in the west and southwest and 5 thus is rich with atmospheric water vapor and it is well covered by a dense network of continous 6 GPS receivers. These conditons make it particularly suitable for atmospheric water vapor 7 studies. Figure 1 shows the topography map of the study area. A set of N=8 ENVISAT ASAR 8 SLC images were acquired over this region for the period between 06 October 2007 to 29 9 November 2008. The image was acquired during descending passes, Track 170, with the 10 average look angle θ = 22.6°. Actually, the value of look angle θ varies over the SAR scene 11 from near range to far range between 16.5° to 23.2°. Accuracy may improve, if local look angle 12 of every pixels within interferogram is considered when calculating the mapping function. We The analysis centers provide tropospheric data products including zenith atmospheric delay are 30 archived at the UNAVCO Data Center and are openly and freely available 31 (http://www.unavco.org/data/data.html). The availability of GPS measurements also allowed 32 us to separate possible surface deformation from the atmospheric signals in differential 1 interferograms. The red triangles in Fig. 1 represent the locations of GPS stations. 2 3 The ERA-Interim reanalysis from European Centre for M edium-Range Weather Forecasts 4 (ECM WF) is used to produce maps of hydrostatic delay and water vapor conversion factor. 5 ERA-Interim is a global atmospheric model which was instigated to address some of the 6 problems seen in ERA-40 (Dee et al., 2011). It is based on 4-dimensional variational 7 assimilation of global surface and satellite meteorological data. The outputs of ERA-Interim 8 used in our study are estimates of temperature, specific humidity, and geopotential height, 9 defined at 37 pressure levels (1000-1 hPa), and a spatial resolution of 0.75° (~75 km). The 10 black crosses in in Fig. 1 shows the distribution of ERA-Interim model grid nodes used in this 11 study. The M ERIS is located together with ASAR sensor on board of ENVISAT satellite 12 (Bennartz and Fischer, 2001), thus the simultaneous water vapor measurements from M ERIS 13 were used as a reference data for comparison and evaluation. 14 15 3 Estimating PWV from InSAR 16 Here, we present the methods for obtaining zenith wet delay from SAR interferogram and 17 converting it to PWV. In Sect.3.1, the retrieval of zenith wet delay from SAR interferogram is 18 described. Section 3.2 describes the method for computing the conversion factor required to 19 map the zenith wet delay onto PWV by using ERA-Interim reanalysis. In Sect.3.2, the approach 20 for calibrating the PWV estimated from InSAR using GPS observations is discussed. 21

Atmospheric delay in InSAR 22
The unwrapped interferometric phase for each pixel in an interferogram is given by the 23 superposition of several components including topography, Earth surface displacement, and 24 atmosphere. It can be written as: 25 where ∅ is the phase contribution from land topography, ∅ represents the ground 27 deformation between the acquisitions, ∅ counts for the phase caused by inaccurate satellite 28 orbit, ∅ indicates the atmospheric state variations during SAR acquisitions and ∅ 29 denotes the noise component including system thermal noise, decorrelation noise, co-registration noise and processing noise. The contribution of topography is compensated for by 1 using an external DEM (30 m SRTM DEM are used in this study, Fig. 2a) Angeles Basin areas was observed in a number of interferograms although those interferograms 7 possess short temporal baselines. The rapidly subsiding displacement area in the Los Angeles 8 Basin region was masked out from all interferograms to avoid mixing the atmospheric signal 9 with surface deformation. After subtracting the topographic phase and orbital ramp, the residual 10 phase in the unwrapped interferograms only result from the atmospheric delay, which can be 11 split into hydrostatic, wet, liquid, and ionospheric components. In this study, we focus only on 12 the hydrostatic and wet components in the troposphere, as the delay induced by liquid water is 13 expected to be small under usual conditions, and the ionspheric components are assumend to 14 be small for C-band SAR signal (Hanssen, 2001). Thus, leading to the troposhperic phase delay 15 ∅ as 16 where 18 (3) 19 The hydrostatic delay ∅ ℎ is calculated using the specific gas constant for hydrostatic air , 21 gravity acceleration at ground level 0 , and air pressure . The wet delay ∅ is computed 22 using the partial pressure of water vapor , water vapor specific gas constant , and 23 temperature . 0 is the ground level and represents a reference height (30 km used in 24 this study) above which the delay is assumed to be nearly unchanged with time. The 25 atmospheric refractivity constants 1 , 2 and 3 are determined in (Smith and Weintraub, 26 is the radar wavelength and 27 − 4 is a scale factor to convert the delay in millimeter into phase in radian. is the radar 28 incidence angle and the factor be up to 15 mm in our study area. Therefore, in order to accurately derive the wet delay, the 10 hydrostatic delay must be precisely estimated and subtracted from the total tropospheric delay. 11 This delay can be calculated if the atmospheric pressure is known along the signal propagation 12 path or along the zenith direction. In this work, we used the vertical profiles of atmospheric 13 pressure provided by ERA-Interim reanalysis products to predict this component of hydrostatic 14 delay. We interpolated the atmospheric pressure onto altitude profiles at each ERA-Interim 15 model grids using a spline interpolation and calculated the hydrostatic delay using Eq. (3). The 16 resulting vertical profiles of hydrostatic delay were horizontally interpolated to the resolution 17 of SAR interferogram using a bilinear interpolation. We also used the outputs of temperature 18 and relative humidity from ERA-Interim to produce the maps of water vapor conversion factor 19 using the same interpolation strategy; this will be discussed in next subsection. The map of 20 hydrostatic delay is displayed in Fig. 2c, this delay represents a long-wavelength signal and is 21 smooth in space, rose up to 1 cm on the mountain areas. The slant wet delay (Fig. 2d

Conversion of ZWD into PWV 27
The zenith wet delay is considered as a measurement of the water vapor content in the 28 atmosphere. The relationship between the ZWD and PWV can be expressed as (Bevis et al., 29 1994): 30 where κ is the water vapor conversion factor and Π = κ −1 is calculated by the following 2 where is the density of the liquid water (listed in Table 2). is the weighted mean 5 temperature of the atmosphere and it is related to the surface temperature in degrees Kelvin 6 (Bevis et al., 1992). 7 = 70.2 + 0.72 × 8 Using this relationship to estimate will produce approximately 2% error in PWV (Bevis et 9 al., 1992). The most accurate way to compute the mean temperature is to calculate the following 10 integral equation between the ground surface 0 and the reference height , given by 11 (Davis et al., 1985), 12 The value of Π is dimensionless and usually ranges from 6.0 to 6.5 (and could be up to 7.0 at 14 some circumstances) (Bevis et al., 1992). For the purpose of rough conversion between ZWD 15 and PWV, an empirical constant Π = 6.25 (κ = 0.16) was used. However, the actual value of 16 κ changes with water vapor pressure and temperature that minor errors in κ could result in 17 significant biases in PWV. For example, using the constant value κ = 0.16 and assuming the 18 ZWD as 200 mm, the corresponding value of PWV is 32 mm. However, if the value of κ is 19 computed using Eqs. (7) and (9) as 0.15, then the value of PWV will be 30 mm. In fact, the 20 larger the ZWD, the more critical is the value of κ. Rather than using the empirical constant 21 value, we evaluated the conversion factor κ at each pixel of the SAR interferogram using 22 ERA-Interim reanalysis. To compute the weighted mean temperature , the outputs of ERA-23 Interim we used are the vertical profiles of temperature and relative humidity. The relative 24 humidity is converted to partial pressure of water vapor by a mixed Clausius-Clapeyron law 25 Interim model grids (indicated as the black crosses in Fig. 1 From Fig. 3, we observed that the value of Π changes with , and Π is in the range of 6.09 1 to 6.79 in the year of 2007 (Fig. 3a), whereas it varies between 6.17 and 6.74 in the year of 2 2008 (Fig. 3b). The fitted average curves linearly decrease with rates of -0.0214/K and -3 0.0221/K, respectively. As expected that the value of Π is much higher on winter days (low 4 temperature) than summer days (high temperature). On the other hand, since the temperature 5 generally decreases with altitude in the troposphere, the conversion factor is correlated with the 6 elevation. Therefore, using the empirical value of κ = 0.16 is not appropriate for the whole 7 study area; rather its value is calculated using global atmospheric model ERA-Interim. Figure  8 4 shows the spatial distribution map of Π on 16 August 2008 produced by ERA-Interim. It can 9 be seen that the value of Π varies spatially and it has a higher value on the mountainous areas 10 than those areas with a flat terrain. We then averaged the spatial maps of Π at the two 11 interferometric acquisition time to derive the conversion factors for mapping the wet delay onto 12 water vapor. where GPS is the number of GPS receivers, ( ) is the number of InSAR pixels located 2 within the circular area around the ℎ GPS receiver, ∆PWV GPS is the temporal difference 3 of PWV between master and slave dates by GPS, ∆PWV InSAR represents the ∆PWV 4 estimated by InSAR. Finally, the relative map of the ∆PWV in interferograms were calibrated 5 by adding the constant K to ∆PWV InSAR map. 6 4 Results and discussion 7 In this section, we will evaluate and validate the performance of InSAR-based water vapor 8 maps by comparing the calibrated ∆PWV estimated from InSAR to ∆PWV measurements 9 from GPS, as well as measured values from M ERIS. The evaluation was conducted as follows.  study. In summer, high temperature causes water to evaporate from the surface of lakes and 7 oceans, resulting in higher PWV content and more variable, whereas in autumn and winter, the 8 lower and smoother PWV were observed due to dry weather conditions. 9

InSAR ∆ measurements 21
The eight ENVISAT ASAR images are used for interferometric processing. The constraints for 22 normal baseline (< 300 m) and temporal baseline (<105 days) are used in order to minimize the 23 effects of ground deformation and decorrelation noise. Table 1  Adaptive power spectrum filter have been applied to interferograms to reduce phase noise 1 (Goldstein and Werner, 1998). All interferograms were multilooked by 40 looks in azimuth and 2 8 looks in range to enhance the coherence quality and improve the phase unwrapping accuracy. 3 The multilook processing resulted in a reduction of the spatial resolution of the interferograms 4 to 160×160 m. The wrapped phases were unwrapped using a branch cut algorithm (Goldstein 5 et al., 1988) and possible orbital errors were corrected by network de-ramping method. 6 Oscillator drifts induce a systematic phase ramp in the interferogram from ENVISAT satellite 7 (M arinkovic and Larsen, 2015), they were removed by the script provided in the StaM PS 8 software. The local rapid ground subsiding region was masked out. The wet delay differences 9 of InSAR were obtained by subtracting the component of hydrostatic delay predicted from 10 ERA-Interim. The wet delay differences were finally mapped onto ∆PWV maps using the 11 water vapor conversion factor as explained in Sect. 3.2. 12 Due to the fact that the unwrapped processing introduced an arbitrary constant into the phase, 13 all the ∆PWV maps from InSAR were relative measurements. Therefore, we shows the lowest water vapor content. The global negative correlation between ∆PWV and 28 altitude (Fig. 8b) implies that the absolute humidity in the bottom layer of atmosphere was 29 smaller at the acquisition time of the slave image than at the acquisition time of the master 30 image. The quantitative comparison of this interferogram is summarized in Table 3. It can be 31 seen that most of differences are smaller than 2 mm. The M AE of ∆PWV between GPS and 32 between InSAR and GPS at stations CGDM , ECFS and WLSN (indicated by the black arrows 1 in Fig. 8) were observed, especially the largest difference (-2.84 mm) at station WLSN. The 2 standard deviations of InSAR pixels located within the circular area around these three GPS 3 stations also show a high value (the fourth column in Table 3). The three GPS stations are 4 located in mountain areas with an altitude 730, 820, 1700 meters for the CGDM , ECFS and 5 WLSN stations, respectively. This interferogram also show a high value for height ambiguity 6 (290.90 m). Therefore, we can conclude that the large discrepancies between InSAR and GPS 7 for these three stations are most possibly due to the topographic phase error during 8 interferometric processing. 9

10
The comparisons of ∆PWV from the two techniques at each GPS station for the ten 11 interferograms are shown in Fig. 9. A good agreement between the InSAR and GPS was found 12 in the whole data sets. Large differences between InSAR and GPS at stations CGDM , ECFS 13 and WLSN were also found on those interferograms with a high value of height ambiguity 14 (interferograms 1, 2, 4 and 7 in Table 1). The PWV estimates from the two techniques 15 characterized by different sampling properties both in space and time. GPS can provide an 16 absolute value of the PWV every five minutes but refers to the parts of atmosphere observed 17 within a cone whose radius depends on the elevation cutoff angle, whereas InSAR gives a high 18 spatial resolution map of the ∆PWV with a time separation of 35 days or more. The high 19 temporal sampling of GPS and high spatial resolution of InSAR are complementary for 20 numerical weather modeling, which will improve the model resolution and give a better 21 understanding of the structure of atmospheric patterns. 22

Validation using water vapor measurements from MERIS 23
In this section, we will evaluate and analyze the accuracy of time series of the calibrated ∆PWV 24 maps derived from InSAR to confirm the performance of this technique as a tool for 25 constructing PWV maps. We carry out a cross-validation pixel by pixel using cloud-free water 26 vapor pixels by M ERIS acquired simultaneously with the ENVISAT ASAR images. The percentage of cloud-free conditions for M ERIS data we used in this study are larger than 1 90% except for the image acquired on 29 November 2008 having a coverage percentage of 80%. 2 For the sake of comparison, we built differences of PWV maps (∆PWV) from M ERIS. This is 3 performed based on the software package called TRAIN (Toolbox for Reducing Atmospheric 4 InSAR Noise) (Bekaert et al., 2015). In Fig. 10, are shown the calibrated ∆PWV maps derived 5 from the ten interferograms (in Table 1) and the corresponding ∆PWV maps from M ERIS data. 6 The first column shows the ∆PWV derived from InSAR that have been calibrated with GPS 7 observations. The In this paper, we presented the results of the temporal evolution of the PWV over Southern 22 California, USA using SAR interferograms during the period from 06 October 2007 to 29 23 November 2008. Interferograms were spatially averaged and spatial resolution was reduced to 24 160×160 m. In order to improve the quality maps of atmopsheric water vapor, the hydrostatic 25 delay was precisely estimated by using ERA-Interim reanalysis products. We also used the 26 outputs from ERA-Interim to produce maps of the conversion factor for mapping zenith wet 27 delay onto PWV at each pixel in the radar scene. All maps of ∆PWV derived from 28 interferograms were calibrated using a network of 29 continuous GPS stations located in the 29 SAR scene. The PWV estimates from InSAR and M ERIS show strong agreement with the data 30 from GPS. Since the GPS PWV estimates represent the average of the tropospheric effect within 31 a cone above the receiver, InSAR and M ERIS pixels were aggregated to enable a proper 32 comparison. The comparative analysis between InSAR and M ERIS ∆PWV maps demonstrate 1 strong spatial correlation with a less than 2 mm standard deviation for the difference. Our study 2 demonstrates that satellite Synthetic Aperture Radar interferometry can be applied to study the 3 spatial distribution of the PWV with a spatial resolution of 20 m and an accuracy of ~2 mm. 4 This advantage of InSAR provides unsurpassed insights in capturing the small-scale water 5 vapor distribution. This property could be important for numerical weather forecasting models. 6 Furthermore, forecasting models could take advantage of this source of water vapor maps to 7 enhance the accuracy of their assimilation sytems. In turn, the more accurate atmospheric 8 prediction models can be used to correct the tropospheric delay affected by water vapor in the 9 application of geodesy.