A Secchi Depth Algorithm Considering the Residual Error in Satellite Remote Sensing Reﬂectance Data

: A scheme to semi-analytically derive waters’ Secchi depth ( Z sd ) from remote sensing reﬂectance ( R rs ) considering the e ﬀ ects of the residual errors in satellite R rs was developed for the China Eastern Coastal Zone (CECZ). This approach was evaluated and compared against three existing algorithms using ﬁeld measurements. As it was challenging to provide the accurately inherent optical properties data for running the three existing algorithms in the extremely turbid waters, the new developed algorithm worked more e ﬀ ective than the latter. Moreover, with both synthetic and match-up data, the results indicated that the proposed algorithm was able to minimize some residual errors in R rs , and thus could generate inter-mission consistent Z sd results from two ocean color missions. Finally, after application of new model to satellite images, we presented the spatial and temporal variations of Secchi depth and trophic state in the CECZ during 2002–2014. The study led to several ﬁndings: Firstly, the Z sd -based trophic state index (TSI) in the East China Sea ﬁrst increased since 2002, and then gradually dropped during 2008–2014. Secondly, more and more waters within 30–35 m and 20–25 m isobaths were deteriorating from oligotrophic to mesotrophic type and from mesotrophic to eutrophic water, respectively, during 2002–2014. Lastly, the TSI increased on average 0.091 and 0.286 m per year respectively in Bohai Sea and Yellow Sea since 2002, and it might only take 14 and 67 years for Bohai Sea and Yellow Sea to deteriorate from mesotrophic to eutrophic water, following their current yearly deterioration rate and trophic trend. These results highlighted the importance to make some strict regulations for protecting the aquatic environment in the CECZ.


Introduction
The coastal zones are populated zones for human society on the earth, and the productivity, trophic state, and health of these have great impact on the lifestyle of peoples living around them. To support the healthy development of the human society near the coastal zones, it is necessary to know the history of their trophic state. However, the optical properties in coastal waters are always very complex [1,2], so it is not easy to obtain ocean color data from satellite with good quality. The Secchi depth (Z sd ) is a measurement of water transparency, which is a simple and universal method in oceanographic environmental survey. As the water transparency is closely related to optically activity components, Z sd is usually used for estimating the order of concentration of optical substances in waters [3]. Moreover, the Z sd can be cheap and easy to obtain, and there is a large archive of historical Secchi

In Situ Dataset
To evaluate and compare the applicability of BGSD, IOPSD, MKSD and CSSD models in quantifying the Zsd from turbid coastal waters, three independent in situ datasets were collected from the CECZ seas (see Figure 1 and Table 1). The field measurements included the simultaneous measurements of the Rrs and Zsd. The Rrs were measured using a field Spectroradiometer (Spectral Devices, Boulder, CO, ASD) from 10:00 to 16:00 local time. To avoid the interference of ship reflect and shadow with measurements, the instrument was positioned at an angle of 90-135° away from the sun, while the view angle was fixed 30-45° with the aplomb direction [27]. Three measurements were made at each station in order to minimize the effects of random uncertainties on measurements. The first dataset included 321 samples used for model calibration (called as initialization dataset), and it was taken from 101 independent cruises during 2006-2012. The second dataset included 285 samples used for model evaluation (called as the test dataset), and it was taken from 55 independent cruise ranging from 2003 to 2005. Moreover, we collected 41 match-up (±3 hours) data during [2003][2004][2005] for evaluating the accuracy of satellite-derived Zsd data (called as synchronized dataset). These data were measured closely according to the rigorous and community-defined protocols for experiment deployment and data collection [27].
Importantly, many studies have fully reported the BGSD, IOPSD, and MKSD models [16,19]. The empirical coefficients of BGSD model were adjusted according to the field measurements in Table  1. The IOPSD and MKSD models were the global algorithms initialized using the global datasets, and have been proven to be valid for most global oceans. Thus the empirical parameters of these two models were not reinitialized according to the initialization dataset in this study.

In Situ Dataset
To evaluate and compare the applicability of BGSD, IOPSD, MKSD and CSSD models in quantifying the Z sd from turbid coastal waters, three independent in situ datasets were collected from the CECZ seas (see Figure 1 and Table 1). The field measurements included the simultaneous measurements of the R rs and Z sd . The R rs were measured using a field Spectroradiometer (Spectral Devices, Boulder, CO, ASD) from 10:00 to 16:00 local time. To avoid the interference of ship reflect and shadow with measurements, the instrument was positioned at an angle of 90-135 • away from the sun, while the view angle was fixed 30-45 • with the aplomb direction [26]. Three measurements were made at each station in order to minimize the effects of random uncertainties on measurements. The first dataset included 321 samples used for model calibration (called as initialization dataset), and it was taken from 101 independent cruises during 2006-2012. The second dataset included 285 samples used for model evaluation (called as the test dataset), and it was taken from 55 independent cruise ranging from 2003 to 2005. Moreover, we collected 41 match-up (±3 h) data during 2003-2005 for evaluating the accuracy of satellite-derived Z sd data (called as synchronized dataset). These data were measured closely according to the rigorous and community-defined protocols for experiment deployment and data collection [26]. Table 1. Descriptive statistics of field-measured optical properties in the initialization dataset, test dataset, and synchronized dataset. Stdev is the standard deviation. Note that the K d (488) is derived from R rs data using neural network model developed by Chen et al. [20]. Importantly, many studies have fully reported the BGSD, IOPSD, and MKSD models [15,18]. The empirical coefficients of BGSD model were adjusted according to the field measurements in Table 1. The IOPSD and MKSD models were the global algorithms initialized using the global datasets, and have been proven to be valid for most global oceans. Thus the empirical parameters of these two models were not reinitialized according to the initialization dataset in this study.

Synthetic Dataset
Data from field measurements always include some uncertainties, thus a synthetic dataset of 1000 samples provided by International Ocean Colour Coordinating Group (IOCCG) [27], for wavelengths of the MODIS bands generated by Hydrolight version 5.2, were used to diagnose the sensitivity of CSSD model to the residual error in R rs data. After considering all potential error sources in the atmospheric correction, Hu et al. [28] estimated error budget of <±0.0003 sr −1 for R rs (667) ( R rs (667)), which was roughly 50% smaller than the results in Gordon [29] (calculated from <±0.0013 sr −1 for R rs (555) using spectral relationship of residual error in Hu et al. [28]). Thus, to see how CSSD performed for residual error in R rs data, two tests were designed. In test 1 we added <±0.0003 sr −1 random noise to the error-free R rs (667) data, while in the test 2 we increased the added random noise to <±0.0006 sr −1 . The spectral relationship for residual error in Chen et al. [21] was used to extrapolate the random noise from 667 nm to shorter wavelengths and added into the error-free R rs data in these two tests. This added residual error could be viewed as errors to R rs data introduced by imperfect sensor calibration and atmospheric correction procedures [28].

Satellite Image and Match-Up Data
To learn the spatial and temporal variation patterns of Z sd and TSI in the CECZ seas, the MODIS daily images during July 2002-December 2014 with 1 km-resolution were collected from the NASA data center. The near-infrared and shortwave infrared wavelengths combined the atmospheric correction algorithm (NIR-SWIR model) developed by Wang and Shi [30] which was used to derive R rs products from the MODIS radiance products at the top-of-atmosphere. The NIR-SWIR atmospheric correction model has been proven to be valid for most coastal waters including the CECZ seas [30,31]. The Hilbert-Huang transform (HHT) approach developed by Huang [32] was suggested for reconstructing the time series of Z sd and TSI in the CECZ seas, as the HHT method can improve the operational efficiency and robustness compared with empirical mode decomposition method.
The SeaWiFS Bio-optical Active and Storage System (SeaBASS) maintains a local repository of field-measured and satellite-observed optical data to support and sustain regular study such as matchup analyses [33]. To illustrate the performances of the band ratio-based and band difference-based CSSD models in minimizing residual error in satellite R rs data, 908 synchronized pairs (±3 h) including MODIS versus SeaWiFS R rs data have been collected from SeaBASS. This was true for evaluating the sensitivity of CSSD model to the residual error in satellite-derived R rs data, as the inter-mission difference between MODIS and SeaWiFS-derived Z sd data would increase with the increase of sensitivity of model to residual error. The unwanted "low-quality" data were discarded following the standard exclusion criteria proposed by Bailey and Werdell [34].

The Overall Scheme for Derivation of Z sd and TSI
The inherent optical properties (IOPs) data processing system (IDAS) is insensitive to residual error in satellite R rs [21], so it can produce a superior performance to the simple band ratios in IOPs retrievals [35]. In the following analyses, the waters are empirically classified into three different optical types by a turbid index (T d ): low and moderate water type (T d < 0.01 where Z sd >~1.0 m), intermediate type (0.01 ≤ T d < 0.014), and extremely turbid water type (T d ≥ 0.014 where Z sd <~0.5 m), where T d is the band difference between R rs at 667 and 488 nm. Note that many intervals for T d have been set for model development, but the empirical results showed that the interval division using T d values with 0.01 and 0.014 was the optimal scheme (with the largest determination coefficient, R 2 ) for Z sd estimation in the CECZ seas. Specifically: (1) The IDAS model is utilized to calculate the a(λ) and b b (λ) from R rs data. This model is a step-wise algorithm, with absorption and backscattering coefficients first derived at 555 nm, and then the estimations are extended to λ after applying a spectral model of the particle backscattering coefficient [21]. (2) Lee et al. [18] showed that the Z sd is a function of K d (λ), while Lee et al. [36] proposed that the K d (λ) is a function of total absorption coefficient (a(λ)) and backscattering coefficient (b b (λ)). So Z sd should be a function of a(λ) and b b (λ). Thus, the relationship between Z sd and IDAS-derived a(λ) and b b (λ) are analyzed here. Based on the results of analyses, with a(488) and b b (488) as inputs, a new algorithm is developed for estimation of Z sd for both low and moderate turbid waters (denoted as Z sd,tc ). (3) The optical properties are commonly very complicated in the extremely turbid waters [1], so it is not easy to derive a(λ) and b b (λ) from R rs data with high quality in these waters. As a result, the semi-analytical algorithm in step 2 works well for low and moderate turbid waters but might be invalid for extremely turbid waters. Thus, an empirical approach is developed for predicting Z sd for extremely turbid waters (denoted as Z sd,et ). This empirical algorithm denotes Z sd as a function of [R rs (748) − R rs (869)] −1 , as the near-infrared wavelengths is more sensitive to changes of optically activity components than blue and green wavelengths in the extremely turbid waters. Following Hu et al. [28] and Chen et al. [21], band difference approach can absorb some residual error in satellite R rs for Z sd retrieval, as the satellite residual error has good linear relationships between each other at the visible bands. (4) These two algorithms need to be combined to generate smooth Z sd data for low, moderate, and extremely turbid water, as well as waters of intermediate type. Thus, the linear weighting function (W) developed by Wang et al. [36] is applied for bridging two types of models for Z sd retrieval in the intermediate water type (denoted as Z sd,mt ). Calibrated using initialization dataset, the optimal CSSD model can be found as follows (see Figure 2d): Here b bw is the backscattering coefficient of pure water; The parameter b bw (488)/b b (488) was used to minimize the effect of changing scattering agents [36,37].  Table 1.
Importantly, there were many scatterplots (~50 anomalous samples) far away from 1:1 line in IOPSD and MKSD models (see Figures 2b and 2c), which in turn led to the performance of IOPSD and MKSD models in initialization dataset were quite different from those in test dataset (see Figures  3b and 3c). The anomalous samples were mainly collected from 17 independent cruises over the South Yellow Sea on April 2007. If those anomalous samples were excluded, the scatterplots of both IOPSD and MKSD models became tighter around 1:1 line as well as test dataset. The possible reason of this result might be attributed to complex optical properties in the sample regions, but we did not have sufficient field measurements to support this statement. However, the impacts of those anomalous samples on CSSD model are very limited. Those results further confirmed the effective of CSSD model in deriving Zsd from optical complex coastal waters.

Model Evaluation and Comparisons
The model evaluation was assessed by comparison of the model-derived Zsd with the field-measured Zsd in the test dataset in Table 1. Figure 3 showed the comparison of the predicted versus measured estimations of Zsd by BGSD, IOPSD, MKSD and CSSD models. Similar to the results obtained from the initialization dataset in Figure 2a, the BGSD model's Zsd estimations were accurate for the higher end (Zsd>10 m), but produced significant underestimation below 6 m. The bias between field measurements and BGSD model predictions increased with the decreasing Zsd, but there was no statistically significant relationship between them. The IOPSD model's Zsd estimations were found to be accurate for waters with Zsd>1m, even though some underestimation could still be observed in this Zsd range (see Figure 3b). However, for the lower end (Zsd<1m), the Zsd was systematically overestimated using IOPSD model. In comparison, the scatter points of the field-measured and model-predicted Zsd for the MKSD model were much tighter around 1:1 line than both BGSD and IOPSD models (see Figure 3c), especially for Zsd lower than 1 m, indicating that the MKSD model was more effective than both BGSD and IOPSD models in Zsd estimation in waters with Zsd<1 m, at least for this dataset. Specifically, the MSPD value of MKSD model was 55.98%, which was 31.77% and 2.35% respective lower than the BGSD and IOPSD models. However, for lower end (Zsd <0.5 m), the MSPD and RMSElog values of MKSD model were 164% and 0.263 m, respectively. When we ignored  Table 1. A multi-parameter trophic index is limited in its usefulness because of the number of parameters including chlorophyll-a, phosphorus, nitrogen, and permanganate concentration must be known, but most of them are hard to obtain using satellite technology. The Z sd -based TSI proposed by Carlson [23] is simple and easy to use, so the Z sd -based TSI is suggested for identifying the trophic state of the CECZ seas. This trophic index is

Statistical Criterions for Product Evaluation
The root-mean-square of percent difference (MSPD) and root-mean-square of log transformed data (RMSE log ) are used to assess the accuracy of the Z sd value retrieval model. These statistical criterions are: Where x md,i is the model-derived Z sd of the i th element; x od,i is the field-measured Z sd of the i th element, and n is the total number of elements. Although many statistical metrics are available for assessing the performance of satellite data products [38], MSPD and RMSE log can be able to provide an appropriate metric for validation exercises due to their common usage.

BGSD, IOPSD, MKSD, and CSSD Models Calibration
Figure 2a-c showed the results of the scatterplots of field-measured against model-derived Z sd for BGSD, IOPSD, and MKSD models based on the initialization dataset, showing that the R 2 varied among models from 0.71 to 0.83. Overall, all of these three existing models could derive accurate Z sd data at the higher end (Z sd > 6 m), but significantly overestimated at the lower end (Z sd < 3 m), indicating that BGSD, IOPSD and MKSD models in waters with high Z sd worked better than waters with low Z sd , at least for this dataset. Judging by the slope, bias, and R 2 values, the BGSD model has a superior performance to the MKSD model, but slightly poorer than the IOPSD model in the CECZ seas.
Using the IDAS model-derived a(488), b b (488) and R rs data as inputs, the model shown in Figure 2d and Equation (1) was suggested as the optimal CSSD model for estimations of Z sd from the CECZ seas, indicating that the CSSD model was an effective algorithm for Z sd estimations with the MODIS spectral bands from optically complex coastal waters, whose R 2 was 0.87. In other words, with CSSD model could account for 87% variations of the Z sd in the initialization dataset. By comparing the results shown in Figure 2a-c, it was found that when judging by R 2 , the CSSD model worked better than the BGSD, IOPSD, and MKSD models in predicting Z sd from the CECZ seas, especially at the lower end (Z sd < 3 m), where the scatter plots of the in situ Z sd versus model Z sd for CSSD model were much tighter around 1:1 line than either of BGSD, IOPSD, and MKSD models, at least for this dataset.
Importantly, there were many scatterplots (~50 anomalous samples) far away from 1:1 line in IOPSD and MKSD models (see Figure 2b,c), which in turn led to the performance of IOPSD and MKSD models in initialization dataset were quite different from those in test dataset (see Figure 3b,c). The anomalous samples were mainly collected from 17 independent cruises over the South Yellow Sea on April 2007. If those anomalous samples were excluded, the scatterplots of both IOPSD and MKSD models became tighter around 1:1 line as well as test dataset. The possible reason of this result might be attributed to complex optical properties in the sample regions, but we did not have sufficient field measurements to support this statement. However, the impacts of those anomalous samples on CSSD model are very limited. Those results further confirmed the effective of CSSD model in deriving Z sd from optical complex coastal waters.   Table 1.
Based on 285 data points in Table 1, Figure 3d showed that the CSSD model-derived Zsd agreed well with the in situ data. Specifically, the slope and R 2 of linear relationship between model Zsd and in situ Zsd were 0.89 and 0.88, respectively, indicating that the CSSD model could account for 88% variations of the Zsd in the test dataset. By comparison, the CSSD model works more effective than the BGSD, IOPSD and MKSD models in deriving Zsd from the optically complex coastal waters like  Table 1.

Model Evaluation and Comparisons
The model evaluation was assessed by comparison of the model-derived Z sd with the field-measured Z sd in the test dataset in Table 1. Figure 3 showed the comparison of the predicted versus measured estimations of Z sd by BGSD, IOPSD, MKSD and CSSD models. Similar to the results obtained from the initialization dataset in Figure 2a, the BGSD model's Z sd estimations were accurate for the higher end (Z sd > 10 m), but produced significant underestimation below 6 m. The bias between field measurements and BGSD model predictions increased with the decreasing Z sd , but there was no statistically significant relationship between them. The IOPSD model's Z sd estimations were found to be accurate for waters with Z sd > 1 m, even though some underestimation could still be observed in this Z sd range (see Figure 3b). However, for the lower end (Z sd < 1 m), the Z sd was systematically overestimated using IOPSD model. In comparison, the scatter points of the field-measured and model-predicted Z sd for the MKSD model were much tighter around 1:1 line than both BGSD and IOPSD models (see Figure 3c), especially for Z sd lower than 1 m, indicating that the MKSD model was more effective than both BGSD and IOPSD models in Z sd estimation in waters with Z sd < 1 m, at least for this dataset. Specifically, the MSPD value of MKSD model was 55.98%, which was 31.77% and 2.35% respective lower than the BGSD and IOPSD models. However, for lower end (Z sd < 0.5 m), the MSPD and RMSE log values of MKSD model were 164% and 0.263 m, respectively. When we ignored points with Z sd < 0.5 m (including 20 samples), the MSPD value of MKSD model decreased significantly, from 56% to 39%.
Based on 285 data points in Table 1, Figure 3d showed that the CSSD model-derived Z sd agreed well with the in situ data. Specifically, the slope and R 2 of linear relationship between model Z sd and in situ Z sd were 0.89 and 0.88, respectively, indicating that the CSSD model could account for 88% variations of the Z sd in the test dataset. By comparison, the CSSD model works more effective than the BGSD, IOPSD and MKSD models in deriving Z sd from the optically complex coastal waters like the CECZ seas. Utilizing the CSSD model produced 35.73% MSPD and 0.135 m RMSE log values, which decreased by >20% MSPD and >0.035 m RMSE log values from the BGSD, IOPSD and MKSD models. Specifically, for waters with Z sd < 1.0 m, the MSPD value of CSSD model was 39.23% which was almost comparable with the MSPD value of MKSD model for Z sd exceeding 0.5 m. These results indicated that the advantage of CSSD model was primarily attributed to the improvement of Z sd retrieval accuracy in waters with Z sd < 3.0 m, where the scatter plots of field-measured and model-derived Z sd for CSSD model were tighter around 1:1 line, but the BGSD, IOPSD, and MKSD models exhibited significant overestimation or underestimation of Z sd . Those results further confirmed that the CSSD model outperformed the other algorithms for these turbid coastal waters, and could be used to quantify the Z sd even in the extremely turbid water, at least for this dataset.

Match-Up Analysis
The model-based satellite-derived Z sd could be retrieved from the MODIS data after the atmospheric correction using NIR-SWIR model. The accuracy of the satellite-predicted Z sd was evaluated by a comparison between the field-measured and synchronized satellite-quantified Z sd (±3 h) from 11 MODIS imageries in the CECZ seas. The approach proposed by Bailey and Werdell [34] was used to generate the Z sd data for matchup analysis. Based on 41 samples (see Table 1), Figure 4 showed the performance of BGSD, IOPSD, MKSD, and CSSD models in the CECZ seas. It was found that the satellite-predicted Z sd agreed well with the in situ data. Specifically, R 2 of the linear relationship between the field-measured and satellite-predicted Z sd did not lower than 0.84. By comparison, the CSSD model exhibited more effective than others, whose MSPD and RMSE log values were 33.43% and 0.115. Using the CSSD model for retrieving Z sd in the CECZ seas decreased by 28.29% MSPD value from BGSD model, 47.52% MSPD value from IOPSD model, and 7.15% MSPD value from MKSD model, respectively. It appeared that only MKSD and CSSD models can be used to retrieve Z sd from turbid coastal waters, but in these waters the CSSD model outperformed the MKSD model. These results indicated that when an atmospheric correction procedure for ocean color data was available, the Remote Sens. 2019, 11,1948 9 of 18 synoptic coverage ocean color data could be used for retrieving Z sd from the optical complex coastal waters using CSSD model. Remote Sens. 2019, 11, x FOR PEER REVIEW 9 of 18 were 33.43% and 0.115. Using the CSSD model for retrieving Zsd in the CECZ seas decreased by 28.29% MSPD value from BGSD model, 47.52% MSPD value from IOPSD model, and 7.15% MSPD value from MKSD model, respectively. It appeared that only MKSD and CSSD models can be used to retrieve Zsd from turbid coastal waters, but in these waters the CSSD model outperformed the MKSD model. These results indicated that when an atmospheric correction procedure for ocean color data was available, the synoptic coverage ocean color data could be used for retrieving Zsd from the optical complex coastal waters using CSSD model.  Table 1.
It is well known that the MODIS-recorded signs were easily saturated in the extremely turbid waters due to the strong water-leaving contributions, which in turn led to the loss of image coverage in the extremely turbid waters. Table 1 showed that the minimum Zsd in synchronized data was 0.5 m, which was much larger than the initialization and test datasets. When the Zsd was very small while the model had substantial errors, the MSPD value would become extremely large. Figures 2 and 3 showed that the scatterplots at the higher ender were much tighter than the lower end. This was to say that the four models in higher end were more effective than the lower end. Due to those reasons, the performance of matchup analysis was slightly better than the test results. The above match-up analysis using synchronized data provided only general guidance on the satellite-derived Zsd products, as the match-up data cannot cover all the variations of the natural waters in the CECZ seas, especially for extremely turbid water due to the signal saturation of MODIS sensor.

Noise Tolerance: IDAS-based Versus QAA-based CSSD Model
To overview the influences of residual error in satellite-derived Rrs data on CSSD model in Zsd retrievals, random noises were added into the Rrs data at all bands. Figures 5a-b showed the sensitivity analysis results of test 1, and Figures 5c-d represented the results of test 2. As expected, compared with the error-free data, the random <±0.0003 sr −1 and <±0.0006 sr −1 for △Rrs(667) could result in, respectively, ~10% and ~21% MSPD values for Rrs(488) data. Fortunately, the MSPD values for CSSD model-derived Zsd data were only ~9% and ~25%, respectively, for two tests. Therefore, the CSSD model was a quite noise tolerant model for addressing satellite-derived Rrs data even for turbid coastal waters. Note that the above sensitivity analysis using synthetic error-included data provided  Table 1.
It is well known that the MODIS-recorded signs were easily saturated in the extremely turbid waters due to the strong water-leaving contributions, which in turn led to the loss of image coverage in the extremely turbid waters. Table 1 showed that the minimum Z sd in synchronized data was 0.5 m, which was much larger than the initialization and test datasets. When the Z sd was very small while the model had substantial errors, the MSPD value would become extremely large. Figures 2 and 3 showed that the scatterplots at the higher ender were much tighter than the lower end. This was to say that the four models in higher end were more effective than the lower end. Due to those reasons, the performance of matchup analysis was slightly better than the test results. The above match-up analysis using synchronized data provided only general guidance on the satellite-derived Z sd products, as the match-up data cannot cover all the variations of the natural waters in the CECZ seas, especially for extremely turbid water due to the signal saturation of MODIS sensor.

Noise Tolerance: IDAS-based Versus QAA-based CSSD Model
To overview the influences of residual error in satellite-derived R rs data on CSSD model in Z sd retrievals, random noises were added into the R rs data at all bands. Figure 5a,b showed the sensitivity analysis results of test 1, and Figure 5c,d represented the results of test 2. As expected, compared with the error-free data, the random <±0.0003 sr −1 and <±0.0006 sr −1 for ∆R rs (667) could result in, respectively,~10% and~21% MSPD values for R rs (488) data. Fortunately, the MSPD values for CSSD model-derived Z sd data were only~9% and~25%, respectively, for two tests. Therefore, the CSSD model was a quite noise tolerant model for addressing satellite-derived R rs data even for turbid coastal waters. Note that the above sensitivity analysis using synthetic error-included data provided only a general guidance on the stability of CSSD model, as the simulated data cannot cover all the variations of the natural waters.  Many band ratio-based models including empirical, quasi-analytical, or semi-analytical methods could be effectively used to retrieve IOPs from both ocean and coastal waters [41,42], but we limited the exercise and show the performance of CSSD model with QAA version 5 [36] model-derived IOPs data because it is simple and easy to use. Figure 6a showed the performance of CSSD model with QAA model's IOPs in predicting Zsd data from the CECZ seas. It was found that even though the QAA model was quite different from IDAS model [22,36], the QAA-based CSSD model (MSPD = 31.64%) produced a slightly better performance than the IDAS-based results from in situ measurements (MSPD = 35.73% in Figure 3d). However, the band ratio parameters were only insensitive to the scaling type error, but sensitive to the addition type error such as residual error in satellite-derived Rrs data [22,23], which in turn impacted on the outputs of QAA-based CSSD model. As a result, based on the 1000 synthetic dataset in test 1, Figure 6b showed that the MSPD of QAA-based results was almost comparable with the IDAS-based result for test 2, but the residual error in Rrs data for test 2 is two times larger than test 1. Specifically, for the 7 points at the end (Zsd<0.2), the Zsd values were significantly overestimated by QAA-based CSSD model. Many band ratio-based models including empirical, quasi-analytical, or semi-analytical methods could be effectively used to retrieve IOPs from both ocean and coastal waters [39,40], but we limited the exercise and show the performance of CSSD model with QAA version 5 [35] model-derived IOPs data because it is simple and easy to use. Figure 6a showed the performance of CSSD model with QAA model's IOPs in predicting Z sd data from the CECZ seas. It was found that even though the QAA model was quite different from IDAS model [21,35], the QAA-based CSSD model (MSPD = 31.64%) produced a slightly better performance than the IDAS-based results from in situ measurements (MSPD = 35.73% in Figure 3d). However, the band ratio parameters were only insensitive to the scaling type error, but sensitive to the addition type error such as residual error in satellite-derived R rs data [21,22], which in turn impacted on the outputs of QAA-based CSSD model. As a result, based on the 1000 synthetic dataset in test 1, Figure 6b showed that the MSPD of QAA-based results was almost comparable with the IDAS-based result for test 2, but the residual error in R rs data for test 2 is two times larger than test 1. Specifically, for the 7 points at the end (Z sd < 0.2), the Z sd values were significantly overestimated by QAA-based CSSD model.
As there was no available Z sd data in the synchronized dataset provided by SeaBASS, the MSPD and RMSE log values of CSSD model were calculated by comparison of MODIS versus SeaWiFS-derived Z sd data. This was true for evaluating the sensitivity of CSSD model to the residual error in satellite-derived R rs data, as the inter-mission difference between MODIS and SeaWiFS-derived Z sd data would increase with the increase of sensitivity of model due to the existing inter-mission difference in R rs data. The fact that the IDAS-based CSSD model can minimize the residual error in satellite R rs data was demonstrated further with synchronized R rs data of MODIS versus with SeaWiFS sensor provided by SeaBASS. When a Z sd model was insensitive to the residual error in satellite data, this model could decrease the biases in Z sd value between one pair synchronized observations using two different satellites. Specifically, the MSPD and RMSE log of IDAS-based CSSD model between synchronized MODIS and SeaWiFS-derived Z sd data were 18% and 0.072 m, respectively, which were~11% and 0.022 m lower than the QAA-based model (see Figure 7). These results indicated that the IDAS model's noise tolerant property was one of the important factors to make the CSSD model be an accurate and reliable Z sd estimator for noise polluted ocean color data.
even though the QAA model was quite different from IDAS model [22,36], the QAA-based CSSD model (MSPD = 31.64%) produced a slightly better performance than the IDAS-based results from in situ measurements (MSPD = 35.73% in Figure 3d). However, the band ratio parameters were only insensitive to the scaling type error, but sensitive to the addition type error such as residual error in satellite-derived Rrs data [22,23], which in turn impacted on the outputs of QAA-based CSSD model. As a result, based on the 1000 synthetic dataset in test 1, Figure 6b showed that the MSPD of QAA-based results was almost comparable with the IDAS-based result for test 2, but the residual error in Rrs data for test 2 is two times larger than test 1. Specifically, for the 7 points at the end (Zsd<0.2), the Zsd values were significantly overestimated by QAA-based CSSD model.  Figure 6b were calculated by comparison of Zsd derived from error-included versus error-free Rrs data.
As there was no available Zsd data in the synchronized dataset provided by SeaBASS, the MSPD and RMSElog values of CSSD model were calculated by comparison of MODIS versus SeaWiFS-derived Zsd data. This was true for evaluating the sensitivity of CSSD model to the residual error in satellite-derived Rrs data, as the inter-mission difference between MODIS and SeaWiFS-derived Zsd data would increase with the increase of sensitivity of model due to the existing inter-mission difference in Rrs data. The fact that the IDAS-based CSSD model can minimize the residual error in satellite Rrs data was demonstrated further with synchronized Rrs data of MODIS versus with SeaWiFS sensor provided by SeaBASS. When a Zsd model was insensitive to the residual error in satellite data, this model could decrease the biases in Zsd value between one pair synchronized observations using two different satellites. Specifically, the MSPD and RMSElog of IDAS-based CSSD model between synchronized MODIS and SeaWiFS-derived Zsd data were 18% and 0.072 m, respectively, which were ~11% and 0.022 m lower than the QAA-based model (see Figure 7). These results indicated that the IDAS model's noise tolerant property was one of the important factors to make the CSSD model be an accurate and reliable Zsd estimator for noise polluted ocean color data.  Figure 8 showed the climatologically monthly Zsd data derived from MODIS daily Rrs data from July 2002 to September, 2014 following the NASA default data processing habits [43]. As expected, the Zsd exhibited a wide varying range (0.01-39 m) in the CECZ seas. The higher values, such as those exceeding 10 m, were found in the southeast CECZ seas, which was a great distance away from the shoreline. Meanwhile, the lower values (<1 m) exhibited in the coastal zones nearby the coastline. These low Zsd regions were mainly due to the rivers-discharged suspended, and the tidal currents and winds-suspended sedimentations, erosion sedimentations from coastline and seabed, as well as other factors. As a branch of Kuroshio current, the northwestwardly Yellow Sea Warm Current pumped the clear and warm sea waters into the noticeably cool and turbid waters in the Yellow Sea and Bohai Sea, while the southeastwardly CECZ Coastal Currents pumped the turbid water from Bohai Sea to East China Sea, so the Zsd value gradually decreased from southeast to northwest in the central CECZ seas. Due to the mixing between turbid and clear waters and the large amount of river-discharged sediments pouring into the coastal zones, many tongue-shaped plumes could be found in coastal zones from Zsd maps (see Figure 8).  Figure 8 showed the climatologically monthly Z sd data derived from MODIS daily R rs data from July 2002 to September, 2014 following the NASA default data processing habits [41]. As expected, the Z sd exhibited a wide varying range (0.01-39 m) in the CECZ seas. The higher values, such as those exceeding 10 m, were found in the southeast CECZ seas, which was a great distance away from the shoreline. Meanwhile, the lower values (<1 m) exhibited in the coastal zones nearby the coastline. These low Z sd regions were mainly due to the rivers-discharged suspended, and the tidal currents and winds-suspended sedimentations, erosion sedimentations from coastline and seabed, as well as other factors. As a branch of Kuroshio current, the northwestwardly Yellow Sea Warm Current pumped the clear and warm sea waters into the noticeably cool and turbid waters in the Yellow Sea and Bohai Sea, while the southeastwardly CECZ Coastal Currents pumped the turbid water from Bohai Sea to East China Sea, so the Z sd value gradually decreased from southeast to northwest in the central CECZ seas. Due to the mixing between turbid and clear waters and the large amount of river-discharged sediments pouring into the coastal zones, many tongue-shaped plumes could be found in coastal zones from Z sd maps (see Figure 8). Remote Sens. 2019, 11, x FOR PEER REVIEW 12 of 18 Figure 8. Climatology monthly and spatial change patterns of Zsd in CECZ seas derived from the MODIS daily Rrs data during July 2002 to September 2014 following the NASA default data processing habits [43]. White color represents clouds or lands. All data within Rrs(555)>0.0 sr −1 are included in the analysis.

Spatial and Monthly Patterns of Secchi Depth
By comparison, the exhibited Zsd value in the Bohai Sea was lower than both the Yellow and East China Seas. This unique spatial and temporal distribution characteristics for Zsd was mostly controlled by both the riverine and coastal currents, which included the river discharges and the mixing of fresh and salt waters, as well as others. Specifically, there were at least one small coastal gyres in the Bohai Bay, Laizhou Bay, and Liaodong Bay [44]. These gyres reduced the waters exchange rate between bay and coastal current, so the Bohai Bay, Laizhou Bay, and Liaodong Bay were perennially filled with high sediment and pollutant waters. Moreover, the suspended sediment was mostly distributed around the southern Laizhou Bay, particularly around the Yellow River estuary in the summer season. It was impacted by the river plume, yet confined to a very limited area around the mouth of the Yellow River, whereas in the winter-spring seasons the sediment concentration became much higher as a result of the active coastal re-suspension, which was induced by the energetic wave actions in the shallow water [45]. Seasonally, the lowest mean Zsd in the Bohai Sea (roughly around 1.63 m) was found in February, while the highest (roughly around 3.54 m) was observed in August, which is consistent with the variations of monthly averaged suspended sediment concentration observed in the south Bohai Strait [46]. The Bohai strait is roughly divided into two parts by Huangcheng Island. The relatively clear northern Yellow Sea water discharges into the Bohai Sea through the northern Bohai strait, and the turbid Bohai Sea waters flow into the Yellow Sea through the southern section [46]. Due to the physical mixing process of the turbid Bohai Sea waters and the clear Yellow Sea waters, a clear tongue-shaped front was found around the Bohai strait in Zsd maps. Figure 8. Climatology monthly and spatial change patterns of Z sd in CECZ seas derived from the MODIS daily R rs data during July 2002 to September 2014 following the NASA default data processing habits [41]. White color represents clouds or lands. All data within R rs (555) > 0.0 sr −1 are included in the analysis.
By comparison, the exhibited Z sd value in the Bohai Sea was lower than both the Yellow and East China Seas. This unique spatial and temporal distribution characteristics for Z sd was mostly controlled by both the riverine and coastal currents, which included the river discharges and the mixing of fresh and salt waters, as well as others. Specifically, there were at least one small coastal gyres in the Bohai Bay, Laizhou Bay, and Liaodong Bay [42]. These gyres reduced the waters exchange rate between bay and coastal current, so the Bohai Bay, Laizhou Bay, and Liaodong Bay were perennially filled with high sediment and pollutant waters. Moreover, the suspended sediment was mostly distributed around the southern Laizhou Bay, particularly around the Yellow River estuary in the summer season. It was impacted by the river plume, yet confined to a very limited area around the mouth of the Yellow River, whereas in the winter-spring seasons the sediment concentration became much higher as a result of the active coastal re-suspension, which was induced by the energetic wave actions in the shallow water [25]. Seasonally, the lowest mean Z sd in the Bohai Sea (roughly around 1.63 m) was found in February, while the highest (roughly around 3.54 m) was observed in August, which is consistent with the variations of monthly averaged suspended sediment concentration observed in the south Bohai Strait [43]. The Bohai strait is roughly divided into two parts by Huangcheng Island. The relatively clear northern Yellow Sea water discharges into the Bohai Sea through the northern Bohai strait, and the turbid Bohai Sea waters flow into the Yellow Sea through the southern section [43]. Due to the physical mixing process of the turbid Bohai Sea waters and the clear Yellow Sea waters, a clear tongue-shaped front was found around the Bohai strait in Z sd maps.
Runoff from the Changjiang River and the upwelling of Kuroshio subsurface water are two major sources of materials and nutrients that sustain the primary production and optical properties in the East China Sea [44]. The Changjiang River's discharge contained large amounts of nutrients, sedimentation, and pollutants, so the Changjiang diluted waters were always more fertile and turbid than the adjacent sea waters. As a result, a clear tongue-shaped structure could be observed in the Changjiang River estuary from Z sd maps (see Figure 8). However, the transport path of Changjiang diluted waters was not only depended on the seasonal changed river's discharge amounts [45], but also related to the seasonal monsoon and Kuroshio current in the East China Sea [46], so the coverage domain and shape of frontal edge associated with Changjiang diluted waters exhibited a pronounced seasonal pattern. The spring and summer discharge of the Changjiang River is large enough to form a surface plume of relatively fresh water which can extend at least 300 km out over the shelf, making the tongue-shaped structure extending toward the southeast to the location of 31 • N and 127 • E. In this period, the frontiers of the tongue-shape edge gradually became clearer, but the coverage domain slowly shrunk. From summer to autumn, the tongue-shaped structure's tip gradually changed from southeastward to northeastward direction, and reached the location near Jeju Island in August. During the winter season, the much weaker Changjiang River discharge is generally confined to flow primarily southward in a narrow band along the Chinese coast. As a result, the tongue-shaped structure's front became blurry. Also in this season, as the winter monsoon become stronger [47], the tongue-shaped structure's front gradually changed from northeastward direction to southeastward direction, and the coverage domains slowly increased with time. Figure 9 showed the monthly changing phases of TSI in Bohai Sea, Yellow Sea, and East China Sea during 2002-2014, derived from the daily global composite R rs data. It was found that Bohai Sea, Yellow Sea, and East China Sea exhibited similar seasonal pattern: the maximal TSI value was found around February, while the minimum TSI value was found around August. However, the inter-annual phases in these three seas were quite different from each other. Specifically, the 13-year average of TSI values were 48.72, 30.82, and 25.57 m, respectively, for Bohai Sea, Yellow Sea, and East China Sea. Following the classification method for eutrophication level proposed by Carlson [23] in Table 2, the water quality in Bohai Sea, Yellow Sea, and East China Sea fall into mesotrophic, mesotrophic, and oligotrophic types respectively. This complicated phase change might be caused by many factors including the Three Gorges Dam, global climate changes, river discharges, wind speed, erosion rate of sedimentation in the intertidal zones, and others [47,48], but this is beyond the scope of this study. Table 2. Classification of water quality using Z sd -based trophic state index (TSI) proposed by Carlson [24].   Figure 10 showed the spatial and temporal changing patterns of water quality in the CECZ seas, indicating that the water qualities in summer season were much better than the spring season. Generally, the eutrophic waters exhibited within 35 m isobaths, while the mesotrophic waters distributed within 85 m isobaths. Little mesotrophic water can be found from zones deeper than 30 m during 2002-2005. However, the coverage range of mesotrophic waters has extended toward to 35 m isobaths since 2005 with the increasing pace (up to 1 m per year during 2010-2014). The eutrophic waters mainly distributed within 5 and 25 m isobaths in summer and winter seasons, respectively. Similar to mesotrophic waters, the coverage domain of eutrophic waters has also extended toward to the outer shelf, but the extension rate was slightly slower than the mesotrophic waters. Overall, more and more deep-waters have the risk of becoming mesotrophic and eutrophic waters in the CECZ seas.  Figure 10 showed the spatial and temporal changing patterns of water quality in the CECZ seas, indicating that the water qualities in summer season were much better than the spring season. Generally, the eutrophic waters exhibited within 35 m isobaths, while the mesotrophic waters distributed within 85 m isobaths. Little mesotrophic water can be found from zones deeper than 30 m during 2002-2005. However, the coverage range of mesotrophic waters has extended toward to 35 m isobaths since 2005 with the increasing pace (up to 1 m per year during 2010-2014). The eutrophic waters mainly distributed within 5 and 25 m isobaths in summer and winter seasons, respectively. Similar to mesotrophic waters, the coverage domain of eutrophic waters has also extended toward to the outer shelf, but the extension rate was slightly slower than the mesotrophic waters. Overall, more and more deep-waters have the risk of becoming mesotrophic and eutrophic waters in the CECZ seas.

Discussion
The models based on the blue/green Rrs ratios usually produced a poor performance in the productive and optically complex coastal waters, despite the fact that these models may produce a good performance in quantifying the Zsd in open ocean. The reason for this may be attributed not only to the fundamental algorithm structure, but also to the variability in the relationships between optically active constituents in different optical water types [1,53]. Specifically, due to the strong absorption of gelbstoff and detritus in the turbid coastal waters, the Rrs at blue band is too low to apply for ocean color remote sensing [1]. However, the Rrs at red and near-infrared bands are very sensitive to the changes of detritus in the turbid coastal waters [1], so the Rrs at longer bands are more suitable than the shorter bands for Zsd retrievals in the turbid coastal waters. More importantly, the impact of atmospheric correction on the water-leaving radiance or Rrs was not spectrally uniform, so the estimation of the Zsd derived by the blue-green ratios can be significant [29]. Fundamentally, ~90% of the top-of-atmosphere radiance is contributed by the atmosphere, and consequently a small error in the atmospheric correction would translate to a much larger relative error in the Rrs [22]. As a result, the BGSD model, which was based on the blue/green Rrs, may not be able to accurately estimate the Zsd in the CECZ.
Semi-analytical models could be deemed as an approximated solution for the radiative transfer equation in the waters based on a few simplified assumptions. As stated by Odermatt et al. [54], the advantage of semi-analytical over empirical models was a matter of adaptivity and justification rather than accuracy. Application of semi-analytical models implied a physical sound procedure with defined Rrs from field measurements or atmospheric correction procedure. To semi-analytically evaluate the desired Zsd from satellite-derived Rrs, it was necessary to know IOPs and Kd(λ) from satellite data. However, due to the imperfect IOPs and Kd(λ) retrieval models, the accuracy and stability of IOPSD and MKSD models were also questionable. Especially for extremely turbid waters, it was hard to accurately determine the IOPs and Kd(λ) using ocean color technology, which in turn led to the unreliable results of Zsd. Thus, we classified the waters into three different groups, and the corresponding Zsd were respective retrieved using three different models. The results showed that the CSSD model exhibits more effective than both MKSD and IOPSD models in deriving Zsd from the CECZ seas.
There were only two samples in the test dataset collected from clear waters whose Zsd values were around 20 m, even though the MSPD values of CSSD model for those two samples did not exceed 28%. In the open ocean, Zsd were suggested for retrieving from IDAS IOPs data in CSSD model. Figures 2-5 showed that the CSSD model was able to provide Zsd data with high quality as long as the data quality of model inputs was not too poor. The IDAS model was originally designed

Discussion
The models based on the blue/green R rs ratios usually produced a poor performance in the productive and optically complex coastal waters, despite the fact that these models may produce a good performance in quantifying the Z sd in open ocean. The reason for this may be attributed not only to the fundamental algorithm structure, but also to the variability in the relationships between optically active constituents in different optical water types [49,50]. Specifically, due to the strong absorption of gelbstoff and detritus in the turbid coastal waters, the R rs at blue band is too low to apply for ocean color remote sensing [1]. However, the R rs at red and near-infrared bands are very sensitive to the changes of detritus in the turbid coastal waters [1], so the R rs at longer bands are more suitable than the shorter bands for Z sd retrievals in the turbid coastal waters. More importantly, the impact of atmospheric correction on the water-leaving radiance or R rs was not spectrally uniform, so the estimation of the Z sd derived by the blue-green ratios can be significant [28]. Fundamentally,~90% of the top-of-atmosphere radiance is contributed by the atmosphere, and consequently a small error in the atmospheric correction would translate to a much larger relative error in the R rs [21]. As a result, the BGSD model, which was based on the blue/green R rs , may not be able to accurately estimate the Z sd in the CECZ.
Semi-analytical models could be deemed as an approximated solution for the radiative transfer equation in the waters based on a few simplified assumptions. As stated by Odermatt et al. [51], the advantage of semi-analytical over empirical models was a matter of adaptivity and justification rather than accuracy. Application of semi-analytical models implied a physical sound procedure with defined R rs from field measurements or atmospheric correction procedure. To semi-analytically evaluate the desired Z sd from satellite-derived R rs , it was necessary to know IOPs and K d (λ) from satellite data. However, due to the imperfect IOPs and K d (λ) retrieval models, the accuracy and stability of IOPSD and MKSD models were also questionable. Especially for extremely turbid waters, it was hard to accurately determine the IOPs and K d (λ) using ocean color technology, which in turn led to the unreliable results of Z sd . Thus, we classified the waters into three different groups, and the corresponding Z sd were respective retrieved using three different models. The results showed that the CSSD model exhibits more effective than both MKSD and IOPSD models in deriving Z sd from the CECZ seas.
There were only two samples in the test dataset collected from clear waters whose Z sd values were around 20 m, even though the MSPD values of CSSD model for those two samples did not exceed 28%. In the open ocean, Z sd were suggested for retrieving from IDAS IOPs data in CSSD model. Figures 2-5 showed that the CSSD model was able to provide Z sd data with high quality as long as the data quality of model inputs was not too poor. The IDAS model was originally designed for deriving IOPs and residual error from satellite R rs in the open oceans [20,21]. Evaluated with satellite observations and field measurements, Chen et al. [21] showed that the IDAS model could provide IOPs products with high quality. Thus, the CSSD model would be robust for Z sd retrievals in the open oceans, even though further evaluations with different datasets are still required in the global open oceans.

Conclusion
It is challenging to accurately estimate Z sd from turbid coastal waters using ocean color technology. In this study, we develop and evaluate a semi-analytical model for estimations of the ocean transparency from space with ocean color data. The band difference method is coupled into the Z sd model to absorb the residual error in satellite R rs data. The model is evaluated by in situ measured dataset, synthetic dataset, and synchronized dataset. The results show that the CSSD model can accurately quality Secchi depth from R rs data, and also can minimize the effects of residual errors in satellite R rs on model's output.
The CSSD model is proposed for retrieving Z sd and TSI from the climatologically monthly R rs data. As expected, low Z sd values are found around the river estuaries, while the high values are found in the far offshore of the CECZ seas. According to the corresponding TSI products, the water quality in the Bohai Sea and Yellow Sea has become more eutrophic during past 13 years, and might continue to deteriorate in the near future. If the current trophic state continues to evolve, it only needs 14 and 67 years for the Bohai Sea and Yellow Sea to fall into the eutrophic waters, respectively. Moreover, our results also show that the trophic state in the deep-waters has been significantly impacted by the trophic state in the shallow regions. These findings have warned us that tighter regulations should be enforced by China Government for protecting the aquatic environment in the CECZ seas.