A Framework for Estimating Clear-Sky Atmospheric Total Precipitable Water (TPW) from VIIRS/S-NPP

Atmospheric water vapor content or total precipitable water (TPW) is a highly variable atmospheric constituent, yet it remains one of the meteorological parameters that is most difficult to characterize accurately. We develop a framework for estimating atmospheric TPW from Visible Infrared Imaging Radiometer Suite (VIIRS) data in this study. First, TPW is retrieved from VIIRS top-of-atmosphere (TOA) radiance of channels 15 and 16 using the refined split-window covariance-variance ratio (SWCVR) method. Then, the VIIRS TPW is blended with the microwave integrated retrieval system (MIRS) derived TPW via Bayesian model averaging (BMA) to improve the accuracy of VIIRS TPW. Three years (2014–2017) of ground measurements collected from SuomiNet sites over North America are used to validate the VIIRS TPW and blended TPW. The mean bias error (MBE) and root mean square error (RMSE) of the VIIRS TPW are 0.21 g/cm2 and 0.73 g/cm2, respectively, and the accuracy of the VIIRS TPW in daytime is much better than at night time. The MBE and RMSE of BMA integrated TPW are 0.06 g/cm2 and 0.35 g/cm2, and the accuracy difference between daytime and nighttime is also removed. The global radiosonde measurements are also collected to validate the BMA integrated VIIRS TPW. The MBE and RMSE of the BMA integrated TPW are 0.09 g/cm2 and 0.44 g/cm2 compared to the radiosonde measurements. This accuracy is also superior to the VIIRS TPW. Therefore, it is concluded that the developed framework can be used to derive accurate clear-sky TPW for VIIRS. This is the first time that we can obtain high accuracy TPW from VIIRS. This study will certainly benefit the study of atmospheric processes and climate change.


Introduction
Total atmospheric water content is the mass of water vapor contained in a vertical atmospheric column from the surface to the top of the atmosphere, also known as total precipitable water (TPW). Atmospheric water vapor is one of the most important atmospheric absorbers of greenhouse gases, accounting for approximately 70% of the total atmospheric absorption and 60% of the total natural greenhouse effects under clear-sky conditions [1], and it plays an important role in the energy exchange between the atmosphere and ground surface, the Earth's hydrological cycle [2], and climate change [2][3][4]. Water vapor is a highly variable atmospheric constituent, yet it remains one of the meteorological parameters that is most difficult to characterize accurately. It is important to improve our knowledge of atmospheric water vapor for a variety of atmospheric research and applications.
There are three types of methods for obtaining TPW: ground-based observations, radiosonde techniques, and satellite remote sensing. It is widely acknowledged that the radiosonde technique The description of the method and the validation results are respectively provided in Sections 3 and 4. A brief discussion and conclusion are provided in Sections 5 and 6.

SuomiNet
"SuomiNet" is a network of GPS receivers, which provides real-time measurements of TPW for atmospheric research and education [8,35]. The network, named after Verner Suomi, a weather satellite pioneer, uses the ground-based GPS receiver to conduct thousands of accurate measurements of the upper and lower atmosphere each day. The GPS sends the radio signal to the global terrestrial GPS receiver. These signals are delayed by the ionosphere and troposphere. The ionospheric effect can be eliminated by using a linear combination of two GPS frequencies. The total tropospheric delay along the zenith is called the zenith tropospheric delay (ZTD), which can be divided into two parts: the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD). ZHD is mainly a function of the surface pressure of the GPS receiver, and ZWD strongly depends on TPW [36,37]. ZWD delay can be converted into TPW along each GPS ray path. SuomiNet can obtain continuous, accurate, all-weather, real-time TPW data, which will help promote mesoscale modeling and data assimilation, bad weather, precipitation, and cloud dynamics, regional climate, and hydrology research.
SuomiNet provides continuous TPW estimates from 828 sites distributed across the United States with accuracy better than 2 mm (~0.2 g/cm 2 ) [8,38]. In this article, three years (2014-2017) of SuomiNet data are downloaded for the validation of the retrieved TPW. The geolocations of SuomiNet GPS stations are displayed in Figure 1.

SuomiNet
"SuomiNet" is a network of GPS receivers, which provides real-time measurements of TPW for atmospheric research and education [8,35]. The network, named after Verner Suomi, a weather satellite pioneer, uses the ground-based GPS receiver to conduct thousands of accurate measurements of the upper and lower atmosphere each day. The GPS sends the radio signal to the global terrestrial GPS receiver. These signals are delayed by the ionosphere and troposphere. The ionospheric effect can be eliminated by using a linear combination of two GPS frequencies. The total tropospheric delay along the zenith is called the zenith tropospheric delay (ZTD), which can be divided into two parts: the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD). ZHD is mainly a function of the surface pressure of the GPS receiver, and ZWD strongly depends on TPW [36,37]. ZWD delay can be converted into TPW along each GPS ray path. SuomiNet can obtain continuous, accurate, all-weather, real-time TPW data, which will help promote mesoscale modeling and data assimilation, bad weather, precipitation, and cloud dynamics, regional climate, and hydrology research.
SuomiNet provides continuous TPW estimates from 828 sites distributed across the United States with accuracy better than 2 mm (~0.2g/cm 2 ) [8,38]. In this article, three years (2014-2017) of SuomiNet data are downloaded for the validation of the retrieved TPW. The geolocations of SuomiNet GPS stations are displayed in Figure 1.

The Joint Polar Satellite System (JPSS) Microwave Integrated Retrieval System (MIRS) Precipitation Product
The microwave integrated retrieval system (MIRS) was developed by the NOAA/National Environmental Satellite, Data, and Information Service (NESDIS) Center for Satellite Application and Research (STAR) and has been put into use in the NOAA/NESDIS Satellite and Product Operations Office (OSPO), which aims to upgrade the current operational microwave surface and precipitation products system (MSPPS) [39,40]. The operational MIRS can produce advanced near-real-time surface and precipitation products in all weather and over all surface conditions using brightness temperatures from the microwave instruments, such as the advanced technology microwave sounder (ATMS) onboard the S-NPP satellite. ATMS is a cross-track scanner with 22 channels in spectral bands from 23 GHz through 183 GHz [41]. The operational products of MIRS include vertical profiles of temperature and TPW, rainfall rate, cloud liquid water, snow cover, snow water equivalent, sea ice

The Joint Polar Satellite System (JPSS) Microwave Integrated Retrieval System (MIRS) Precipitation Product
The microwave integrated retrieval system (MIRS) was developed by the NOAA/National Environmental Satellite, Data, and Information Service (NESDIS) Center for Satellite Application and Research (STAR) and has been put into use in the NOAA/NESDIS Satellite and Product Operations Office (OSPO), which aims to upgrade the current operational microwave surface and precipitation products system (MSPPS) [39,40]. The operational MIRS can produce advanced near-real-time surface and precipitation products in all weather and over all surface conditions using brightness temperatures from the microwave instruments, such as the advanced technology microwave sounder (ATMS) onboard the S-NPP satellite. ATMS is a cross-track scanner with 22 channels in spectral bands from 23 GHz through 183 GHz [41]. The operational products of MIRS include vertical profiles of temperature and TPW, rainfall rate, cloud liquid water, snow cover, snow water equivalent, sea ice concentration, ice water path, and land surface temperature. The spatial resolution of the MIRS TPW product derived from ATMS is 15 km at Nadir. Suomi NPP orbits the Earth 14.1875 times each day in a Sun-synchronous, polar orbit that allows ATMS to observe nearly the entire global atmosphere twice daily.

Visible Infrared Imaging Radiometer Suite (VIIRS)
The Visible Infrared Imaging Radiometer Suite (VIIRS), aboard the joint NASA/NOAA Suomi National Polar-Orbiting Partnership (S-NPP) satellite, was developed based on the heritage of legacy instruments, including the advanced very high resolution radiometer (AVHRR) and MODIS. VIIRS has 22 spectral channels with wavelengths covering from 0.41 µm to 12.5 µm, which can be used for environmental monitoring and numerical weather forecasting. The results from the on-orbit verification in the postlaunch check-out, intensive calibration, and validation have shown that VIIRS is performing very well [42]. More than 20 environmental data records have been produced operationally from VIIRS data, including clouds, land surface temperature, sea surface temperature, aerosol optical thickness, ocean color, polar wind, vegetation fraction, aerosol, fire, snow and ice, vegetation index, etc. However, to our knowledge, VIIRS does not provide TPW product.
The VIIRS data utilized in this study include VIIRS sensor data records (SDR) and VIIRS cloud mask intermediate product (VCM). The VIIRS SDR contains the day-night channel, imagery channel, moderate resolution channel, and geolocation data. Two thermal-infrared channels of VIIRS, M15 (10.75 µm) and M16 (12.01 µm), are selected. The VIIRS VCM product is used to identify whether a pixel is clear or cloudy [43]. The pixels with a confident clear flag are identified as clear-sky pixels. In addition, geolocation and sensor view zenith angle data are also utilized.

Seebor Atmosphere Profile
This Seebor global profiles database (named SeeBor V5.0) consists of 15,704 profiles distributed all over the world [44]. The temperature profile of Seebor has 101 levels, and the corresponding atmospheric pressure varies from 1100 hPa to 0.005 hPa. The atmosphere profiles are collected from NOAA-88, a European Centre for Medium-Range Weather Forecasts(ECMWF) 60 L training set, TIGR-3, ozonesondes from 8 NOAA Climate Monitoring and Diagnostics Laboratory (CMDL) sites, and radiosondes from 2004 in the Sahara Desert. All the profiles are checked with the following saturation criteria to eliminate the cloudy atmosphere profile, i.e., the relative humidity value for each profile at each level below the 250 hPa pressure level must be less than 90%.

Split-Window Covariance-Variance Ratio Method
Assuming that the atmosphere state is invariant within a certain area, and the variation of surface emissivity is minimal, the split-window covariance-variance ratio method (SWCVR) [22,23] retrieves the TPW using the channel transmission ratio of two split-window channels (approximately 11 µm and 12 µm). The transmittance ratio can be obtained from the TOA brightness temperature of the two split-window channels. The SWCVR is expressed as: where TPW represents total precipitable water, k represents pixel k, the subscripts i and j represent the channels 15 and 16 of the VIIRS,τ i and ε i represent transmittance and emissivity of VIIRS channel i respectively, T i,k represents TOA brightness temperature of VIIRS channel i for pixel k, T i represents mean TOA brightness temperature of VIIRS channel i over N neighboring pixels, and r 2 denotes the correlation coefficient of T i and T j over the N neighboring pixels, and is used to determine whether the previous assumption is valid in the derivation process of the transmittance ratio. The transmission ratio can be retrieved using Equation (2) in the case that the emissivity ratio of two split-window channels is known. R ji is calculated from the atmospheric top brightness temperature of the two split-window channels using Equation (3). The Advanced Spaceborne Thermal Emission and Reflection Radiometer(ASTER) emissivity library [45] and moderate resolution imaging spectroradiometer (MODIS) University of California, Santa Barbara(UCSB) emissivity library (https://icess.eri.ucsb.edu/modis/EMIS/html/em.html) is used to calculate the channel emissivity of VIIRS split-window channels. Figure 2 shows the variation of ε i /ε j for different natural surface types including rocks, vegetation, soils, water, and snow. The emissivity ratios of channels 11 µm and 12 µm for most vegetation, water, snow, and soil are between 0.98 and 1.02. Considering that most pixels at a scale of 0.75 km are usually covered by soils, water, vegetation, snow, or their mixtures, the emissivity ratio can be approximated as unity in practice. Thus, the channel transmission ratio can be calculated as: where TPW represents total precipitable water, k represents pixel k, the subscripts i and j represent the channels 15 and 16 of the VIIRS, i τ and i ε represent transmittance and emissivity of VIIRS channel i respectively, , i k T represents TOA brightness temperature of VIIRS channel i for pixel k, i T represents mean TOA brightness temperature of VIIRS channel i over N neighboring pixels, and 2 r denotes the correlation coefficient of i T and j T over the N neighboring pixels, and is used to determine whether the previous assumption is valid in the derivation process of the transmittance ratio.
The transmission ratio can be retrieved using Equation (2) in the case that the emissivity ratio of two split-window channels is known. ji R is calculated from the atmospheric top brightness temperature of the two split-window channels using Equation (3). The Advanced Spaceborne Thermal Emission and Reflection Radiometer(ASTER) emissivity library [45] and moderate resolution imaging spectroradiometer (MODIS) University of California, Santa Barbara(UCSB) emissivity library (https://icess.eri.ucsb.edu/modis/EMIS/html/em.html) is used to calculate the channel emissivity of VIIRS split-window channels. Figure 2 shows the variation of i j ε ε for different natural surface types including rocks, vegetation, soils, water, and snow. The emissivity ratios of channels 11 μm and 12 μm for most vegetation, water, snow, and soil are between 0.98 and 1.02. Considering that most pixels at a scale of 0.75 km are usually covered by soils, water, vegetation, snow, or their mixtures, the emissivity ratio can be approximated as unity in practice. Thus, the channel transmission ratio can be calculated as: The number of neighboring pixels is a key parameter for the SWCVR method. The window size N (N × N neighboring pixels) should not be too large or too small, in case the atmospheric conditions may not be constant or cannot guarantee a high correlation between two channels' TOA brightness The number of neighboring pixels is a key parameter for the SWCVR method. The window size N (N × N neighboring pixels) should not be too large or too small, in case the atmospheric conditions may not be constant or cannot guarantee a high correlation between two channels' TOA brightness temperatures [28]. A simple simulation experiment is conducted to investigate the effect of window size on the accuracy of the SWCVR method using the MODIS emissivity product MOD11B1. Seebor atmosphere profiles with a standard deviation of 2.5 K. The TOA brightness temperatures of VIIRS channels 15 and 16 are simulated using the moderate resolution atmospheric radiative transfer code (MODTRAN 5.0) for different window sizes. The TPW is retrieved using the SWCVR method from simulated VIIRS TOA brightness temperatures for different window sizes, and the root mean square error (RMSE) which are displayed in Figure 3. Obviously, the RMSE decreases as the window size increases, and the retrieval accuracy does not improve significantly with the increasing number of the neighboring pixels when the window is larger than 18. Therefore, the window size is set as 18 in this study.
temperatures [28]. A simple simulation experiment is conducted to investigate the effect of window size on the accuracy of the SWCVR method using the MODIS emissivity product MOD11B1. The emissivities of MODIS channels 31 and 32 are used as proxies of VIIRS channels M15 and M16, and the standard deviations of emissivities, approximately 0.01, are collected for different window sizes. The land surface temperature (LST) varies around the bottom layer temperature of different Seebor atmosphere profiles with a standard deviation of 2.5 K. The TOA brightness temperatures of VIIRS channels 15 and 16 are simulated using the moderate resolution atmospheric radiative transfer code (MODTRAN 5.0) for different window sizes. The TPW is retrieved using the SWCVR method from simulated VIIRS TOA brightness temperatures for different window sizes, and the root mean square error (RMSE) which are displayed in Figure 3. Obviously, the RMSE decreases as the window size increases, and the retrieval accuracy does not improve significantly with the increasing number of the neighboring pixels when the window is larger than 18. Therefore, the window size is set as 18 in this study. T are used instead of the mean values, as the median is a more robust estimator than the mean value. 2) Remove the cloudy pixels. The VIIRS cloud mask product is used to identify the cloudy pixels, and only confidently clear pixels are used to calculate ji R .
3) Eliminate the invalid pixels. Since j τ (at 12 μm) is smaller than i τ (at 11 μm) and the ratio of emissivities is very close to unity for most conditions, j i τ τ should be smaller than unity. The absolute value of 15 15 ( ) T T − must be larger than the absolute value of  Eliminating the invalid pixel can effectively avoid invalid inversion of transmittance ratios (i.e., a transmittance ratio lager than 1). For more details, please refer to Li et al. [23].  The SWCVR algorithm is implemented in five steps:

4) Estimate
(1) Calculate the mean values of TOA brightness temperature T i and T j . Herein, the median values T i and T j are used instead of the mean values, as the median is a more robust estimator than the mean value. (2) Remove the cloudy pixels. The VIIRS cloud mask product is used to identify the cloudy pixels, and only confidently clear pixels are used to calculate R ji . (3) Eliminate the invalid pixels. Since τ j (at 12 µm) is smaller than τ i (at 11 µm) and the ratio of emissivities is very close to unity for most conditions, τ j /τ i should be smaller than unity. The absolute value of (T 15 − T 15 ) must be larger than the absolute value of (T 16 − T 16 ) and the values of (T 15 − T 15 ) and (T 16 − T 16 ) should have the same sign. Eliminating the invalid pixel can effectively avoid invalid inversion of transmittance ratios (i.e., a transmittance ratio lager than 1). For more details, please refer to Li et al. [23]. (4) Estimate r 2 and reject R ji if r 2 is less than 0.95. (5) Calculate the TPW for different view zenith angles. The relationship between TPW and the transmission ratio is established at several specific angles. When the view zenith angle is not equal to the specific angles, TPW is linearly interpolated from TPWs predicted by the model with an adjacent view zenith angle. View zenith angles exceeding 75 • are not considered.  Table 1. The fitting accuracy decreases with the increasing view angle. The determination coefficient ranges from 0.976 to 0.983, and RMSE ranges from 0.191 to 0.223 g/cm 2 . predicted by the model with an adjacent view zenith angle. View zenith angles exceeding 75° are not considered. Figure 4 shows scatterplots of TPW versus the transmission ratio for VIIRS under different view angles (0°, 15°, 30°, 45°, 60°, 75°). The TPW and transmission ratio are simulated by MODTRAN 5 using Seebor atmosphere profiles and VIIRS filter response functions. The cubic polynomial is used to fit the relationship between matchup of the TPW of each Seebor atmosphere profile and corresponding VIIRS transmission ratio. The fitting results are provided in Table 1. The fitting accuracy decreases with the increasing view angle. The determination coefficient ranges from 0.976 to 0.983, and RMSE ranges from 0.191 to 0.223 g/cm 2 .   is unsure which one is the best [46][47][48]. BMA does not depend on a single "best" model, but fuses multiple models to obtain optimal results by considering the errors of different models, so the results of BMA are more stable than those of a single model. The BMA predictive probability density function (PDF) is a weighted average of PDFs centered on each of the models considered, where the weights of each model are equal to posterior probabilities of these models. The weights can reflect relative contributions of each model to the final output of the BMA model over the training period.
BMA can be used to obtain a more accurate and stable estimate of TPW by blending the TPW retrieved from the thermal infrared data and that retrieved from passive microwave data. For convenience, we denote the TPW value to be forecast from two data sources as y and the ground measurement at a certain time to be tpw. f 1 , f 2 , f 3 . . . . f N are forecast results from each model. According to the total probability formula, the forecast PDF is given by where the p(y f k ) is the forecast PDF based on f k alone and p( f k tpw) is the posterior probability of f k being correct given the training data. The posterior model probabilities add up to one, N k=1 p( f k tpw) = 1, such that they can be viewed as weights. The BMA predictive model is then where w k is the posterior probability of forecast k. It is assumed the conditional PDF of y is normally distributed with mean a k + b k f k and standard deviation σ, which can be written as: In this case, the BMA predictive mean is the conditional expectation of y given the forecasts, namely a k and b k can be estimated by simple linear regression using training data, which can be viewed as a bias correction process. w k and σ can be estimated by maximum likelihood from the training data [49]. The expectation-maximization algorithm is used to maximize the likelihood function in this study [50].

Validating VIIRS TPW and MIRS-Derived TPW
Two indicators are selected to depict the accuracy of the derived TPW, the mean bias error (MBE) and root mean square error (RMSE). These indicators are defined as follows: where N is the number of samples; E i,true is the true value of TPW, derived from atmosphere profile or measured by the ground-based GPS network; E i,derived is the derived TPW.
In this study, three years (2014-2017) of VIIRS data are used to derive the TPW using the refined SWCVR. The MIRS-derived TPW product synchronized in time and space with VIIRS is also extracted. The TPW measured by SuomiNet GPS receivers is used to validate the above two TPWs, Figures 5  and 6 display the validation results. The MBE and RMSE of VIIRS TPW are 0.21 g/cm 2 and 0.73 g/cm 2 , whereas MBE and RMSE of MIRS-derived TPW are 0.11 g/cm 2 and 0.38 g/cm 2 . Clearly, the accuracy of MIRS-derived TPW is better than that of VIIRS TPW. Regarding the validation results for daytime and nighttime, no significant difference is found for MIRS-derived TPW. The MBE and RMSE are 0.12 g/cm 2 and 0.39 g/cm 2 for daytime, and 0.10 g/cm 2 and 0.36 g/cm 2 for nighttime. However, the accuracy of VIIRS TPW in daytime is better than that at nighttime. The MBE and RMSE are 0.05 g/cm 2 and 0.61 g/cm 2 for daytime, and 0.35 g/cm 2 and 0.83 g/cm 2 for nighttime. There may be two reasons for this result. First, the SWCVR may achieve good results when the surface emissivity varies little and the surface temperature varies greatly. However, the LST is more homogeneous in the nighttime than that in the daytime [51]. Second, cloud detection at night is worse than that in the daytime because there is no shortwave observation at night [52]. The VIIRS cloud mask cannot filter all cloudy pixels, and cannot entirely identify whether the pixels around the surface GPS station are confidently clear. This can partly explain the underestimation of VIIRS TPW under the high TPW condition in the nighttime. In addition, an obvious overestimation of VIIRS TPW can be found when TPW is low for both daytime and nighttime. The detailed reasons are discussed in a later section.
where N is the number of samples; i,true E is the true value of TPW, derived from atmosphere profile or measured by the ground-based GPS network; ,derived i E is the derived TPW.
In this study, three years (2014-2017) of VIIRS data are used to derive the TPW using the refined SWCVR. The MIRS-derived TPW product synchronized in time and space with VIIRS is also extracted. The TPW measured by SuomiNet GPS receivers is used to validate the above two TPWs, Figures 5 and 6 display the validation results. The MBE and RMSE of VIIRS TPW are 0.21 g/cm 2 and 0.73 g/cm 2 , whereas MBE and RMSE of MIRS-derived TPW are 0.11 g/cm 2 and 0.38 g/cm 2 . Clearly, the accuracy of MIRS-derived TPW is better than that of VIIRS TPW. Regarding the validation results for daytime and nighttime, no significant difference is found for MIRS-derived TPW. The MBE and RMSE are 0.12 g/cm 2 and 0.39 g/cm 2 for daytime, and 0.10 g/cm 2 and 0.36 g/cm 2 for nighttime. However, the accuracy of VIIRS TPW in daytime is better than that at nighttime. The MBE and RMSE are 0.05 g/cm 2 and 0.61 g/cm 2 for daytime, and 0.35 g/cm 2 and 0.83 g/cm 2 for nighttime. There may be two reasons for this result. First, the SWCVR may achieve good results when the surface emissivity varies little and the surface temperature varies greatly. However, the LST is more homogeneous in the nighttime than that in the daytime [51]. Second, cloud detection at night is worse than that in the daytime because there is no shortwave observation at night [52]. The VIIRS cloud mask cannot filter all cloudy pixels, and cannot entirely identify whether the pixels around the surface GPS station are confidently clear. This can partly explain the underestimation of VIIRS TPW under the high TPW condition in the nighttime. In addition, an obvious overestimation of VIIRS TPW can be found when TPW is low for both daytime and nighttime. The detailed reasons are discussed in a later section.

Blending TPWs Using BMA
In this section, VIIRS TPW and clear-sky MIRS-derived TPW are integrated by BMA. The pixel spatial resolution difference between VIIRS TPW and MIRS-derived TPW is ignored during the process of integration. Two-thirds of the VIIRS TPW data and MIRS TPW data in the daytime are used for BMA model training, the remaining data (including one-third of the data in the daytime and all data in the nighttime) are used for validation. The blending results are displayed in Figure 7. The MBE and RMSE for BMA-derived TPW are 0.06 g/cm 2 and 0.35 g/cm 2 , respectively, which is lower than that of MIRS-derived TPW (Figure 6a). Obviously, the accuracy of blended TPW is better than that of VIIRS TPW and MIRS-derived TPW on the whole, which demonstrates the good performance of BMA in integrating multi-source TPW data.

Blending TPWs Using BMA
In this section, VIIRS TPW and clear-sky MIRS-derived TPW are integrated by BMA. The pixel spatial resolution difference between VIIRS TPW and MIRS-derived TPW is ignored during the process of integration. Two-thirds of the VIIRS TPW data and MIRS TPW data in the daytime are used for BMA model training, the remaining data (including one-third of the data in the daytime and all data in the nighttime) are used for validation. The blending results are displayed in Figure 7. The MBE and RMSE for BMA-derived TPW are 0.06 g/cm 2 and 0.35 g/cm 2 , respectively, which is lower than that of MIRS-derived TPW (Figure 6a). Obviously, the accuracy of blended TPW is better than that of VIIRS TPW and MIRS-derived TPW on the whole, which demonstrates the good performance of BMA in integrating multi-source TPW data.
To further compare the three TPW data sets, we divided the validation results into four sub-regions according to the range of water vapor. The statistical results are provided in Tables 2 and 3. The accuracy of MIRS-derived TPW and BMA-derived TPW are both better than the accuracy of VIIRS TPW. When the TPW is less than 3 g/cm 2 , the accuracy of BMA-derived TPW is better than that of MIRS-derived TPW. When the TPW is greater than 3 g/cm 2 , the accuracy of BMA-derived TPW is slightly worse than that of MIRS-derived TPW in the daytime and nighttime. This accuracy reducing may be caused by the undetected cloud pixels used in the SWCVR algorithm. With the upgrade of the VIIRS cloud detection algorithm, this problem may be alleviated to some extent. As a whole, blending VIIRS and MIRS-derived TPW using BMA improves the accuracy of TPW retrieval. To further compare the three TPW data sets, we divided the validation results into four subregions according to the range of water vapor. The statistical results are provided in Tables 2 and 3. The accuracy of MIRS-derived TPW and BMA-derived TPW are both better than the accuracy of VIIRS TPW. When the TPW is less than 3 g/cm 2 , the accuracy of BMA-derived TPW is better than that of MIRS-derived TPW. When the TPW is greater than 3 g/cm 2 , the accuracy of BMA-derived TPW is slightly worse than that of MIRS-derived TPW in the daytime and nighttime. This accuracy reducing may be caused by the undetected cloud pixels used in the SWCVR algorithm. With the upgrade of the VIIRS cloud detection algorithm, this problem may be alleviated to some extent. As a whole, blending VIIRS and MIRS-derived TPW using BMA improves the accuracy of TPW retrieval.

Production of VIIRS TPW
To examine the possibility of producing VIIRS TPW using the developed framework, we select a VIIRS granule image with relatively less cloud cover and the corresponding MIRS-derived TPW product beyond the acquisition time of validation data. The overpass time of the selected VIIRS granule image is at 20:18 (UTC), 25 September, 2018, and the MIRS-derived TPW product is acquired from 20:10 to 20:20 (UTC) on the same day. The VIIRS TPW is retrieved by the SWCVR method. The MIRS-derived TPW is interpolated to the VIIRS pixels using the nearest neighbor method. The resampled MIRS TPW is blended with VIIRS TPW using BMA. Figure 8 shows the produced BMA TPW, as well as VIIRS TPW and MIRS-derived TPW. The spatial pattern of VIIRS TPW is similar to that of MIRS-derived TPW, but the VIIRS TPW exhibits more details than the interpolated MIRS-derived TPW. The difference between TPW values in the VIIRS TPW and MIRS-derived TPW can be easily observed. The weight of MIRS-derived TPW is larger than that of VIIRS TPW because MIRS-derived TPW is more accurate than that of VIIRS TPW. Thus, the spatial pattern of BMA TPW is more similar to that of MIRS-derived TPW. We also extract the SuomiNet TPW in the coverage of the VIIRS granule and validate three TPWs. The validation results are shown in Figure 9. The MBE and RMSE of BMA TPW are −0.024 and 0.314 g/cm 2 respectively, which is superior to VIIRS TPW and comparable to MIRS-derived TPW whose MBE and RMSE are −0.238 and 0.535 g/cm 2 , and −0.029 and 0.337g/cm 2 , respectively. This validation agrees with that presented in Section 4.2. We also extract the SuomiNet TPW in the coverage of the VIIRS granule and validate three TPWs. The validation results are shown in Figure 9. The MBE and RMSE of BMA TPW are −0.024 and 0.314 g/cm 2 respectively, which is superior to VIIRS TPW and comparable to MIRS-derived TPW whose MBE and RMSE are −0.238 and 0.535 g/cm 2 , and −0.029 and 0.337 g/cm 2 , respectively. This validation agrees with that presented in Section 4.2. We also extract the SuomiNet TPW in the coverage of the VIIRS granule and validate three TPWs. The validation results are shown in Figure 9. The MBE and RMSE of BMA TPW are −0.024 and 0.314 g/cm 2 respectively, which is superior to VIIRS TPW and comparable to MIRS-derived TPW whose MBE and RMSE are −0.238 and 0.535 g/cm 2 , and −0.029 and 0.337g/cm 2 , respectively. This validation agrees with that presented in Section 4.2.

Validating BMA Integrated TPW Using Radiosonde Measurements
To verify the applicability of the BMA method in other areas, two years (2014-2015) of global climate observing system (GCOS) reference upper-air network (GRUAN) radiosonde data are collected. The GRUAN measurements provide appropriate data for studying atmospheric processes. Details about the GRUAN radiosonde measurements can be found at https://www.gruan.org. The TPW of radiosonde measurement is estimated by vertically integrating the specific humidity: where q is specific humidity, g is acceleration of gravity, and p is the atmospheric pressure.
The VIIRS granule image is considered when the difference between the VIIRS overpass time and the balloon launch time is less than ten minutes. The VIIRS TPW is calculated using the SWCVR method, and then is blended with MIRS-derived TPW using the BMA method. According to Figure  10, the MBE and RMSE of BMA TPW are 0.09 g/cm 2 and 0.44 g/cm 2 when compared to the radiosonde measurements. The result is more accurate than both VIIRS TPW and MIRS-derived TPW, whose MBE and RMSE are −0.27 g/cm 2 and 0.82 g/cm 2 , and 0.28 g/cm 2 and 0.52 g/cm 2 , respectively. The results show that the BMA developed in North America can also be applied at a global scale.

Validating BMA Integrated TPW Using Radiosonde Measurements
To verify the applicability of the BMA method in other areas, two years (2014-2015) of global climate observing system (GCOS) reference upper-air network (GRUAN) radiosonde data are collected. The GRUAN measurements provide appropriate data for studying atmospheric processes. Details about the GRUAN radiosonde measurements can be found at https://www.gruan.org. The TPW of radiosonde measurement is estimated by vertically integrating the specific humidity: where q is specific humidity, g is acceleration of gravity, and p is the atmospheric pressure. The VIIRS granule image is considered when the difference between the VIIRS overpass time and the balloon launch time is less than ten minutes. The VIIRS TPW is calculated using the SWCVR method, and then is blended with MIRS-derived TPW using the BMA method. According to Figure 10, the MBE and RMSE of BMA TPW are 0.09 g/cm 2 and 0.44 g/cm 2 when compared to the radiosonde measurements. The result is more accurate than both VIIRS TPW and MIRS-derived TPW, whose MBE and RMSE are −0.27 g/cm 2 and 0.82 g/cm 2 , and 0.28 g/cm 2 and 0.52 g/cm 2 , respectively. The results show that the BMA developed in North America can also be applied at a global scale.
The VIIRS granule image is considered when the difference between the VIIRS overpass time and the balloon launch time is less than ten minutes. The VIIRS TPW is calculated using the SWCVR method, and then is blended with MIRS-derived TPW using the BMA method. According to Figure  10, the MBE and RMSE of BMA TPW are 0.09 g/cm 2 and 0.44 g/cm 2 when compared to the radiosonde measurements. The result is more accurate than both VIIRS TPW and MIRS-derived TPW, whose MBE and RMSE are −0.27 g/cm 2 and 0.82 g/cm 2 , and 0.28 g/cm 2 and 0.52 g/cm 2 , respectively. The results show that the BMA developed in North America can also be applied at a global scale.

SWCVR
In Section 3, we explained that most of the emissivity ratios of two split-window channels for vegetation, water, snow, and soil are between 0.98 and 1.02, then that the transmission ratio of the split-window channel can be calculated using Equation (3). Since the transmittance in the channel at 12 μm is smaller than that in the channel at 11 μm, j i τ τ is smaller than unity. To avoid bad calculations of j i τ τ , the pixel filtering criteria proposed by Li et al. [23] is adopted, i.e., the absolute value of 15 15 ( ) T T − must be larger than the absolute value of directly. The emissivity standard deviation for each window is also calculated. The surface temperature of these pixels is set as a Gaussian distribution, with the mean value equal to the bottom layer temperature of each Seebor atmosphere profile. The standard deviations are 2.5 K and 1 K for different conditions. The atmosphere transmission and TOA bright temperature of VIIRS are simulated using MODTRAN 5, and the atmosphere transmission ratio is calculated using Equation (2) with the simulated VIIRS TOA bright temperature. Detailed settings are provided in Table 4. Six sets of simulated data are separately used to calculate transmission ratios respectively. Then the calculated transmission ratio is compared to the true transmission ratio of each atmosphere profile. The comparison results are displayed in Figure 11. As shown in Figure 11, it is apparent that pixel filtering criteria can actually avoid unreasonable transmission ratios when using SWCVR. However, the criteria may cause underestimation of transmission ratios when the transmission ratio is greater

SWCVR
In Section 3, we explained that most of the emissivity ratios of two split-window channels for vegetation, water, snow, and soil are between 0.98 and 1.02, then that the transmission ratio of the split-window channel can be calculated using Equation (3). Since the transmittance in the channel at 12 µm is smaller than that in the channel at 11 µm, τ j /τ i is smaller than unity. To avoid bad calculations of τ j /τ i , the pixel filtering criteria proposed by Li et al. [23] is adopted, i.e., the absolute value of (T 15 − T 15 ) must be larger than the absolute value of (T 16 − T 16 ), moreover, values of (T 15 − T 15 ) and (T 16 − T 16 ) should have the same sign. This criterion can ensure the rationality of the calculated transmission ratio and improve the accuracy of water vapor retrievals, which has been demonstrated in Li et al. [23]. However, it may cause underestimation of the transmission ratio under certain circumstances.
To exploit the effect of pixel filtering criteria in refined SWCVR, we conduct a simulation using MOD11B1 and Seebor atmosphere profiles. One year (2014) of MOD11B1 data covering North America is downloaded. An 18 × 18 window is used to extract the emissivities of channels 31 and 32. The spectral response functions of VIIRS channels 15 and 16 are similar to that of MODIS channel 31 and 32. The emissivities of VIIRS channels 15 and 16 are represented by MODIS channels 31 and 32 directly. The emissivity standard deviation for each window is also calculated. The surface temperature of these pixels is set as a Gaussian distribution, with the mean value equal to the bottom layer temperature of each Seebor atmosphere profile. The standard deviations are 2.5 K and 1 K for different conditions. The atmosphere transmission and TOA bright temperature of VIIRS are simulated using MODTRAN 5, and the atmosphere transmission ratio is calculated using Equation (2) with the simulated VIIRS TOA bright temperature. Detailed settings are provided in Table 4. Six sets of simulated data are separately used to calculate transmission ratios respectively. Then the calculated transmission ratio is compared to the true transmission ratio of each atmosphere profile. The comparison results are displayed in Figure 11. As shown in Figure 11, it is apparent that pixel filtering criteria can actually avoid unreasonable transmission ratios when using SWCVR. However, the criteria may cause underestimation of transmission ratios when the transmission ratio is greater than 0.85. Additionally, with the increase in the standard deviation of surface emissivity and the decrease in the standard deviation of surface temperature, the underestimation of transmission ratios becomes more significant under low TPW conditions (i.e., a transmission ratio lager than 0.85). Since the transmittance ratio is inversely proportional to the TPW, this can account for the overestimation of VIIRS TPW when the TPW is less than 1.5 g/cm 2 .

The Limitations of this Study
A scale mismatch exists during the BMA integration of VIIRS TPW and MIRS-derived TPW. If both VIIRS and ATMS TPW are clear-sky pixels, it is reasonable to ignore the scale mismatch because the atmosphere is uniform. However, we cannot integrate VIIRS and ATMS TPW directly if the ATMS TPW pixel is partially cloudy, due to the fact that the 15 km ATMS pixel represents the overall water vapor state while VIIRS only represents the 750 m clear-sky condition. In order to deal with this condition, ATMS pixels need to be decomposed into clear sky and cloudy sky portions. Then the ATMS clear-sky TPW can be used to fuse with VIIRS TPW, and the cloudy sky TPW can be used to fill the gap of VIIRS TPW considering the scale effect. We will work on this in the future.

Conclusions
Atmospheric water content plays a pivotal role in atmospheric processes and the Earth's climate. To our knowledge, VIIRS does not provide an operational TPW product. We developed a framework for estimating atmospheric TPW from Visible Infrared Imaging Radiometer Suite (VIIRS) data in this study. First, the refined SWCVR method was used to retrieve TPW from the VIIRS TOA brightness  According to the above analysis, incorporating surface emissivity or normalized difference vegetation index (NDVI) data to select the relatively uniform part in the window size for TPW calculation in the SWCVR algorithm may improve the accuracy of the SWCVR method. This will be attempted in a future study.

The Limitations of this Study
A scale mismatch exists during the BMA integration of VIIRS TPW and MIRS-derived TPW. If both VIIRS and ATMS TPW are clear-sky pixels, it is reasonable to ignore the scale mismatch because the atmosphere is uniform. However, we cannot integrate VIIRS and ATMS TPW directly if the ATMS TPW pixel is partially cloudy, due to the fact that the 15 km ATMS pixel represents the overall water vapor state while VIIRS only represents the 750 m clear-sky condition. In order to deal with this condition, ATMS pixels need to be decomposed into clear sky and cloudy sky portions. Then the ATMS clear-sky TPW can be used to fuse with VIIRS TPW, and the cloudy sky TPW can be used to fill the gap of VIIRS TPW considering the scale effect. We will work on this in the future.

Conclusions
Atmospheric water content plays a pivotal role in atmospheric processes and the Earth's climate. To our knowledge, VIIRS does not provide an operational TPW product. We developed a framework for estimating atmospheric TPW from Visible Infrared Imaging Radiometer Suite (VIIRS) data in this study. First, the refined SWCVR method was used to retrieve TPW from the VIIRS TOA brightness temperature of channels M15 and M16. Then a more accurate and stable TPW product was derived by blending the VIIRS TPW and MIRS-derived TPW using BMA.
The empirical relationships between TPW and channel transmission ratio were established at 0 • , 15 • , 30 • , 45 • , 60 • , and 75 • viewing angles using a huge amount of representative samples produced by extensive radiate transfer modeling. The established formulae accounted for more than 97.6% of the variation in the simulation database, the MBE was 0, and the RMSEs ranged from 0.19 g/cm 2 to 0.22 g/cm 2 . Three years (2014-2017) of VIIRS data were used to retrieve TPW using the refined SWCVR. The retrieval results were validated using the ground measurements collected from SuomiNet sites over North America. The MBE and RMSE of VIIRS TPW were 0.21 g/cm 2 and 0.73 g/cm 2 , respectively. The accuracy of VIIRS TPW in daytime was much better than at night time. The VIIRS TPW was blended with the MIRS-derived TPW product via BMA. The accuracy of MIRS-derived TPW was better than that of VIIRS TPW, with MBE and RMSE values of 0.11 g/cm 2 and 0.38 g/cm 2 when validated by the same data used for VIIRS TPW. The MBE and RMSE of BMA integrated TPW were 0.06 g/cm 2 and 0.35 g/cm 2 , respectively, for which the accuracy difference in daytime and night time was removed. The global radiosonde measurements were also collected at the global scale to validate the BMA integrated VIIRS TPW. The MBE and RMSE of the BMA integrated TPW were 0.09 g/cm 2 and 0.44 g/cm 2 respectively compared to the radiosonde measurements, which was also more accurate than the VIIRS TPW. It was concluded that the developed framework can be used to derive accurate clear-sky TPW for VIIRS. This is the first time that we can obtain high accuracy TPW from VIIRS. This study will certainly benefit the study of atmospheric processes and climate change.