Evaluation of Remote-Sensing Reﬂectance Products from Multiple Ocean Color Missions in Highly Turbid Water (Hangzhou Bay)

: Validation of remote-sensing reﬂectance (Rrs) products is necessary for the quantitative application of ocean color satellite data. While validation of Rrs products has been performed in low to moderate turbidity waters, their performance in highly turbid water remains poorly known. Here, we used in situ Rrs data from Hangzhou Bay (HZB), one of the world’s most turbid estuaries, to evaluate agency-distributed Rrs products for multiple ocean color sensors, including the Geostationary Ocean Color Imager (GOCI), Chinese Ocean Color and Temperature Scanner aboard HaiYang-1C (COCTS/HY1C), Ocean and Land Color Instrument aboard Sentinel-3A and Sentinel-3B, respectively (OLCI/S3A and OLCI/S3B), Second-Generation Global Imager aboard Global Change Observation Mission-Climate (SGLI/GCOM-C), and Visible Infrared Imaging Radiometer Suite aboard the Suomi National Polar-orbiting Partnership satellite (VIIRS/SNPP). Results showed that GOCI and SGLI/GCOM-C had almost no effective Rrs products in the HZB. Among the others four sensors (COCTS/HY1C, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP), VIIRS/SNPP obtained the largest correlation coefﬁcient (R) with a value of 0.7, while OLCI/S3A obtained the best mean percentage differences (PD) with a value of − 13.30%. The average absolute percentage difference (APD) values of the four remote sensors are close, all around 45%. In situ Rrs data from the AERONET-OC ARIAKE site were also used to evaluate the satellite-derived Rrs products in moderately turbid coastal water for comparison. Compared with the validation results at HZB, the performances of Rrs from GOCI, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP were much better at the ARIAKE site with the smallest R (0.77) and largest APD (35.38%) for GOCI, and the worst PD for these four sensors was only − 13.15%, indicating that the satellite-retrieved Rrs exhibited better performance. In contrast, Rrs from COCTS/HY1C and SGLI/GCOM-C at ARIAKE site was still signiﬁcantly underestimated, and the R values of the two satellites were not greater than 0.7, and the APD values were greater than 50%. Therefore, the performance of satellite Rrs products degrades signiﬁcantly in highly turbid waters and needs to be improved for further retrieval of ocean color components.


Introduction
Ocean color remote sensing provides abundant data for ocean-based research, but the quality of remote sensing data can be affected by many factors (radiometric accuracy of sensors, vicarious calibration, atmospheric correction, etc.) [1,2], and atmospheric correction from the GOCI based on in situ data of the East China Sea, Japan/East Sea, and Yellow Sea. Feng et al. (2020) [61] and Du et al. (2021) [62] used neural network algorithms and empirical algorithms based on UV bands to retrieve GOCI's turbidity, respectively, and GOCI-derived turbidity products were in good agreement with the in situ data. GOCI applies a new algorithm based on the near-infrared band, and the retrieved Rrs products are in good agreement with in situ data in the Western Pacific region [63]. Song et al. (2019) evaluated the alternative calibration coefficient estimation model of COCTS and achieved good results in short-term data verification [64]. Chen et al. (2020) evaluated the Rrs products of COCTS/HY1C using in situ data from the South Pacific Gyre region with good spatial homogeneity and high temporal stability and found that the uncertainty of Rrs products in the blue band was less than 5% [65]. Zibordi et al. (2018) evaluated the Lwn products of OLCI/S3A using in situ data from AERONET-OC and found systematic underestimation of the satellite Lwn products in the red and blue bands [66]. Tilstone et al. (2021) evaluated the performance of OLCI/S3A, OLCI/S3B, and VIIRS/SNPP in the Atlantic Ocean and found that both OLCI/S3A and OLCI/S3B underestimated chlorophylla (Chl-a), and the trend was greater for OLCI/S3A compared to OLCI/S3B [67]. Research on the consistency of OLCI products and in situ data [68][69][70][71][72][73][74][75][76][77][78], the performance comparison with other sensor's products [79][80][81][82][83][84][85][86][87][88], and the performance of OLCI products with different retrieval algorithms [80,83,86,[89][90][91][92] all reveal the great potential of OLCI. In addition, research on SGLI/GCOM-C has mainly concentrated on cloud and aerosols [93][94][95][96][97][98][99][100], with level-2 ocean color remote sensing products poorly studied [101,102]. In situ hyperspectral radiometric data measured at the Marine Optical Buoy (MOBY) [22] were routinely used to monitor VIIRS ocean color products [36,37,[103][104][105][106], and AERONET-OC data were also used for VIIRS ocean color products validation [34][35][36]87,88,107]. Other remote sensor products (such as MODIS, OLCI) were also used to evaluate the performance of VIIRS products. Overall, these validations highlight the potential of satellite ocean color products in the open ocean and shelf waters [108,109]. However, the performance of ocean color products in highly turbid waters remains unknown.
In this study, we used time-series in situ Rrs data from the highly turbid Hangzhou Bay (HZB) to evaluate the performance of Rrs products from the GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP. In Section 2, we provide detailed information on the in situ Rrs data measured in HZB and AERONET-OC station (ARIAKE site) and the satellite-derived Rrs products from multiple ocean color sensors. In Section 3, we provide the validation results of the Rrs products from multiple ocean color sensors in highly turbid waters. In Section 4, we compare and discuss the performance of the Rrs products at HZB and ARIAKE site and in highly and lowly-moderately turbid waters.

In Situ Rrs Data Measured in HZB
Hangzhou Bay, which narrows from east to west, is part of the Qiantang River estuary, adjacent to the north side of the Yangtze River estuary, as shown in Figure 1. In total, 4.8 × 10 8 t·y −1 of sediment is transported by the Yangtze River each year, some of which enters HZB [110] and 7.6 × 10 6 t·y −1 of sediment is transport by the Qiantang River, which is directly discharged into HZB [111]. The average water depth of HZB is 10 m, and the submarine topography rises from east to west, with an overall gradient of~0.06 m/km [112]. Narrowing of the bay and uplift of the seabed topography together increase the tidal range as the tide propagates to the top of the bay [112]. Under the influence of sediment input from the two rivers and the strong tides, the concentration of total suspended sediment (TSS) in HZB can reach 5000mg/L, and daily variation amplitude can reach 1000mg/L [110,111,113]. Therefore, HZB is an ideal region to examine the performance of ocean color products in highly turbid waters. where ρ is the sea surface reflecting coefficient, with value of 0.028 according to Mobley (1999) [114] and Dai et al., (2015) [115]. During the system operation, we carry out routine maintenance once a month, including data backup, sensor cleaning, and power supply inspection. When the system is running abnormally or there is extreme weather such as typhoons, we carry out special maintenance. As water spectral measurements can be influenced by platform shadow and sun glint at certain times, quality control of the data was first carried out. As per the methods proposed by Zibordi et al. (1999Zibordi et al. ( , 2002Zibordi et al. ( , 2004 [24,26,27,29] and Dai et al. (2015) [115], several conditions were adopted to select high-quality samples. Firstly, when measured Rrs is affected by platform shadow, the spectral shape changes rapidly [115], thus abnormally low Rrs spectra were removed. Secondly, we eliminated sun glint coefficients greater than 0.005 [116,117]. Thirdly, Rrs recorded under cloudy days were removed based on atmospheric diffuse transmittance. Actual atmospheric diffuse transmittance is denoted as and calculated as cos ⁄ , where θ is the solar zenith angle and is the extraterrestrial solar irradiance. Atmospheric diffuse transmittance under an ideal clear sky (assuming aerosol optical thickness of 0.3 corresponding to high aerosol load case) is denoted as λ and calculated by Gordon's approximate model [48]. Here, if 750 nm was less than 750 nm , it was considered as non-clear sky and the sample was removed [115]. Fourthly, He et al. (2012) found that A high-frequency water spectrum observation system at Hai-Tian-Yi-Zhou (HTYZ) site (121.12528 • E, 30.46278 • N) was constructed at the middle of the HZB Bridge, which is 16 and 19 km from the south and north coasts, respectively. Sea surface upward radiance (Lwater), downward sky radiance (Lsky), and downward irradiance (E s ) were measured by two hyperspectral radiance sensors and one hyperspectral irradiance sensor (Trios RAMSES), respectively. The azimuth angle for two radiance sensors was fixed about 139.64 • (referring to northward), and the observation zenith angles of Lsky and Lwater were 40 • and 140 • (referring to upward), respectively ( Figure 1). The sensors were calibrated before they were put into operation. The operating time of the system was from 07:00 a.m. to 17:00 p.m. local time with a step of 15-min. The measured Spectral data covered wavelengths from 320 to 950 nm with a spectral resolution of 3.3 nm. Finally, Rrs was estimated by: where ρ is the sea surface reflecting coefficient, with value of 0.028 according to Mobley (1999) [114] and Dai et al., (2015) [115]. During the system operation, we carry out routine maintenance once a month, including data backup, sensor cleaning, and power supply inspection. When the system is running abnormally or there is extreme weather such as typhoons, we carry out special maintenance. As water spectral measurements can be influenced by platform shadow and sun glint at certain times, quality control of the data was first carried out. As per the methods proposed by Zibordi et al. (1999Zibordi et al. ( , 2002Zibordi et al. ( , 2004 [24,26,27,29] and Dai et al. (2015) [115], several conditions were adopted to select high-quality samples. Firstly, when measured Rrs is affected by platform shadow, the spectral shape changes rapidly [115], thus abnormally low Rrs spectra were removed. Secondly, we eliminated sun glint coefficients greater than 0.005 [116,117]. Thirdly, Rrs recorded under cloudy days were removed based on atmospheric diffuse transmittance. Actual atmospheric diffuse transmittance is denoted as t(λ) and calculated as t(λ) = E s (λ)/F 0 (λ) cos θ, where θ is the solar zenith angle and F 0 is the extraterrestrial solar irradiance. Atmospheric diffuse transmittance under an ideal clear sky (assuming aerosol optical thickness of 0.3 corresponding to high aerosol load case) is denoted as t(λ) and calculated by Gordon's approximate model [48]. Here, if t(750 nm) was less than t(750 nm) , it was considered as non-clear sky and the sample was removed [115]. Fourthly, He et al. (2012) found that Lwn in highly turbid water is low at the ultraviolet band [13]; thus, here, if Lwn (350 nm) was more than 3 mW/ cm 2 ·µm·sr the sample was removed [115]. Fifthly, to eliminate samples under weak illumination, we also excluded Rrs data when the sun zenith angle was greater than 70 • . We used in situ Rrs data from the HTYZ site from July 2019 to July 2020, except for 20 days in December 2019 when effective data were missing due to abnormal power supply.

In Situ Rrs Data at ARIAKE Site
AERONET-OC can support satellite ocean color investigations and guarantees the consistency and accuracy of long-term data across sites through unified measurement, single calibration, and uniform processing code [31]. AERONET-OC sampling occurs from 08.00 a.m. to 16.00 p.m. local time, with a sampling interval of half an hour. The AERONET-OC products were converted from Lwn to Rrs using the formula Rrs = Lwn/F 0 , and the value of F 0 was taken from Thuiller et al. (2003) [118].
In this study, in situ OC data retrieved at specific AERONET-OC sites, namely ARIAKE TOWER (ARIAKE site), were used to evaluate the multi-source sensor data. The ARIAKE site (130.27195 • E, 33.10362 • N) is located in the Ariake Sea of Japan at an elevation of 15 m, about 5 km from the coast of Saga. The average depth of the Ariake Sea is 15 m [119]. At this depth the contribution of sea bottom reflection to the Rrs of the water surface is negligible [120]. The station started to provide data in February 2018. It is managed and operated by Saga University. The average annual total suspended sediment (TSS) value at the ARIAKE site is 3.92 mg/L (calculated using GOCI's TSS product from July 2019 to July 2020), which means that the optical properties of the water near the ARIAKE site are complex. In this study, Version 3 Level 1.5 Lwn_IOP data from ARIAKE from July 2019 to July 2020 were used.

Satellite Data
According to the operating time of the two in situ sites (HTYZ and ARIAKE), we obtained the Level-2A Rrs products from multiple satellite ocean color sensors. Specifically, GOCI Level-2A Rrs products at three noontime observations (11:28, 12:28, and 13:28 Beijing time) were obtained from the Korea Ocean Satellite Center, with six visible light wavelength bands (412, 443, 490, 555, 660, and 680 nm) and a spatial resolution of~500 m. The COCTS/HY1C Level-2A Rrs products were obtained from the National Satellite Oceanic Application Center of China, with six visible light wavelength bands (412, 443, 490, 520, 565, and 670 nm) and a spatial resolution of~1.1 km. The full resolution (FR, 300 m) Level-2 Rrs products of OLCI/S3A and OLCI/S3B were obtained from the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT), with 10 visible light wavelength bands (400, 412, 442, 490, 510, 560, 620, 665, 673, and 681 nm). OLCI uses the Baseline Atmospheric Correction (BAC) algorithm inherited from MERIS in Case-1 water and the Alternative Atmospheric Correction (AAC) algorithm for turbid and highly absorbing Case-2 water [20,21,[121][122][123]. We obtained the L2-NWLR products of SGLI/GCOM-C from the Japan Aerospace Exploration Agency (JAXA), with seven visible light wavelength bands (380, 412, 443, 490, 530, 565, and 670 nm) and a spatial resolution of 250 m. The Level-2 products of VIIRS/SNPP were obtained from the National Oceanic and Atmospheric Administration (NOAA), which uses the Multi-Sensor Level-1 to Level-2 processing system (MSL12) as an atmospheric correction algorithm. The VIIRS/SNPP has five visible light wavelength bands (410, 443, 486, 551, and 671 nm) and a spatial resolution of 750 m.

Matchup between Satellite and In Situ Rrs Data
A temporal window of ±30 minutes was used to match the satellite and in situ Rrs data at the HTYZ site due to the strong hydrodynamics at this site [110]. At the ARIAKE site, we adopted a ±1 h temporal window due to the lower hydrodynamics at this site [107].
For spatial matchup between the satellite and in situ Rrs data, average satellite Rrs in a 5 × 5 pixels box centered on the site was used to match the in situ Rrs, except for OLCI/S3A, OLCI/S3B, and SGLI/GCOM-C at HTYZ site. Because the two OLCI sensors and SGLI/GCOM-C have a high spatial resolution, their Rrs products at the HTYZ site could be significantly affected by the platform, thus we shifted the matchup box by four Remote Sens. 2021, 13, 4267 6 of 23 pixels to the northeast. We used several criteria to determine the representativeness of the matchups [27,31]. Specifically, if the number of pixels with effective Rrs in the box was less than half of the total pixel number in the box, then the sample was discarded. Moreover, the coefficient of variation (CV = σ/µ, where σ and µ are the standard deviation and mean values of Rrs in the box, respectively) was adopted to examine spatial homogeneity in the box [124]. If the CV was larger than 0.4, then the sample was discarded. If a satellite Rrs spectrum with all visible bands met the above spatial quality control conditions, it was considered to be an effective spectrum.
In addition to the spatiotemporal matchups, spectral matchups were also considered. The HTYZ site provides hyperspectral Rrs data from 320-950 nm, which matched all visible light wavelength bands of all six satellite sensors. However, the ARIAKE site only provides in situ Rrs data at eight bands (400-667 nm), and thus did not match all bands of the satellite sensors. We ignored wavelength differences if the center wavelengths between the satellite and in situ bands were less than 5 nm. For cases where the wavelength difference was greater than 5 nm and the wavelength of the satellite sensor band was within the in situ wavelength range (400-667 nm), we estimated the Rrs for the satellite sensor band by interpolating the in situ Rrs, otherwise there was no matchup (resulting in no effective matchups for SGLI/GCOM-C 380 nm, OLCI 681 nm, and GOCI 680 nm at the ARIAKE site). The number of all effective spectra and matchups by each remote sensor at HTYZ site and ARIAKE site are shown in Table 1.

Evaluation Methods
Based on the matchup samples, statistical analyses were carried out to evaluate the performance of the satellite Rrs using liner regression analysis. Specifically, average percentage difference (PD) and absolute percentage difference (APD) were calculated using: where x i and y i represent the in situ and satellite Rrs data, respectively, for a specific matchup sample. Clearly, PD represents bias between satellite and in situ data. A positive PD represents overestimation of the satellite Rrs data compared with the in situ data; a negative PD represents underestimation. In this paper, two standard deviation filtering procedures were applied to ensure that the statistical results obtained from the overall matchups were not skewed by a few extreme abnormal cases. Specifically, any matchups with APD i values greater than µ APD + 2σ APD , where µ APD and σ APD represent initial average and standard deviation of all matchups, respectively, were excluded from matchup comparison analysis [27,107]. SGLI/GCOM-C, and VIIRS/SNPP, respectively. The average curve of the in situ Rrs increased with the increase in wavelength from 380-670 nm. After a peak at 670 nm, Rrs decreased slightly with the increase in wavelength. Overall, the spectral shape of the satellite-derived Rrs was consistent with the in situ average Rrs, which peaked at 670 nm, except for the average Rrs of GOCI with a peak value at 555 nm. In terms of magnitude, the Rrs from in situ, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP at HTYZ ranged from 0.005 − 0.09 sr −1 , 0.023 − 0.050 sr −1 , 0.03 − 0.058 sr −1 , 0.02 − 0.068 sr −1 , 0.015 − 0.037 sr −1 , and 0.02 − 0.05 sr −1 at 670 nm, respectively. It should be noted that the Rrs from OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP had negative values at the blue band (412 nm), which was likely caused by the overestimation of aerosol scattering radiance in atmospheric correction. Figure 2 shows the spectral comparisons among all effective satellite-derived Rrs from GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, and SGLI/GCOM-C and in situ Rrs at the HTYZ site. In total, we identified 1509 effective samples of Rrs at HTYZ, and 10, 38, 72, 76, 7, and 80 effective Rrs samples for GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP, respectively. The average curve of the in situ Rrs increased with the increase in wavelength from 380-670 nm. After a peak at 670 nm, Rrs decreased slightly with the increase in wavelength. Overall, the spectral shape of the satellite-derived Rrs was consistent with the in situ average Rrs, which peaked at 670 nm, except for the average Rrs of GOCI with a peak value at 555 nm. In terms of magnitude, the Rrs from in situ, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP at HTYZ ranged from 0.005 0.09 sr , 0.023 0.050 sr , 0.03 0.058 sr , 0.02 0.068 sr , 0.015 0.037 sr , and 0.02 0.05 sr at 670 nm, respectively. It should be noted that the Rrs from OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP had negative values at the blue band (412 nm), which was likely caused by the overestimation of aerosol scattering radiance in atmospheric correction. Further comparisons of satellite and in situ Rrs data for each sensor are shown in Figure 3. When the wavelength was less than 555 nm, the average Rrs of the GOCI and in situ data showed similar increasing trends as wavelength increased, though their differences also increased significantly as wavelength increased. However, when the Further comparisons of satellite and in situ Rrs data for each sensor are shown in Figure 3. When the wavelength was less than 555 nm, the average Rrs of the GOCI and in situ data showed similar increasing trends as wavelength increased, though their differences also increased significantly as wavelength increased. However, when the wavelength was greater than 555 nm, average GOCI Rrs was inconsistent with the in situ Rrs spectrum, and the Rrs values by GOCI were significantly underestimated. The spectral shape of the average COCTS/HY1C Rrs was consistent with the in situ data, however COCTS/HY1C systematically underestimated Rrs. The average Rrs spectra of OLCI/S3A and OLCI/S3B were slightly lower than that of the in situ data. It should be noted that the performances of OLCI/S3A and OLCI/S3B were not identical. The σ values of OLCI/S3A decreased with increasing wavelength, while the σ values of OLCI/S3B hardly changed. As seen in Figures 2 and 3, the differences in average Rrs between OLCI/S3A and OLCI/S3B were most obvious at 681 nm. The average Rrs values of VIIRS/SNPP and SGLI/GCOM-C were lower than the in situ data in all bands. Although the average Rrs values of VIIRS/SNPP and SGLI/GCOM-C increased with wavelength, like the in situ data, the growth rate of the satellite-derived Rrs decreased significantly after 565 nm. Therefore, Remote Sens. 2021, 13, 4267 8 of 23 in the long-wavelength bands, the differences between the satellite (VIIRS/SNPP and SGLI/GCOM-C) and in situ Rrs values were greater.

Analysis of Spectral Consistency
OLCI/S3A and OLCI/S3B were slightly lower than that of the in situ data. It should be noted that the performances of OLCI/S3A and OLCI/S3B were not identical. The σ values of OLCI/S3A decreased with increasing wavelength, while the σ values of OLCI/S3B hardly changed. As seen in Figures 2 and 3, the differences in average Rrs between OLCI/S3A and OLCI/S3B were most obvious at 681 nm. The average Rrs values of VIIRS/SNPP and SGLI/GCOM-C were lower than the in situ data in all bands. Although the average Rrs values of VIIRS/SNPP and SGLI/GCOM-C increased with wavelength, like the in situ data, the growth rate of the satellite-derived Rrs decreased significantly after 565 nm. Therefore, in the long-wavelength bands, the differences between the satellite (VIIRS/SNPP and SGLI/GCOM-C) and in situ Rrs values were greater.

Time-Series Comparisons
As band settings of the different satellites are not completely consistent, we found several differences in the bands of different satellites in time-series analysis. At the HTYZ site with ultra-high turbidity, the Rrs value increased with increasing wavelength, so the minimum and maximum values were obtained at violet and red bands, respectively, and the increasing trend in Rrs changed ~560 nm. Thus, assessment of seasonal variations in the Rrs(λ) data can be determined using these bands. Therefore, we selected remote sensor bands at 412, 560, and 670 nm at the HTYZ sites. Figure 4 shows the time-series comparison of Rrs values between the satellitederived and in situ (HTYZ) data. Due to the high turbidity at HTYZ, only 10 effective Rrs records were available from GOCI in a year and had only seven effective Rrs records were available from SGLI/GCOM-C in a year (effective Rrs spectra were obtained in the

Time-Series Comparisons
As band settings of the different satellites are not completely consistent, we found several differences in the bands of different satellites in time-series analysis. At the HTYZ site with ultra-high turbidity, the Rrs value increased with increasing wavelength, so the minimum and maximum values were obtained at violet and red bands, respectively, and the increasing trend in Rrs changed~560 nm. Thus, assessment of seasonal variations in the Rrs(λ) data can be determined using these bands. Therefore, we selected remote sensor bands at 412, 560, and 670 nm at the HTYZ sites. Figure 4 shows the time-series comparison of Rrs values between the satellite-derived and in situ (HTYZ) data. Due to the high turbidity at HTYZ, only 10 effective Rrs records were available from GOCI in a year and had only seven effective Rrs records were available from SGLI/GCOM-C in a year (effective Rrs spectra were obtained in the summer when water was relatively clean). The remaining four sensors obtained better Rrs records at HTYZ, but the satellite-derived Rrs data were significantly lower than the in situ data. In addition, the in situ Rrs data showed significant seasonal changes, with high values in winter and low values in summer. However, the satellite-derived Rrs data showed no such variation. in situ data at HTYZ, and thus comparisons were not applicable for these two sensors. However, we found 20, 26, 32, and 51 matchups for COCTS/HY1C, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP, respectively (Table 2). Overall, satellite-derived Rrs values showed a significant underestimation with the increase in in situ Rrs values.

Scatter Plot Comparisons at HTYZ Site
Remote Sens. 2021, 13, x FOR PEER REVIEW 9 of 24 summer when water was relatively clean). The remaining four sensors obtained better Rrs records at HTYZ, but the satellite-derived Rrs data were significantly lower than the in situ data. In addition, the in situ Rrs data showed significant seasonal changes, with high values in winter and low values in summer. However, the satellite-derived Rrs data showed no such variation.  Figure 5 shows scatter-plot comparisons between satellite-derived and HTYZ in situ Rrs data for the six sensors (GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP) with all bands together. We found only one matchup between the GOCI and in situ data, and only two matchups between the SGLI/GCOM-C and in situ data at HTYZ, and thus comparisons were not applicable for these two sensors. However, we found 20, 26, 32, and 51 matchups for COCTS/HY1C, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP, respectively (Table 2). Overall, satellite-derived Rrs values showed a significant underestimation with the increase in in situ Rrs values.

Scatter Plot Comparisons at HTYZ Site
The regression lines of COCTS/HY1C, OLCI/S3A, and OLCI/S3B also showed a significant underestimation compared to the in situ data. With the increase in in situ Rrs values, the underestimation of satellite-derived Rrs increased significantly. VIIRS/SNPP also showed systematic underestimation compared to the in situ data, but the increase in underestimation was not obvious.
The correlation coefficients (R) of COCTS/HY1C, OLCI/S3A, OLCI/S3B, VIIRS/SNPP, and HTYZ (in situ) data were 0.57, 0.49, 0.45, and 0.70 for all bands, respectively ( Figure  5). The VIIRS/SNPP Rrs showed the highest R value. All single-band R values for COCTS/HY1C were negative ( Table 2). All single-band R values for the two OLCI sensors were negative, though the values of OLCI/S3A were greater than that of OLCI/S3B for all bands less than 620 nm, while the R values of OLCI/S3B were greater than that of OLCI/S3A for all bands greater than 620 nm ( Table 2). The single-band R values of VIIRS/SNPP were positive at 671 nm (but <0.1) but negative at all other bands (Table 2). 490 nm, the single-band APD of OLCI/S3A was greater than that of OLCI/S3B, and at wavelengths 510 nm, the single-band APD of OLCI/S3B was greater than that of OLCI/S3A, indicating that the performance of the two OLCI sensors was inconsistent.    The regression lines of COCTS/HY1C, OLCI/S3A, and OLCI/S3B also showed a significant underestimation compared to the in situ data. With the increase in in situ Rrs values, the underestimation of satellite-derived Rrs increased significantly. VIIRS/SNPP also showed systematic underestimation compared to the in situ data, but the increase in underestimation was not obvious.
The correlation coefficients (R) of COCTS/HY1C, OLCI/S3A, OLCI/S3B, VIIRS/SNPP, and HTYZ (in situ) data were 0.57, 0.49, 0.45, and 0.70 for all bands, respectively ( Figure 5). The VIIRS/SNPP Rrs showed the highest R value. All single-band R values for COCTS/HY1C were negative ( Table 2). All single-band R values for the two OLCI sensors were negative, though the values of OLCI/S3A were greater than that of OLCI/S3B for all bands less than 620 nm, while the R values of OLCI/S3B were greater than that of OLCI/S3A for all bands greater than 620 nm ( Table 2). The single-band R values of VIIRS/SNPP were positive at 671 nm (but <0.1) but negative at all other bands ( Table 2). These negative and low positive R values indicate poor performance of the satellite-derived Rrs in highly turbid waters.

Variation in APD with Water Turbidity
In this study, in situ Rrs values in the red-light wavelength band (670 nm) were used to characterize water turbidity and examine the performance of multi-sensor Rrs products under different turbidities, as shown in Figure 6. Overall, the APD values at 412, 443, 490, and 670 nm from VIIRS/SNPP and COCTS/HY1C increased rapidly with increasing water turbidity, indicating the degradation of Rrs product performance for these two sensors with increasing turbidity. Moreover, VIIRS/SNPP did not obtain effective data when Rrs (670 nm) was larger than 0.055 sr −1 , whereas COCTS/HY1C, OLCI/S3A, and OLCI/S3B obtained effective Rrs data, even Rrs values (670 nm) up to 0.08 sr −1 .   When the in situ Rrs (670 nm) was less than 0.055 sr −1 , all four remote sensors obtained effective matchups (Figure 7). Compared with Figure 5, the overall statistical results for the four common bands among sensors showed that VIIRS/SNPP had systematic underestimation with PD value of −36.00%, while there was overestimation with PD values of 4.55% and 26.72% for the OLCI/S3A and OLCI/S3B, respectively. The APD values were comparable with values of 50.26%, 56.73%, and 48.92% for the VIIRS/SNPP, OLCI/S3A, and OLCI/S3B, respectively. Due to only having eight matchups, it is hard to objectively evaluate the performance of COCTS/HYIC when Rrs (670 nm) is less than 0.055 sr −1 at HTYZ site.    Figure  5, the underestimation of all the three sensors has increased. Specifically, the underestimation of COCTS/HY1C has increased significantly; in contrast, the underestimation of the two OLCIs has increased slightly. Overall, none of the four sensors can accurately retrieve the water spectrum in the HZB.

Different Performances in High and Moderate Turbidity Waters
As a comparison, we used the in situ Rrs data from the ARIAKE site to examine the performance of the Rrs products from the multiple sensors under moderately turbid water. Figure 9 shows the average Rrs spectrum comparison between the satellite-derived and in situ Rrs at the ARIAKE site. The numbers of the effective Rrs samples from in situ measurement, GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP were 282, 392, 58, 72, 77, 82, and 163, respectively. Overall, the spectrum shape of Rrs derived by all six sensors was consistent with that of the in situ data. The OLCI/S3A, OLCI/S3B, and VIIRS/SNPP data were the closest to the in situ data in both spectral shape and magnitude. Of note, GOCI underestimated Rrs at 555 nm, and SGLI/GCOM-C and COCTS/HY1C underestimated Rrs at all bands. In addition, COCTS/HY1C showed slight deviation in spectral shape, with two peaks at 490 and 565 nm. For the 412, 443, and 490 nm bands, with increasing water turbidity, the APD values for OLCI/S3A and OLCI/S3B increased at first, then showed a decreasing trend to a minimum value, and then increased again. Considering that neural network atmospheric correction algorithm is a non-linear algorithm, the change trend of OLCI's APD may be caused by atmospheric correction. The neural network atmospheric correction model of OLCI is based on the simulated data generated by HYDROLIGHT radiative transfer code plus a bio-optical model relating scattering and absorption coefficients to material concentrations [20,50,51]. The data of water optical properties are mainly from cruises in the North Sea, partly in the Baltic Sea, Mediterranean Sea, and North Atlantic [20,50,51]. These optical properties may not fully cover the water optical types in the Hangzhou Bay with extreme turbidity. Thus, the Rrs products of OLCI/S3A and OLCI/S3B may have higher accuracy under certain water turbidity. At the 670 nm band, the APD values of all four sensors increased with increasing water turbidity.

Different Performances in High and Moderate Turbidity Waters
As a comparison, we used the in situ Rrs data from the ARIAKE site to examine the performance of the Rrs products from the multiple sensors under moderately turbid water. Figure 9 shows the average Rrs spectrum comparison between the satellite-derived and in situ Rrs at the ARIAKE site. The numbers of the effective Rrs samples from in situ measurement, GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP were 282, 392, 58, 72, 77, 82, and 163, respectively. Overall, the spectrum shape of Rrs derived by all six sensors was consistent with that of the in situ data. The OLCI/S3A, OLCI/S3B, and VIIRS/SNPP data were the closest to the in situ data in both spectral shape and magnitude. Of note, GOCI underestimated Rrs at 555 nm, and SGLI/GCOM-C and COCTS/HY1C underestimated Rrs at all bands. In addition, COCTS/HY1C showed slight deviation in spectral shape, with two peaks at 490 and 565 nm.

Different Performances in High and Moderate Turbidity Waters
As a comparison, we used the in situ Rrs data from the ARIAKE site to examine the performance of the Rrs products from the multiple sensors under moderately turbid water. Figure 9 shows the average Rrs spectrum comparison between the satellite-derived and in situ Rrs at the ARIAKE site. The numbers of the effective Rrs samples from in situ measurement, GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP were 282, 392, 58, 72, 77, 82, and 163, respectively. Overall, the spectrum shape of Rrs derived by all six sensors was consistent with that of the in situ data. The OLCI/S3A, OLCI/S3B, and VIIRS/SNPP data were the closest to the in situ data in both spectral shape and magnitude. Of note, GOCI underestimated Rrs at 555 nm, and SGLI/GCOM-C and COCTS/HY1C underestimated Rrs at all bands. In addition, COCTS/HY1C showed slight deviation in spectral shape, with two peaks at 490 and 565 nm. In contrast to obtaining almost no effective Rrs data at the highly turbid HTYZ site, GOCI and SGLI/GCOM-C obtained more effective Rrs data at the ARIAKE station, despite significant underestimation. COCTS/HY1C also showed significant underestimation at HTYZ and ARIAKE, and the differences between the COCTS/HY1C and in situ data at ARIAKE were relatively larger at the short-wavelength bands (<565 nm). Similar to the results for HTYZ site, the performances of OLCI/S3A and OLCI/S3B at ARIAKE were also slightly different. The OLCI/S3A data were slightly lower than the in situ data, while the OLCI/S3B data were closer to the in situ data. Overall, the difference in σ between OLCI/S3A and OLCI/S3B at ARIAKE was significantly smaller than the difference at HTYZ. The average Rrs value of SGLI/GCOM-C was markedly lower than the in situ value, even under the moderately turbid waters at ARIAKE. Compared to the significant underestimation at HTYZ, the average Rrs value of VIIRS/SNPP at ARIAKE matched the in situ data quite well. Figure 10 shows the comparisons between satellite-derived and in situ Rrs data for each sensor at ARIAKE. The number of matchups between satellite-derived and in situ Rrs data at ARIAKE was 154, 27,22,29,30, and 48 for the GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP, respectively (Table 3). Overall, when the in situ Rrs was less than 0.007 sr −1 , GOCI overestimated Rrs values, or else, it underestimated Rrs values. COCTS/HY1C and SGLI/GCOM-C were systematically lower than the in situ data. In contrast, the Rrs values of OLCI/S3A, OLCI/S3B, and VIIRS/SNPP were relatively consistent with the in situ data, thus showing much better performance than that at the HTYZ site.
OLCI/S3B data were closer to the in situ data. Overall, the difference in σ between OLCI/S3A and OLCI/S3B at ARIAKE was significantly smaller than the difference at HTYZ. The average Rrs value of SGLI/GCOM-C was markedly lower than the in situ value, even under the moderately turbid waters at ARIAKE. Compared to the significant underestimation at HTYZ, the average Rrs value of VIIRS/SNPP at ARIAKE matched the in situ data quite well. Figure 10 shows the comparisons between satellite-derived and in situ Rrs data for each sensor at ARIAKE. The number of matchups between satellite-derived and in situ Rrs data at ARIAKE was 154, 27,22,29,30, and 48 for the GOCI, COCTS/HY1C, OLCI/S3A, OLCI/S3B, SGLI/GCOM-C, and VIIRS/SNPP, respectively (Table 3). Overall, when the in situ Rrs was less than 0.007 sr , GOCI overestimated Rrs values, or else, it underestimated Rrs values. COCTS/HY1C and SGLI/GCOM-C were systematically lower than the in situ data. In contrast, the Rrs values of OLCI/S3A, OLCI/S3B, and VIIRS/SNPP were relatively consistent with the in situ data, thus showing much better performance than that at the HTYZ site.  Table 3 shows the performances of the multi-sensors at the ARIAKE site. The R values for OLCI/S3A, OLCI/ S3B, and VIIRS/SNPP (all bands) were greater than 0.80 and  Table 3 shows the performances of the multi-sensors at the ARIAKE site. The R values for OLCI/S3A, OLCI/S3B, and VIIRS/SNPP (all bands) were greater than 0.80 and much higher than the values at HTYZ (0.49, 0.45, and 0.70, respectively). Similarly, the R value for COCTS/HY1C was 0.70 at ARIAKE, which was much higher than that at HTYZ.
At HTYZ, the R values for each band were mostly negative. In contrast, the R values for each band at ARIAKE were all positive. Except for COCTS/HY1C and SGLI/GCOM-C, the R values of the other four sensors in the red band at ARIAKE were greater than 0.7 ( Table 3).
The comparison of the scatter points of the four common bands shared by the multisource remote sensor was shown in Figure 11. Compared with the total PD and APD of all bands as shown in Figure 10, the PD and APD of all four common bands slightly increased. The statistical results showed that the accuracy of OLCI/S3A, OLCI/S3B, and VIIRS/SNPP were close (APD are 32.33%, 31.76%, and 31.77%, respectively). The deviation of OLCI/S3B (PD = −4.82%) was the smallest, and the deviation of OLCI/S3A (PD = −20.07%) was slightly higher than that of VIIRS/SNPP (PD = −16.43%). At ARIAKE site, the statistical results of four common bands were slightly different from that of all bands. However, at the HTYZ site, the total statistical values of the four common bands were quite different from the statistical values of all bands. This means that in moderately turbid water, the various bands of the multi-source sensors can capture the water spectrum well. VIIRS/SNPP were close (APD are 32.33%, 31.76%, and 31.77%, respectively). The deviation of OLCI/S3B (PD = −4.82%) was the smallest, and the deviation of OLCI/S3A (PD = −20.07%) was slightly higher than that of VIIRS/SNPP (PD = −16.43%). At ARIAKE site, the statistical results of four common bands were slightly different from that of all bands. However, at the HTYZ site, the total statistical values of the four common bands were quite different from the statistical values of all bands. This means that in moderately turbid water, the various bands of the multi-source sensors can capture the water spectrum well. In summary, at ARIAKE site, GOCI significantly overestimated Rrs at 410 and 443 nm, and even though it underestimated Rrs at the long-wavelength (>443 nm) bands, the overall PD value indicated overestimation. In addition, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP showed relatively low R values at the short-wavelength (<486 nm) bands, but overall performance was good. COCTS/HY1C and SGLI/GCOM-C significantly underestimated Rrs at all bands. OLCI/S3A performed better than OLCI/S3B at HTYZ, while OLCI/S3B had a smaller deviation and showed better performance than OLCI/S3A at ARIAKE.

Analysis of Possible Reasons for the Performance of Multi-Source Remote Sensing Products at HTYZ Site
Most of the GOCI Rrs products released by KOSC are masked in the Hangzhou Bay. The atmospheric correction method for GOCI in turbid waters uses adjacent nonturbid water to obtain water-leaving reflectance in the red-light band [125], and then uses empirical relationships between NIR band and red-light band to obtain water-leaving reflectance at the NIR bands as [126]: where p w (745), p w (660), and p w (865) are the water-leaving reflectance at 745, 660, 865 nm respectively, and j n and k n are empirical coefficients. This empirical relationship is only applicable to waters with TSS less than 200 mg/L [126], while the TSS concentration in HZB is as higher as 5000 mg/L [127], which may cause the failure of GOCI Rrs products in Hangzhou Bay. Similar to GOCI, COCTS/HY1C has no SWIR band and its atmospheric correction adopts the framework of Gordon and Wang (1994) [48,108,128], and the overestimation of the aerosol contribution of NIR bands may be the reason for its poor performance at HTYZ site.
The neural network atmospheric correction model of OLCI in turbid waters is based on the simulated data generated by HYDROLIGHT radiative transfer code plus a bio-optical model relating scattering and absorption coefficients of material concentrations [20,50,51]. The data of water optical properties are mainly from cruises in the North Sea, partly in the Baltic Sea, Mediterranean Sea, and North Atlantic [20,50,51]. These optical properties may not fully cover the water optical types in the Hangzhou Bay with extreme turbidity.
It is well known that atmospheric correction based on SWIR band can improve the performance of remote sensing products in turbid waters [40]. SGLI/GCOM-C's SWIRalgorithm uses one NIR band (865 nm) and one SWIR band (1630 nm) to estimate the aerosols scattering radiance, but the water-leaving reflectance, in HZB, at 865 nm of is still not negligible [129]. As a result, the aerosol contribution may be overestimated, which led to an underestimation of water-leaving reflectance. VIIRS/SNPP uses SWIR bands of 1238 and 1601 nm to fulfil atmospheric correction in turbid waters [130]. Similarly, in HZB, due to the non-negligible water-leaving radiance at the 1238 nm [129], the aerosol contribution can also be overestimated. Another possible reason is that the use of SWIR bands for atmospheric correction increases noise errors [63,131]. This phenomenon is more obvious in the blue band that is farther from the SWIR band.
Therefore, atmospheric correction algorithms for highly turbid waters should be developed further. In highly turbid waters, water-leaving radiance at ultraviolet (UV) band might much less than that at the visible-light (VIS) bands or even the NIR bands [13,110]. Therefore, the UV band can be used to estimate aerosol scattering radiance in such cases [13]. Although GOCI and COCTS/HY1C do not have UV band, it is still feasible to use the 412 nm band instead of the UV band [13,110]. Singh et al. (2019) proposed an atmospheric correction algorithm based on a combination of UV band (387 nm), blue-light band (443 nm), and NIR bands (named as UVNIR-ex), and obtained a good performance in the Yangtze River Estuary [14]. Therefore, in extremely turbid water, for remote sensors without SWIR bands, such as GOCI and COCTS/HY1C, an atmospheric correction algorithm based on the UV band is an option. Hieronymi et al. (2017) proposed an algorithm consisting of several blended neural networks that were specialized for 13 different optical water classes, and the maximum concentration of inorganic suspended matter was up to 1500 mg/L [91]. If the in situ data of water types like HZB are also added to the neural network model training, the atmospheric correction in extremely turbid water may achieve better performance. In addition, in high turbid water, since water-leaving radiance at 1240 nm is still not negligible [129], using longer wavelength SWIR bands combination, such as 1635-2209 nm and 1601-2257 nm may improve the performance of atmospheric correction for SGLI/GCOM-C and VIIRS/SNPP. However, considering that the extrapolation from SWIR band to VIS band is too far, resulting in increased noise errors in the blue band [63,131,132], the atmospheric correction by combining UV and SWIR bands may provide a choice for extremely turbid waters.

Conclusions
In this study, we evaluated the Rrs product performance of multiple satellite ocean color sensors in the highly turbid waters of HZB in comparison to in situ data. Results showed that the Rrs products of GOCI and SGLI/GCOM-C contained limited effective Rrs data at the HTYZ site. COCTS/HY1C, OLCI/S3A, OLCI/S3B, and VIIRS/SNPP significantly underestimated Rrs at the HTYZ site and showed poor consistency with the in situ data, with PD values of -35.19%, -20.16%, -25.30%, and -34.87%, respectively, and APD values of 44.91%, 47.36%, 46.12%, and 44.48%, respectively. Among these four sensors, VIIRS/SNPP obtained the largest R value, while OLCI/S3A obtained the best PD value.