Validation of Sentinel-3 OLCI Integrated Water Vapor Products Using Regional GNSS Measurements in Crete, Greece

Water vapor is one of the essential variables in monitoring the Earth’s climate. The Ocean and Land Color Instrument (OLCI) on-board the Copernicus Sentinel-3 missions measures the Integrated Water Vapor (IWV) column over land and ocean surfaces. Post-launch calibration and validation of satellite measurements constitutes a key process in the operational phase of Earth observation satellites. This work presents the external and independent validation of OLCI-A IWV product using the regional network of continuously operating Global Navigation Satellite System (GNSS) comprised 10 stations distributed over the island of Crete in the eastern Mediterranean. The Sentinel-3A/-3B OLCI imagery that captures in a single scene the entire area of Crete has been examined. For each OLCI image, the IWV value of cloud-free pixels containing the GNSS stations have been derived and compared against simultaneous GNSS-derived measurements. The absolute as well as the relative bias between OLCI-A and OLCI-B IWV measurements have been determined. There is a good agreement between OLCI and GNSS with a bias of −0.57 mm ± 2.90 mm for OLCI(A) and +2.42 ± 3.41 mm for OLCI(B). The results of this regional validation activity are compared against other studies and the regular validation carried out at the Sentinel-3 Mission Performance Center. This work concludes that the accuracy of the OLCI IWV products is within its design requirements. The potential synergy between Sentinel-2 and Sentinel-3 IWV products is also discussed.


Introduction
Water vapor is one of the essential climate variables [1] which seems to have a strong impact on the Earth's hydrological cycle (evaporation, condensation and precipitation) and correlates to the planet's energy balance. The rise of the Earth's temperature increases evaporation from ocean and inland waters and subsequently builds up water vapor in the atmosphere. Although water vapor does not cause by itself global warming, it contributes to higher concentrations of water vapor in the atmosphere which in turn amplify the effect of temperature rise and eventually produce instability in the climate system.

Methodology and Dataset
The aim of this work is to compare the OLCI products for the IWV of Sentinel-3 against values derived from permanent GNSS stations on the ground. This external validation of OLCI products for IWV data is in agreement with the requirements for FRM redundancy. This section presents the methodology followed to determine the IWV values from OLCI and GNSS observations, the regional GNSS network, and the OLCI dataset implemented to perform this work.

Estimation of the IWV from OLCI
The OLCI multi-spectral imaging spectrometer on board the Copernicus Sentinel-3 family (Sentinels-3A and -3B are already operational whereas Sentinels-3C/-3D are to replace them in the near future) is based on the heritage of the ENVISAT's medium resolution imaging spectrometer (MERIS). Its main difference with its predecessor is that OLCI has 21 spectral bands, unlike the 15 bands on MERIS, the design of OLCI's five cameras mitigate sun-glint contamination, and it provides a resolution of 300 m over all target surfaces [21]. Each OLCI camera has its own charge couple device (CCD) composed of 384,800 imaging pixels (740 columns × 520 lines), resulting in 3700 ground pixels (740 pixels × 5 cameras) and covering a swath of 1270 km on the Earth's surface over a spectral range from 390 nm to 1040 nm [22].
Retrieval of water vapor from OLCI measurements relies upon the differential absorption technique [23] which allots measured radiance at a non-absorbing spectral band (O19 spectral band at 900 nm) and the reference water vapor absorption band (O18 spectral band at 885 nm). Although the measured radiance per spectral band is a function of solar irradiance, atmospheric transmittance and surface albedo (ignoring diffusion effects), following the differential absorption technique, it can be shown that IWV [kg/m 2 or mm] is given as [24]: where L 18 , L 19 is the measured spectral radiance at the O18 and O19 spectral bands, respectively, and k 19 is the mass extinction coefficient of water vapor at the O19 spectral band corresponding to the presumed vertical profile of temperature, pressure and water vapor over that wavelength. The algorithm used by Sentinel-3 to retrieve water vapor from OLCI measurements is described in [25].

Estimation of the IWV from GNSS Measurements
The signals transmitted by the GNSS satellites are refracted by the neutral atmosphere (troposphere) before reaching a ground GNSS receiver. This refraction introduces a signal propagation delay in the troposphere, called zenith total delay (ZTD), and estimated through processing (relative positioning, precise point positioning, etc.) of GNSS observations [26]. This ZTD value corresponds to the integrated signal delay caused by the zenith hydrostatic delay (ZHD) and the zenith wet delay (ZWD). The temporal and spatial variations of ZHD are small, contrary to the ZWD which change rapidly both in space and time [27].
The IWV may be calculated using GNSS-derived ZWD measurements using Equation (2) [28]: where IWV is the water vapor content in kg m −2 , ρ is the mass density of water (=1 gr cm −3 ), R V is the gas constant for water vapor (=0.4615 N·m·K −1 ·gr −1 ), C 1 (=3.776 × 10 5 K 2 ·hPa −1 ) and C 2 (=22.10 K hPa −1 ) are refractivity coefficients, T m is the weighted temperature of the atmosphere in • K, and ZWD is the delay caused by the wet troposphere in mm. The weighted temperature of the atmosphere is approximated with a 2% uncertainty for all weather condition by Equation (3) [16]: where T s is the surface temperature in • K.

The Regional GNSS Monitoring Network
A network of 15 permanent, continuously operating GNSS reference stations has been established and maintained by the Technical University of Crete, in western Crete for monitoring tectonic deformation and for supporting the operations of the ESA's Permanent Facility for Altimetry Calibration in Crete [29]. In this work, GNSS data from 10 reference stations have been processed ( Figure 1) using relative positioning implemented via the GAMIT scientific software (version 10.71) developed by MIT, USA [30]. Table A1 in Appendix A presents the positioning results (coordinates and velocities) for each GNSS station as computed using the latest (2014) realization of the International Terrestrial Reference Frame [31]. The time series of the daily coordinates of the CDN0 GNSS station, operating on a mountain in West Crete, is given in Figure 2.
Remote Sens. 2020, 12, x FOR PEER REVIEW  4 of 23 where is the surface temperature in °K.

The Regional GNSS Monitoring Network
A network of 15 permanent, continuously operating GNSS reference stations has been established and maintained by the Technical University of Crete, in western Crete for monitoring tectonic deformation and for supporting the operations of the ESA's Permanent Facility for Altimetry Calibration in Crete [29]. In this work, GNSS data from 10 reference stations have been processed ( Figure 1) using relative positioning implemented via the GAMIT scientific software (version 10.71) developed by MIT, USA [30]. Table A1 in Appendix A presents the positioning results (coordinates and velocities) for each GNSS station as computed using the latest (2014) realization of the International Terrestrial Reference Frame [31]. The time series of the daily coordinates of the CDN0 GNSS station, operating on a mountain in West Crete, is given in Figure 2.   The International Terrestrial Reference Frame 2014 (ITRF 2014) [31] has been here implemented and the GAMIT time series have been analyzed using the Quasi-Observation Combination Analysis software [32]. In addition to this regional GNSS network, extra observations from reference stations which are part of the International GNSS Service as well as the European Reference Frame networks, have been included in the database for processing. The GNSS observations have been carried out using precise satellite orbits, a 10 • elevation cut-off angle at a rate of 30 s. Solid-earth, polar motion and oceanic loading tides have been also applied [33,34].
In GNSS processing, the zenith wet delay is given as the difference between the zenith tropospheric delay and the zenith hydrostatic delay. The Saastamoinen model [35] has been used to provide an accurate estimate for the zenith hydrostatic delay as a function of the station atmospheric pressure (P) and latitude (ϕ) and orthometric height (H): Measurements of the atmospheric pressure have been made by collocated barometers and/or provided by meteorological models. The CDN0, CRS1 and TUC2 stations are equipped with meteorological sensors, while for the remaining GNSS stations the atmospheric pressure is estimated Remote Sens. 2020, 12, 2606 5 of 22 using atmospheric pressure loading models [36] and the Vienna Mapping Function 1 (VMF1) model [37]. The VMF1 model is applied in case of malfunction of an in-situ operating barometer.   Different input (i.e., VMF1 model and in-situ sensors in our case) impact the produced water vapor product. In order to quantify this contribution, the zenith wet delay (1-min rate records) using both the VMF1 model and the collocated meteorological sensors at the CDN0, CRS1 and TUC2 GNSS stations has been determined ( Figure 3). Apparently, this comparison has been carried out when the meteorological sensors at the respective sites have been operating. The Pearson correlation coefficient has been calculated to be 0.96, 0.98 and 0.98 for CDN0, CRS1 and TUC2 respectively, and the mean bias is at the sub-millimeter level. These correlation results along with the values of the R 2 coefficient (>0.92, statistically significant) justify the integration of the two atmospheric pressure data sources.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 23 The International Terrestrial Reference Frame 2014 (ITRF 2014) [31] has been here implemented and the GAMIT time series have been analyzed using the Quasi-Observation Combination Analysis software [32]. In addition to this regional GNSS network, extra observations from reference stations which are part of the International GNSS Service as well as the European Reference Frame networks, have been included in the database for processing. The GNSS observations have been carried out using precise satellite orbits, a 10° elevation cut-off angle at a rate of 30 s. Solid-earth, polar motion and oceanic loading tides have been also applied [33,34].
In GNSS processing, the zenith wet delay is given as the difference between the zenith tropospheric delay and the zenith hydrostatic delay. The Saastamoinen model [35] has been used to provide an accurate estimate for the zenith hydrostatic delay as a function of the station atmospheric pressure (P) and latitude (φ) and orthometric height (H): Measurements of the atmospheric pressure have been made by collocated barometers and/or provided by meteorological models. The CDN0, CRS1 and TUC2 stations are equipped with meteorological sensors, while for the remaining GNSS stations the atmospheric pressure is estimated using atmospheric pressure loading models [36] and the Vienna Mapping Function 1 (VMF1) model [37]. The VMF1 model is applied in case of malfunction of an in-situ operating barometer.
Different input (i.e., VMF1 model and in-situ sensors in our case) impact the produced water vapor product. In order to quantify this contribution, the zenith wet delay (1-min rate records) using both the VMF1 model and the collocated meteorological sensors at the CDN0, CRS1 and TUC2 GNSS stations has been determined ( Figure 3). Apparently, this comparison has been carried out when the meteorological sensors at the respective sites have been operating. The Pearson correlation coefficient has been calculated to be 0.96, 0.98 and 0.98 for CDN0, CRS1 and TUC2 respectively, and the mean bias is at the sub-millimeter level. These correlation results along with the values of the R 2 coefficient (>0.92, statistically significant) justify the integration of the two atmospheric pressure data sources. Figure 3. Scatterplots of the IWV (Precipitable Water) at the CDN0, CRS1 and TUC2 GNSS stations as determined using the VFM1 model (IWV_model) and in-situ measurements from meteorological sensors (IWV_meteo). It can be seen that the slope of the linear fit to the data is close to unity while its intercept is at the sub-mm level, with R 2 > 0.92 at the worst case.

The OLCI Dataset
The Sentinel-3A mission was launched on 16 February 2016 and Sentinel-3B was launched on 25 April 2018. Initially, Sentinel-3B was operating in tandem with Sentinel-3A and eventually reached its nominal orbit in December 2018. This work applies the full range of Sentinel-3A products, covering Figure 3. Scatterplots of the IWV (Precipitable Water) at the CDN0, CRS1 and TUC2 GNSS stations as determined using the VFM1 model (IWV_model) and in-situ measurements from meteorological sensors (IWV_meteo). It can be seen that the slope of the linear fit to the data is close to unity while its intercept is at the sub-mm level, with R 2 > 0.92 at the worst case.

The OLCI Dataset
The Sentinel-3A mission was launched on 16 February 2016 and Sentinel-3B was launched on 25 April 2018. Initially, Sentinel-3B was operating in tandem with Sentinel-3A and eventually reached its nominal orbit in December 2018. This work applies the full range of Sentinel-3A products, covering a period of about four years of operations (February 2016 to February 2020). For Sentinel-3B, only OLCI products covering the full year of 2019 (1 January to 31 December) have been evaluated at this stage.
The repeat cycle of Sentinel-3A/-3B satellites is 27 days. Nonetheless, OLCI permits global coverage in 2-3 days because of its large swath of 1270 km. This revisit time of OLCI is further reduced to about daily coverage with the two satellites in orbit.
Only OLCI images containing the entire island of Crete in a single scene have been taken into consideration for the present work. Specifically, twelve Sentinel-3A and Sentinel-3B Passes The repeat cycle of Sentinel-3A/-3B satellites is 27 days. Nonetheless, OLCI permits global coverage in 2-3 days because of its large swath of 1270 km. This revisit time of OLCI is further reduced to about daily coverage with the two satellites in orbit.
Only OLCI images containing the entire island of Crete in a single scene have been taken into consideration for the present work. Specifically, twelve Sentinel-3A and Sentinel-3B Passes (No. 7, 21, 64, 78, 121, 135, 178, 235, 278, 292, 335, and 349 of each mission) fulfill this coverage requirement ( Figure 4). The Sentinel-3 OLCI Level-2 Land Full Resolution (OL_2_LFR) product contains the IWV band. Full resolution stands for a 300-m pixel size on OLCI IWV imagery. The Copernicus Open Access Application Programing Interface Hub (API) and a custom-built script written in Python have been used to download the OLCI imagery, and in the end, to determine the IWV values.
Initially, the Sentinelsat module (https://pypi.org/project/sentinelsat/) has been used to execute a spatial query through the Copernicus API per relative pass and to identify the Sentinel products that fulfill the following criteria: (1) Mission: Sentinel-3, (2) Product: OL_2_LFR, and (3) Region: Crete as defined in the script with a polygon in GeoJson (https://geojson.org/) format. The identified Sentinel-3 products are automatically downloaded and archived as separate folders. Then, the IWV value, its error, geographical coordinates, date & time stamps, contained within each Sentinel-3 product folder, are examined.
In the sequel, for each GNSS station, the script estimates the nearest to it pixel in the OLCI image. Two approaches are followed for the determination of the OLCI-derived IWV values: (1) single pixel The Sentinel-3 OLCI Level-2 Land Full Resolution (OL_2_LFR) product contains the IWV band. Full resolution stands for a 300-m pixel size on OLCI IWV imagery. The Copernicus Open Access Application Programing Interface Hub (API) and a custom-built script written in Python have been used to download the OLCI imagery, and in the end, to determine the IWV values.
Initially, the Sentinelsat module (https://pypi.org/project/sentinelsat/) has been used to execute a spatial query through the Copernicus API per relative pass and to identify the Sentinel products that fulfill the following criteria: (1) Mission: Sentinel-3, (2) Product: OL_2_LFR, and (3) Region: Crete as defined in the script with a polygon in GeoJson (https://geojson.org/) format. The identified Sentinel-3 products are automatically downloaded and archived as separate folders. Then, the IWV Remote Sens. 2020, 12, 2606 7 of 22 value, its error, geographical coordinates, date & time stamps, contained within each Sentinel-3 product folder, are examined.
In the sequel, for each GNSS station, the script estimates the nearest to it pixel in the OLCI image. Two approaches are followed for the determination of the OLCI-derived IWV values: (1) single pixel and (2) average of pixels within an area of influence.
In the first approach, the IWV derived from GNSS is directly compared against the OLCI pixel value nearest to the location of the GNSS station. This approach has been recently followed to match the GNSS IWV observations with those derived from the IASI, MIRS, MODIS and MODIS-FUB satellites [38].
The second approach takes into account a GNSS area of influence, as the zenith total delay represents the average of all slant signal paths for the visible satellites at the troposphere layer height [39]. The size of this area of influence varies and depends on local conditions (i.e., visibility to GNSS satellites, variation in GNSS constellation geometry, etc.). In this work, a rectangular area of 31 × 31 pixels of OLCI (about 9.3 km 2 ) has been selected to represent the effective area of influence for the GNSS-derived IWV values. The center of this rectangle is located near in space or on the GNSS station ( Figure 5). This influence area is similar in size to that applied in [40]. The OLCI IWV value to be compared against the GNSS-derived IWV is the average of 961 pixels in the vicinity of the central pixel (GNSS station location). Only OLCI IWV values that are flagged as "LAND" (land surface and clear sky) have been taken into consideration for this averaging.
Remote Sens. 2020, 12, x FOR PEER REVIEW 7 of 23 the GNSS IWV observations with those derived from the IASI, MIRS, MODIS and MODIS-FUB satellites [38]. The second approach takes into account a GNSS area of influence, as the zenith total delay represents the average of all slant signal paths for the visible satellites at the troposphere layer height [39]. The size of this area of influence varies and depends on local conditions (i.e., visibility to GNSS satellites, variation in GNSS constellation geometry, etc.). In this work, a rectangular area of 31 × 31 pixels of OLCI (about 9.3 km 2 ) has been selected to represent the effective area of influence for the GNSS-derived IWV values. The center of this rectangle is located near in space or on the GNSS station ( Figure 5). This influence area is similar in size to that applied in [40]. The OLCI IWV value to be compared against the GNSS-derived IWV is the average of 961 pixels in the vicinity of the central pixel (GNSS station location). Only OLCI IWV values that are flagged as "LAND" (land surface and clear sky) have been taken into consideration for this averaging. Graphical illustration of the two approaches implemented to derive the OLCI IWV value per image per GNSS station. The image above shows such an example for the Sentinel-3B OLCI image acquired on 21 June 2020 and the RDK1 GNSS station used for validation. In the "single-point" approach, the IWV value of the pixel which contains the GNSS station (RDK1 in image) is used. In the "area-of-influence" approach, the average IWV value of an OLCI window of 30 × 30 pixels is calculated to be compared against the GNSS-derived IWV.
For the OLCI pixel nearest to GNSS station as well as for those pixels within the area of influence, the value and error of the IWV, the sensing date and time are from the related OLCI band records. All those are archived for further analysis. Finally, a file is created per Sentinel-3 orbit (Table 1), containing values and parameters, along with the satellite cycle and relative orbit number. For each Sentinel-3 OLCI IWV image, the IWV values and associated errors are obtained for that pixel (300-m size) which contains a GNSS station along with the 30 × 30 pixels in its vicinity. Table 2 presents the dataset (cycles) used per satellite orbit and mission. The cycles that contained no IWV values at the GNSS stations (pixel covered by clouds for example) are also identified. It can be seen Figure 5. Graphical illustration of the two approaches implemented to derive the OLCI IWV value per image per GNSS station. The image above shows such an example for the Sentinel-3B OLCI image acquired on 21 June 2020 and the RDK1 GNSS station used for validation. In the "single-point" approach, the IWV value of the pixel which contains the GNSS station (RDK1 in image) is used. In the "area-of-influence" approach, the average IWV value of an OLCI window of 30 × 30 pixels is calculated to be compared against the GNSS-derived IWV.
For the OLCI pixel nearest to GNSS station as well as for those pixels within the area of influence, the value and error of the IWV, the sensing date and time are from the related OLCI band records. All those are archived for further analysis. Finally, a file is created per Sentinel-3 orbit (Table 1), containing values and parameters, along with the satellite cycle and relative orbit number.
For each Sentinel-3 OLCI IWV image, the IWV values and associated errors are obtained for that pixel (300-m size) which contains a GNSS station along with the 30 × 30 pixels in its vicinity. Table 2 presents the dataset (cycles) used per satellite orbit and mission. The cycles that contained no IWV values at the GNSS stations (pixel covered by clouds for example) are also identified. It can be seen that at least one comparison point (GNSS station) exists in all Sentinel-3A and Sentinel-3B cycles.
This section briefly presented the way the integrated water vapor using OLCI and GNSS measurements are determined. The satellite and GNSS dataset used to perform the external validation of the OLCI-derived IWV values have been also described.

Results
This section compares the IWV values measured by the OLCI instrument with the respective values derived from the GNSS reference stations on the ground using both the "single point" and the "area-of-influence" approaches. Absolute calibration as well as relative validation of Sentinel-3A/-3B OLCI IWV products has been carried out. Absolute calibration is performed using the GNSS-derived IWV as reference, whereas relative validation refers to comparison of simultaneously captured Sentinel-3A and Sentinel-3B OLCI IWV products.

Absolute Bias of Sentinel-3 OLCI IWV Products
The differences of IWV values per GNSS station, as measured by OLCI and derived by GNSS stations using the "single point" approach have been calculated and used for the OLCI validation. Figure 6 presents the histograms of these differences for Sentinel-3A. Differences of IWV values have also been obtained using the "area-of-influence" approach.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 23 This section briefly presented the way the integrated water vapor using OLCI and GNSS measurements are determined. The satellite and GNSS dataset used to perform the external validation of the OLCI-derived IWV values have been also described.

Results
This section compares the IWV values measured by the OLCI instrument with the respective values derived from the GNSS reference stations on the ground using both the "single point" and the "area-of-influence" approaches. Absolute calibration as well as relative validation of Sentinel-3A/-3B OLCI IWV products has been carried out. Absolute calibration is performed using the GNSS-derived IWV as reference, whereas relative validation refers to comparison of simultaneously captured Sentinel-3A and Sentinel-3B OLCI IWV products.

Absolute Bias of Sentinel-3 OLCI IWV Products
The differences of IWV values per GNSS station, as measured by OLCI and derived by GNSS stations using the "single point" approach have been calculated and used for the OLCI validation. Figure 6 presents the histograms of these differences for Sentinel-3A. Differences of IWV values have also been obtained using the "area-of-influence" approach.  At this preliminary stage of analysis, histograms have been created to help us examine the shape of the statistical distribution of these differences. The shape of the histograms seems to follow approximately normal distribution, although measurements of physical variables, such as IWV, are not expected to give an ideal "bell-like" shape.
The "single point" approach could not be implemented at the RDK1 site, as the corresponding IWV error in the OLCI image pixel was high (>25 mm, whereas the expected IWV error is about 1-2 mm). This occurred because the pixel is flagged as "WATER", implying that the OLCI IWV retrieval algorithm considers that over a water surface. The RDK1 proximity to the sea seems to be irrelevant for this outcome as other coastal stations (i.e., CRS1, PALC) are properly flagged as land surfaces, thus allowing the OLCI algorithms to operate well over these stations (Figure 7). Remote Sens. 2020, 12, x FOR PEER REVIEW 10 of 23 At this preliminary stage of analysis, histograms have been created to help us examine the shape of the statistical distribution of these differences. The shape of the histograms seems to follow approximately normal distribution, although measurements of physical variables, such as IWV, are not expected to give an ideal "bell-like" shape.
The "single point" approach could not be implemented at the RDK1 site, as the corresponding IWV error in the OLCI image pixel was high (>25 mm, whereas the expected IWV error is about 1-2 mm). This occurred because the pixel is flagged as "WATER", implying that the OLCI IWV retrieval algorithm considers that over a water surface. The RDK1 proximity to the sea seems to be irrelevant for this outcome as other coastal stations (i.e., CRS1, PALC) are properly flagged as land surfaces, thus allowing the OLCI algorithms to operate well over these stations (Figure 7). (right) is flagged as "WATER" and "LAND" respectively. This introduces large errors in IWV determination at the RDK1 station thus they were excluded from further analysis when the "single point" approach was followed. Table 3 shows some statistical results for the determined differences of IWV in Sentinel-3A using both the "single point" and the "area-of-influence" approaches. The mean of the differences for the OLCI minus GNSS values for all validation sites lie within the range of −1.09 to +1.94 mm for the "single point" and −1.59 to +2.31 mm for the "area of influence" approach. The mean bias of Sentinel-3A OLCI IWV is calculated to be +0.57 mm ± 4.02 mm and −0.23 mm ± 3.16 mm under the "single point" and "area of influence" approaches, respectively. It can be seen that these different approaches have generated similar results for Sentinel-3 bias determination. The magnitude of this bias lies within the specification of the OLCI instrument, as well as within the accuracy of the IWV determination using GNSS processing. The OLCI image pixel corresponding to the coastal GNSS stations of RDK1 (left) and CRS1 (right) is flagged as "WATER" and "LAND" respectively. This introduces large errors in IWV determination at the RDK1 station thus they were excluded from further analysis when the "single point" approach was followed. Table 3 shows some statistical results for the determined differences of IWV in Sentinel-3A using both the "single point" and the "area-of-influence" approaches. The mean of the differences for the OLCI minus GNSS values for all validation sites lie within the range of −1.09 to +1.94 mm for the "single point" and −1.59 to +2.31 mm for the "area of influence" approach. The mean bias of Sentinel-3A OLCI IWV is calculated to be +0.57 mm ± 4.02 mm and −0.23 mm ± 3.16 mm under the "single point" and "area of influence" approaches, respectively. It can be seen that these different approaches have generated similar results for Sentinel-3 bias determination. The magnitude of this bias lies within the specification of the OLCI instrument, as well as within the accuracy of the IWV determination using GNSS processing. Table 3. The mean and standard deviation (SD) for the IWV values derived from Sentinel-3A OLCI using the "single point" (IWV sgl ) and the "area-of-influence" (IWV avg ) approaches, and from the GNSS, as well as for their differences at each GNSS station. All values are in [mm]. The overall mean and SD per IWV determination are given in bold. In order to investigate the existence of a correlation between the two independent techniques for the IWV estimation (OLCI and GNSS), scatterplots of Sentinel-3A OLCI and the GNSS-derived IWV values have been created for both the "single point" and "area of influence" approaches. For illustration purposes only, the "area of influence" scatterplots are shown in Figure 8. The "single point" scatterplots are given in Figure A1 in Appendix A. These figures also present a line fitted with regression on the results for each station.

Site
Only OLCI pixels flagged as "LAND "are used to calculate the average IWV value per OLCI image. Furthermore, it has been observed that the vast majority of IWV error of valid OLCI pixels range between 0.3-2.0 mm. Pixels with IWV error larger than 2.0 mm have been excluded for the remaining analysis. The offset and trend of the linear regression lines for both OLCI retrieval approaches are presented in Table 4. The Pearson linear correlation coefficient and its significance value (p value) is also given. It seems that significant (p-value less than the strict arbitrary significance level of 0.003) strong (Pearson coefficient >0.80) correlation between the two datasets exists.  In order to investigate the existence of a correlation between the two independent techniques for the IWV estimation (OLCI and GNSS), scatterplots of Sentinel-3A OLCI and the GNSS-derived IWV values have been created for both the "single point" and "area of influence" approaches. For illustration purposes only, the "area of influence" scatterplots are shown in Figure 8. The "single point" scatterplots are given in Figure A1 in Appendix A. These figures also present a line fitted with regression on the results for each station. Only OLCI pixels flagged as "LAND "are used to calculate the average IWV value per OLCI image. Furthermore, it has been observed that the vast majority of IWV error of valid OLCI pixels range between 0.3-2.0 mm. Pixels with IWV error larger than 2.0 mm have been excluded for the remaining analysis. The offset and trend of the linear regression lines for both OLCI retrieval approaches are presented in Table 4. The Pearson linear correlation coefficient and its significance Figure 8. Scatterplots of simultaneous IWV measurements obtained by the Sentinel-3A OLCI instrument as derived by the "area of influence" approach and the GNSS station at the respective station location. Table 5 presents the statistical analysis of the differences between the Sentinel-3B OLCI IWV measurements and the IWV values derived from GNSS observations. The mean bias for Sentinel-3B OLCI is +1.07 mm ± 4.80 mm using the "single point" approach and +0.24 mm ± 3.93 mm under the "area of influence" approach. Table 5. The mean and standard deviation (SD) for the IWI derived from Sentinel-3B OLCI, and GNSS as well as for their differences at the GNSS stations. All values are in [mm]. The overall mean and SD per IWV determination approach are also given in bold.  Figure 9 provides a graphical representation of the mean IWV values of OLCI(A) and OLCI(B) instruments per validation site and per OLCI retrieval approach along with the corresponding GNSS-derived IWV measurements. It can be clearly concluded that each independent technique (OLCI and GSNS processing) produces consistent results per validation site. Their difference is also, in most of the cases, similar although differentiation of the magnitude of the OLCI(A/B) absolute bias (−0.57 mm for OLCI(A) versus +2.42 mm for OLCI(B)) may be attributed to the different lengths of the available datasets.

Site
Remote Sens. 2020, 12, x FOR PEER REVIEW 14 of 23 Table 5. The mean and standard deviation (SD) for the IWI derived from Sentinel-3B OLCI, and GNSS as well as for their differences at the GNSS stations. All values are in [mm]. The overall mean and SD per IWV determination approach are also given in bold.  Figure 9 provides a graphical representation of the mean IWV values of OLCI(A) and OLCI(B) instruments per validation site and per OLCI retrieval approach along with the corresponding GNSSderived IWV measurements. It can be clearly concluded that each independent technique (OLCI and GSNS processing) produces consistent results per validation site. Their difference is also, in most of the cases, similar although differentiation of the magnitude of the OLCI(A/B) absolute bias (−0.57 mm for OLCI(A) versus +2.42 mm for OLCI(B)) may be attributed to the different lengths of the available datasets.   Figure 10 presents the scatterplots of Sentinel-3B OLCI IWV values against the GNSS observations per validation site. A near-perfect linear correlation between OLCI-B and GNSS IWV values is clearly visible from these scatterplots. Only the "area of influence" OLCI IWV values are presented as similar results exist for the "single point" approach. Table 6 provides the correlation between Sentinel-3B and GNSS IWV values.  Figure 10 presents the scatterplots of Sentinel-3B OLCI IWV values against the GNSS observations per validation site. A near-perfect linear correlation between OLCI-B and GNSS IWV values is clearly visible from these scatterplots. Only the "area of influence" OLCI IWV values are presented as similar results exist for the "single point" approach. Table 6 provides the correlation between Sentinel-3B and GNSS IWV values.  Figure 10 presents the scatterplots of Sentinel-3B OLCI IWV values against the GNSS observations per validation site. A near-perfect linear correlation between OLCI-B and GNSS IWV values is clearly visible from these scatterplots. Only the "area of influence" OLCI IWV values are presented as similar results exist for the "single point" approach. Table 6 provides the correlation between Sentinel-3B and GNSS IWV values.   The main conclusions of this first assessment for the integrated water values of the OLCI instrument in Sentinel-3A/-3B and with respect to GNSS-derived values from a Crete network, are:

IWVsgl
• The mean bias of OLCI(A) is determined as +0.57 mm and −0.23 mm implementing the "single point" and the "area of influence" approach respectively and using data covering a 4-year operational period; • The mean bias of OLCI(B) is determined as +1.07 mm and +0.24 mm implementing the "single point" and the "area of influence" approach respectively, over a one-year period of validation data; • The magnitude of OLCI(A) and OLCI(B) bias is within the operational capabilities of the satellite instrument but also the accuracy of GNSS-derived IWV values that have been used as reference in this study; • The OLCI and GNSS IWV values are significantly correlated irrespectively of validation sites and approach used.

Relative Bias of Sentinel-3 OLCI IWV Products
The relative bias between Sentinel-3A and Sentinel-3B can be estimated over simultaneous measurements made by both satellites and over the same target area on the ground. More specifically, the available dataset for Sentinel-3A, Sentinel-3B has been investigated to identify specific times at which the two OLCI instruments observe Crete GNSS sites together. Besides the simultaneity date, other criteria have been taken into consideration: (1) both OLCI(A) and OLCI(B) should produce valid IWV measurements, and (2) GNSS observations are available for both satellite overpasses. These criteria for concurrent OLCI(A) and OLCI(B) measurements have been met thirty-two times over the Crete GNSS reference stations.
In these cases, Sentinel-3A OLCI imagery has been captured at 08:12 UTC, while Sentinel-3B at 08:51 UTC. To compensate for any variation taken place in between these 39 min in the wet troposphere, the GNSS delays have been examined at the two specific times of satellite passes. The difference between the two GNSS-derived delays has been used to adjust the IWV value of OLCI(B) at 08:51 UTC, as if it had been measured at the same time of Sentinel-3A (08:12 UTC) (pseudo-simultaneity). Then, the value of OLCI(B) minus OLCI(A) has been reported as the relative bias in the integrated water vapor.
Following the analysis of this rather small sample of IWV values, the relative bias is estimated to be of the order of +3 mm. In other words, when compared to OLCI(A), the OLCI(B) overestimates integrated water vapor by +3 mm. Nonetheless, the derived results are based on a small sample and cannot be completely statistically significant.

Discussion
This work attempts to validate the integrated water vapor observed by OLCI in Sentinel-3 against GNSS measurements conducted through a permanent network. This is the first time that a regional GNSS network is implemented to provide an independent IWV source for assessing the OLCI performance. In order to cross examine the validation results previously presented, the official OLCI performance report as published by the Sentinel-3 Mission Performance Center (S3MPC) [41] has been considered.
The OLCI versus GNSS comparison at the S3MPC is carried out using the SUOMI-NET GNSS network operating mainly in North and Central America [41]. In the latest S3MPC report, the bias for OLCI(A) and OLCI(B) is reported as +1.24 and +1.02 mm, respectively. This S3MPC corresponds to Cycles 56 and 37 of Sentinel-3A and Sentinel-3B covering approximately the same period as the one implemented in the present investigation. Thus, it is possible to directly compare these S3MPC values for the OLCI(A) and OLCI(B) against those provided here, using the regional GNSS network in Crete ( Table 7). The difference of the two processing and analysis strategies is 0.78 to 1.47 mm, and is within the reported uncertainties. Table 7. Bias of the OLCI(A) and OLCI(B) integrated water vapor measurements as determined in this investigation, using the regional Crete GNSS network and the S3MPC analysis based on the SUOMI-NET GNSS network.

Site
Sentinel The S3MPC also reports an overestimation of the OLCI-derived IWV by about 11-12% with respect to the SUOMI-NET GNSS network. The present Tables 3 and 5 result in an overestimation of less than 10% over the specific geographical area. Finally, the S3MPC conclusion that both instruments of OLCI(A) and OLCI(B) present an equal behavior [41], is in agreement with the qualitative interpretation of Figure 7.
Further to the comparison of this study results with the S3MPC quality assessment, there are also some key findings that need to be described:

•
The OLCI pixel that contains the RDK1 station is flagged as "WATER". Thus, it is expected that the OLCI IWV retrieval algorithm will not perform well at this location. Even though implementation of the "area of influence" approach resolved this problem and made it possible to obtain valid OLCI IWV values in the vicinity of the RDK1 station, the RDK1 site presents the largest deviation among all examined GNSS stations. A few explanations for this deviation are possible: The geometry of the GPS satellites is worse (weak) compared to other GNSS sites, as high mountains in the North (see Figure 7) block the signal reception from satellites in that direction. This may be resolved if a multi-constellation (i.e., GPS, GLONASS and Galileo) receiver is operated at this location. Another possible explanation is the absence of in-situ meteorological sensor. That entails that the applied VMF1 model does not perform well in the RDK1 area and cannot sense the particularities that the specific geomorphology causes to the water vapor. • Two alternative techniques have been employed to determine the OLCI IWV value at a GNSS station location. Through the "area of influence" approach, it was possible to derive valid OLCI IWV values, even when the single OLCI pixel nearest to the GNSS station produced invalid IWV values. This is a valuable outcome when monitoring the OLCI IWV performance in coastal regions. The results presented in Tables 3 and 5 cannot help us arrive at a definite conclusion on the performance of each approach. In the cases of CDN0 and TUC2 (land only stations), the "area of influence" approach produces largest IWV values than the "single point" approach. The opposite stands true for the rest of the GNSS stations. The difference between the IWV values is within ±3 mm by these two approaches.

•
At present, there is no "golden rule" for defining a representative size for the area of influence in the IWV value. The actual size of the "area of influence" varies both spatially and temporally and depends on local conditions. In this work, we have chosen to use an area of 9.3 km 2 , which is close to the size of 10 km 2 applied in [41]. This size of 9.3 km 2 has been selected in our case to compare both results directly from the two studies. Nevertheless, a more detailed and independent investigation on the effective area of influence of GNSS water vapor products is needed to be carried out in the future.

•
The number of validation points to be used for the comparison analysis affects the reliability and the conclusions drawn out of a statistical analysis. Even though the results presented in this study are statistically significant (good p-value, adequate number of sampling points), a compromise to the quality of the IWV values may have been made. For example, when the "area of influence" approach was employed, the percentage of valid pixels within each window of the OLCI image has not been taken into consideration. This was made intentionally to enlarge the number of valid OLCI-GNSS sampling points for validation. Using stricter criteria (i.e., exclude an OLCI image window if it has less than 50% of its pixels flagged as "LAND", within the area of influence) would significantly reduce the number of validation points. Note that the global analysis carried out in [41] used only 5% of the total available products for their validation.

•
The continuous operation of the GNSS stations is a vital requirement for carrying out an effective validation. Because it has to be ensured that whenever a valid OLCI image is captured, there will be always ground measurements available to compare. In some cases, improper operation of GNSS stations limited comparison affecting the total number of validation points. Also, collocation of the GNSS stations with meteorological stations reduces the possibility of reverting to global models for the estimation of the atmospheric pressure.

Conclusions
This work presented for the first time an external validation of the integrated water vapor products of the OLCI instrument in Sentinel-3 against GNSS-derived results obtained by a regional GNSS network in Crete, Greece. The selection of the island of Crete has been made for three main reasons: It hosts a network of permanent facilities that support calibration and validation of satellite altimeters, such as, and not only, the SRAL instrument of Sentinel-3A and Sentinel-3B. Thus, this work promotes the synergy between instruments operating at the same satellite platform, as for the first time the same ground infrastructure is used for the calibration and validation of the microwave (altimeter) as well as the optical observations of Sentinel-3. The second reason is that Crete is covered by a relatively dense network of GNSS stations, continuously operating and for a long time (some more than 20 years) making it possible to calibrate the same OLCI image with several GNSS stations. The third reason is that several of these GNSS stations operate next to the sea, making them ideal for evaluating the performance of OLCI IWV retrieval algorithm in the sea-land boundary. This is a first attempt to estimate the absolute and relative bias of the Sentinel-3 OLCI instrument integrated water vapor measurements. It seems that the uncertainty of the produced IWV values is within OLCI specifications.
In summary, the absolute bias of OLCI(A) in integrated water vapor is closer to zero than that of the OLCI(B) [−0.57 mm for OLCI(A) versus +2.42 for OLCI(B)], but both are highly correlated with GNSS. This may be attributed to the relatively small dataset particularly for Sentinel-3B.
Furthermore, there is a need for accurate determination of the hydrostatic delays at the GNSS stations and the execution of an uncertainty budget analysis for GNSS-derived IWV values. This is a requirement posed by the FRM strategy for ground measurements involved in the calibration and validation of Copernicus satellite missions' products. The work presented in [27] constitutes a good starting point for such an uncertainty analysis. Also, the GNSS processing may be carried out using precise point positioning, which provides better temporal resolution of the wet tropospheric delays.
Future plans include continuation of OLCI-sensed IWV validation against the regional GNSS network including also images which capture part of Crete island, utilization of alternative water vapor estimation sources such as measurements obtained by the microwave radiometer that is collocated with the CDN0 GNSS station at the permanent facility for altimetry calibration and/or numerical weather models, and the comparison against Sentinel-2 water vapor products [41].

Conflicts of Interest:
The authors declare no conflict of interest.  Figure A1. Scatterplots of simultaneous IWV measurements obtained by the Sentinel-3A OLCI instrument as derived by the "single point" approach and the GNSS station at the respective station location.