Evaluation of Four Atmospheric Correction Algorithms for GOCI Images over the Yellow Sea

Atmospheric correction (AC) for coastal waters is an important issue in ocean color remote sensing. AC performance is fundamental in retrieving reliable water-leaving radiances and then bio-optical parameters. Unlike polar-orbiting satellites, geostationary ocean color sensors allow high-frequency (15–60 min) monitoring of ocean color over the same area. The first geostationary ocean color sensor, i.e., the Geostationary Ocean Color Imager (GOCI), was launched in 2010. Using GOCI data acquired over the Yellow Sea in summer 2017 at three principal overpass times (02:16, 03:16, 04:16 UTC) with ±1 and ±3 h match-up times, this study compared four GOCI AC algorithms: (1) the standard near infrared (NIR) algorithm of NASA (NASA-STD), (2) the Korea Ocean Satellite Center (KOSC) standard algorithm for GOCI (KOSC-STD), (3) the diffuse attenuation coefficient at 490 nm Kd (490)-based NIR correction algorithm (Kd-based), and (4) the Management Unit of the North Sea Mathematical Models (MUMM). The GOCI-estimated remote sensing reflectance (Rrs), aerosol parameters [aerosol optical thickness (AOT), Angström Exponent (AE)], and chlorophyll-a (Chla) were validated using in situ data. For Rrs, AOT, AE, and Chla, GOCI-retrieved results performed well within the ±1 h temporal window, but the number of match-ups was extended within the ±3 h match-up window. For ±3 h GOCI-derived Rrs, all algorithms had an absolute percentage difference (APD) at 490 and 555 nm of <40%, while other bands showed larger differences (APD > 60%). Compared with in situ values, the APD of the Rrs(490)/Rrs(555) band ratio was <20% for all ACs. For AOT and AE, the APD was >40% and >200%, respectively. Of the four algorithms, the KOSC-STD algorithm demonstrated satisfactory performance in deriving Rrs for the region of interest (Rrs APD: 22.23%–73.95%) in the visible bands. The Kd-based algorithm worked well obtaining Ocean Color 3 GOCI Chla because Rrs(443) is more accurate than the KOSC-STD. The poorest Rrs retrievals were achieved using the NASA-STD and the MUMM algorithms. Statistical analysis indicated that all methods had optimal performance at 04:16 UTC.


Introduction
Ocean color remote sensing with daily-hourly sampling frequency and broad spatial coverage plays a critical role in the investigation of the bio-optical properties and the biogeochemical parameters of nearshore coastal waters. Although they represent only 7% of the total ocean surface, coastal and inland waters produce up to 40% of marine and freshwater biomasses inventoried today and 85% of marine and freshwater resources exploited by humans. Moreover, 60% of the world's population lives within 100 km of a coast, whilst inland waters provide key ecosystem services with direct linkages to human health [1]. These waters, which are very often optically complex, are generally identified as The four AC algorithms used in this study are briefly described in the following subsections. They each extend the GW-AC approach to derive Lw from TOA reflectance for turbid waters.

NASA Standard Algorithm (NASA-STD)
The NASA-STD AC algorithm was initially developed by Gordon and Wang [6], extended for application to turbid waters by Stumpf et al. [19], and revised by Bailey et al. and Ahmad et al. [28,34]. The latter revision is used by default in the present SeaDAS package version 7.X. The algorithm is based on an iterative process that accounts for the non-zero ρ w in the NIR bands. This method uses 80 aerosol models built from AERONET observations and vector radiative transfer code for the ocean-atmosphere system [34]. First, the black-pixel assumption is adopted using GW-AC to retrieve ρ w at 443 and 555 nm. Next, these two ρ w values are used as input for a bio-optical model (standard OBPG Chla OC3 algorithm [35], formula (1) in Figure 1) to obtain initial estimates of the chlorophyll-a (Chla) concentration, which then makes further efforts to determine particulate and CDOM absorption in the red band (formula (2) in Figure 1), a(660). Then, a(660) and ρ w (660) are used to compute the particulate backscattering in the red band, b bp (660), following which, b bp (NIR) can be estimated in accordance with a power exponent function [36] (formula (3) in Figure 1). Therefore, values of ρ w in the NIR bands (or equivalently Rrs) are generated on the basis of these three relationships, and they are removed from ρ rc (NIR). Ultimately, the procedure is repeated until convergence is reached. NASA uses this algorithm to generate its official GOCI L2 products.

Kd-Based NIR Correction Algorithm (Kd-Based)
Wang et al. [13] developed an algorithm specifically for processing GOCI ocean color data acquired in the highly turbid western Pacific region, including the Bohai Sea, the Yellow Sea, and the East China Sea. It is based on the regional relationship between nLw(745), nLw(865), and the diffuse attenuation coefficient at 490 nm [Kd(490)], which is derived from long-term MODIS-Aqua measurements (2002-2009) using NIR-based ocean color data processing (formula (1) in Figure 1) [38]: Wang et al. [13] developed an algorithm specifically for processing GOCI ocean color data acquired in the highly turbid western Pacific region, including the Bohai Sea, the Yellow Sea, and the East China Sea. It is based on the regional relationship between nLw(745), nLw(865), and the diffuse attenuation coefficient at 490 nm [Kd(490)], which is derived from long-term MODIS-Aqua measurements (2002-2009) using NIR-based ocean color data processing (formula 1 in Figure 1) [38]: where c1 = 0.465, c2 = −0.385, c3 = 0.152, and c4 = −0.0121. Similarly, nLw(865) can be formulated as (formula 2 in Figure 1): where b1 = 0.368 and b2 = 0.040. Then, an iterative process is conducted to calculate nLw(745) and nLw(865) with the inputs of Kd(490) (Figure 1). The fourth algorithm uses an analytical method that is referred to as the Management Unit of the North Sea Mathematical Models (MUMM) [30]. This algorithm replaces the invalid black-pixel The fourth algorithm uses an analytical method that is referred to as the Management Unit of the North Sea Mathematical Models (MUMM) [30]. This algorithm replaces the invalid black-pixel assumption in the NIR bands over coastal waters with two alternative assumptions-one concerns the water optical properties, and the other concerns the atmosphere.
The first assumption originates from the fact that the shape of the NIR water spectrum is dominated primarily by pure water absorption and hence is invariant, while the magnitude of the signal is approximately proportional to the backscatter coefficient [30,39], which can be regarded as constant in the region of interest. The ratio of any two values of NIR water-leaving reflectance, named α, is defined as (formula (1) in Figure 1): The second assumption arises from the fact that the aerosol concentration does not usually vary over spatial scales of about 100 km. Therefore, the ratio of multiple-scattering aerosol reflectance ρ a (λ) + ρ ra (λ), named ε, can be considered constant over the region of interest. For turbid waters with clear waters nearby, ε can also be calculated using the values of ρ rc (λ) in the clear waters as follows (formula (2) in Figure 1): One of the key points for the MUMM algorithm is to determine the values of α and ε. Here, α is set to 1.932 in accordance with Ruddick et al. [30]. To estimate ε(745,865) for each image, we produced scatterplots of ρ rc (745) versus ρ rc (865) for the region of interest, and we calculated the slope of the line for values of ρ rc (865) < 0.015. Using the set value of α and the estimated values of ε, the following equations are defined: The values of ρ A (745) and ρ A (865) are estimated using Equations (9) and (10) to select appropriate aerosol models. Finally, the determined aerosol models are reentered into the GW-AC scheme.

Chla Retrievals
The accuracy of AC algorithms determines the precision of the extraction of ocean color parameters (e.g., Chla). To retrieve Chla, different types of algorithm (e.g., OC2, OC3, OC4, and GSM01) were adopted in previous studies [40][41][42][43]. In this study, considering the specific spectral bands of the GOCI sensor, Chla concentrations were derived from all atmospherically corrected images using the OC3G algorithm (Ocean Color 3 GOCI) [37]. OC3G is an empirically derived algorithm developed as an extension of OC4 and OC2 [43]. The general form of the OC3G algorithm can be expressed as: R = log 10 max(Rrs(443), Rrs(490)) Rrs(555) with

Study Area
The Yellow Sea (YS) of China is the largest marginal sea of the northwestern Pacific Ocean ( Figure 2). The YS lies within a shallow basin (average water depth:~44 m) that is bounded by the Chinese mainland to the west and the Korean Peninsula to the east [44]. It covers an area of 417,000 km 2 [45]. One of the main characteristics of the YS is that the sea temperature and the salinity have prominent spatiotemporal diurnal variations [46,47]. The hydrologic and the circulation processes are governed by monsoon wind systems and the Kuroshio Current [30]. Because of the effects of these dominant processes, the Cold Water Mass, which is an important phenomenon within the YS, is formed with characteristics of low temperature (5-12 • C) and high salinity (31.5-32.5 psu) beneath a strong pycnocline [48], which is more prevalent during summer. Anthropogenic inputs such as pollution, eutrophic materials, and substantial sediment transport from the Yangtze River mean the YS region is well known for its high turbidity [49]. are governed by monsoon wind systems and the Kuroshio Current [30]. Because of the effects of these dominant processes, the Cold Water Mass, which is an important phenomenon within the YS, is formed with characteristics of low temperature (5-12 °C) and high salinity (31.5-32.5 psu) beneath a strong pycnocline [48], which is more prevalent during summer. Anthropogenic inputs such as pollution, eutrophic materials, and substantial sediment transport from the Yangtze River mean the YS region is well known for its high turbidity [49].

In Situ Measurements
In situ bio-optical and atmospheric optical parameters were acquired during a sea cruise in the YS during 8-26 August 2017, which was conducted by the National Ocean Technology Center and the National Satellite Ocean Application Service. The data fields are exhibited in Figure 2. Overall, measurements were taken at 92 stations. Observations of remote sensing reflectance, atmospheric optical properties, and significant water constituents were acquired during the field experiment.

In Situ Measurements
In situ bio-optical and atmospheric optical parameters were acquired during a sea cruise in the YS during 8-26 August 2017, which was conducted by the National Ocean Technology Center and the National Satellite Ocean Application Service. The data fields are exhibited in Figure 2. Overall, measurements were taken at 92 stations. Observations of remote sensing reflectance, atmospheric optical properties, and significant water constituents were acquired during the field experiment.

Measurement of Remote Sensing Reflectance
At stations with suitable conditions of solar illumination (generally between 09:00 and 15:00 local time), above-water optical measurements were conducted using a field spectrophotometer (ASD FieldSpec®3, full range: 350-2500 nm). The total upwelling radiance from the surface water (L t ), downwelling irradiance above the water surface (E s ), and downward sky radiance (L sky ) were measured. To avoid sunglint contamination, the zenith and the azimuth angles used to observe Lt were about 40 • and 135 • (referring to the solar plane), respectively. Furthermore, we selected the optimal orientation to minimize the influence of ship shading and whitecaps. Then, the remote sensing reflectance (Rrs, sr −1 ) was calculated as follows [50]: where L w+ is the Lw just above the sea surface, and ρ s is the Fresnel reflectance at the air-water interface, which depends on viewing and illumination geometry, wind speed, cloud, and wavelength [51,52]. For this study, ρ s (λ) was retrieved using the nonlinear spectral optimization method and a bio-optical model for each station [53,54].

Measurement of Aerosol Optical Properties
Aerosol optical thickness (AOT) is defined as the integrated extinction coefficient over a vertical column of unit cross-section. It is a proxy for the concentration of aerosols within the air column. To validate the results of aerosol optical properties derived from the four AC algorithms, shipborne AOT was measured using a hand-held Microtops II sun photometer (SolarLight, USA) at five central wavelengths (i.e., 440, 500, 675, 870, and 1020 nm), denoted as AOT(440), AOT(500), AOT(675), AOT(870), and AOT(1020), respectively ( Figure 3). For convenience of evaluation between satellite and field data, spline interpolation of AOT(500) and AOT(675) was applied to obtain AOT(555). The method for the retrieval of AOT from direct solar irradiance measurements is described by [50,55]. For details regarding cloud-screening and quality control procedures, the reader is referred to [56]. Additionally, the Angström exponent, AE, is commonly used to provide basic information on the particle size distribution and the type of aerosols. In this study, in situ AE between 440 and 870 nm, i.e., AE(440,870), was determined using linear regression with log-transformed in situ AOT(440) and AOT(870) as follows: Additionally, the Angström exponent, AE, is commonly used to provide basic information on the particle size distribution and the type of aerosols. In this study, in situ AE between 440 and 870 nm, i.e., AE(440,870), was determined using linear regression with log-transformed in situ AOT(440) and AOT(870) as follows:

Measurement of Chla Data
In situ surface Chla concentrations were used for validation of the imagery-derived Chla products of the four AC algorithms. Chla samples (n = 92) were collected at the surface using a continuous shipboard laboratory pump. The surface water samples were filtered using 0.7 μm Whatman GF/F glass fiber filters following Ocean Optics protocols [57]. Chla samples were stored at -20 °C until the samples were processed. Chla pigment concentrations were extracted using a high-performance liquid chromatography (HPLC) system. The derived mean values and the standard deviations of in situ Chla were 0.78 and 0.51, respectively.

GOCI Data
The Geostationary Ocean Color Imager (GOCI), the world's first geostationary ocean color spaceborne instrument, was launched in June 2010 and provides eight images per day. It is one of the three payloads onboard the Communication, Ocean and Meteorological Satellite (COMS).
The GOCI instrument was designed to provide hourly data in eight bands in the visible and the NIR parts of the spectrum (412-865 nm) with 500 m spatial resolution. The region of GOCI observations, which covers an area of 2500 km × 2500 km (21.54°-46.99°N, 116.41°-148.67°E, centered at 36°N, 130°E), includes the coast of East China, the Korean Peninsula, and Japan [ Figure  2a] [58].

Measurement of Chla Data
In situ surface Chla concentrations were used for validation of the imagery-derived Chla products of the four AC algorithms. Chla samples (n = 92) were collected at the surface using a continuous shipboard laboratory pump. The surface water samples were filtered using 0.7 µm Whatman GF/F glass fiber filters following Ocean Optics protocols [57]. Chla samples were stored at −20 • C until the samples were processed. Chla pigment concentrations were extracted using a high-performance liquid chromatography (HPLC) system. The derived mean values and the standard deviations of in situ Chla were 0.78 and 0.51, respectively.

GOCI Data
The Geostationary Ocean Color Imager (GOCI), the world's first geostationary ocean color spaceborne instrument, was launched in June 2010 and provides eight images per day. It is one of the three payloads onboard the Communication, Ocean and Meteorological Satellite (COMS).
The GOCI instrument was designed to provide hourly data in eight bands in the visible and the NIR parts of the spectrum (412-865 nm) with 500 m spatial resolution. The region of GOCI observations, which covers an area of 2500 km × 2500 km ( Figure 2a] [58]. The solar and the satellite viewing zenith angles of GOCI are reasonably small (<30 • ) between 02:16 and 04:16 UTC; therefore, available daily GOCI L1B images during this period and corresponding to the days of the sea cruise were downloaded from KOSC (http://kosc.kiost.ac.kr/eng/p30/kosc_p34.html). The GOCI images were processed from L1B to L2 for the KOSC algorithms using GDPSv2.0 and using the SeaDAS software package version 7.5 (SeaDASv7.5, OBPG http://oceancolor.gsfc.nasa.gov/) for NASA and MUMM algorithms.

Match-Ups Procedures
Match-ups between the in situ and the GOCI-retrieved AOT and Rrs were selected based on locations and overpass times. The slight difference in the wavelengths between the in situ and the satellite-retrieved values was ignored. The match-up criteria adopted were similar to the approach described in Bailey and Werdell [59]. Match-up time-windows of ±1 h and ±3 h were compared for AOT, AE, Rrs, and Chla. First, 3 × 3 pixel boxes were extracted from the GOCI images centered on the measurement sites. Second, a coefficient of variation (CV; standard deviation divided by mean values) was calculated for each band to account for the spatial homogeneity of the pixels within each 3 × 3 box. Match-ups with CV values >0.2 in the 3 × 3 pixel boxes for Rrs(555) were excluded. We required that at least five valid pixels within the 3 × 3 pixel box be valid. Finally, the mean value of the remaining pixels was calculated.

Statistical Metrics
Statistical metrics were used to evaluate the four algorithms performance between in situ Microtops II aerosol products versus GOCI aerosol products, in situ Rrs versus GOCI Rrs products, and in situ Chla versus GOCI Chla products. Statistical parameters included the absolute percentage difference (APD, Equation (14)), the root mean square error (RMSE, Equation (15)), the Bias (Equation (16)), the correlation coefficient (R 2 ), and the slope and the intercept of the linear regression: where X i is the i-th in situ observation, Y i is the i-th GOCI observation, and N is the number of match-ups between the in situ measurements and the GOCI-retrieved values. These metrics represent unbiased statistics that reflect how accurately the GOCI-derived values agree with the field data. Both APD and RMSE are sensitive to outliers, Bias indicates the deviation level of the GOCI-retrieved values from the in situ data, and the R 2 value reflects the linear consistency between the in situ data and the GOCI measurements, which is related to data distributions. For this study, statistical significance was defined at the 95% confidence level.

Results
The criterion for GOCI cloud masking is different in GDPS than in SeaDAS [60]. Minor differences in flags between GDPS and SeaDAS mean the number of match-up pairs varies among the AC algorithms. We analyzed time-windows of ±1 and ±3 h for Rrs, AOT, AE, and Chla between the shipborne data and the GOCI overpass to choose the best time-window. Within a ±3 h match-up window, of the 78 available in situ Rrs measurements, 30, 29, 31, and 32 match-ups for Rrs(555) were available for the NASA-STD, the KOSC-STD, the Kd-based, and the MUMM algorithms, respectively.
The number of match-ups for each algorithm was greater with a time-window of ±3 h than with a time-window of ±1 h (i.e., 10, 10, 10, and 11 match-ups, respectively). The number of match-ups for AOT with ±3 and ±1 h time-windows was similar to Rrs(555).

Comparison of Rrs
Statistical results for GOCI-derived Rrs values from the four AC methods versus in situ data are given in Tables 1-4, and corresponding scatterplots are shown in Figure 4. The method providing the least accurate Rrs is MUMM, which has the highest uncertainties for all bands. The NASA-STD method performs slightly better than the MUMM algorithm, as reflected in the improved APD and the RMSE ranging from 27.58%-94.55% and 0.0003-0.0028 sr −1 , respectively, in the VIS bands. The Kd-based algorithm is similar to the KOSC-STD method. Satisfactory results are obtained for the KOSC-STD algorithm, which shows the most accurate performance for a time-window of ±3 h with values of APD and Bias of 22.08%-73.95% and 0.0008-0.0019 sr −1 , respectively (for a time-window of ±1 h: APD is 12.37%-100.53% and Bias is −2 × 10 −6 to 1 × 10 −3 sr −1 ), and slopes closer to the 1:1 relationship. Overall, the spectrum shapes and the magnitudes at all sites for the KOSC-STD algorithm are much more consistent with in situ Rrs (not shown here).
Rrs(443) 7 (12) Rrs(745) 6 (11) 22.08 (17.16) 69.79 68.92 (100.53) 103.27 (95.65) 140.56 (125.78)   The difference between the GOCI-derived and the field-based Rrs data as a function of wavelength among the three GOCI overpass times is shown in Figure 5. The APD at 04:16 UTC is lower than at 02: 16   The difference between the GOCI-derived and the field-based Rrs data as a function of wavelength among the three GOCI overpass times is shown in Figure 5. The APD at 04:16 UTC is lower than at 02:16 and 03:16 UTC. The retrievals at Rrs(555) are the most accurate for the three overpass times for all four AC methods. The values of RMSE and Bias generally decrease with increasing wavelength. The lowest RMSE at 04:16 UTC occurs with the KOSC-STD method. Bias decreases quickly as wavelength increases from 412 to 490 nm, although it is nearly stable at wavelengths beyond 490 nm for all overpass times; this is especially true at 04:16 UTC.   To develop bio-optical algorithms, band ratios of remote-sensing reflectance, i.e., Rrs(443)/Rrs(555) and Rrs(490)/Rrs(555), are commonly used to retrieve biogeochemical parameters [46,47]. The band ratios for each algorithm were evaluated against in situ band ratios (Table 5 and Figure 6). The distribution of the ratio Rrs(490)/Rrs(555) is much more concentrated than Rrs(443)/Rrs(555), and the former is closer to the 1:1 line (Figure 6). Within a ±3 h match-up time, the ratio Rrs(490)/Rrs(555) exhibits a relatively better performance (APD: 11%-19%, RMSE: 0.272-0.309 sr −1 ) in comparison with Rrs(443)/Rrs(555) for all four methods. The ratio Rrs(490)/Rrs(555) retrieved by the Kd-based algorithm has the best consistency with the in situ data, whereas Rrs(443)/Rrs(555) has the poorest performance for the KOSC-STD algorithm. To develop bio-optical algorithms, band ratios of remote-sensing reflectance, i.e., Rrs(443)/Rrs(555) and Rrs(490)/Rrs(555), are commonly used to retrieve biogeochemical parameters [46,47]. The band ratios for each algorithm were evaluated against in situ band ratios (Table 5 and Figure 6). The distribution of the ratio Rrs(490)/Rrs(555) is much more concentrated than Rrs(443)/Rrs(555), and the former is closer to the 1:1 line (Figure 6). Within a ±3 h match-up time, the ratio Rrs(490)/Rrs(555) exhibits a relatively better performance (APD: 11%-19%, RMSE: 0.272-0.309 sr −1 ) in comparison with Rrs(443)/Rrs(555) for all four methods. The ratio Rrs(490)/Rrs(555) retrieved by the Kd-based algorithm has the best consistency with the in situ data, whereas Rrs(443)/Rrs(555) has the poorest performance for the KOSC-STD algorithm.
In conclusion, the ratio Rrs(490)/Rrs(555) estimated by all four AC methods results in better retrievals, which means this band ratio is more suitable for retrieving the biogeochemical properties of the YS region in summer. In comparison with the three single band values, i.e., Rrs(443), Rrs(490), and Rrs(555) (Tables 1-4), the band ratios reduce the systematic uncertainty of the AC procedure. Furthermore, the band ratios also improve the strength of the correlations relative to the single bands.    In conclusion, the ratio Rrs(490)/Rrs(555) estimated by all four AC methods results in better retrievals, which means this band ratio is more suitable for retrieving the biogeochemical properties of the YS region in summer. In comparison with the three single band values, i.e., Rrs(443), Rrs(490), and Rrs(555) (Tables 1-4), the band ratios reduce the systematic uncertainty of the AC procedure. Furthermore, the band ratios also improve the strength of the correlations relative to the single bands.

Comparison of AOT and AE
Aerosol optical products are by-products derived from the AC algorithms. To understand the uncertainties of aerosol products for AC schemes, it is necessary to analyze the features of aerosols. Many studies have investigated the wavelength dependence of AOT and the AE. The accuracy of the retrievals of AOT(443), AOT(555), AOT(680), AOT(865), and AE(443,865) of our current study is presented in Tables 6-9 and Figure 7. Summary statistics for GOCI versus field data aerosol products of ±3 and ±1 h are also shown in Table 3       To analyze the error between the retrieved and the in situ AOT (at 443, 555, 680, and 865 nm) at the three GOCI overpass times, Figure 8 illustrates the APD, the RMSE, and the Bias per time and per method between the satellite retrieval values and the field-based AOT. The spectral behavior of APD for AOT is similar to that of Rrs ( Figure 5, upper row). The performance of each algorithm is best at 04:16 UTC and worst at 02:16 UTC. This is consistent with the results obtained for Rrs. For all approaches, the APD in the visible bands is stable over the spectrum with a peak at 865 nm (around 50% at 04:16 UTC and near 100% at 02:16 UTC, Figure 8a-d). The sign of the bias is the same for the four AC methods, i.e., negative values in blue bands and positive values in red and NIR bands. The  The results of the AOT data show that the four methods exhibit underestimation at AOT bands. AOT(443) and AOT(865) retrievals present the highest RMSEs and relative errors, while AOT(555) performs best in the YS. All retrieved AE values from the four algorithms are <2.0, with the majority in the range 1.2-1.6. This indicates a size distribution dominated by fine particle sizes, which are considered absorbing aerosol types (coastal and pollution aerosols) [34,[61][62][63]. The distribution of the retrieved values around the 1:1 line are more scattered for AE(443,865) than AOT, which is consistent with previous published articles [13,14].
To analyze the error between the retrieved and the in situ AOT (at 443, 555, 680, and 865 nm) at the three GOCI overpass times, Figure 8 illustrates the APD, the RMSE, and the Bias per time and per method between the satellite retrieval values and the field-based AOT. The spectral behavior of APD for AOT is similar to that of Rrs ( Figure 5, upper row). The performance of each algorithm is best at 04:16 UTC and worst at 02:16 UTC. This is consistent with the results obtained for Rrs. For all approaches, the APD in the visible bands is stable over the spectrum with a peak at 865 nm (around 50% at 04:16 UTC and near 100% at 02:16 UTC, Figure 8a-d). The sign of the bias is the same for the four AC methods, i.e., negative values in blue bands and positive values in red and NIR bands. The time of image acquisition affects the magnitude of the errors but not the shape. The aerosol optical properties generated from the four methods can be generalized as follows: (1) the MUMM algorithm produces the largest uncertainties in retrievals for both AOT and AE(443,865), the NASA-STD performs slightly better than MUMM, the performance of KOSC-STD is similar to the Kd-based algorithm, and the Kd-based algorithm obtains the most accurate values; (2)

Comparison of OC3G Chla Retrievals
GOCI-OC3G Chla retrievals from the four AC methods versus in situ Chla concentrations were compared for the two different match-up windows (±1 and ±3 h) (Table 10)

Comparison of OC3G Chla Retrievals
GOCI-OC3G Chla retrievals from the four AC methods versus in situ Chla concentrations were compared for the two different match-up windows (±1 and ±3 h) (Table 10)

Discussion
This research compared four AC algorithms (NASA-STD, KOSC-STD, Kd-based, and MUMM) using in situ aerosol optical properties (AOT and AE) as well as above-water Rrs and Chla concentration. These four algorithms extend the GW-AC algorithm to Case-2 waters for which the black-pixel assumption cannot be considered. The features of GOCI-derived Rrs, aerosol optical

Discussion
This research compared four AC algorithms (NASA-STD, KOSC-STD, Kd-based, and MUMM) using in situ aerosol optical properties (AOT and AE) as well as above-water Rrs and Chla concentration.
These four algorithms extend the GW-AC algorithm to Case-2 waters for which the black-pixel assumption cannot be considered. The features of GOCI-derived Rrs, aerosol optical products, and Chla concentrations derived using the four algorithms can be summarized as follows. (1) Match-up time-windows of ±1 and ±3 h were compared for AOT, Rrs, and Chla (Tables 1-4) to obtain the best balance in achieving the highest number of match-ups and reducing the water mass and particle dynamics. Overall, all methods had better performance with a ±1 h time-window. However, the number of match-ups was too low, thus the match-up time-window was extended to ±3 h to ensure a suitable number of match-ups (Figures 3-8); (2) all methods overestimated Rrs and underestimated AOT data in comparison with in situ values; (3) for visible bands, all algorithms performed well for Rrs and AOT at 555 nm, but larger uncertainties appeared at other wavelengths (especially 412 nm) primarily because of incorrect estimation of the NIR ocean contributions [4]. For NIR bands, the error is larger than for the visible bands; (4) GOCI-derived Rrs and AOT were sensitive to the image acquisition time. The errors were lower for images acquired at 04:16 UTC in comparison with the two other times (02:16 and 03:16 UTC) because the solar and the satellite angles were lower at 04:16 UTC; (5) for GOCI-OC3G Chla retrievals, the Kd-based and the KOSC-STD algorithms had lower uncertainties for Rrs, AOT, and AE than the other two algorithms. However, the strengths and the weaknesses of each approach affected the accuracy of the AC procedure. The discrepancy could possibly have been caused by uncertainties in the NIR correction models, the aerosol models, or GOCI slot effects.
The essential distinction between the four AC algorithms is the NIR correction model adopted. In our study and with a limited dataset, the KOSC-STD approach appears the method most appropriate for Rrs retrieval, whereas the Kd-based algorithm is slightly better than KOSC-STD for estimation of Chla over the YS region. For the KOSC-STD algorithm, the Lw at 660 nm can be calculated, then the GW-AC algorithm iteratively updates the fields of Lw in the NIR bands. This method has been validated in turbid coastal regions such as the Korean Peninsula [3,27,64]. However, the KOSC-STD algorithm has a further limitation, i.e., the polynomial equation of Lw between the red-NIR bands is variable depending on the concentrations of suspended particulate matter, CDOM, and Chla [62,65]. Similarly, for the Kd-based method, which is a regional model for the GOCI overpass area, Kd(490) plays an important role in the iterative computation of Rrs in the NIR bands. On the one hand, the relationship between Kd(490) and NIR nLw might not fit for turbid waters in the short term. On the other hand, nLw(660) might be insensitive to changes of nLw(745) in extremely turbid waters [66]. In other words, nLw(660) would become constant, meaning that Kd(490) derived from nLw(660) could not be used to estimate NIR nLw for turbid waters with Kd(490) > 5.0 m −1 [13]. The NASA-STD algorithm establishes the relations between Rrs and Chla indirectly, and the YS region is characterized by high concentrations of total suspended matter and CDOM. Therefore, the iterative NIR correction might not be a complete fit for the YS region. The MUMM method that assumes spatial homogeneity of the ratio of Rrs and ρ A might not always be suitable for the YS region in summer [9]. Some studies have shown that the ratio of Lw is not usually constant in the NIR bands for turbid coastal waters [65]. This ratio can vary in the range 1.0-2.0 with a change in water turbidity [67], which could be described as a function of the absorption coefficient of pure water in the NIR bands and the suspended particulate matter [16,68].
The aerosol models selected by the four algorithms are different. Aerosol concentrations and components are difficult to determine because of their high spatiotemporal variation [62]. The intermediate products for AC, i.e., aerosol products, can be substantially different with different aerosol models. The four AC algorithms use the two NIR bands to select the two most suitable models, but the aerosol products generated from SeaDASv7.5 (NASA-STD and MUMM algorithms) consider 80 aerosol models with eight relative humidity values (30%, 50%, 70%, 75%, 80%, 85%, 90%, and 95%) and 10 aerosol particle size distributions for each value of relative humidity [12]. Conversely, the KOSC-STD method uses the simplified nine aerosol models from GW94: Oceanic with relative humidity of 99% (O99), Maritime (M50, M70, M90, and M99), Coastal (C50 and C70), and Tropospheric (T50 and T90) [69]. The hypothesis used in SeaDASv7.5, which is based on using single scattering to determine the aerosol models (weighting factor) and to apply the weighting factor directly to all bands, neglects the nonlinear relationship not only between suspended sediments and multiple-scattering aerosol reflectance but also between the interband relations of multiple scattering. To avoid the unsuitability of this condition, the spectral relationships in the aerosol multiple-scattering reflectance contributions are updated in GDPSv2.0 [70], which are improved by 1.1% in comparison with suspended sediments. Moreover, aerosol models used in SeaDASv7.5 and GPDSv2.0 are non-or weakly-absorbing aerosols that might not be completely representative of the aerosols over the YS region [observed values of AOT(870) are >0.3 with AE values of around 1.5].
Slot effects are evident features of GOCI data. An L1B GOCI image consists of 16 slot images. The slot border stray-light effect is a particular artifact of the GOCI optical system that creates clear inconsistency between adjacent slots. It primarily has an impact on the AE [29]. To minimize its influence, an image-based method called the minimum noise fraction transform has been proposed for the partial removal of the slot border [71]. In the future, this method will be applied in GDPS Software. However, NASA's Ocean Biology Processing Group has still not considered using the method to mitigate the problem. Based on our study, it is considered that slot effects have obvious influence on the images ( Figure 10). Slot effects are evident features of GOCI data. An L1B GOCI image consists of 16 slot images. The slot border stray-light effect is a particular artifact of the GOCI optical system that creates clear inconsistency between adjacent slots. It primarily has an impact on the AE [29]. To minimize its influence, an image-based method called the minimum noise fraction transform has been proposed for the partial removal of the slot border [71]. In the future, this method will be applied in GDPS Software. However, NASA's Ocean Biology Processing Group has still not considered using the method to mitigate the problem. Based on our study, it is considered that slot effects have obvious influence on the images ( Figure 10). Vicarious calibration (VC) of the GOCI system is conducted to ensure that a relative radiometric calibration to GOCI minimizes uncertainties in the derived remote sensing reflectance products. The system VC gains for GOCI are slightly different between NASA and KIOST/KOSC. For NASA, it is applied as a one-time gain without any in-depth or long-term evaluation [72], whereas the default VC gains processed by SeaDASv7.5 are 0.9726, 0.9520, 0.9258, 0.8974, 0.9007, 0.8719, 0.9430, and 1.0 for each GOCI band [72]. In the present GDPSv2.0, the VC gains for 412-865 nm are 1.0118, 0.9954, 0.9715, 0.93431, 0.9596, 0.9669, 0.96125, and 1.0 [73,74]. The influence of VC on AC between KIOST/KOSC and NASA will be evaluated in future research. The MUMM and the Kd-based algorithms were not calibrated vicariously in this study (thus, it is beyond the scope of this article).

Conclusions
This study evaluated the performances of the NASA-STD, the KOSC-STD, the Kd-based, and the MUMM algorithms for GOCI over the YS at three observation times (02:16, 03:16, and 04:16 UTC) with temporal windows of ±1 and ±3 h using in situ aerosol optical measurements (AOT and AE), Rrs, and Chla data from shipborne observations. Our quantitative results indicate that the KOSC-STD method was the algorithm most appropriate for deriving Rrs and aerosol indices in Vicarious calibration (VC) of the GOCI system is conducted to ensure that a relative radiometric calibration to GOCI minimizes uncertainties in the derived remote sensing reflectance products. The system VC gains for GOCI are slightly different between NASA and KIOST/KOSC. For NASA, it is applied as a one-time gain without any in-depth or long-term evaluation [72], whereas the default VC gains processed by SeaDASv7.5 are 0.9726, 0.9520, 0.9258, 0.8974, 0.9007, 0.8719, 0.9430, and 1.0 for each GOCI band [72]. In the present GDPSv2.0, the VC gains for 412-865 nm are 1.0118, 0.9954, 0.9715, 0.93431, 0.9596, 0.9669, 0.96125, and 1.0 [73,74]. The influence of VC on AC between KIOST/KOSC and NASA will be evaluated in future research. The MUMM and the Kd-based algorithms were not calibrated vicariously in this study (thus, it is beyond the scope of this article).

Conclusions
This study evaluated the performances of the NASA-STD, the KOSC-STD, the Kd-based, and the MUMM algorithms for GOCI over the YS at three observation times (02:16, 03:16, and 04: 16 UTC) with temporal windows of ±1 and ±3 h using in situ aerosol optical measurements (AOT and AE), Rrs, and Chla data from shipborne observations. Our quantitative results indicate that the KOSC-STD method was the algorithm most appropriate for deriving Rrs and aerosol indices in summer over the YS region with ±1 and ±3 h windows. The Kd-based NIR correction algorithm was the second most accurate algorithm, and the MUMM method was the least accurate. To interpret the temporal changes of GOCI-derived Rrs and AOT products, the results suggest that Rrs and AOT products were closer to in situ values at 04:16 UTC.
In terms of Rrs, all algorithms showed better performance at 490 and 555 nm for all stations at all GOCI observing times and low accuracies in the blue (412 nm) and the NIR bands (745-865 nm). For the band ratio of Rrs commonly used to estimate Chla, the use of Rrs(490)/Rrs(555) is recommended for the YS region rather than Rrs(443)/Rrs(555). For Chla retrievals, the Kd-based and the KOSC-STD algorithms produced Chla concentration estimations that were most accurate; the NASA-STD and the MUMM methods did not perform well in this process in our study.
For AOT, the uncertainties of satellite-retrieved values were greatest at 443 and 865 nm. AOT(555) showed reasonably satisfactory agreement with in situ AOT. However, none of the algorithms properly estimated AE, which means further work is needed to expand consideration of aerosol type and to deal with the influence of moderately and strongly absorbing aerosols.
Further improvements and optimization for AC in coastal turbid waters could be undertaken as follows. (1) The choice of NIR correction models has a crucial effect on AC accuracy because such models are often based on several hypotheses or relationships between interband Lw. The classification of water type in coastal oceans [65,75,76] will help to refine the bio-optical coefficients in NIR correction models; however, some studies have classified China's coastal waters into several types. (2) The AC process requires a considerable amount of in situ aerosol data (e.g., aerosol particle size distribution and index of refraction) [77] to validate and develop regional AC algorithms. A large network of ground stations such as AERONET [55], which relies on CE-318 automated sun photometers, could provide wide coverage; however, there are few such stations distributed in Chinese coastal waters at present. Therefore, it is urgently required that an aerosol observation network be established over China's coast. (3) The latest demands regarding sensors have been for the addition of UV bands because such bands are sensitive to aerosol absorption [78]. Consequently, GOCI-II will have UV bands [79]. (4) We also need realistic moderately and strongly absorbing aerosol models for generating realistic aerosol look-up tables. (5) One of the most important future tasks is to tackle the limited number of Rrs match-ups to assess the performance of AC over coastal waters with several different optical constituents. Considerable efforts are being made to enable the instruments to obtain continuous Rrs data near Chinese coastal waters. Additionally, the present algorithms should be validated with long-term observations covering different seasons and different coastal areas.  Figure A1 represent the spectra of Rrs over five subregions. According to these figures, the peaks of the spectra for PSD, YT, and TZ are around 555 nm and around 490 nm for WH. For CWM, the spectra of Rrs decrease rapidly when the wavelength increases in most of the stations.

Appendix A
Remote Sens. 2019, 11, x FOR PEER REVIEW 23 of 27 Figure A1 represent the spectra of Rrs over five subregions. According to these figures, the peaks of the spectra for PSD, YT, and TZ are around 555 nm and around 490 nm for WH. For CWM, the spectra of Rrs decrease rapidly when the wavelength increases in most of the stations.