Evaluation and Calibration of MODIS Near-Infrared Precipitable Water Vapor over China Using GNSS Observations and ERA-5 Reanalysis Dataset

: Water vapor is one of the most important parameters in climatic studies. MODerate-resolution Imaging Spectroradiometer (MODIS) is a key instrument and can provide spatially continuous precipitable water vapor (PWV) products. This study was focused on the performance evaluation of the MODIS near-infrared PWV product (MOD-NIR-PWV) over China. For a comprehensive assessment of the performance of MOD-NIR-PWV, PWV retrieved from the measurements at the global navigation satellite systems (GNSS) stations (i.e., GNSS-PWV) and the ERA5 reanalysis dataset (ERA-PWV) from 2013 to 2018 were used as the reference. To investigate the suitability of using ERA-PWV as the reference for the evaluation, ERA-PWV was compared to the high-accuracy GNSS-PWV at 246 GNSS stations and PWV retrieved from radiosonde observations (RS-PWV) at 78 radiosonde stations over China. The results showed that the mean bias and mean root-mean-square (RMS) of the differences between ERA-PWV and GNSS-PWV across all the stations were 0.5 and 1.7 mm, respectively, and the mean correlation coefﬁcient of the two datasets was above 0.96. The values were 0.4 and 1.9 mm and 0.97, respectively, for the differences between ERA-PWV and RS-PWV. This suggests the suitability of ERA-PWV as the reference for the evaluation of MOD-NIR-PWV. In addition, MOD-NIR-PWV was compared with both GNSS-PWV and ERA-PWV, and their mean bias and mean RMS were 2.9 and 3.8 mm (compared to GNSS-PWV) and 2.1 and 3.0 mm (compared to ERA-PWV), respectively. The positive bias values and the non-normal distribution of the differences between MOD-NIR-PWV and both reference datasets imply that a considerable systematic overestimation of MOD-NIR-PWV over China may exist. To mitigate the systematic bias, ERA-PWV was utilized as the sample data due to its spatial continuities, and a grid-based calibration model was developed based on the annual and semiannual periodicities in the differences between MOD-NIR-PWV and ERA-PWV at each grid point. After applying the calibration model to correct MOD-NIR-PWV, the calibrated MOD-NIR-PWV was compared with ERA-PWV and GNSS-PWV for precision and accuracy analysis, respectively. The comparison showed that the model could signiﬁcantly improve the precision by 94% and accuracy by 53%, which manifested the effectiveness of the calibration model in improving the performance of MOD-NIR-PWV over China.


Introduction
Atmospheric water vapor plays a key role in the Earth's radiative balance, hydrological process, energy circulation throughout surface evaporation, transportation, condensation, Liu Z. et al. [13] assessed the accuracy of MOD-IR-PWV and MOD-NIR-PWV over a period of several years in the Hong Kong region.Liu H. et al. [54] compared MOD-IR-PWV and MOD-NIR-PWV over China against RS-PWV over 83 radiosonde stations from the IGRA in 2012.Additionally, several comparisons for the performance assessment of MODIS PWV products in the Tibetan Plateau were also conducted [39,51,55,56].However, since the reference data used in the above studies were mostly from a small number of single-point observations, no confirmative conclusions were made on the performance of MODIS PWV products over China.To address this issue, new reference data with good spatial coverage are preferred for the validation of MODIS PWV products.This is the main reason for this study to use both GNSS-PWV and ERA-PWV as references to evaluate the accuracy of MODIS PWV products over China.
Due to some systematic errors contained in MODIS PWV products, it is desirable to perform calibration for correcting or mitigating such errors before the products are used for climatic applications, e.g., the construction of accurate PWV maps based on multi-source data, if a calibration model is available [54,[57][58][59].The accuracy of MODIS PWV products is dependent upon the performance of the calibration model.The model is usually developed based on the differences between MODIS PWV products and selected reference PWV at co-located ground-based sites, e.g., GNSS stations and/or radiosonde stations, which are the so-called sample data of the modeling.Khaniani et al. [60] used the differences between MOD-NIR-PWV and GNSS-PWV at 38 ground-based GNSS stations over Iran, and a fitting function of height was developed for the calibration model, which effectively reduced the error in MOD-NIR-PWV due to the uniform distribution of the 38 GNSS stations and the slight variation of PWV over Iran.In China, a linear fitting model based on the linear relationship of the differences between MOD-NIR-PWV and reference PWV obtained from other technologies over the whole continent of China is the most common calibration model [54,55,58].The main problem of using the model is that it performs differently in different regions; hence, Liu B. et al. [61] divided China into five regions, and five linear regional models were proposed for linear PWV fitting.However, the performance of the regional models was still not satisfactory due to the low coverage of GNSS stations and the characteristics of strong seasonal variations of the PWV, e.g., large differences between MOD-NIR-PWV and reference PWV in summer and small differences in winter.To consider these two factors, in this study, the ERA5 datasets were adopted as the sample data for better spatial coverage, and a harmonic model fitting the seasonal characteristic of the PWV difference at each ERA5 grid point was established for the calibration model of the grid point.A 0.25 • × 0.25 grid-based calibration model was eventually developed over China.
The outline of this paper is as follows: Various PWV products used in this study are introduced in Section 2. Section 3 presents comprehensive comparisons among PWV products obtained from different techniques over China, followed by the establishment of the grid-based calibration model proposed in this study for the mitigation of the systematic bias in MOD-NIR-PWV in Section 4. Summary and conclusions are given in the last section.

Datasets 2.1.1. MODIS
MODIS is a remote scanning spectroradiometer with 36 discrete spectral bands between 0.645 and 14.235 µm.PWV from the surface to the upper troposphere is obtained based on the difference of transmittance between the absorption bands (near 0.905, 0.936, and 0.94 µm) and the non-absorption bands (near 0.865 and 1.24 µm) of water vapor [35].Since MODIS observations are sensitive to clouds in the atmosphere, the cloud-free PWV data in MODIS Level-2 water vapor product from Terra and Aqua (i.e., MOD05 and MYD05) is used in this study with the aid of the MODIS cloud mask product.The MOD05 and MYD05 products consist of two types of PWV data: MOD-NIR-PWV and MOD-IR-PWV.The former, with 1 km spatial resolution, is available during daytime with a typical error level of 5-10%.The latter, with 5 km spatial resolution, is available during both daytime and nighttime, but its accuracy is poorer than the former.Considering the higher spatial resolution and better accuracy, MOD-NIR-PWV during a six-year period from 2013 to 2018 over China, collected from Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC), was used in this study.
It is worth mentioning that different height systems are used for different types of PWV data, e.g., MOD-NIR-PWV refers to the orthometric height, while ERA5 data and radiosonde observations refer to the geopotential, and GNSS observations refer to ellipsoidal heights.In this study, the orthometric height system used in MODIS was treated as the standard.Ellipsoidal heights of GNSS stations, geopotentials of ERA5 atmospheric data and geopotential heights of radiosonde observations were converted into orthometric heights in the following sections.

ERA5
ERA5 is the latest generation of the reanalysis dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF), which provides hourly atmospheric profile data at 37 pressure levels from 1000 hPa to 1 hPa over 0.25 • × 0.25 • latitude-longitude grids.To adapt to the MOD-NIR-PWV height range, the ERA-PWV used in this study is also referred to the content of PWV from the height of the MODIS pixel to the pressure level of 1 hPa.
ERA-PWV can be calculated by integrating atmospheric variables along all pressure levels as defined by [24]: where ρ w is the density of liquid water, g is the gravity acceleration (in unit of m/s 2 ), dP is the increment of pressure (in unit of hPa) between two adjacent pressure levels, subscript i denotes the ith pressure level, N is the number of all pressure levels, and q is the specific humidity, which can be obtained from [62]: T C −29.65 (2) where P v and P s are the partial pressure of water vapor and saturation water vapor pressure (in hPa), respectively, RH is the relative humidity, and T C is the Celsius temperature.
To mitigate the error caused by the height difference between an ERA5 grid point and its corresponding MODIS pixel, the geopotential of each pressure level at the ERA5 grid point was transformed into the orthometric height to unify the two height systems using the following formulas [63,64]: where ϕ is the latitude of the grid point, H, EH, and H gp are orthometric height (in km), ellipsoidal height (in km), and geopotential (in 10 −3 m 2 /s 2 ), respectively, H geoid is the geoid height (in km) obtained from Earth Gravitational Model 2008 (EGM2008), and g ϕ and R ϕ are the geoid gravity (in m/s 2 ) and the curvature radius (in km) of the Earth at the latitude of ϕ, respectively.
Upon the unification of the height systems, atmospheric variables including relative humidity, temperature, and pressure at the height of the MODIS pixel were calculated for the final integration.If the MODIS pixel was above the lowest pressure level, the atmospheric variable values for the MODIS pixel were obtained by simple linear interpolation based on the variable values at the two adjacent pressure levels that contain the pixel.If the MODIS pixel was below the lowest pressure level, different strategies for different variables were used: relative humidity was assigned to the constant value, which was the same as that of the lowest pressure level; the temperature was linearly extrapolated based on the temperature value at the lowest pressure level and a constant temperature lapse rate of 6.5 • C/km; the pressure was calculated from the following hydrostatic and ideal gas equation [65]: where R d = 8.31462 is the ideal gas constant (in mol −1 K −1 ), H M and H L denote the orthometric heights of the MODIS pixel and the lowest pressure level of ERA5, respectively, P M and T M are the pressure (in hPa) and temperature (in K), respectively, at the MODIS pixel, P L and T L are the pressure (in hPa) and temperature (in K), respectively, at the lowest pressure level.After this data preprocess was completed, the atmospheric variable values were used in the integral of Equation ( 1) to obtain the corresponding ERA-PWV from the height of the MODIS pixel to the pressure level of 1 hPa.

Radiosonde
Radiosonde observations over 78 radiosonde stations were used in this study from the IGRA Version 2 dataset.The geographical distribution of these selected 78 radiosonde stations is shown in Figure 1.The dataset provides various parameters with a temporal resolution of twice a day over China, including pressure, temperature, relative humidity, etc.These derived RS-PWV were calculated based on the same integration method as described in Equation ( 1).The height match process was also conducted as described in the previous section.After the process, these RS-PWV were used in the performance evaluation of ERA-PWV.It should be noted that RS-PWV was not used in the performance evaluation of MOD-NIR-PWV because of the sparse distribution of radiosonde stations and the large time difference (more than 2 h) between the MODIS overpass time and the radiosonde launch epoch.

GNSS
GNSS-PWV used in this study were retrieved from hourly zenith tropospheric delay (ZTD) during the period from 2013 to 2018 over 246 primary ground-based GNSS stations in the Crustal Movement Observation Network of China (CMONOC) (at http://www.cgps.ac.cn, (accessed on 20 July 2019).The geographical distribution of these 246 stations in Figure 1 indicates a good spatial coverage over most regions in China, i.e., sufficiently dense stations in the central region, but sparse stations in the northeastern region, southeastern region, and the Tibetan Plateau.
In GNSS data processing, for the estimation of unknown parameters, mainly including position and atmospheric errors, etc., the slant tropospheric delay errors contained in those GNSS observations from all the satellites at the same station are projected or mapped onto the station's zenith direction (using a mapping function) to reduce the number of unknown parameters to be solved for.The estimated tropospheric delay error at each station is the zenith tropospheric delay (ZTD), which can be decomposed into zenith hydrologic delay (ZHD) and zenith wet delay (ZWD): The ZHD can be obtained accurately from a standard model that is a function of the surface meteorological variables of the site.The most commonly used model is the Saastamoinen model as expressed below [66]: where P 0 is the surface pressure (in hPa), h 0 is the height above geoid (in km), and f (ϕ, h 0 ) is the function given by: f (ϕ, h 0 ) = 1 − 0.00266 cos(2ϕ) − 0.0028h 0 (7) After subtracting the ZHD from the ZTD, what remains is the ZWD, which can be converted to PWV using the following equation [29,67]: where ρ w is the density of liquid water, R w is the specific gas constant of water vapor, the two atmospheric refraction constants k 2 = 16.52 k 2 /hPa and k 3 = 3.776 × 10 5 k 2 /hPa, T m is the atmospheric weighted mean temperature, and its definition (by the following integral) as well as approximation (the summation) is: where T is the atmospheric temperature (in Kelvin) along the vertical direction of the site, dz is the incremental step ∆z is the difference between the heights of the two adjacent pressure levels (in meters).Equations ( 6)- (9) indicate that meteorological data are essential in the conversion of GNSS-ZTD into PWV.However, some GNSS stations are neither equipped with meteorological sensors nor co-located with a radiosonde station that can provide such measurements.In this case, the only way is to use assimilated values from the reanalysis dataset from NWP models (e.g., ERA-Interim and ERA5) for the conversion.Considering the consistency of these atmospheric data, ERA5 reanalysis data rather than meteorological measurements were used in this study.After the atmospheric variable values at the position of the GNSS antenna are calculated from the procedure introduced in Section 2.1.2,T m and ZHD can be obtained and then converted into PWV.This GNSS-PWV, denoted by PWV GPS,0 , is then corrected to the height of each MODIS pixel, denoted by PWV GPS,M , using the vertical correction function below: where PWV GPS,0 and PWV GPS,M are in millimeters, EH and H M are the ellipsoidal height and orthometric height of the GNSS antenna and MODIS pixel, respectively, and the constant λ = 0.439 is the water vapor lapse rate given by [68].

Statistical Metrics
Three statistical metrics were used to evaluate the performance of MOD-NIR-PWV in this study.They were bias, root-mean-squared (RMS), and correlation coefficient (r) and their definitions are listed below: In Equations ( 11)-( 13), the subscriptions e and r denote "evaluated" and "reference", respectively, PWV e and PWV r are the mean values of PWV e and PWV r , respectively, of all the samples, and n is the number of PWV samples.

Evaluation Results
In this section, the accuracy of ERA-PWV over China is firstly evaluated using GNSS-PWV and RS-PWV as reference data; then, the performance and the systematic discrepancies of MOD-NIR-PWV relative to both GNSS-PWV and ERA-PWV are evaluated.Finally, a grid-based calibration model is proposed based on the differences between MOD-NIR-PWV and ERA-PWV.

Geographical Distribution of PWV over China
Figure 2 shows the geographical distribution of the mean PWVs calculated from the ERA-PWV in a period from 2013 to 2018 over China.PWV values over China present a relationship with the geographical location: a clear decline from the southeastern coastal region to the northwestern inland region.The variation tendency is consistent with the result in Gui [49] and Zhang [46].
The maximum value was found in the southeastern coastal region with a mean PWV greater than 40 mm, which was mainly due to its moist condition.The mean PWVs over central China varied from 15 to 30 mm.In the southwestern regions, the mean PWVs were affected by the topography, in particular, the Tibetan Plateau with the highest elevation presented the lowest mean PWVs (about 5 mm).In the northwestern inland region with dry conditions, PWV values remained at a low level in the range from 5 to 15 mm.The PWVs in the northeastern regions also presented low values ranging from 10 to 15 mm due to the high latitude and the effect of the winter monsoon.greater than 40 mm, which was mainly due to its moist condition.The mean PWVs over central China varied from 15 to 30 mm.In the southwestern regions, the mean PWVs were affected by the topography, in particular, the Tibetan Plateau with the highest elevation presented the lowest mean PWVs (about 5 mm).In the northwestern inland region with dry conditions, PWV values remained at a low level in the range from 5 to 15 mm.The PWVs in the northeastern regions also presented low values ranging from 10 to 15 mm due to the high latitude and the effect of the winter monsoon.

Evaluation of ERA-PWV 3.2.1. Comparison between GNSS-PWV and ERA-PWV
The hourly GNSS-PWV at the aforementioned 246 GNSS stations (see Figure 1) is in a different form from the hourly ERA-PWV, i.e., single-point and two-dimensional gridded observations.The comparisons between the two types of data needed to be made over the same points.Therefore, ERA-PWV at a given GNSS station was obtained from an inverse distance weighting interpolation of the ERA-PWV values at the four grids surrounding the GNSS station in this study.
The statistical result of the differences between GNSS-PWV and ERA-PWV in the six-year period from 2013 to 2018 over each of the selected 246 GNSS stations is shown in Figure 3.The biases at all GNSS stations were in the range from −2.0 to 3.2 mm with a mean value of 0.5 mm.The positive bias values at 167 stations imply that ERA-PWV were slightly overestimated.The correlation coefficient values were in the range of 0.959-0.997,which were considerably high, meaning that the temporal variation between GNSS-PWV and ERA-PWV was consistent.As for the RMSs, the minimum, maximum, and mean of the RMSs at all stations were 0.8, 4.1, and 1.7 mm, respectively.The largest RMS was at the SCSM station (latitude 29.2 • N, longitude 102.4 • E) in Sichuan province close to the Tibetan Plateau.In addition, the RMS values seemed to be geographical location-dependent: the values gradually decreased from the southeastern coastal region to the northwest inland region.This decreasing trend was related to the nature of the regional climate, i.e., the tropical monsoon climate in the southeastern coastal region, the humid continental climate in the northwestern inland region, and the plateau climate in the Tibetan Plateau.
Figure 4 shows the distribution, represented by the frequency (in percentage) of data falling into each of the predefined ranges, of the differences between two of the three sets of PWV over the same six-year period.It is noted that GNSS-PWV and RS-PWV are based on single-point results, while ERA-PWV is based on the entire gridded points.The total frequency accumulated in the first five ranges, which are negative data, was 37% for the difference between GNSS-PWV and ERA-PWV, and about 86% of the differences concentrate in the range from −2 to 3 mm.The distribution means a slight overestimation of ERA-PWV.A normal distribution shown by all the blue bars implies that the discrepancies between the two sets of PWV were not systematic biases, which also means a good agreement between the two sets of PWV.
the SCSM station (latitude 29.2° N, longitude 102.4°E) in Sichuan province close to the Tibetan Plateau.In addition, the RMS values seemed to be geographical location-dependent: the values gradually decreased from the southeastern coastal region to the northwest inland region.This decreasing trend was related to the nature of the regional climate, i.e., the tropical monsoon climate in the southeastern coastal region, the humid continental climate in the northwestern inland region, and the plateau climate in the Tibetan Plateau. Figure 4 shows the distribution, represented by the frequency (in percentage) of data falling into each of the predefined ranges, of the differences between two of the three sets of PWV over the same six-year period.It is noted that GNSS-PWV and RS-PWV are based on single-point results, while ERA-PWV is based on the entire gridded points.The total frequency accumulated in the first five ranges, which are negative data, was 37% for the difference between GNSS-PWV and ERA-PWV, and about 86% of the differences concentrate in the range from −2 to 3 mm.The distribution means a slight overestimation of ERA-PWV.A normal distribution shown by all the blue bars implies that the discrepancies between the two sets of PWV were not systematic biases, which also means a good agreement between the two sets of PWV.

Comparison between RS-PWV and ERA-PWV
The RS-PWV at 78 radiosonde stations in the six-year period from 2013 to 2018 was used as another reference to evaluate the performance of the ERA-PWV over China.The inverse distance weighting interpolation was also conducted to release the non-co-location problem in the spatial domain between the single-point RS-PWV and two-dimensional ERA-PWV.
The statistical result of the differences between GNSS-PWV and ERA-PWV over each radiosonde station is shown in Figure 5.The minimum, maximum, and mean values of the biases at these radiosonde stations were −3.0, 2.2 and 0.4 mm, respectively.The RMSs varied from 0.5 to 4.2 mm with a mean value of 1.9 mm.The biases and RMSs were similar with those referring to GNSS-PWV, indicating the equivalent accuracy between RS-PWV and GNSS-PWV over China.The correlation coefficients between RS-PWV and ERA-PWV were in the range from 0.897 to 0.990, and the mean correlation coefficient was 0.970.The extremely low correlation coefficient values were distributed in northeastern regions of the Tibetan Plateau, which may be attributed to the poor reliability of ERA-PWV over complex terrain.In addition, about 34% of the differences between ERA-PWV and RS-PWV were accumulated in the negative ranges (see Figure 4).The apparent normal distribution of the differences in Figure 4 indicated no obvious systematic discrepancy between ERA and RS.
These results referring to GNSS-PWV and RS-PWV agree well with each other, implying the perfect accuracy of ERA-PWV over China.Therefore, ERA-PWV can be used as a new reference for the evaluation of MOD-NIR-PWV due to its better spatial continuity and coverage.

Comparison between RS-PWV and ERA-PWV
The RS-PWV at 78 radiosonde stations in the six-year period from 2013 to 2018 was used as another reference to evaluate the performance of the ERA-PWV over China.The inverse distance weighting interpolation was also conducted to release the non-co-location problem in the spatial domain between the single-point RS-PWV and two-dimensional ERA-PWV.
The statistical result of the differences between GNSS-PWV and ERA-PWV over each radiosonde station is shown in Figure 5.The minimum, maximum, and mean values of the biases at these radiosonde stations were −3.0, 2.2 and 0.4 mm, respectively.The RMSs varied from 0.5 to 4.2 mm with a mean value of 1.9 mm.The biases and RMSs were similar with those referring to GNSS-PWV, indicating the equivalent accuracy between RS-PWV and GNSS-PWV over China.The correlation coefficients between RS-PWV and ERA-PWV were in the range from 0.897 to 0.990, and the mean correlation coefficient was 0.970.The extremely low correlation coefficient values were distributed in northeastern regions of the Tibetan Plateau, which may be attributed to the poor reliability of ERA-PWV over complex terrain.In addition, about 34% of the differences between ERA-PWV and RS-PWV were accumulated in the negative ranges (see Figure 4).The apparent normal distribution of the differences in Figure 4 indicated no obvious systematic discrepancy between ERA and RS.

Comparison between RS-PWV and ERA-PWV
The RS-PWV at 78 radiosonde stations in the six-year period from 2013 to 2018 was used as another reference to evaluate the performance of the ERA-PWV over China.The inverse distance weighting interpolation was also conducted to release the non-co-location problem in the spatial domain between the single-point RS-PWV and two-dimensional ERA-PWV.
The statistical result of the differences between GNSS-PWV and ERA-PWV over each radiosonde station is shown in Figure 5.The minimum, maximum, and mean values of the biases at these radiosonde stations were −3.0, 2.2 and 0.4 mm, respectively.The RMSs varied from 0.5 to 4.2 mm with a mean value of 1.9 mm.The biases and RMSs were similar with those referring to GNSS-PWV, indicating the equivalent accuracy between RS-PWV and GNSS-PWV over China.The correlation coefficients between RS-PWV and ERA-PWV were in the range from 0.897 to 0.990, and the mean correlation coefficient was 0.970.The extremely low correlation coefficient values were distributed in northeastern regions of the Tibetan Plateau, which may be attributed to the poor reliability of ERA-PWV over complex terrain.In addition, about 34% of the differences between ERA-PWV and RS-PWV were accumulated in the negative ranges (see Figure 4).The apparent normal distribution of the differences in Figure 4 indicated no obvious systematic discrepancy between ERA and RS.
These results referring to GNSS-PWV and RS-PWV agree well with each other, implying the perfect accuracy of ERA-PWV over China.Therefore, ERA-PWV can be used as a new reference for the evaluation of MOD-NIR-PWV due to its better spatial continuity and coverage.These results referring to GNSS-PWV and RS-PWV agree well with each other, implying the perfect accuracy of ERA-PWV over China.Therefore, ERA-PWV can be used as a new reference for the evaluation of MOD-NIR-PWV due to its better spatial continuity and coverage.Since MOD-NIR-PWV and GNSS-PWV are not spatially co-located, a spatial interpolation process is needed to make the former co-located with the latter.Both sets of PWV are also not simultaneous observations in the temporal domain, i.e., the overpass time of the MODIS over a GNSS station differs from the epoch of any hourly GNSS-PWV, thus, preprocessing for the selection of the GNSS-PWV that adapts to the overpass time of the MODIS needs to be performed.In this study, the mean of the GNSS-PWV at four consecutive epochs, two of which were immediately before the MODIS overpass and the other two were immediately after the overpass, was used in the later comparison for addressing the simultaneity problem.In the spatial domain, a vertical correction for GNSS-PWV and a horizontal interpolation for MOD-NIR-PWV were conducted to reduce the effects caused by the height difference and non-co-location, respectively, between the GNSS station and MODIS pixel.More specifically, for a given GNSS station, the optimal height of the corresponding MODIS image was determined from the inverse distance interpolation of the heights of nine MODIS pixels in a 3 × 3 pixel window centered at the GNSS station.Then, to match the optimal MODIS height, the GNSS-PWV over the station obtained from the mean value of consecutive four epochs was corrected based on Equation (10) and named as corrected GNSS-PWV.For MOD-NIR-PWV, i.e., in the horizontal domain, the inverse distance weighted interpolation of MOD-NIR-PWV data of the same 3 × 3 pixel window was performed for its final result, which was compared to the above corrected GNSS-PWV.It is noted that 25 GNSS stations were discarded because their MOD-NIR-PWV time series were less than 20% of the six-year period studied.The statistical metrics for the differences between the two sets of PWV over each of the 221 GNSS stations are shown in Figure 6.
lation process is needed to make the former co-located with the latter.Both sets of PWV are also not simultaneous observations in the temporal domain, i.e., the overpass time of the MODIS over a GNSS station differs from the epoch of any hourly GNSS-PWV, thus, preprocessing for the selection of the GNSS-PWV that adapts to the overpass time of the MODIS needs to be performed.In this study, the mean of the GNSS-PWV at four consecutive epochs, two of which were immediately before the MODIS overpass and the other two were immediately after the overpass, was used in the later comparison for addressing the simultaneity problem.In the spatial domain, a vertical correction for GNSS-PWV and a horizontal interpolation for MOD-NIR-PWV were conducted to reduce the effects caused by the height difference and non-co-location, respectively, between the GNSS station and MODIS pixel.More specifically, for a given GNSS station, the optimal height of the corresponding MODIS image was determined from the inverse distance interpolation of the heights of nine MODIS pixels in a 3 × 3 pixel window centered at the GNSS station.Then, to match the optimal MODIS height, the GNSS-PWV over the station obtained from the mean value of consecutive four epochs was corrected based on Equation (10) and named as corrected GNSS-PWV.For MOD-NIR-PWV, i.e., in the horizontal domain, the inverse distance weighted interpolation of MOD-NIR-PWV data of the same 3 × 3 pixel window was performed for its final result, which was compared to the above corrected GNSS-PWV.It is noted that 25 GNSS stations were discarded because their MOD-NIR-PWV time series were less than 20% of the six-year period studied.The statistical metrics for the differences between the two sets of PWV over each of the 221 GNSS stations are shown in Figure 6.
It can be seen that both the biases and RMSs were considerably large, while the correlation coefficient values were small, especially in the southeastern coastal region.Overall, the bias of each GNSS station varied from −1.2 to 8.3 mm with a mean value of 2.9 mm, and biases at 215 stations were positive.The RMSs were in the range from 0.9 to 9.4 mm with a mean value of 3.8 mm.The biases and RMSs indicate the poor performance of MOD-NIR-PWV over China.Moreover, some relatively higher biases (>5 mm) and RMSs (>6 mm) were observed in the southeastern coastal region, which suggests the much poorer performance of MOD-NIR-PWV.A similar decreasing trend of the biases and RMSs from the southeast to the northwest was also noticeable.Additionally, the minimum, maximum, and mean of the correlation coefficient values were 0.954, 0.998, and 0.991, respectively, implying similar temporal variation trends between the two sets of PWV.This indicates the potential to improve the precision of MOD-NIR-PWV throughout some calibrations.It can be seen that both the biases and RMSs were considerably large, while the correlation coefficient values were small, especially in the southeastern coastal region.Overall, the bias of each GNSS station varied from −1.2 to 8.3 mm with a mean value of 2.9 mm, and biases at 215 stations were positive.The RMSs were in the range from 0.9 to 9.4 mm with a mean value of 3.8 mm.The biases and RMSs indicate the poor performance of MOD-NIR-PWV over China.Moreover, some relatively higher biases (>5 mm) and RMSs (>6 mm) were observed in the southeastern coastal region, which suggests the much poorer performance of MOD-NIR-PWV.A similar decreasing trend of the biases and RMSs from the southeast to the northwest was also noticeable.Additionally, the minimum, maximum, and mean of the correlation coefficient values were 0.954, 0.998, and 0.991, respectively, implying similar temporal variation trends between the two sets of PWV.This indicates the potential to improve the precision of MOD-NIR-PWV throughout some calibrations.
The non-normal distribution shown by the red bars in Figure 7 indicates that the differences between MOD-NIR-PWV and GNSS-PWV were systematic biases.Moreover, the total frequency in the five ranges, i.e., the frequency of negative differences, was less than 11%; about 66% of the differences were concentrated in the range between 0 and 4 mm, implying that the large systematic biases were due to the significant overestimation of MOD-NIR-PWV.This result is consistent with the findings from the comparisons with radiosonde data in Hongkong by Liu Z. et al. [13] and also in mainland China by Liu H. et al. [54].However, a more confirmative conclusion needs to be made based on more comparisons with spatially continuous reference data such as ERA-PWV to evaluate the performance of MOD-NIR-PWV over China.This will be discussed in the next section.
Remote Sens. 2021, 13, x FOR PEER REVIEW 12 of 21 The non-normal distribution shown by the red bars in Figure 7 indicates that the differences between MOD-NIR-PWV and GNSS-PWV were systematic biases.Moreover, the total frequency in the five ranges, i.e., the frequency of negative differences, was less than 11%; about 66% of the differences were concentrated in the range between 0 and 4 mm, implying that the large systematic biases were due to the significant overestimation of MOD-NIR-PWV.This result is consistent with the findings from the comparisons with radiosonde data in Hongkong by Liu Z. et al. [13] and also in mainland China by Liu H. et al. [54].However, a more confirmative conclusion needs to be made based on more comparisons with spatially continuous reference data such as ERA-PWV to evaluate the performance of MOD-NIR-PWV over China.This will be discussed in the next section.

Comparison between ERA-PWV and MOD-NIR-PWV
The high accuracy of ERA-PWV has been validated using GNSS-PWV as the reference over China in Section 3.2, hence, ERA-PWV was used as the reference of MOD-NIR-PWV in this section.The two sets of data also needed to be preprocessed for a valid comparison, i.e., the two types of datasets needed to correspond to the same position and also the same time.In this study, the 3 × 3 pixel window algorithm was adopted to determine the MOD-NIR-PWV values for the position of the ERA-5 grid points in the MODIS image, and the four-consecutive-epoch interpolation was used to solve the problem of the mismatch in the temporal domain, as described in the previous section.Figure 8 shows the statistical results of the comparison between ERA-PWV and MOD-NIR-PWV in the same six-year period studied.Note that a few grid points show white color, meaning lack of values, the reasons for which were (1) MODIS sensors being sensitive to sunlight (only MOD-NIR-PWV in the bright land is available) and ( 2) the lengths of the MOD-NIR-PWV time series at the grid points being not sufficient for a reasonable statistic result.

Comparison between ERA-PWV and MOD-NIR-PWV
The high accuracy of ERA-PWV has been validated using GNSS-PWV as the reference over China in Section 3.2, hence, ERA-PWV was used as the reference of MOD-NIR-PWV in this section.The two sets of data also needed to be preprocessed for a valid comparison, i.e., the two types of datasets needed to correspond to the same position and also the same time.In this study, the 3 × 3 pixel window algorithm was adopted to determine the MOD-NIR-PWV values for the position of the ERA-5 grid points in the MODIS image, and the four-consecutive-epoch interpolation was used to solve the problem of the mismatch in the temporal domain, as described in the previous section.Figure 8 shows the statistical results of the comparison between ERA-PWV and MOD-NIR-PWV in the same six-year period studied.Note that a few grid points show white color, meaning lack of values, the reasons for which were (1) MODIS sensors being sensitive to sunlight (only MOD-NIR-PWV in the bright land is available) and ( 2) the lengths of the MOD-NIR-PWV time series at the grid points being not sufficient for a reasonable statistic result.The non-normal distribution shown by the red bars in Figure 7 indicates that the differences between MOD-NIR-PWV and GNSS-PWV were systematic biases.Moreover, the total frequency in the five ranges, i.e., the frequency of negative differences, was less than 11%; about 66% of the differences were concentrated in the range between 0 and 4 mm, implying that the large systematic biases were due to the significant overestimation of MOD-NIR-PWV.This result is consistent with the findings from the comparisons with radiosonde data in Hongkong by Liu Z. et al. [13] and also in mainland China by Liu H. et al. [54].However, a more confirmative conclusion needs to be made based on more comparisons with spatially continuous reference data such as ERA-PWV to evaluate the performance of MOD-NIR-PWV over China.This will be discussed in the next section.

Comparison between ERA-PWV and MOD-NIR-PWV
The high accuracy of ERA-PWV has been validated using GNSS-PWV as the reference over China in Section 3.2, hence, ERA-PWV was used as the reference of MOD-NIR-PWV in this section.The two sets of data also needed to be preprocessed for a valid comparison, i.e., the two types of datasets needed to correspond to the same position and also the same time.In this study, the 3 × 3 pixel window algorithm was adopted to determine the MOD-NIR-PWV values for the position of the ERA-5 grid points in the MODIS image, and the four-consecutive-epoch interpolation was used to solve the problem of the mismatch in the temporal domain, as described in the previous section.Figure 8 shows the statistical results of the comparison between ERA-PWV and MOD-NIR-PWV in the same six-year period studied.Note that a few grid points show white color, meaning lack of values, the reasons for which were (1) MODIS sensors being sensitive to sunlight (only MOD-NIR-PWV in the bright land is available) and ( 2) the lengths of the MOD-NIR-PWV time series at the grid points being not sufficient for a reasonable statistic result.As shown in Figure 8a, the most positive bias values mean the general overestimation of the MOD-NIR-PWV in comparison with the ERA-PWV.The biases were mostly in the range from −0.6 to 8.2 mm with a mean value of 2.1 mm.Most RMS values (see Figure 8b) were in the range from 0.7 to 10.0 mm with a mean value of 3.0 mm.It is clear that poor results were from the regions under a relatively humid climate condition, such as mid-latitude and southeastern coastal regions, especially in Hainan Province, the island in southeastern China where the maximum biases and RMSs were above 7 and 9 mm, respectively.Moreover, both biases and RMSs had similar variation trends: they increased with increasing amounts of PWV from the southeastern coastal region to the northwestern inland region.The overall correlation coefficient values in Figure 8c varied largely, from 0.616 to 0.998 with a mean of 0.975, which was different from the comparison with GNSS-PWV discussed in the above section.In contrast to the high variation in biases and RMSs, the correlation coefficient values showed slight variation in most regions, but abnormally low values (<0.8) occurred on some grid points in the Tibetan Plateau and southwestern region.
The frequency distribution shown in yellow bars in Figure 7 also indicates the nonnormal feature of the differences between ERA-PWV and MOD-NIR-PWV.The frequency of negative differences (about 17%) was significantly lower than that in all of the other ranges, which accounted for over 83% of all the differences.Moreover, about 76% of the differences were concentrated in the range from −1 to 3 mm.These results confirm the considerable systematic overestimation of MOD-NIR-PWV over China.
The distributions of seasonal biases and RMSs were also compared; see Figures 9 and 10, which show little seasonal differences in the Tibetan Plateau, thus, the Tibetan Plateau was discarded in the following investigation.Apart from the Tibetan Plateau, all the other regions showed distinctive seasonal differences: the lowest biases and RMSs were all found in winter in all regions with means of 0.7 and 1.2 mm, respectively, while the highest values occurred in summer or autumn in different regions.In most of these regions, except for the southwestern region, the highest values were found in the summer, mainly due to the dramatic wet bias of MODIS sensors in monsoon seasons.In summer, a large amount of water vapor is continuously transported into inland regions from the western Pacific Ocean by the East Asian monsoon, leading to an increase in PWV in eastern China.These continuous wet conditions result in relatively high wet bias and bad statistical results.In winter, dry and cool conditions caused by the air-flow pattern of the winter monsoons lead to the low biases in MOD-NIR-PWV.Therefore, relatively better agreement is found in winter in these regions.This comparison confirms that the wet conditions in summer are the main reason for the significant systematic bias between MOD-NIR-PWV and the references (including GNSS-PWV and ERA-PWV).In the southwestern region, the highest value occurred in autumn.This phenomenon may be related to few synoptic stations deployed in these regions, which further affects the accuracy of the ERA5 dataset.

Grid-Based Calibration Modeling for MOD-NIR-PWV
One of the main aims of this study was to mitigate systematic biases in the MODIS water vapor products over China.GNSS-PWV and ERA-PWV are the two types of reference data to calibrate these biases.Several studies used the difference between MOD-NIR-PWV and GNSS-PWV as sample data to achieve it [53,54,57,61].However, the relationship between GNSS PWV and MOD-NIR-PWV in different geographical locations is different, thus, the performance of these calibration models may perform differently in different regions, especially in the regions where only a few GNSS stations are deployed (see Figure 6c).To have consistent calibration performance, the differences between ERA-PWV and MOD-NIR-PWV in the six-year period at each grid point were used as the sample data to develop a grid-based calibration model due to the better spatial coverage and continuity of ERA-PWV in comparison with those of GNSS-PWV.As a result, the difference between

Grid-Based Calibration Modeling for MOD-NIR-PWV
One of the main aims of this study was to mitigate systematic biases in the MODIS water vapor products over China.GNSS-PWV and ERA-PWV are the two types of reference data to calibrate these biases.Several studies used the difference between MOD-NIR-PWV and GNSS-PWV as sample data to achieve it [53,54,57,61].However, the relationship between GNSS PWV and MOD-NIR-PWV in different geographical locations is different, thus, the performance of these calibration models may perform differently in different regions, especially in the regions where only a few GNSS stations are deployed (see Figure 6c).To have consistent calibration performance, the differences between ERA-PWV and MOD-NIR-PWV in the six-year period at each grid point were used as the sample data to develop a grid-based calibration model due to the better spatial coverage and continuity of ERA-PWV in comparison with those of GNSS-PWV.As a result, the difference between

Grid-Based Calibration Modeling for MOD-NIR-PWV
One of the main aims of this study was to mitigate systematic biases in the MODIS water vapor products over China.GNSS-PWV and ERA-PWV are the two types of reference data to calibrate these biases.Several studies used the difference between MOD-NIR-PWV and GNSS-PWV as sample data to achieve it [53,54,57,61].However, the relationship between GNSS PWV and MOD-NIR-PWV in different geographical locations is different, thus, the performance of these calibration models may perform differently in different regions, especially in the regions where only a few GNSS stations are deployed (see Figure 6c).To have consistent calibration performance, the differences between ERA-PWV and MOD-NIR-PWV in the six-year period at each grid point were used as the sample data to develop a grid-based calibration model due to the better spatial coverage and continuity of ERA-PWV in comparison with those of GNSS-PWV.As a result, the difference between the original MOD-NIR-PWV and the calibrated (or corrected) MOD-NIR-PWV was in fact the residual of the calibration model which represents the internal agreement (i.e., precision) of the samples used in the modeling process.In addition, the accuracy of the model was also evaluated by comparing the calibrated MOD-NIR-PWV against the reference of the GNSS-PWV at all the 221 stations in the same six-year period, which were out-of-sample data.This can be considered as an indication of external agreement (i.e., accuracy).The reason for the use of the same six-year period (2013-2018) GNSS-PWV as the test data, rather than the follow-up of one or two years (2019-2020) was that the amount of available MOD-NIR-PWV in the two years was insufficient for a sound statistical result.
Although the seasonal biases and RMSs at all grids have been shown in Figures 9 and 10, for clearer comparisons, the seasonal differences at four randomly selected grid points are extracted and shown in Table 1.The periodic characteristics of the differences in the six-year period studied at the four points were also analyzed using the Lomb-Scargle periodogram method [69], and results are shown in Figure 11.The differences shown in the left column in Figure 11 indicate clear annual and semiannual periodicities at all four grid points, and the periodograms in the right column also illustrate the most notable components in the differences.Based on this characteristic, a harmonic regression model that contained annual and semiannual periodicities and fitted the differences of the two sets of PWV at a given grid point was constructed for the point.The model was: where y is the difference, i.e., MOD-NIR-PWV-ERA-PWV (in millimeters) at the point, t is the epoch (in unit of year), y 0 is the constant intercept, v is the slope, c 1 and s 1 are the sine and cosine coefficients for the annual period, respectively, c 2 and s 2 are the sine and cosine coefficients for the semiannual period, respectively, and ε is the fitting residual of the model.All unknown parameters were estimated from the least-squares method and the sample difference data at the grid point.Different grid points had different fitting models as different sample data were used.These models for all 0.25 • × 0.25 grids formed a grid-based model, which was used to predict the difference value at any target point as a correction of the MOD-NIR-PWV on the point.It should be noted that if the target site was one of the grid points, then the predicted difference (or correction) value could be directly obtained from the model of the grid point.However, for any non-grid point, its predicted correction value was obtained from an interpolation of the model-predicted values at the four grid points surrounding the target point.
To evaluate the performance of the constructed calibration model, all the original MOD-NIR-PWV during the six-year period studied were calibrated at first, then, the differences of the calibrated MOD-NIR-PWV compared to both ERA-PWV and GNSS-PWV were used to measure the precision and accuracy of the model, respectively.Their statistical results are shown in Figures 12 and 13, and the frequency distribution of the differences between the calibrated MOD-NIR-PWV and the two sets of reference data over the same six-year period studied is shown in Figure 14.It is worth mentioning that Figure 12 shows spatially continuous data (except for a few white or blank points), while Figure 13 shows discrete points, because ERA-PWV and GNSS-PWV are from regular grid points and non-regular (i.e., non-grid) points, respectively.
iodogram method [69], and results are shown in Figure 11.The differences shown in the left column in Figure 11 indicate clear annual and semiannual periodicities at all four grid points, and the periodograms in the right column also illustrate the most notable components in the differences.Based on this characteristic, a harmonic regression model that contained annual and semiannual periodicities and fitted the differences of the two sets of PWV at a given grid point was constructed for the point.The model was: where y is the difference, i.e., MOD-NIR-PWV-ERA-PWV (in millimeters) at the point, In Figure 13, the bias values were in the range from −3.4 to 2.2 mm, and the minimum bias of −3.4 mm was further away from zero, compared to the original minimum bias of −1.2 mm, and the frequency accumulated in negative ranges was higher than that between GNSS-PWV and original MOD-NIR-PWV.However, the distribution of the differences was more consistent with a normal distribution (see Figure 14), and most differences were concentrated in the range from −1 to 2 mm, resulting in smaller systematic biases of the calibrated MOD-NIR-PWV.Moreover, the RMSs were decreased by 53% with a smaller range from 0.6 to 4.3 mm, and the correlation coefficient values showed slight variation due to good temporal agreement between the original MOD-NIR-PWV and GNSS-PWV at the GNSS stations.This suggests that the grid-based calibration model, which was based on the differences between ERA-PWV and MOD-NIR-PWV, could effectively improve the accuracy of MODIS water vapor products over China.In Figure 13, the bias values were in the range from −3.4 to 2.2 mm, and the minimum bias of −3.4 mm was further away from zero, compared to the original minimum bias of −1.2 mm, and the frequency accumulated in negative ranges was higher than that between GNSS-PWV and original MOD-NIR-PWV.However, the distribution of the differences was more consistent with a normal distribution (see Figure 14), and most differences were concentrated in the range from −1 to 2 mm, resulting in smaller systematic biases of the calibrated MOD-NIR-PWV.Moreover, the RMSs were decreased by 53% with a smaller range from 0.6 to 4.3 mm, and the correlation coefficient values showed slight variation due to good temporal agreement between the original MOD-NIR-PWV and GNSS-PWV at the GNSS stations.This suggests that the grid-based calibration model, which was based on the differences between ERA-PWV and MOD-NIR-PWV, could effectively improve the accuracy of MODIS water vapor products over China.

Conclusions
In this study, the accuracy of MODIS near-infrared PWV products collected from the Terra and Aqua platform during the six-year period from 2013 to 2018 over China was evaluated by comparing them against the two reference datasets: single-point PWV retrieved from 246 ground-based GNSS tracking stations and the two-dimensional PWV retrieved from ERA5 reanalysis datasets over China.To validate the suitability of using ERA-PWV as the reference for the grid-based calibration model developed for China, ERA-PWV was compared with GNSS-PWV and RS-PWV.The statistical results showed good agreement: the means of biases and RMSs were 0.5 and 1.7 mm compared to GNSS-PWV, 0.4 and 1.9 mm compared to RS-PWV, and their correlation coefficients were above 0.959 and 0.894, respectively.
After the validation, GNSS-PWV and ERA-PWV were used as the reference to evaluate the performance of MOD-NIR-PWV, since most of the large differences between the MODIS overpass time and radiosonde launch time were more than two hours.From the comparison of MOD-NIR-PWV to both GNSS-PWV and ERA-PWV in the six-year period studied, the non-normal distribution of the differences between the MOD-NIR-PWV and the two sets of reference data suggested that there exists a considerable systematic discrepancy in MOD-NIR-PWV over China.More specifically, the mean values of biases and RMSs were 2.9 and 3.8 mm compared to GNSS-PWV, and the values were 2.1 and 3.0 mm compared to ERA-PWV, which implies the overestimation of MOD-NIR-PWV over From Figure 12, the calibrated MOD-NIR-PWV had a better agreement with ERA-PWV than the original MOD-NIR-PWV over China, which was reasonable, since the calibration model at each grid point was based on the differences of ERA-PWV and MOD-NIR-PWV at the point.After the calibration, the minimum, maximum, and mean biases were −0.2, 0.5 and 0.1 mm, respectively; the corresponding values of RMSs were 0.2, 4.6 and 1.0 mm, respectively.The biases and RMSs were reduced by 94% and 71%, respectively, compared to those of the original MOD-NIR-PWV, especially in the southeastern coastal region.The calibration also had significant positive effects on the correlation coefficients in the southwestern region, but little effects in the other regions.In addition, the normal distribution of the yellow bars in Figure 14 also indicates the effectiveness in the reduction of the systematic biases of MOD-NIR-PWV at the GNSS stations.Therefore, the constructed calibration model had good precision.
In Figure 13, the bias values were in the range from −3.4 to 2.2 mm, and the minimum bias of −3.4 mm was further away from zero, compared to the original minimum bias of −1.2 mm, and the frequency accumulated in negative ranges was higher than that between GNSS-PWV and original MOD-NIR-PWV.However, the distribution of the differences was more consistent with a normal distribution (see Figure 14), and most differences were concentrated in the range from −1 to 2 mm, resulting in smaller systematic biases of the calibrated MOD-NIR-PWV.Moreover, the RMSs were decreased by 53% with a smaller range from 0.6 to 4.3 mm, and the correlation coefficient values showed slight variation due to good temporal agreement between the original MOD-NIR-PWV and GNSS-PWV at the GNSS stations.This suggests that the grid-based calibration model, which was based on the differences between ERA-PWV and MOD-NIR-PWV, could effectively improve the accuracy of MODIS water vapor products over China.

Conclusions
In this study, the accuracy of MODIS near-infrared PWV products collected from the Terra and Aqua platform during the six-year period from 2013 to 2018 over China was evaluated by comparing them against the two reference datasets: single-point PWV retrieved from 246 ground-based GNSS tracking stations and the two-dimensional PWV retrieved from ERA5 reanalysis datasets over China.To validate the suitability of using ERA-PWV as the reference for the grid-based calibration model developed for China, ERA-PWV was compared with GNSS-PWV and RS-PWV.The statistical results showed good agreement: the means of biases and RMSs were 0.5 and 1.7 mm compared to GNSS-PWV, 0.4 and 1.9 mm compared to RS-PWV, and their correlation coefficients were above 0.959 and 0.894, respectively.
After the validation, GNSS-PWV and ERA-PWV were used as the reference to evaluate the performance of MOD-NIR-PWV, since most of the large differences between the MODIS overpass time and radiosonde launch time were more than two hours.From the comparison of MOD-NIR-PWV to both GNSS-PWV and ERA-PWV in the six-year period studied, the non-normal distribution of the differences between the MOD-NIR-PWV and the two sets of reference data suggested that there exists a considerable systematic discrepancy in MOD-NIR-PWV over China.More specifically, the mean values of biases and RMSs were 2.9 and 3.8 mm compared to GNSS-PWV, and the values were 2.1 and 3.0 mm compared to ERA-PWV, which implies the overestimation of MOD-NIR-PWV over China.In addition, the high correlation coefficient values between the MOD-NIR-PWV and both sets of reference data in most regions indicated the similar temporal trends of the three sets of PWV.
To mitigate systematic biases in the MODIS water vapor products over China, the differences between ERA-PWV and MOD-NIR-PWV from 2013 to 2018 at each grid point were used as the sample data to develop a 0.25 • × 0.25 grid-based calibration model of MOD-NIR-PWV based on a harmonic model with annual and semiannual periods.This was achieved through a periodic characteristic analysis of the sample data.To evaluate the performance of the developed calibration model, the calibrated MOD-NIR-PWV in the above-mentioned six-year period was compared to ERA-PWV and GNSS-PWV for the evaluation of precision and accuracy, respectively, of the model.The results show that the calibration model could not only significantly improve precision by 94% and accuracy by 53%, but also weaken the non-normal distribution features of the differences between MOD-NIR-PWV and both sets of reference PWV.These results suggest that the developed grid-based calibration model can be used to improve the accuracy of MOD-NIR-PWV over China.
The main shortcoming of this study is that the test data (i.e., GNSS-PWV) used to evaluate the accuracy of the calibration model were from the same period as the sample data used to develop the model (2013-2018), due to the limited amount of available MOD-NIR-PWV in 2019 and 2020, which is desirable for a solid statistical result.Our future work will focus on using test data of multiple years later than those of the sample data (e.g., after 2018) to evaluate the performance of the grid-based calibration model.

21 Figure 1 .
Figure 1.Geographical distribution of 78 radiosonde stations (blue triangle) and 246 ground-based GNSS stations in the CMONOC (orange circle).2.1.4.GNSS GNSS-PWV used in this study were retrieved from hourly zenith tropospheric delay (ZTD) during the period from 2013 to 2018 over 246 primary ground-based GNSS stations

Figure 2 .
Figure 2. Mean PWV (millimeters) calculated from the ERA-PWV in the six-year period from 2013 to 2018 over China.

Figure 2 .
Figure 2. Mean PWV (millimeters) calculated from the ERA-PWV in the six-year period from 2013 to 2018 over China.

Figure 3 .
Figure 3. (a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and ERA-PWV in the six-year period studied from 2013 to 2018 over each of the 246 GNSS stations.

Figure 3 .
Figure 3. (a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and ERA-PWV in the six-year period studied from 2013 to 2018 over each of the 246 GNSS stations.Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 21

Figure 4 .
Figure 4. Distribution of the differences (in millimeters) between GNSS-PWV, RS-PWV, and ERA-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 4 .
Figure 4. Distribution of the differences (in millimeters) between GNSS-PWV, RS-PWV, and ERA-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 4 .
Figure 4. Distribution of the differences (in millimeters) between GNSS-PWV, RS-PWV, and ERA-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 5 .
Figure 5. (a) Bias, (b) RMS, and (c) correlation coefficient between RS-PWV and ERA-PWV in the six-year period studied from 2013 to 2018 over each of the 78 radiosonde stations.

Figure 5 .
Figure 5. (a) Bias, (b) RMS, and (c) correlation coefficient between RS-PWV and ERA-PWV in the six-year period studied from 2013 to 2018 over each of the 78 radiosonde stations.

Figure 6 .
Figure 6.(a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and MOD-NIR-PWV in the six-year period studied from 2013 to 2018 over each of the 221 GNSS stations.

Figure 6 .
Figure 6.(a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and MOD-NIR-PWV in the six-year period studied from 2013 to 2018 over each of the 221 GNSS stations.

Figure 7 .
Figure 7. Distribution of the differences (in millimeters) between GNSS-PWV, ERA-PWV, and MOD-NIR-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 8 .
Figure 8.(a) Bias, (b) RMS, and (c) correlation coefficient between ERA-PWV and MOD-NIR-PWV in the six-year period from 2013 to 2018 over each grid point.

Figure 7 .
Figure 7. Distribution of the differences (in millimeters) between GNSS-PWV, ERA-PWV, and MOD-NIR-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 7 .
Figure 7. Distribution of the differences (in millimeters) between GNSS-PWV, ERA-PWV, and MOD-NIR-PWV in the six-year period studied from 2013 to 2018 over China.

Figure 8 .
Figure 8.(a) Bias, (b) RMS, and (c) correlation coefficient between ERA-PWV and MOD-NIR-PWV in the six-year period from 2013 to 2018 over each grid point.

Figure 8 .
Figure 8.(a) Bias, (b) RMS, and (c) correlation coefficient between ERA-PWV and MOD-NIR-PWV in the six-year period from 2013 to 2018 over each grid point.

Figure 9 .
Figure 9. Distribution of seasonal biases between ERA-PWV and MOD-NIR-PWV in the six-year period studied over each grid point in (a) spring, (b) summer, (c) autumn, and (d) winter.

Figure 10 .
Figure 10.Distribution of seasonal RMSs between ERA-PWV and MOD-NIR-PWV in the six-year period studied over each grid point in (a) spring, (b) summer, (c) autumn, and (d) winter.

Figure 9 . 21 Figure 9 .
Figure 9. Distribution of seasonal biases between ERA-PWV and MOD-NIR-PWV in the six-year period studied over each grid point in (a) spring, (b) summer, (c) autumn, and (d) winter.

Figure 10 .
Figure 10.Distribution of seasonal RMSs between ERA-PWV and MOD-NIR-PWV in the six-year period studied over each grid point in (a) spring, (b) summer, (c) autumn, and (d) winter.

Figure 10 .
Figure 10.Distribution of seasonal RMSs between ERA-PWV and MOD-NIR-PWV in the six-year period studied over each grid point in (a) spring, (b) summer, (c) autumn, and (d) winter.

Figure 13 .
Figure 13.(a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and calibrated MOD-NIR-PWV in the sixyear period studied over each site of the 221 GNSS stations.

Figure 12 .Figure 12 .
Figure 12.(a) Bias, (b) RMS, and (c) correlation coefficient between ERA-PWV and calibrated MOD-NIR-PWV in the six-year period studied over each grid point.

Figure 13 .
Figure 13.(a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and calibrated MOD-NIR-PWV in the sixyear period studied over each site of the 221 GNSS stations.

Figure 13 . 21 Figure 14 .
Figure 13.(a) Bias, (b) RMS, and (c) correlation coefficient between GNSS-PWV and calibrated MOD-NIR-PWV in the six-year period studied over each site of the 221 GNSS stations.Remote Sens. 2021, 13, x FOR PEER REVIEW 18 of 21

Figure 14 .
Figure 14.Distribution of the differences (in millimeters) between GNSS-PWV, ERA-PWV, and calibrated MOD-NIR-PWV in the six-year period studied over China.

Table 1 .
Geographical location, seasonal bias, and RMS in the six-year period studied at each of four randomly selected grid points.

Table 1 .
Geographical location, seasonal bias, and RMS in the six-year period studied at each of four randomly selected grid points.

Bias (mm) RMS (mm) Spring Summer Autumn Winter Spring Summer Autumn Winter
Remote Sens. 2021, 13, x FOR PEER REVIEW 16 of 21