Estimation of Spatio-temporal Distribution of Precipitable Water Using Modis and Avhrr Data: a Case Study for Cyprus

In this paper, the atmospheric precipitable water (PW) over the area of Cyprus was estimated by means of Advanced Very High Resolution Radiometer (AVHRR) thermal channels brightness temperature difference (T). The AVHRR derived T was calculated in a grid of 5 × 5 km cells; the corresponding PW value in each grid cell was extracted from Moderate Resolution Imaging Spectroradiome-ter (MODIS) Level 2 product (near-infrared algorithm). Once the PW – T relationship coefficients corresponding to the area of Cyprus were calculated, the relationship was applied to AVHRR data for one month period. Radiosonde derived PW values, as well as MODIS independent PW values were used to validate the estimations and a good agreement was noted.


Introduction
Atmospheric water vapor is the principal contributor to the greenhouse effect and plays a key role in our understanding of the Earth's climate. Precipitable water (PW) is the amount of vertically integrated water vapor and can be expressed in g cm −2 , or as the height of an equivalent column of liquid water in cm. PW is an important component of the hydrological cycle and it is adopted as an input variable in global climatological studies. Moreover, it has the potential to support hydrological, biospheric and atmospheric modeling efforts, both at the local and regional scales, since it is widely used in energy budget and evapotranspiration studies. PW is also of an essential requirement in atmospheric correction of high spatial resolution satellite data and it also necessary Correspondence to: D. Hadjimitsis (d.hadjimitsis@cut.ac.cy) for the enhancement of the precision of land surface temperature estimates obtained from satellite data. The satellite derived high-resolution PW retrievals can be useful in anticipating the distribution of precipitation patterns. Thus, it is of particular importance to monitor PW variability on regional scales in order to monitor drought conditions and desertification processes.
Retrieving PW using Moderate Resolution Imaging Spectroradiometer (MODIS) near-infrared records is based on detecting the water vapor absorption of the reflected solar radiation, after it has been transmitted down through the atmosphere, reflected at the surface and subsequently transmitted up through the atmosphere to the sensor. The equivalent total vertical amount of water vapor is derived from a comparison between the reflected solar radiation in the absorption channel and the reflected solar radiation in nearby nonabsorption channels . Three nearinfrared channels were implemented on MODIS for water vapor sensing. PW is estimated over clear land areas, clouds and oceanic areas with sun glint (King et al., 2003). Typical errors for MODIS derived PW over land range between 5% and 10% .
A number of algorithms have been proposed to derive PW from Advanced Very High Resolution Radiometer (AVHRR) thermal infrared observations and can be classified into four main categories: simple split-window of thermal channels (Dalu, 1986), variance ratio (Jedlovec, 1990), regression slope (Goward et al., 1994), covariance -variance method and look-up table approach (Czajkowski et al., 2002). It should be noted that emissivity over land surfaces presents a significant complication in interpreting split window temperature differences, however accuracies of better than 0.5 cm water can be achieved (Eck and Holben, 1994 Theoretical analysis proposes an almost linear relationship between PW and T , with slope and intercept being determined by the surface emissivity. Spectral variation of emissivity is found to be the major factor that causes deviation of the intercept from a near-zero value . Assuming that the surface is Lambertian and the landatmosphere system is spatially homogeneous, a split window algorithm can be applied to AVHRR thermal infrared channels (Choudhury and Di Girolamo, 1995). As is has been shown in past studies (Dalu, 1986;Eck and Holben, 1994;Choudhury and DiGirolamo, 1995;Czajkowski et al., 2002), if the radiative transfer equation is linearized, the PW is proportional to the difference between the two AVHRR thermal channels brightness temperature. This split-window technique is not a universal solution to PW retrieval due to the dependence of the difference T on air temperature and surface characteristics. Linear relationships characterized by different coefficients are expected to hold in different locations. The coefficients of the linear relationship between PW and T are site specific and can be derived as output from a radiative transfer model or as observational data at each location.
Let T 4 and T 5 denote the brightness temperatures of AVHRR channels 4 and 5, respectively. Analysis of the split window temperature differences T = T 4 −T 5 , for land surfaces in several different climatic regions with differing soil compositions and vegetative covers, revealed a linear relationship between T and PW (Eck and Holben, 1994;Chrysoulakis et al., 2008). The accuracy of the split window technique is uncertain in arid environment during daytime due to strong near-surface air temperature gradient (Choudhury and Di Girolano, 1995). Provided that the channels 4 and 5 radiances are accurately calibrated, the AVHRR split window technique of PW estimation has the potential to provide consistent estimates of PW. This technique can be referenced to a given time period of independent data or analyzed PW fields and thus provide consistent estimates for each location (Chrysoulakis et al., 2008). As discussed by , spatial and temporal changes of T and PW can also be associated with the changes of land surface characteristics and hence the emissivity. In this paper, the MODIS precipitable water Level 2 product is related to the AVHRR derived T . All spatial calculations were performed in a grid of 5 × 5 km cells over the study area. Data examination confirmed a linear relationship between adjusted PW and T . The PW -T relationship was applied to the AVHRR data to provide a PW spatial distribution. The distribution was validated against radiosonde derived PW time series at one synoptic station and against independent MODIS PW measurements. AVHRR data acquired during July 2006 were used. In a previous study, Chrysoulakis et al. (2008) have calculated a T -PW relationship, corresponding to a different geographic area. AVHRR PW values were estimated using the relationship derived from their study for comparison. 4 were spatially aggregated in 5 x 5 km cells.
In particular, ΔΤ was estimated at AVHRR pixel level and the mean value in each grid cell was computed. Likewise, the PW values extracted at 1.1 x 1.1 km level from the MODIS Level 2 product were aggregated in 5 x 5 km cells. To calibrate the PW -ΔT relationship, 16 MODIS and AVHRR acquisitions were used. Fig. 1. The study area. The MODIS derived PW and the AVHRR derived T were spatially averaged in 5 × 5 km. For this reason a grid of 5 × 5 km cells was developed (Projection: Geographic/WGS84). The location of the Athalassa synoptic station is shown.

The study area
The study area is the island of Cyprus. The area is shown in Fig. 1 and a grid comprising 5 × 5 km cells is superimposed. The PW -T relationship was calibrated at 5 × 5 km cell level; for this reason, the 1 × 1 km AVHRR and MODIS derived data were spatially aggregated in 5 × 5 km cells.
In particular, T was estimated at AVHRR pixel level and the mean value in each grid cell was computed. Likewise, the PW values extracted at 1.1 × 1.1 km level from the MODIS Level 2 product were aggregated in 5 × 5 km cells. To calibrate the PW -T relationship, 16 MODIS and AVHRR acquisitions were used.
The time difference between the NOAA 17/AVHRR or NOAA 18/AVHRR and Aqua/MODIS or Terra/MODIS satellites was found to be between 7 and 37 min, for each one of the selected 16 days. The time lag between the two satellites is small enough, allowing us to consider the satellites overpasses as concurrent. Near nadir, AVHRR acquisitions over the study area were used, therefore, the effect of satellite zenith angle was ignored. Furthermore, based on MODIS Level 2 aerosol product, it was concluded that the atmospheric aerosol content was low for the selected days and thus the effect of aerosols in thermal infrared radiation absorption was ignored.

The MODIS data
The MODIS radiometer has 36 channels that cover the spectral region between 0.4 and 15 µm. channels with decreasing absorption coefficients. The strong absorption channel at 0.936 µm is most useful for dry conditions, while the weak absorption channel at 0.905 mm is most useful for very humid conditions, or low solar elevation . The MODIS Level 2 precipitable water product consists of column water vapor amounts. The results of the nearinfrared retrieval algorithm over land were used in this study; this algorithm  is applied over clear land areas of the globe and above clouds, over both land and ocean, during daytime. Over clear ocean areas, water vapor estimates are provided over the extended glint area. The Level 2 data are generated at the 1-km spatial resolution of MODIS using the near-infrared algorithm during the day. In this study, the Aqua and Terra MODIS PW retrievals for the study area and for 16 cases in July 2006 (Table 1) were selected; the selection criteria are explained in Sect. 3.1.

The AVHRR data
Records of AVHRR onboard NOAA 17 and NOAA 18 acquired over the study area for July 2006 were used to derive the spatio-temporal distribution of PW. The AVHRR is a five channel instrument, with three of the spectral channels located in the visible, near-infrared and mid infrared regions of the spectrum, while the remaining two are located in the thermal infrared with effective wavelengths centered around 10.8 µm (channel 4), and 12 µm (channel 5). It has the same spatial characteristics with previous versions of AVHRR: spatial resolution of 1.1 km at nadir and swath coverage of 2700 km. For the present study, the AVHRR data were available from both NOAA Satellite Active Archive and the HRPT ground receiving station of the Foundation for Research and Technology -Hellas (Chrysoulakis et al., 2008). For the PW -T relationship calibration, data corresponding to the 16 cases displayed in Table 1 were used.

Overview
In this study, PW is estimated using MODIS as an independent dataset. Therefore, the MODIS derived PW can be related to the AVHRR derived T by: where, PW(x,y) is the total column precipitable water in the (x,y) cell of the grid x,y are the coordinates of the center of the cell, a and b are coefficients that need to be estimated for the study area, T = T 4 − T 5 is the brightness temperatures for channels 4 and 5 in the (x,y) cell. The data for AVHRR channels 4 and 5 are given in digital numbers (Level 1b format), which are converted to radiance at the sensor and finally to brightness temperatures to the allowable accuracy of 0.5 K (Chrysoualkis et al., 2008). In this study, the calculation of the coefficients a and b was based on MODIS and AVHRR data acquired during the selected days of June 2006 (Table 1). The selection criteria for these days were: (a) as much cloud free land area as possible and b) near nadir satellite acquisitions over the study area.
The PW -T relationship was applied to the AVHRR dataset to provide PW spatial distributions for a summer month (July 2006) and was validated using radiosonde derived PW values from the synoptic station of Athalassa (35.15 • E 33.40 • N), as well as MODIS PW retrievals.

MODIS derived precipitable water
As noted before, MODIS PW retrievals for July 2006 were used. The data were available as images (PW scenes from collection 5) in Hierarchical Data Format (HDF). A Matlab script was used to read the HDF files, to extract the PW values for each pixel and to assign a flag value to all the cloudy pixels of each PW scene.
Information about cloud cover was included in HDF format. A cloud mask in bit-level representation provided criteria that guide one to decide if a pixel of the PW scene was cloudy or not. Pixels that considered being less than 5% cloudy were used.
Navigation parameters were also included in HDF format; these parameters were used to correct the panoramic distortion and to project each PW scene in a common projection system (Geographic/WGS84). The nearest neighbor resampling method was used to geometrically correct the PW scenes. Subsequently, the MODIS derived PW was transformed to geometrically corrected scenes at 1 km spatial resolution. Real PW values were assigned to non cloudy pixels whereas the cloudy pixels were masked. Finally, all scenes were spatially averaged to the 5 × 5 km cell grid covering the study area. It should be noted that the spatial averaging was performed if all pixels corresponding to each grid cell were non-cloudy. If one or more cloudy pixels were included, the corresponding cell was not taken into account in calculations.

AVHRR derived brightness temperature difference
Standard pre-processing steps were employed to convert the raw AVHRR records to calibrated values. The calibration procedure is based on the conversion of digital numbers of each image to radiance at the sensor values for visible and near-infrared channels (channels 1 and 2) and to brightness temperature values for the thermal infrared channels (channels 3, 4 and 5) with the use of the non-linear conversion equations given in the NOAA KLM Polar Orbiter Data Users Guide (Chrysoulakis et al., 2008). The intercept and gain for each AVHRR channel, as well as the parameters related to the central wave number and the constants needed to parameterize the calibration equations were extracted from the NOAA Level 1b header of each image file. The platform flight altitude was also extracted from the same source and it was used for the geometric correction and the correction of the panoramic distortion in calibrated AVHRR images. The standard NOAA/AVHRR projection system (Geographic/WGS72) was used. The images were finally reprojected to Geographic/WGS84 to be combined with the projected MODIS PW products. The nearest neighborhood resampling method was used in all cases.
To derive the spatial distribution of T at 5 × 5 km level, a cloud mask (Chrysoulakis and Cartalis, 2003) was used to detect cloudy pixels. Following, T was calculated at AVHRR pixel level (1.1 km) for non-cloudy pixels. T values were spatially averaged at the 5 × 5 km level assigning a mean T value in each grid cell. As for the MODIS data, the spatial averaging was performed if all pixels corresponding to each grid cell were non-cloudy. The above steps were repeated for all NOAA 17 and 18 acquisitions of July 2006. Cloudy acquisitions, as well as acquisitions of significant time difference with Terra or Aqua over the study area, were considered worthless and were not included in the analysis.

Cross-calibration, application and validation
The coefficients in Eq. (1) were derived by applying a linear regression procedure. The derived PW -T relationship was applied to produce PW spatial distributions over the study area using AVHRR thermal infrared measurements. As noted previously, T was estimated for cloud free pixels of geometrically corrected scenes acquired over the study area for July 2006. Daytime acquisitions were used. Radiosonde derived PW values at Athalassa synoptic station were used (35.15 • N, 33.40 • E) for validation. The radiosonde data (pressure, temperature and humidity measurements at several levels at 12:00 UTC) for the study period (July 2006) were obtained from the UK Meteorological Office (available online at: http://badc.nerc.ac.uk/data/ radiosglobe). The radiosonde measurements provided a set of time series that was used as an independent dataset. Radiosondes are expected to produce PW with an uncertainty of 0.1-0.2 cm, which is considered to be the accuracy standard of PW for meteorological applications. The PW was estimated using the approach described by Cartalis and Chrysoulakis (1997). Missing records at some heights in the radiosonde data were estimated using a linear regression (Chrysoulakis et al., 2008).
It was assumed that the radiosonde measurements at the synoptic station were representative for the 5 × 5 km cell in which the station was included. Therefore, the radiosonde derived PW corresponding to the station location, was assigned to the respective grid cell. In this way, the radiosonde derived PW and the AVHRR derived PW can be compared at cell level. The Root Mean Squared Error (RMSE) error was calculated and scatterplot diagrams were produced for the visualization of the validation results.
To validate further the PW estimates, these were compared to a MODIS PW independent dataset for July 2006. To measure how close PW estimates were to MODIS PW values, the RMSE was calculated for the cloud free cells. Validation was performed for 20 independent cases in July 2006. Moreover, to visualize biases in estimation, a scatterplot diagram of the estimated PW values versus the MODIS PW values was used. in the present study gave more accurate results. This confirms the assertion that PW -ΔΤ relationship is site-specific and depends on land cover and emissivity.

Fig. 3. Scatterplots of radiosonde derived PW versus AVHRR derived PW, estimated
for cloud free cells corresponding to the synoptic station of Athalassa, a) using the relationship formulated for the study area and b) using the relationship formulated for the broader area of Greece. Additional comparisons were performed to compare the results to a previous study referring to a different geographic area. Chrysoulakis et al. (2008) have estimated AVHRR T -PW linear relationship coefficients for the broader area of Greece. These coefficients were also estimated for a 5 × 5 km cell grid. Cyprus is geographically close to Greece and has similar topographic characteristics. The equation calculated for the area of Greece should result in relatively good PW estimations for the study area. To test this assumption, PW values were estimated for June 2006, using the relationship of Chrysoulakis et al. (2008).
It is worth mentioning the attempt to perform all the calculations in this study at a lower scale: the methodology described above was applied to the initial available data but all calculations were performed in a grid of 3 × 3 km cells; however, the PW -T values referring to 3 × 3 km cells were not well correlated, so the 5 × 5 km grid was finally selected.

Results and discussion
The PW -T relationship was calibrated using cloud free cells from the 16 cases in Table 1. A significant correlation was identified between T and PW (Pearson r = 0.83, p < 0.01). The regression sample consisted of 587 PW -T pairs. A simple linear regression method was used for the derivation of the coefficients in Eq. (1). The results of the regression analysis are shown in Table 2. The t-statistic together with the corresponding p-values indicate that all estimates are statistically significant. A PW -T scatterplot depicting the regression analysis is shown in Fig. 2. The derived PW -T relationship was applied to the dataset for July 2006 and PW distributions derived from AVHRR records. AVHRR acquisitions were chosen to be between 10:00 and 12:00 UTC, because the radiosonde measurements were taken at 11:00 UTC. The resulted AVHRR derived PW values were validated using radiosonde measurements for the Athalassa synoptic station. As mentioned before, a number of radiosonde measurements at some heights were missing and linear interpolation methods were used to fill in the missing values. This may have resulted in not so accurate radiosonde derived PW values.
A relatively good agreement between radiosonde and AVHRR derived PW values was observed (RMSE = 0.76 cm). Figure 3a shows the radiosonde derived PW values versus the AVHRR derived PW values. AVHRR derived PW values seem to be overestimated, but this is the result of underestimated radiosonde PW values. This underestimation may have been resulted by the interpolated radiosonde measurements at some heights, as it has been explained above. The RMSE between the radiosonde derived and MODIS derived PW values found to be 0.75 cm.  As mentioned before, additional comparisons were performed to compare the results to the previous study (Chrysoulakis et al., 2008). AVHRR PW was estimated for the study area for July 2006 using the relationship proposed by Chrysoulakis et al. (2008). In Fig. 3b the radiosonde derived PW values versus the AVHRR derived PW values are shown. Comparing Fig. 3a and b, it is observed that the relationship formulated in the present study gave more accurate results. This confirms the assertion that PW -T relationship is site-specific and depends on land cover and emissivity.
Additionally, AVHRR PW estimates, were compared to the MODIS PW independent dataset for July 2006. A good agreement was observed between the AVHRR and the MODIS derived PW value for 20 independent cases in July 2006 (RMSE = 0.71 cm). In Fig. 4, MODIS derived PW values versus the AVHRR derived PW values are shown in a scatterplot of for the 20 independent cases in July 2006. Moreover, Fig. 5 shows the scatterplot of the MODIS derived PW values versus the AVHRR PW values derived using the relationship formulated for the broader area of Greece (Chrysoulakis et al., 2008) for the same 20 independent cases in July 2006.
Finally, in Fig. 6 the spatial distribution of AVHRR derived PW is shown. Mean PW values were estimated by averaging the daily PW values for the whole month (July 2006), for cloud free cells. The effect of topography is obvious. Lower PW values noted at Tilliria, Marathasa and Troodos mountains, (mean PW min = 1.61 cm), while higher PW values can be observed in north-eastern area of the island (mean PW max = 3.07 cm).

Conclusions
The relationship between PW and T was calibrated by simple regression using the MODIS derived PW and validated using radiosonde measurements and independent MODIS derived PW values. The application area was the island of Cyprus and all spatial calculations were performed for a 5 × 5 km cell grid. The PW -T relationship was calibrated using cloud free cells from 16 cases in July 2006 (Table 1). This equation was applied to daily T data for July 2006. The resulted PW values were validated using radiosonde measurements and independent MODIS derived PW values. A good agreement between radiosonde and AVHRR derived PW values was observed, as well as a good agreement between MODIS derived PW values and AVHRR derived PW values (RMSE = 0.76 cm). The time lag between MODIS and AVHRR observations was low enough to consider the satellite overpasses concurrent, however the above RMSE value would be smaller if this time lag was lower. This RMSE value can be therefore improved if the above time difference is compensated, using PW rates of change, according to Chrysoulakis et al. (2008). A correspondence of the spatial patterns of MODIS derived PW and AVHRR derived T was observed. The effect of topography was evident since lower PW and T values were observed over the mountainous area of Troodos (mean PW min = 1.61 cm) and higher PW values over the area of Famagusta (mean PW max = 3.07 cm).