Atmospheric Correction of GOCI Using Quasi ‐ Synchronous VIIRS Data in Highly Turbid Coastal Waters

: The Geostationary Ocean Color Imager (GOCI) sensor, with high temporal and spatial resolution (eight images per day at an interval of 1 hour, 500 m), is the world’s first geostationary ocean color satellite sensor. GOCI provides good data for ocean color remote sensing in the Western Pacific, among the most turbid waters in the world. However, GOCI has no shortwave infrared (SWIR) bands making atmospheric correction (AC) challenging in highly turbid coastal regions. In this paper, we have developed a new AC algorithm for GOCI in turbid coastal waters by using quasi ‐ synchronous Visible Infrared Imaging Radiometer Suite (VIIRS) data. This new algorithm estimates and removes the aerosol scattering reflectance according to the contributing aerosol models and the aerosol optical thickness estimated by VIIRS’s near ‐ infrared (NIR) and SWIR bands. Comparisons with other AC algorithms showed that the new algorithm provides a simple, effective, AC approach for GOCI to obtain reasonable results in highly turbid coastal waters.


Introduction
The total radiance measured by an ocean color sensor is primarily composed of the water-leaving radiance, sea surface radiance, Rayleigh scattering radiance caused by air molecules, and aerosol scattering radiance (which includes aerosol single-scattering radiance and interactive scattering radiance between molecules and aerosols). In open waters, the atmospheric radiance can account for significant portions of the total satellite-measured radiance [1]. Therefore, to determine the properties of the upper ocean, such as colored dissolved organic matter, suspended particulate matter, and chlorophyll-a, accurate atmospheric correction (AC) is required. The Rayleigh scattering radiance can be theoretically computed accurately, owing to the stable distribution of the necessary atmospheric components [2,3]. Therefore, it is key to accurately estimate the aerosol scattering by aerosols and Rayleigh-aerosol interactions in order to determine the water-leaving radiance. For clear waters, the AC algorithm developed by Gordon and Wang [4] (herein named the GW94 algorithm) works quite well. It estimates the aerosol optical properties based on the black pixel assumption, according to which the water-leaving radiance at near-infrared (NIR) bands is assumed to be zero because the water can strongly absorb the light in these bands. However, in turbid coastal

Study Sites
The Bohai Sea, Yellow Sea, and East China Sea are parts of the Western Pacific marginal sea. Three highly turbid coastal regions are outlined in the boxes in this region in Figure 1a: (I) Bohai Sea, (II) the southwest coast of Korea, and (III) Changjiang Estuary.
The Bohai Sea is a shallow semi-enclosed sea with an average water depth of 18 m and a maximum depth of ~70 m. It deepens gradually from the coastal bays to the Central Bohai Sea and is characterized by a basin shape [20]. The dominant sediment source in this region is the sediment delivered by the Huanghe [21,22]. The high suspended sediment concentration is mainly distributed in the south of Bohai Sea, especially in the area around the Huanghe Delta [21].
The water depth of the area along the Southwestern Korean coast is less than 50 m. Numerous islands and vast tidal flats are located here and the coastlines in this area are complicated. The suspended sediment concentration of this area is relatively high (>20 g/m 3 ). The highest values (>100 g/m 3 ) occur in the winter owing to the stronger northwestern monsoon and shallow water depth, which can induce a resuspension of bottom sediments [23].
The Changjiang Estuary is located along the central-eastern coast of China, with the Hangzhou Bay to the south and the Subei Shallow to the north. The highly suspended sediment concentrations were observed here year-round [24]. The Changjiang River discharges about 390 × 10 6 tons of sediment into the East China Sea annually [25]. The sediment transportation inside the Hangzhou Bay is considerably affected by the secondary Changjiang plume [26]. Another major sediment source for the East China Sea is the resuspension of sediment at the Subei Shallow [27]. Wind-induced vertical mixing and bottom stress tend to resuspend a large amount of sediment, leading to high suspended sediment concentrations, especially in the winter [28,29].
One cloud-free example is chosen from each of these regions. The RGB pictures of the three examples are composed of the total radiance at 490 nm (B), 555 nm (G), and 680 nm (R) bands of GOCI and are shown in Figure 1b

Method
Reflectance ρ is defined as: where λ denotes the wavelength, F0 is the extraterrestrial solar irradiance [30], θs is the solar zenith angle, and L is the radiance. It is more convenient to work with ρ because it is dimensionless. The total reflectance ρt(λ) measured by a sensor can be written as: where ρr denotes Rayleigh scattering reflectance, ρa denotes aerosol multiple scattering reflectance, ρw denotes water-leaving reflectance, and tv denotes the Rayleigh-aerosol diffuse transmittance from the sea surface to the satellite. It should be noted that the reflectance contributions from the sun glint and whitecaps are ignored [31,32]. Determining the remote sensing reflectance (Rrs, unit: sr -1 ) and normalized water-leaving reflectance (ρwn) is the ultimate goal of AC because they are fundamental parameters that are widely used for deriving the water constituent concentrations and water quality parameters. These can be calculated by: where ts is the diffuse transmittance from the sun to the sea surface.
In this section, we briefly review the previous AC algorithms that were adopted in the main GOCI processing software: SeaDAS (SeaWiFS Data Analysis System) and GDPS (GOCI data process system). Then, the development of the new algorithm is described.

Previous AC Algorithms
The key problem of AC over turbid waters is the removal of the water-leaving reflectance at NIR bands. Numerous red/NIR modeling approaches have been investigated to deal with non-zero ρw(red/NIR) values within the AC process. Three kinds of AC approaches are briefly reviewed herein-the B2010 AC algorithm proposed by Bailey et al. [33] is currently adopted in SeaDAS, while the A2012 algorithm proposed by Ahn et al. [34] and the L2013 algorithm proposed by Lee et al. [35] are implemented in GDPS.
The B2010 algorithm uses an iterative solution to separate ρw(NIR) and ρa(NIR). The ρw(NIR) is first assumed to be zero to complete the GW94 AC process. This gives an initial estimate of Rrs(λ). Next, the concentration of chlorophyll-a is preliminarily estimated by the initial value of Rrs(λ). Then, the absorption coefficient in the red band is obtained via a chlorophyll-a-based empirical relationship. The absorption coefficient and Rrs in the red band are used to solve the backscattering coefficient, which is used to extrapolate the ρw at the NIR band. Finally, the new values of ρw(NIR) are used to correct the non-zero water contribution, and the GW94 AC step is repeated for a new iteration.
The A2012 algorithm is the GOCI standard AC algorithm. It is theoretically based on the GW94 algorithm with partial modifications. The initial ρw(NIR) values are also assumed to be zero to execute the GW94 algorithm. The newly estimated ρwn in the red band is used to correct the ρwn(NIR) by an empirical polynomial relationship model [36]: where jn and kn are known polynomial coefficients. The L2013 algorithm is from the Management Unit of the North Sea Mathematical Models (MUMM) [37] algorithm, with some modified steps for extremely turbid water. In the original MUMM, the parameter α, representing the ratio of ρw at the NIR wavelengths, is assumed to have spatial homogeneity; however, the value of α changes with the concentration of suspended particles in extremely turbid waters, which would cause an underestimation of the water-leaving radiances [38]. The L2013 algorithm calculates the α value adaptively using an NIR water-leaving reflectance model: where the cn values are known polynomial coefficients. The algorithms mentioned above can improve the data quality in turbid waters, but the empirical relationships are highly reliant on in situ datasets; even the empirical relationship used in L2013 is derived from the nearest non-turbid AC algorithm [39]. These algorithms are, therefore, only regionally applicable [40]. For example, the B2010 method only works properly in low to moderately turbid waters [41,42], as the chlorophyll-a-based relationship used in the bio-optical model might not be appropriate to extrapolate the ρw(NIR) in waters whose optical properties are dominated by non-algal particles [41,42]. The A2012 method is suitable for sediment-dominated waters but also fails for extremely turbid waters because the ρwn(660 nm) can become optically saturated [43].

The Development of the New AC Algorithm
In order to distinguish the two sensors, the superscripts V and G represent VIIRS and GOCI, respectively. The procedure for the new algorithm can be divided into four parts:  Extracting the ρa(NIR V ) of VIIRS;  Estimating the aerosol properties by ρa(NIR V );  Calculating the ρa(λ G ) at GOCI observing geometries according to the aerosol properties; and  Removing the ρa(λ G ) and completing the GOCI AC process.
The removal of the ρw(NIR V ) of VIIRS follows the Shortwave Infrared Exponential (SWIRE) algorithm [12], in which the Rayleigh-corrected reflectance values ρrc (ρrc = ρt − ρr) of clear waters were assumed to be an exponential function of wavelength at the NIR/SWIR bands: where ρef is the extrapolated Rayleigh-corrected reflectance, and a and b are the fitting coefficients.
The ρrc values at the three VIIRS SWIR bands (1238, 1601, and 2257 nm) are used to calculate the coefficients a and b because the ρrc at the NIR bands is strongly influenced by suspended sediments scattering in turbid waters. The newly estimated ρef(NIR) can be considered to be the ρrc of the optical equivalent clear water, influenced only by aerosol scattering, and can be used for the GW94 AC of the visible bands. The SWIR-based methods improve the ocean color products in the turbid coastal waters, but they require higher signal-to-noise ratios (SNR) for SWIR bands because of the longer extrapolated distances and lower signals [44]. Wang and Shi [14] used two SWIR bands to derive the aerosol scattering radiance, but their method did not achieve better results than the NIR-based methods in non-turbid waters because significant noise occurred in the derived products [45,46]. Generally, the more SWIR bands used, the more accurately can the aerosol scattering reflectance be estimated, and the requirements of SNR can probably be reduced. As a result, the SWIRE method represented an improvement over the SWIR-based methods [47].
The newly estimated ρef(NIR V ) can be considered as the ρa(NIR V ) and can be applied to evaluate the aerosol properties of the observing area, including the contributing aerosol models and aerosol optical thickness τa. The contributing aerosol models are selected by the look-up-tables (LUTs) for 80 aerosol tables constructed by Ahmad et al. [48]. The parameters used to derive the aerosol single-scattering albedo (ωa), aerosol scattering phase function (Pa), and extinction coefficient (βa) are stored in the LUTs for each aerosol model and various geometries. The LUTs also contain the coefficients relating single to multiple scattering and the coefficients (A and B) used to calculate the Rayleigh-aerosol, diffuse transmittance in the following form [49]: The aerosol optical thickness is the normalized extinction coefficient due to absorption and scattering from the direct beam and is a key parameter used to define the optical state of the atmosphere [50,51].
The estimation of the aerosol properties follows the GW94 method [4], which is one of the AC algorithms widely used for clear waters. In this algorithm, ρa(NIR V ) is first converted into ρas(NIR V ) [52] by solving the quadratic equation: where the pn values are the quadratic coefficients stored in the LUTs. The AC parameter is defined as: where NIRS is the short-wavelength NIR band, and NIRL is the long-wavelength NIR band. For the VIIRS sensor, the NIRS and NIRL are 745 nm and 862 nm, respectively. Theoretically, each candidate aerosol model has its own ε value (ε M ) for a specific geometry: Two aerosol models with different weighting factors whose ε(NIRS, NIRL) values are the closest to the ε M (NIRS, NIRL) averaged over all candidate models, are chosen to represent the aerosol type over the pixel. The selection of the closest aerosol models can dominate the estimation of the water-leaving reflectance, especially over turbid coastal waters [53]. The τa(745 nm) values of the two contributing models can be retrieved by [4]: The VIIRS sensor is located in a sun-synchronous orbit and provides one image per day, whereas the GOCI sensor is located in a geostationary orbit and provides eight images per day. With the hourly observations from the GOCI, the VIIRS image has a corresponding GOCI image with a time difference of less than half an hour, which means the atmospheric conditions are nearly invariant. Therefore, the aerosol scattering reflectance of the GOCI observing geometries can be derived according to the known contributing aerosol models, weighting factor, and the τa(745 nm). First, calculate the τa(λ G ) of each contributing aerosol model by: Second, calculate the ρas(λ G ) by Equation (13). Next, convert ρas(λ G ) to ρam(λ G ) via Equation (10). Finally, the effects of the two contributing aerosol models are superimposed by weighting factors to estimate the total aerosol scattering reflectance of GOCI. In order to maintain the higher spatial resolution of GOCI, the aerosol parameters estimated from VIIRS are resampled to the GOCI data grid by the nearest neighbor method.

Results and Discussion
In the absence of in situ data, the results of the new algorithm (represented by QSV) are compared with those processed by the B2010 [33], A2012 [34], and L2013 [35] algorithms. VIIRS and GOCI Level-1 data were obtained online (https://oceancolor.gsfc.nasa.gov). For the QSV algorithm, the Rayleigh-corrected reflectance values were obtained by using the SeaDAS 7.4 l2gen processor, and these values were used as the inputs of the aerosol scattering correction procedure, which was performed by using Interactive Data Language (IDL). With default settings, the B2010 algorithm was realized by the SeaDAS 7.4 software, while the A2012 and L2013 algorithms were realized by the GDPS v1.4.1 software. The three examples listed in Table 1 are processed. Figures 2-4 represent the Rrs values retrieved by four different algorithms centered at the Bohai Sea, the southwest coast of Korea, and the Changjiang Estuary, respectively. Panels (a-d), (e-h), and (i-l) are results at wavelengths of 490, 555, and 680 nm, respectively. The regions in which the algorithms fail to make an AC because of highly turbid waters or due to the presence of clouds and land, are shown with RGB pictures. It can be observed that the Rrs distributions are similar for all four algorithms, but the B2010, A2012, and L2013 algorithms result in varying degrees of failures, while the QSV algorithm can successfully execute the AC in the highly turbid coastal waters. The unmasked abnormal low estimations from L2013 also can be observed at 35 °N, 126 °E in Figure 3b, f, j, highlighted by the red circles.  To further demonstrate the differences, we also extracted the Rrs(λ) values for pixels along the red arrow in Figure 4a and plotted their profiles in Figure 5. A distance of 0 along the x-axis indicates the location closest to the coast. In general, the Rrs data are low in the open ocean and high near the coast. The mean absolute difference in the Rrs values between the A2012 and QSV algorithms is the smallest, with values of 0.00334, 0.00138, 0.00100, 0.00091, 0.00110, and 0.00182 sr -1 at wavelengths of 412, 443, 490, 555, 660, and 680 nm, respectively. At around 0 to 110 pixels along the transect, the B2010 algorithm is invalid. The results of the L2013 are relatively higher than those of the A2012 and QSV algorithms, at distances larger than 50 pixels along the transect. Abnormal fluctuations of the L2013 curves can be observed at around 25-50 pixels. The A2012 and L2013 methods both fail at around 0-25 pixels. The QSV is the only algorithm that can retrieve continuous and reasonable values of Rrs along the total transect.  Three sets of Rrs spectra derived from four different AC algorithms in highly turbid waters at station (a), moderately turbid water at station (b), and clear water at station (c) are compared and shown in Figure 6. At station (c), the curves of the B2010, A2012, and QSV algorithms tend to be coincident, while the Rrs values retrieved by L2013 are relatively high. At station (b), the B2010 algorithm shows slightly lower estimations. At station (a), the B2010 algorithm has no valid results, and the difference between L2013 and A2012/QSV is more pronounced. Since these four algorithms are developed for turbid waters, they can all improve the data quality in turbid waters. The B2010, A2012, and L2013 algorithms are based on empirical models and are highly regionally dependent. Significant errors or failures can be induced when these methods are applied to the ocean regions they are not suitable for. Even though the QSV algorithm can be applied to only one scene of the GOCI data, it can obtain reasonable results in extremely turbid coastal waters. The use of three SWIR bands can improve the data quality of the AC retrievals in clear waters without obvious noise. In general, the QSV algorithm provides a simple, effective, new AC approach for GOCI to obtain reasonable results in highly turbid waters and provides an alternative method to validating other AC methods.

Conclusion
An alternative AC algorithm using quasi-synchronous VIIRS data for GOCI in highly turbid waters is presented in this paper. GOCI covers highly turbid waters in the Western Pacific region; AC has been a challenge in these turbid ocean regions because of the lack of SWIR bands. The new algorithm shows its superiority in highly turbid coastal waters, and the results could be at least be used for the validation of other AC methods. However, there are still some limitations that should be noted. The major disadvantage of this algorithm is that only one VIIRS image per day is not sufficient for GOCI, and more in situ data are needed for the validation. In the future, GOCI-2 will have full-disk coverage with higher resolution and five more bands than GOCI does, but it will still not have any SWIR bands [54]. Including SWIR bands on a sensor is expensive, therefore, further research is necessary on how to best take advantage of SWIR bands from other ocean color sensors.