Biases and Abrupt Shifts of Monthly Precipitable Water from Terra MODIS

: Monthly atmospheric precipitable water (PW) from the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Terra satellite was assessed over land at 60 ◦ S–60 ◦ N. MODIS provides two PW products by using infrared (IR) and near-IR (NIR) algorithms, respectively. An assessment was performed for both MODIS PW data from 2000 to 2014, comparing them with the measurements at international stations of the global positioning systems and with a reanalysis to detect abrupt changes through monthly variations. It is noted that MODIS IR systematically underestimated PW in over 75% of stations, and that PW estimation declines with time. MODIS NIR signiﬁcantly overestimated PW for tropical land and experienced two abrupt shifts. These data defects result in large spurious decreasing trends in MODIS IR and increasing trends in MODIS NIR. The two MODIS PW products are currently not suitable for a climatic-trend analysis, highlighting the need for data reprocessing and calibration.


Introduction
Atmospheric precipitable water (PW), a measure of total column water vapor, plays an important role in climate change through water-vapor feedback [1,2]. To investigate the feedback, understanding the trends and variability in PW is fundamental, which is usually gained based on monthly data series [3][4][5][6][7]. Therefore, the quality of monthly PW is vital for ascertaining trends and variability.
Since March 2000, the Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Terra satellite has been operational to provide two PW products on an instantaneous, daily, eight-day, and monthly mean basis. The first PW dataset (MOD07_L2, Collection 6, SDS: Water_Vapor, [8]) was the integration of moisture profiles that were estimated using MODIS infrared (IR) observations by an empirical multiple regression [9]. The regression coefficients were generated from simulated MODIS radiances and historical moisture radiosonde profiles [10].
The second PW dataset (MOD05_L2, Collection 6, SDS: Water_Vapor_Near_Infrared, [11]) was derived using MODIS near-IR (NIR) observations of water-vapor attenuation under clear conditions or above clouds. The ratios of the water-vapor absorption channels with atmospheric window channels were used to quantify total water-vapor transmittances, which are used to theoretically calculate PW [12]. Since the method relies on NIR observations, derivations were made only over reflective surfaces during the daytime, resulting in an inevitable lack of global data ( Figure 1).
Instantaneous MODIS PW data were evaluated against ground measurements from Global Positioning Systems (GPS) and radiosondes over many regions [13][14][15][16][17]. Those evaluations indicate that instantaneous MODIS PW correlates well with ground measurements, but a good correlation does Remote Sens. 2019, 11, 1315 2 of 8 not mean that the monthly variations of PW can be captured by the MODIS data due to the wet or dry bias. Although it is shown that, under clear conditions, the PW data of MODIS NIR agree well with the monthly variations of GPS-measured PW over the Tibet Plateau [18], the monthly PW data of MODIS IR and NIR have not been assessed well globally.
In this paper, the monthly PW data of MODIS IR and NIR are evaluated by exploiting the PW measurements from GPS stations of the International Global Navigation Satellite Systems Service (IGS), and assessing their monthly variations with the comparison of reanalyzed PW from the Modern-Era Retrospective analysis for Research and Applications (MERRA, Version 2) for global land.

Global Distribution of MODIS PW
The mean PW of MODIS IR, MODIS NIR, and their differences for 2000-2014 is illustrated in Figure 1

Global Distribution of MODIS PW
The mean PW of MODIS IR, MODIS NIR, and their differences for 2000-2014 is illustrated in Figure 1 for four seasons (MAM, JJA, SON, and DJF), respectively. Here, PW derived by MODIS near-infrared observations under clear conditions is referred to as MODIS NIR. Ocean surfaces are dark and reflect very little light that MODIS NIR cannot use to measure PW. Over clear ocean areas, the PW data of MODIS NIR were estimated only over the extended areas with sun glint. Therefore, there are no PW retrievals of MODIS NIR over large areas of the ocean. Both types of MODIS PW data show climatological convergence zones and the monsoon pattern of water-vapor transport in the Indian Ocean and South China Sea. Large systematic differences were revealed, especially for the tropical belt (20°S-20°N), with a higher PW in MODIS IR from 120°E to 180°E, but a lower PW in Brazil, central Africa, and the Pacific and Indian Oceans. 180 • E, but a lower PW in Brazil, central Africa, and the Pacific and Indian Oceans. MODIS IR shows a global coverage of PW but generally a lower PW over tropical lands. The PW of MODIS NIR covers most global land areas and tropical oceans, while the lack of data significantly varies for the Arctic region and the Antarctic continent. Therefore, the study area was designed over land at 60 • S-60 • N, where neither MODIS IR nor NIR had missing PW data, for a better comparison with the IGS and MERRA data.

Comparison of MODIS PW with IGS PW
The IGS PW data were produced by the National Center for Atmospheric Research Earth Observing Laboratory from the zenith path delay (ZPD) derived from GPS measurements, and are available every two hours. The used ZPD data are obtained by combining the solutions from seven IGS data analysis centers. Historically, four centers use an elevation cutoff angle of 15 • , while the others use cutoff angles of 7 • and 10 • [19]. Usually, the combined ZPD is the piecewise continuous weighted mean of the ZPD values from all available IGS data analysis centers. More than 70% of the stations report ZPD solutions from three or more analysis centers. The quality of the combined ZPD is at the level of 4-6 mm, corresponding to ≤1 mm in PW [19,20]. Since no cut-off elevation angles are reported in the IGS PW data, the impact of different cut-off angles on the comparison of MODIS PW with IGS PW is unclear.
A previous validation [21] compared IGS PW at 98 collocated GPS and radiosonde stations around the globe, showing a mean difference of 1.08 mm (drier in the radiosonde data) and a standard deviation of the differences of 2.68 mm. The comparisons of the IGS PW with the global radiosonde and microwave radiometer data showed a good agreement with no systematic errors [22]. There are 127 IGS stations located at 60 • S-60 • N. The monthly IGS PW at each station is calculated by averaging the two-hourly GPS PW data in a calendar month. The spatial collocation is to take the MODIS grid whose center is the closest to the ground-based GPS station. The linear distance between the GPS station and MODIS grid center was constrained to less than 15 km. Since PW values are largely dependent on elevation, the IGS PW should be corrected to the elevation of the MODIS grids before the comparison, in order to verify each other. Since MODIS does not have information on the grid elevation, it was substituted with the surface elevation from the digital elevation model of Shuttle Radar Topography Mission [23]. The equation for the elevation correction is as follows [24]: where PW IGS is the IGS PW data, PW cor represents the corrected IGS PW values after matching to the MODIS grid elevation, ∆h is the elevation difference between the IGS station and MODIS grid, and c is a constant parameter of 0.439, given by Leckner [25]. Figure 2 shows the comparison of IGS PW with MODIS IR and MODIS NIR, respectively. The metrics used to evaluate the comparison were the correlation coefficient (R), mean bias (between MODIS PW and IGS PW), and root mean square error (RMSE) and its normalized value (NRMSE).
Both MODIS PW types had a high correlation with IGS PW, with an R of 0.95 and 0.96, respectively. However, Figure 2a clearly indicates an overall PW underestimation of MODIS IR, with a negative bias of −1.9 mm and an NRMSE of 21.2%. Compared to MODIS IR, the PW data of MODIS NIR in Figure 2b had a better accuracy, with a slight bias of −0.1 mm and an NRMSE of 18.6%. The slope value of 1.1 for the fitted line indicates that MODIS NIR would overestimate when the PW values are high. It is interesting to find that the two independent monthly PW datasets (MODIS NIR and IGS) showed a good agreement.  Figure 3 shows the comparisons at each IGS station for the NRMSE and bias. The details over the regions of Europe and North America can be referred to in Figure S1. The dark-red points on the Middle East and Africa in Figure 3a,b demonstrate that the PW accuracies over those regions are questionable for both types of MODIS data, and MODIS IR underestimated PW at more than 75% of the stations (Figure 3c). The biases between MODIS and IGS PW, along with the latitudinal bins ( Figure S2), show that large negative biases (underestimation) of MODIS IR and large positive biases (overestimation) of MODIS NIR were identified in the land at 30°S-30°N, e.g., Saudi Arabia, southern Africa and Brazil, which corresponds to incompatible and unreasonable drying ( Figure  S4a) and wetting ( Figure S4b) trends in these regions. The monthly differences between MODIS IR and IGS PW, averaged over all GPS stations for 2000-2014, indicate a decreasing tendency ( Figure  S3a) and increasing tendency for the differences between MODIS NIR and IGS PW ( Figure S3b). However, these tendencies cannot help in reaching a solid conclusion. The problem is that there are not many consecutive long-term IGS PW observations available to validate the time series, while they are suitable for diagnosing errors, such as, for instance, the large dry bias for MODIS IR and wet bias for MODIS NIR of the POVE station in Brazil ( Figure S3c).  Figure 3 shows the comparisons at each IGS station for the NRMSE and bias. The details over the regions of Europe and North America can be referred to in Figure S1. The dark-red points on the Middle East and Africa in Figure 3a,b demonstrate that the PW accuracies over those regions are questionable for both types of MODIS data, and MODIS IR underestimated PW at more than 75% of the stations (Figure 3c). The biases between MODIS and IGS PW, along with the latitudinal bins ( Figure S2), show that large negative biases (underestimation) of MODIS IR and large positive biases (overestimation) of MODIS NIR were identified in the land at 30 • S-30 • N, e.g., Saudi Arabia, southern Africa and Brazil, which corresponds to incompatible and unreasonable drying ( Figure S4a) and wetting ( Figure S4b) trends in these regions. The monthly differences between MODIS IR and IGS PW, averaged over all GPS stations for 2000-2014, indicate a decreasing tendency ( Figure S3a) and increasing tendency for the differences between MODIS NIR and IGS PW ( Figure S3b). However, these tendencies cannot help in reaching a solid conclusion. The problem is that there are not many consecutive long-term IGS PW observations available to validate the time series, while they are suitable for diagnosing errors, such as, for instance, the large dry bias for MODIS IR and wet bias for MODIS NIR of the POVE station in Brazil ( Figure S3c).  Figure 3 shows the comparisons at each IGS station for the NRMSE and bias. The details over the regions of Europe and North America can be referred to in Figure S1. The dark-red points on the Middle East and Africa in Figure 3a,b demonstrate that the PW accuracies over those regions are questionable for both types of MODIS data, and MODIS IR underestimated PW at more than 75% of the stations (Figure 3c). The biases between MODIS and IGS PW, along with the latitudinal bins ( Figure S2), show that large negative biases (underestimation) of MODIS IR and large positive biases (overestimation) of MODIS NIR were identified in the land at 30°S-30°N, e.g., Saudi Arabia, southern Africa and Brazil, which corresponds to incompatible and unreasonable drying ( Figure  S4a) and wetting ( Figure S4b) trends in these regions. The monthly differences between MODIS IR and IGS PW, averaged over all GPS stations for 2000-2014, indicate a decreasing tendency ( Figure  S3a) and increasing tendency for the differences between MODIS NIR and IGS PW ( Figure S3b). However, these tendencies cannot help in reaching a solid conclusion. The problem is that there are not many consecutive long-term IGS PW observations available to validate the time series, while they are suitable for diagnosing errors, such as, for instance, the large dry bias for MODIS IR and wet bias for MODIS NIR of the POVE station in Brazil ( Figure S3c).

Comparison of MODIS PW with MERRA PW
The analyzed PW fields of MERRA were processed by an incremental analysis update embodying the influences of observations as well as the biases in model physics [26][27][28], which enables the PW fields to smoothly evolve in time instead of evolving with jumps at times of analysis. Therefore, the monthly series of MERRA PW is well adapted for identifying abrupt changes at the times of the MODIS PW data. The monthly MERRA PW has a spatial resolution of 0.5 • latitude and 0.625 • longitude, which was interpolated to a 1 • latitude/longitude grid to match the spatial resolution of MODIS PW. The interpolated monthly PW data of MERRA were compared with MODIS IR and MODIS NIR for the land at 60 • S-60 • N (Figure 4).

Comparison of MODIS PW with MERRA PW
The analyzed PW fields of MERRA were processed by an incremental analysis update embodying the influences of observations as well as the biases in model physics [26,27,28], which enables the PW fields to smoothly evolve in time instead of evolving with jumps at times of analysis. Therefore, the monthly series of MERRA PW is well adapted for identifying abrupt changes at the times of the MODIS PW data. The monthly MERRA PW has a spatial resolution of 0.5° latitude and 0.625° longitude, which was interpolated to a 1° latitude/longitude grid to match the spatial resolution of MODIS PW. The interpolated monthly PW data of MERRA were compared with MODIS IR and MODIS NIR for the land at 60°S-60°N (Figure 4). All the PW series in Figure 4a demonstrate obvious periodic variations. Figure 4a illustrates the PW of MODIS NIR, and MERRA seems comparable. However, the two MODIS PW datasets show opposite changing statuses: MODIS NIR was increasing while MODIS IR was significantly decreasing. Theoretically, a continuous declination in PW at a large scale over land is not possible, because the consistent increase in temperature physically accompanies ever-increasing water vapor according to the Clausius-Clapeyron relation.
The differences between the PW series further detail temporal variations (Figure 4b). Since July 2001, the differences between MODIS NIR and MERRA have integrally been uplifted. The uplifting tendency was enhanced after July 2010, particularly for July, August, and September, which led to a much larger and unreliable PW trend estimation from MODIS NIR (see Figure S4). Compared to July 2000, the differences between MODIS IR and MERRA in July 2001 began to drop, and were gradually enlarged afterwards, associated with a significant decrease. These data errors not only occur for large expanses of land as an integral, but are also exhibited for small regions. The differences of monthly PW between MODIS IR and MODIS NIR with MERRA for the grid box over All the PW series in Figure 4a demonstrate obvious periodic variations. Figure 4a illustrates the PW of MODIS NIR, and MERRA seems comparable. However, the two MODIS PW datasets show opposite changing statuses: MODIS NIR was increasing while MODIS IR was significantly decreasing. Theoretically, a continuous declination in PW at a large scale over land is not possible, because the consistent increase in temperature physically accompanies ever-increasing water vapor according to the Clausius-Clapeyron relation.
The differences between the PW series further detail temporal variations (Figure 4b). Since July 2001, the differences between MODIS NIR and MERRA have integrally been uplifted. The uplifting tendency was enhanced after July 2010, particularly for July, August, and September, which led to a much larger and unreliable PW trend estimation from MODIS NIR (see Figure S4). Compared to July 2000, the differences between MODIS IR and MERRA in July 2001 began to drop, and were gradually enlarged afterwards, associated with a significant decrease. These data errors not only occur for large expanses of land as an integral, but are also exhibited for small regions. The differences of monthly PW between MODIS IR and MODIS NIR with MERRA for the grid box over Brazil (range between 5 • S-5 • N and 70 • W-60 • W) are also illustrated in Figure S5, showing a similar temporal variation as that over the land at 60 • S-60 • N, with a continuous decrease in PW for MODIS IR and two abrupt shifts for MODIS NIR.

Discussion
It is evident that there is a substantial decline in MODIS IR PW, which is not the natural variation but a data defect. Moreover, two critical time points could be seemingly identified from the differences between MODIS NIR and MERRA, July 2001 and July 2010, because an abrupt shift of the differences occurred at the time (Figure 4b).
These abnormal changes in the MODIS PW data directly affect the assessment of the convincing trends. Previous studies have shown that the average trend for the global oceans from 1988 to 2001 was 0.40 ± 0.09 mm/decade [4], while the trend for the land from 1995 to 2011 was 0.26 ± 0.08 mm/decade [29] and the significant trend for the globe was 0.39 ± 0.15 mm/decade for 1996-2006 [5].
The ± values define the 95% confidence intervals for the trends, and all trends were significant with a p-value < 0.05. Although those trends were derived for different time periods, they provide a general judgment of the trend magnitudes of PW. For MODIS PW, positive trends in MODIS NIR of 1.77 ± 0.08 mm/decade for land at 60 • S-60 • N are obviously too large to be reliable; negative trends in MODIS IR of −1.79 ± 0.09 mm/decade for land at 60 • S-60 • N are completely wrong. The spatial pattern of linear trends in PW for MODIS IR, MODIS NIR, and MERRA is shown in Figure S4. As compared with MERRA, there are large spurious decreasing (increasing) trends of MODIS IR (MODIS NIR) around the tropical belt, especially for Brazil, Middle Africa, and South Asia. For the grid box in Brazil ( Figure S5), the linear trend of MERRA PW is 1.03 ± 0.43 mm/decade, while the trend is −4.66 ± 0.71 mm/decade for MODIS IR PW, and 3.08 ± 0.47 mm/decade for MODIS NIR PW.
For MODIS IR, an underestimation of PW during night-time [15,30] could cause the low monthly average. Meanwhile, radiance biases of MODIS IR bands should be checked, as they are not adjusted in the PW retrieval [10]. The advantage of MODIS IR PW is global coverage. Since the PW differences between MODIS IR and MERRA are linearly exaggerated, one possible way to correct MODIS IR is to use MERRA through the differential-linear-calibration model [30].
For MODIS NIR, the errors are greater for the retrievals from the data collected over dark surfaces or under hazy conditions [12]. As to the shifts in July 2001 and July 2010, there are no official documents stating that the algorithm was changed or adjusted in the meantime. It is claimed that systematic trends caused by a calibration degradation were corrected in the visible and near-infrared bands of the Collection 6 MODIS Terra [31]. However, these calibrations do not involve the near-infrared bands of water-vapor absorption (0.905, 0.936, and 0.940 µm) for MODIS NIR PW, and the thermal infrared bands for MODIS IR PW. A comparison between MODIS NIR and IGS shows that the PW data over land areas within the tropical belt are apparently overestimated, with acutely positive trends that are much higher than the other areas ( Figure S4). Great care should be taken in examining tropical regions. To correct the overestimation of MODIS NIR, Diedrich et al. [16] showed that the wet bias of MODIS NIR PW can be eliminated by introducing empirical correction coefficients for the transmittance calculation of the satellite radiance data.

Conclusions
In this paper, the monthly PW data from MODIS IR and MODIS NIR were assessed over the land at 60 • S-60 • N for 2000-2014. First, the MODIS PW data after an elevation correction were compared with the IGS PW data measured from GPS stations. The comparisons indicate that MODIS NIR PW showed a better agreement with IGS PW than with MODIS IR PW, where PW was underestimated at over 75% of the stations. Over tropical land, there were overall large wet (positive) biases for MODIS NIR and dry (negative) biases for MODIS IR.
Second, the monthly variations of MODIS PW were compared with the PW from the MERRA reanalysis data. For land, MODIS IR PW linearly declined over time, and MODIS NIR PW rose by two significant level shifts. These data defects resulted in unreal trends in MODIS PW: a significant decrease for MODIS IR and a significant increase for MODIS NIR. Compared to the MERRA PW, both the PW from MODIS IR and MODIR NIR had an apparent shift (but in opposite directions) at the July 2001 time point. The question of whether this is due to the instrument system should be explored in the future.
The need for reprocessing the MODIS PW data, including checking and calibrating the water-vapor-absorption and thermal-infrared bands, was highlighted by the issues revealed in this paper, especially with regard to the significant decreasing tendency in MODIS IR PW, and the abrupt shifts of increase in MODIS NIR PW after 2001 and 2010, respectively.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/11/11/1315/s1, Figure S1. Same as Figure 3, but zooming into the North America and Europe regions, Figure S2: Biases between MODIS and IGS PW along with latitudinal bins, Figure S3 Figure S5: Differences of monthly PW between MODIS IR and MODIS NIR with MERRA for grid-box area in Figure S4.