Remote Sensing of Secchi Depth in Highly Turbid Lake Waters and Its Application with MERIS Data

The Secchi disk depth (ZSD, m) has been used globally for many decades to represent water clarity and an index of water quality and eutrophication. In recent studies, a new theory and model were developed for ZSD, which enabled its semi-analytical remote sensing from the measurement of water color. Although excellent performance was reported for measurements in both oceanic and coastal waters, its reliability for highly turbid inland waters is still unknown. In this study, we extend this model and its evaluation to such environments. In particular, because the accuracy of the inherent optical properties (IOPs) derived from remote sensing reflectance (Rrs, sr−1) plays a key role in determining the reliability of estimated ZSD, we first evaluated a few quasi-analytical algorithms (QAA) specifically tuned for turbid inland waters and determined the one (QAATI) that performed the best in such environments. For the absorption coefficient at 443 nm (a(443), m−1) ranging from ~0.2 to 12.5 m−1, it is found that the QAATI-derived absorption coefficients agree well with field measurements (r2 > 0.85, and mean absolute percentage difference (MAPD) smaller than ~39%). Furthermore, with QAATI-derived IOPs, the MAPD was less than 25% between the estimated and field-measured ZSD (r2 > 0.67, ZSD in a range of 0.1–1.7 m). Furthermore, using matchup data between Rrs from the Medium Resolution Imaging Spectrometer (MERIS) and in-situ ZSD, a similar performance in the estimation of ZSD from remote sensing was obtained (r2 = 0.73, MAPD = 37%, ZSD in a range of 0.1–0.9 m). Based on such performances, we are confident to apply the ZSD remote sensing scheme to MERIS measurements to characterize the spatial and temporal variations of ZSD in Lake Taihu during the period of 2003–2011.


Introduction
A white or black-and-white disc (diameter 20 to 30 cm) is commonly lowered into the water to measure water clarity (or transparency) of aquatic environments, and the depth at which it just disappears from a viewer at the surface is called the Secchi disk depth (Z SD , m) [1].The measurement of Z SD originated in the nineteenth century and is the earliest oceanographic measurement of water clarity.Z SD is a useful index of water quality and eutrophication [1][2][3].Changes in Z SD values have been widely used to describe variations in water properties and light availability [4,5], and these changes have been linked to variations in eutrophication and phytoplankton biomass [6].Although different types of sophisticated electro-optical systems are currently employed to measure water clarity or quality, Secchi disks are still regularly used in both limnology and oceanography studies because of their low cost and easy operation [7].
Accurate determination of Z SD can be obtained in-situ, but this approach cannot provide data in a synoptic manner.Instead, satellite remote sensing is widely used to map the distribution of Z SD values over broad areas.In recent decades, many remote-sensing algorithms have been developed to estimate Z SD in various water types [8][9][10][11][12][13][14].These algorithms can be briefly categorized as follows: (1) empirical relationships between Z SD and remote-sensing reflectance (R rs , sr −1 ) [8,9,15] and (2) semi-analytical models that are based on the relationship between Z SD and inherent optical properties (IOPs), with the later derived from R rs [11][12][13].Compared to empirical algorithms, these semi-analytical algorithms [12] are based on the classical theory regarding Secchi disk observation [13] and presumably not restricted to the region of waters.
Recently, a new theoretical model was developed to estimate Z SD [16,17], which was verified with in-situ Z SD measurements in a wide range of water types.In this new model, Z SD is related to the diffuse attenuation of downwelling irradiance (K d , m −1 ) and R rs in water's spectral transparent window in the visible domain (~400-700 nm).K d is a function of IOPs (total absorption coefficient, a (m −1 ), and the backscattering coefficient, b b (m −1 )) [18], and these IOPs can be estimated semi-analytically from R rs [19][20][21], as can Z SD [16].Since the IOPs play a key role in this process, the accuracy of the remote sensing of Z SD mainly depends on the accuracy of the IOPs derived from R rs .
In the scheme of Lee et al. [16], the Quasi-Analytical Algorithm (QAA, version 6) was adopted to obtain the values of IOPs.The current version of the QAA performs very well in most open ocean and coastal waters [22], but the estimated IOPs could be questionable when this algorithm is applied to highly turbid or eutrophic waters [23,24].In recent decades, substantial efforts have been made to tune QAA for turbid eutrophic waters by shifting the reference wavelength to a longer wavelength or re-parameterize the semi-analytical relationships [23,25,26].However, those studies are still limited to small areas [23] with a relatively small range of the total absorption coefficient, (e.g., a(443) from ~0.5 to 8.0 m −1 [25]) or chlorophyll concentrations, (e.g., from ~0.5 to 12.9 µg/L [26]).For highly turbid waters, the IOPs estimated via these QAA still have relatively large uncertainties [24].
Therefore, before application of the new Z SD scheme to such turbid or eutrophic waters, it is imperative to revise the QAA in order to obtain more accurate estimates of IOPs.Specifically, the empirical relationships employed in QAA_v6 are modified first with a dataset of great dynamic range (chlorophyll concentration in a range of 4.9 to 412.7 µg/L) measured in Lake Taihu, China.This revised QAA for turbid inland waters (QAA TI ) is subsequently validated with two independent datasets from Lake Taihu and Lake Erie, USA, which are typical turbid lake waters with a(443) ranging from ~0.2 to 12.5 m −1 .Furthermore, with IOPs derived from QAA TI , values of Z SD of such waters were derived from both field R rs and from the medium resolution imaging spectrometer (MERIS) following the new Z SD scheme.These estimates were then evaluated with in-situ Z SD measurements.Lastly, using a 2003-2011 time series of R rs from MERIS, we characterized the spatial and temporal variations of water clarity in Lake Taihu, and the potential factors associated with such variation are discussed.

Lake Taihu (LT) Dataset
This dataset (LT) was divided into two sub-sets (LT-A and LT-B in Table 1, Figure 1a).One was used to develop QAA TI by tuning QAA_v6, and the other was used for evaluation of the QAA TI .There are 250 sites in this dataset, including concurrent measurements of hyperspectral R rs , a(λ), and Z SD .All the measurements were generally carried out between 09:30 and 16:00 when the sky was cloudless.
A total of five surveys were conducted with 50 sampling sites in each survey.For dataset LT-A, there are 150 sites collected in three surveys (January, July, and October in 2006).For dataset LT-B, there are 100 sites collected in two separate months (January and April in 2007).

Measurements Number Usage
LT-A(Figure 1a) In-situ R rs , a(λ) and Z SD 150 QAA retune LT-B (Figure 1a) In-situ R rs , a(λ) and Z SD 100 Validation Erie (Figure 1b) In-situ R rs , a(λ) and Z SD

Satellite Dataset
Another dataset (SAT) was further used to validate the ZSD derived from MERIS.This dataset includes 74 MERIS images and 85 matchups between MERIS Rrs and in-situ ZSD measurements.A total of 74 Level 1 (L1) MERIS images of Lake Taihu from 2003 to 2011 were downloaded from the European Space Agency (ESA) Earth Online (http://earth.esa.int/).The images were then processed using BEAM 5.0 software for geometric referencing and smile correction.There are three plug-in processors provided by the BEAM 5.0 for atmospheric correction of MERIS level 1 data, namely Case-2 Regional Processer (C2R), Boreal Lakes Processor (BOR), and Eutrophic Lakes Processor (EUL).Both C2R and BOR processors were not developed for use in highly eutrophic lake waters [33,34].Therefore, the atmospheric correction was performed by the artificial neural network scheme included in the EUL processor.Furthermore, a Phycocyanin index (PCI) [35] was used to eliminate the pixels with thick Microcystis bloom.A total of 85 matchups (time difference less than one day) between satellite data and the in-situ ZSD were obtained, where the in-situ ZSD were collected during six field surveys (see Table 2).The total absorption coefficient of the bulk water was calculated as the sum of colored dissolved organic matter (CDOM), particulate matter (phytoplankton and tripton), and pure seawater, where the absorption coefficient of pure seawater (a w (λ), m −1 ) was a constant [27,28].In this study, the absorption coefficient of CDOM (a CDOM , m −1 ) and total particle absorption (a p , m −1 ) were measured using the quantitative filter technique (QFT) following the National Aeronautics and Space Administration (NASA) ocean optical measurement protocol [29].The procedures are described in Liu et al. [30].Moreover, the concentrations of suspended particulate matter (SPM, mg/L) and chlorophyll a (including phaeopigments) (Chla, µg/L) were also determined from water samples, and the procedures were described in Liu et al. [30].In this study, the SPM ranged from ~10.3 to 285.6 mg/L, and the Chla ranged from ~4.9 to 412.7 µg/L.For the absorption coefficients, the a(443) ranged from ~1.4 to 12.5 m −1 .A standard Secchi disk with a 30-cm diameter was used to measure Z SD , which ranged from ~0.1 to 0.9 m.
The measurement of R rs was made following the above-surface approach [31] with an ASD (Analytical Spectral Devices) field spectrometer (Analytical Devices, Inc., Boulder, CO, USA).The radiometer measurements include the total-above surface radiance (L t ), downwelling sky radiance (L sky ) from the reciprocal angle of measuring L t , and radiance leaving a standard grey card (L g ).The R rs was then calculated.The equation is shown below.
where ρ is the reflectance of the gray card, and F is the Fresnel reflectance of the air-sea surface, which was approximated as 2.2%.Note that an extra procedure called second-order correction was conducted for calculating R rs following the procedure of Lee et al. [32].The measurement of R rs in our study covers wavelengths up to 1050 nm, so the value of ∆ was determined by setting the average of R rs (950-1000 nm) to be 0 for each station because the absorption coefficient of pure seawater reaches a local maximum (~50 m −1 ) in this spectral range.

Lake Erie Dataset
One more dataset (Lake Erie, USA) was used to evaluate the QAA TI and the Z SD scheme.The measurements (R rs , a(λ), and Z SD ) were collected in western Lake Erie (Figure 1b), and the data were downloaded from the SeaWiFS Bio-optical Archive and Storage System (SeaBASS).A total of 14 sites was included in this dataset, where Z SD values are generally less than 2 m (0.4-1.7 m) with a(443) ranging from ~0.2 to 3.9 m −1 .
Above-surface radiance measurements were taken with the Field Spec Pro™ VNIR-NIR1 portable spectrometer system from ASD.The R rs values were derived following Equation (1).

Satellite Dataset
Another dataset (SAT) was further used to validate the Z SD derived from MERIS.This dataset includes 74 MERIS images and 85 matchups between MERIS R rs and in-situ Z SD measurements.A total of 74 Level 1 (L1) MERIS images of Lake Taihu from 2003 to 2011 were downloaded from the European Space Agency (ESA) Earth Online (http://earth.esa.int/).The images were then processed using BEAM 5.0 software for geometric referencing and smile correction.There are three plug-in processors provided by the BEAM 5.0 for atmospheric correction of MERIS level 1 data, namely Case-2 Regional Processer (C2R), Boreal Lakes Processor (BOR), and Eutrophic Lakes Processor (EUL).Both C2R and BOR processors were not developed for use in highly eutrophic lake waters [33,34].Therefore, the atmospheric correction was performed by the artificial neural network scheme included in the EUL processor.Furthermore, a Phycocyanin index (PCI) [35] was used to eliminate the pixels with thick Microcystis bloom.A total of 85 matchups (time difference less than one day) between satellite data and the in-situ Z SD were obtained, where the in-situ Z SD were collected during six field surveys (see Table 2).
To understand the impact of wind speed on the variation of Z SD , the long-term wind speed data of Lake Taihu from 2003 to 2011 were downloaded from the China Meteorological Data Sharing Service System (http://data.cma.cn/).The wind data were collected at the site of Huzhou meteorological station (30 • 52 N, 120 • 3 E) beside Lake Taihu.The wind data were recorded every 5 min.The daily-averaged wind speed was calculated from all available wind speed data on that day, and the monthly-averaged data were then calculated from the daily-averaged wind speed.

Accuracy Assessment
The performance of an algorithm was quantified based on the mean absolute percentage difference (MAPD), which is defined below.

MAPD =
100 Remote Sens. 2019, 11, 2226 5 of 19 where x is the evaluated variable, and n is the number of data points.

Semi-Analytical Z SD Scheme
Based on a new law of contrast reduction, Lee et al. [16] proposed an innovative theoretical model to infer Z SD .Unlike the classical model [12,13] that primarily relies on the beam attenuation coefficient, this model only relies on K d within the spectral transparent window, and a simplified equation is shown below.
where Min(K d (λ)) represents the minimum K d for wavelengths between 443 nm and 665 nm, R pc rs is the R rs value that corresponds to the wavelength with the minimum K d , and C r t is the contrast threshold of the human eye in radiance reflectance and is treated as a constant (0.013 sr −1 ).Thus, the estimate of Z SD from remote sensing primarily depends on K d (λ), which can be calculated following the K d model of Lee et al. [18] from a(λ) and b b (λ).In the Z SD scheme of Lee et al. [18], the values of a(λ) and b b (λ) were retrieved from R rs with QAA_v6.
In order to obtain accurate estimates of Z SD in turbid inland water, it is critical to have accurate estimation of IOPs from R rs , but the challenge is that the QAA_v6 and other QAA tuned versions [23,25,26] have difficulty in providing reliable IOPs in such water, which subsequently yields Z SD with high uncertainties (see Section 3.1.1).To extend the applicability of QAA, we provide a revised version described in Section 3.1.2.

Evaluation of Current IOPs Estimation Approaches for Turbid Waters
Some efforts have been made to extend the original QAA to retrieve IOPs in turbid waters.For instance, Le et al. [23] attempted to improve the QAA in Meiliang Bay in Lake Taihu by shifting the reference wavelength from the red region to near-infrared (called QAA.Le hereafter).Huang et al.
[26] developed a regional QAA (QAA.PY) in the turbid Poyang Lake based on hydro light simulations.Moreover, Yang et al. [25] proposed a revised QAA (QAA.Turbid) for turbid waters by assuming that the optical properties of the near-infrared bands were dominated by the absorption of pure water.Using the entire LT (both LT-A and LT-B) dataset and the Lake Erie dataset, the above three algorithms along with QAA_v6 (two difference reference wavelengths 670 and 710 nm were used separately) were evaluated, in order to understand their performance in turbid inland waters.This comparison of algorithms is intended to (a) characterize the performance of the many revised QAAs in turbid waters, and (b) use Z SD data as an independent data source to measure the performance of these algorithms.
Figure 2 shows histograms of the ratio of the estimated a(λ) (λ = 413, 443, 490, 510, 560, 620, and 665 nm) to in-situ measurements.For the entire LT dataset, the three algorithms (QAA_v6-670, QAA_v6-710, and QAA.Turbid) systematically underestimated the a values with a mean ratio of ~0.6.For the regional algorithms (QAA.Le and QAA.PY), although the mean ratio was relatively closer to 1 (~0.9)than the previous three algorithms, there were still some data points with extremely high ratios (>2) indicating overestimated a values.Similarly, when the above algorithms were evaluated with the other independent dataset (Erie), the absorption coefficients derived from all five algorithms were significantly overestimated with the MAPD ranging from 95.3% to 244.1% (Figure 2a-e).
Furthermore, to understand the uncertainties of derived Z SD with the IOPs estimated from the above QAA algorithms, a comparison of derived Z SD and in-situ measurements is shown in Figure 3.For the entire LT dataset, systematical overestimations were obtained for the algorithms QAA_v6-670 and QAA_v6-710, but underestimated Z SD for the algorithms QAA.Le and QAA.PY (Figure 3a-d).
For the Erie dataset, the Z SD derived by all above algorithms were significantly underestimated with the MAPD ranging from 35% to 62% (Figure 3), which is consistent with the systematical overestimation of IOPs (Figure 2a-e).Separately, the QAA.Turbid shows better performance in estimating Z SD than the other four algorithms with a MAPD of 29% and a determination coefficient (r 2 ) of 0.65 for the entire LT dataset (Figure 3e), but for the dataset Erie, the Z SD was systematically underestimated with a MAPD of 39% (Figure 3e).Furthermore, to understand the uncertainties of derived ZSD with the IOPs estimated from the above QAA algorithms, a comparison of derived ZSD and in-situ measurements is shown in Figure 3.For the entire LT dataset, systematical overestimations were obtained for the algorithms QAA_v6-670 and QAA_v6-710, but underestimated ZSD for the algorithms QAA.Le and QAA.PY (Figures 3a-3d).For the Erie dataset, the ZSD derived by all above algorithms were significantly underestimated with the MAPD ranging from 35% to 62% (Figure 3), which is consistent with the systematical overestimation of IOPs (Figures 2a-2e).Separately, the QAA.Turbid shows better performance in estimating ZSD than the other four algorithms with a MAPD of 29% and a determination coefficient (r 2 ) of 0.65 for the entire LT dataset (Figure 3e), but for the dataset Erie, the ZSD was systematically underestimated with a MAPD of 39% (Figure 3e).Furthermore, to understand the uncertainties of derived ZSD with the IOPs estimated from the above QAA algorithms, a comparison of derived ZSD and in-situ measurements is shown in Figure 3.For the entire LT dataset, systematical overestimations were obtained for the algorithms QAA_v6-670 and QAA_v6-710, but underestimated ZSD for the algorithms QAA.Le and QAA.PY (Figures 3a-3d).For the Erie dataset, the ZSD derived by all above algorithms were significantly underestimated with the MAPD ranging from 35% to 62% (Figure 3), which is consistent with the systematical overestimation of IOPs (Figures 2a-2e).Separately, the QAA.Turbid shows better performance in estimating ZSD than the other four algorithms with a MAPD of 29% and a determination coefficient (r 2 ) of 0.65 for the entire LT dataset (Figure 3e), but for the dataset Erie, the ZSD was systematically underestimated with a MAPD of 39% (Figure 3e).

Modification of QAA for Highly Turbid Waters
The QAA [36] estimates IOPs (a and b b ) in seven steps: two are analytical, three are semi-analytical, and two are empirical (QAA_v6, http://ioccg.org/resources/software/). In oceanic waters, the two empirical steps (Step 2 and 4) are considered of second-order importance because their range of variation and influence on the final output were relatively low.However, the optical properties of turbid inland waters are more complex than oceanic waters and can dramatically change even within small spatial scales.Thus, it is important and necessary to provide reliable relationships for these two steps (Step 2 and 4) in the QAA.In this study, the empirical relationships (Step 2 and 4) are revised with the field measurements collected in Lake Taihu (LT-A dataset, n = 150), and this updated QAA is termed as QAA TI herein.
In Step 2, we proposed a new relationship between the absorption coefficient at a reference wavelength (a(λ 0 )) and R rs at 560, 665, and 780 nm, which are aimed at applications of MERIS measurements.A key difference of QAA TI compared with previous algorithms is the use of 780 nm, which reflects that (1) a w dominates the total absorption coefficient at this wavelength, and (2) the R rs (780) value of these waters varies in a large range, from 0.002 to 0.025 sr −1 , mainly due to large variations in SPM (Figure 4b).Additionally, the usage of 665 nm is considered for less turbid waters.The 560 nm was employed as the reference wavelength (Figure 4a) to reflect the significantly high correlation between a(560) and the R rs ratios of the above wavelengths (see Figure 5).

Modification of QAA for Highly Turbid Waters
The QAA [36] estimates IOPs (a and bb) in seven steps: two are analytical, three are semianalytical, and two are empirical (QAA_v6, http://ioccg.org/resources/software/). In oceanic waters, the two empirical steps (Step 2 and 4) are considered of second-order importance because their range of variation and influence on the final output were relatively low.However, the optical properties of turbid inland waters are more complex than oceanic waters and can dramatically change even within small spatial scales.Thus, it is important and necessary to provide reliable relationships for these two steps (Step 2 and 4) in the QAA.In this study, the empirical relationships (Step 2 and 4) are revised with the field measurements collected in Lake Taihu (LT-A dataset, n = 150), and this updated QAA is termed as QAATI herein.
In Step 2, we proposed a new relationship between the absorption coefficient at a reference wavelength (a(λ0)) and Rrs at 560, 665, and 780 nm, which are aimed at applications of MERIS measurements.A key difference of QAATI compared with previous algorithms is the use of 780 nm, which reflects that (1) aw dominates the total absorption coefficient at this wavelength, and (2) the Rrs(780) value of these waters varies in a large range, from 0.002 to 0.025 sr −1 , mainly due to large variations in SPM (Figure 4b).Additionally, the usage of 665 nm is considered for less turbid waters.The 560 nm was employed as the reference wavelength (Figure 4a) to reflect the significantly high correlation between a(560) and the Rrs ratios of the above wavelengths (see Figure 5).
where 0.062 (m −1 ) is the absorption coefficient of pure seawater at 560 nm.
It is necessary to know the power coefficient (η) of b bp (λ) (Step 4 of QAA) in all semi-analytical algorithms [21,36,37], and, here, values of η were first derived from the measured R rs (λ) and a(λ).Specifically, b bp (λ) can be calculated from the following equation.
where u(λ) is the ratio of b b to (a + b b ), which can be easily calculated from r rs (λ) (Lee et al., 2002).b bw (λ) is the spectrum of the backscattering coefficient of pure seawater and considered to be a constant [38].After b bp (λ) is known, η is calculated by linearly fitting the logarithm of b bp (λ) between 400 and 650 nm.
After values of η are known, we used a similar empirical model in QAA_v6 to estimate η but with r rs at 443 nm and 490 nm, and Step 4 of QAA TI is written as the equation below (see Figure 5b).It is necessary to know the power coefficient (η) of bbp(λ) (Step 4 of QAA) in all semi-analytical algorithms [21,36,37], and, here, values of η were first derived from the measured Rrs(λ) and a(λ).Specifically, bbp(λ) can be calculated from the following equation.
where u(λ) is the ratio of bb to (a + bb), which can be easily calculated from rrs(λ) (Lee et al., 2002).bbw(λ) is the spectrum of the backscattering coefficient of pure seawater and considered to be a constant [38].
After bbp(λ) is known, η is calculated by linearly fitting the logarithm of bbp(λ) between 400 and 650 nm.To validate the QAA TI model, we compared the retrieved a with the measured a at the MERIS bands of 413, 443, 490, 510, 560, 620, and 665 nm using the LT-B dataset.It was found that the estimated and measured a(λ) agree with each other very well (Figure 6), where the r 2 is 0.93, and the MAPD is 15.3%.

Validation of the QAATI Model: LT_B Dataset and Erie Dataset
To validate the QAATI model, we compared the retrieved a with the measured a at the MERIS bands of 413, 443, 490, 510, 560, 620, and 665 nm using the LT-B dataset.It was found that the estimated and measured a(λ) agree with each other very well (Figure 6), where the r 2 is 0.93, and the MAPD is 15.3%.Furthermore, with the Rrs in the dataset from Lake Erie, the a(λ) was also derived by the QAATI.Figure 7 shows a comparison of measured a(λ) with those derived from Rrs.The r 2 is 0.85 with a MAPD as 39.7% (see Figure 7), which suggests very consistent a(λ) between measurements and retrievals.As with Lake Taihu, Lake Erie is also an eutrophic turbid lake, but its optical properties differ significantly from those in Lake Taihu, where the optical properties are highly influenced by rivers and rich terrestrial input from nearby factories [39].Even so, QAATI performed very well for these two lakes on two different continents, which suggests strongly applicability of QAATI in other similar turbid lakes.Furthermore, with the R rs in the dataset from Lake Erie, the a(λ) was also derived by the QAA TI .Figure 7 shows a comparison of measured a(λ) with those derived from R rs .The r 2 is 0.85 with a MAPD as 39.7% (see Figure 7), which suggests very consistent a(λ) between measurements and retrievals.As with Lake Taihu, Lake Erie is also an eutrophic turbid lake, but its optical properties differ significantly from those in Lake Taihu, where the optical properties are highly influenced by rivers and rich terrestrial input from nearby factories [39].Even so, QAA TI performed very well for these two lakes on two different continents, which suggests strongly applicability of QAA TI in other similar turbid lakes.

Validation of the QAATI Model: LT_B Dataset and Erie Dataset
To validate the QAATI model, we compared the retrieved a with the measured a at the MERIS bands of 413, 443, 490, 510, 560, 620, and 665 nm using the LT-B dataset.It was found that the estimated and measured a(λ) agree with each other very well (Figure 6), where the r 2 is 0.93, and the MAPD is 15.3%.Furthermore, with the Rrs in the dataset from Lake Erie, the a(λ) was also derived by the QAATI.Figure 7 shows a comparison of measured a(λ) with those derived from Rrs.The r 2 is 0.85 with a MAPD as 39.7% (see Figure 7), which suggests very consistent a(λ) between measurements and retrievals.As with Lake Taihu, Lake Erie is also an eutrophic turbid lake, but its optical properties differ significantly from those in Lake Taihu, where the optical properties are highly influenced by rivers and rich terrestrial input from nearby factories [39].Even so, QAATI performed very well for these two lakes on two different continents, which suggests strongly applicability of QAATI in other similar turbid lakes.

LT-A and LT-B Datasets
The new Z SD scheme was first evaluated with the LT-A and LT-B datasets.Figure 8 shows a comparison of the derived Z SD with in-situ measurements, where the derived Z SD match well with the measured Z SD (r 2 > 0.65 and MAPD < 25%), even though the Z SD values varied within a very narrow range of ~0.1 to 0.9 m.Considering that both in-situ Z SD measurements and derivation from R rs have errors/uncertainties, these results suggest a very good performance of the new Z SD scheme for application in such highly turbid lake waters.

LT-A and LT-B Datasets
The new ZSD scheme was first evaluated with the LT-A and LT-B datasets.Figure 8 shows a comparison of the derived ZSD with in-situ measurements, where the derived ZSD match well with the measured ZSD (r 2 > 0.65 and MAPD < 25%), even though the ZSD values varied within a very narrow range of ~0.1 to 0.9 m.Considering that both in-situ ZSD measurements and derivation from Rrs have errors/uncertainties, these results suggest a very good performance of the new ZSD scheme for application in such highly turbid lake waters.

Erie Dataset
The new ZSD scheme was further applied to the Lake Erie dataset (ZSD in a range of 0.4-1.7 m).As shown in Figure 9, an excellent agreement was also found between the estimated ZSD and in-situ measurements (r 2 = 0.79, MAPD = 18.6%),where the slope of the regression line is close to 1 and the intercept is close to 0. These results further demonstrated very good performance of the innovative ZSD scheme across a range from clear to turbid waters.

Erie Dataset
The new Z SD scheme was further applied to the Lake Erie dataset (Z SD in a range of 0.4-1.7 m).As shown in Figure 9, an excellent agreement was also found between the estimated Z SD and in-situ measurements (r 2 = 0.79, MAPD = 18.6%),where the slope of the regression line is close to 1 and the intercept is close to 0. These results further demonstrated very good performance of the innovative Z SD scheme across a range from clear to turbid waters.

LT-A and LT-B Datasets
The new ZSD scheme was first evaluated with the LT-A and LT-B datasets.Figure 8 shows a comparison of the derived ZSD with in-situ measurements, where the derived ZSD match well with the measured ZSD (r 2 > 0.65 and MAPD < 25%), even though the ZSD values varied within a very narrow range of ~0.1 to 0.9 m.Considering that both in-situ ZSD measurements and derivation from Rrs have errors/uncertainties, these results suggest a very good performance of the new ZSD scheme for application in such highly turbid lake waters.

Erie Dataset
The new ZSD scheme was further applied to the Lake Erie dataset (ZSD in a range of 0.4-1.7 m).As shown in Figure 9, an excellent agreement was also found between the estimated ZSD and in-situ measurements (r 2 = 0.79, MAPD = 18.6%),where the slope of the regression line is close to 1 and the intercept is close to 0. These results further demonstrated very good performance of the innovative ZSD scheme across a range from clear to turbid waters.

SAT Dataset
To validate the applicability of the new Z SD scheme to MERIS measurements, the scheme was applied to 74 MERIS images to obtain the spatial distribution of Z SD in Lake Taihu.A total of 85 matchups were obtained with the Z SD ranging from ~0.15 to 1.0 m.A comparison between the satellite-derived Z SD and in-situ measurements is shown in Figure 10.Statistically, the r 2 value from a linear regression analysis between the estimated and measured Z SD values is 0.73, and the MAPD is 37% (Figure 10).In the previous study of Doron et al. [12], several Z SD algorithms were evaluated using coincident satellite reflectance and in-situ measurements from both oceanic and coastal waters.It was found that the MAPD between MERIS-derived Z SD and in-situ measurements ranged from 60% to 140%.Comparatively, the Z SD scheme in this study showed much better performance when it is applied to MERIS images, especially for turbid waters.

SAT Dataset
To validate the applicability of the new ZSD scheme to MERIS measurements, the scheme was applied to 74 MERIS images to obtain the spatial distribution of ZSD in Lake Taihu.A total of 85 matchups were obtained with the ZSD ranging from ~0.15 to 1.0 m.A comparison between the satellite-derived ZSD and in-situ measurements is shown in Figure 10.Statistically, the r 2 value from a linear regression analysis between the estimated and measured ZSD values is 0.73, and the MAPD is 37% (Figure 10).In the previous study of Doron, et al. [12], several ZSD algorithms were evaluated using coincident satellite reflectance and in-situ measurements from both oceanic and coastal waters.It was found that the MAPD between MERIS-derived ZSD and in-situ measurements ranged from 60% to 140%.Comparatively, the ZSD scheme in this study showed much better performance when it is applied to MERIS images, especially for turbid waters.The monthly average ZSD of the entire lake ranged from ~0.24 m to 0.74 m, with a mean of 0.49 m.These values confirmed that Lake Taihu is highly turbid or opaque.Strong seasonal variability in ZSD exists across the entire lake.The monthly mean ZSD of the entire lake is relatively low from December to March, which ranges from 0.24 to 0.30 m with an average of 0.27 m.January exhibited the lowest mean ZSD throughout the year (0.24 m).In April and May, the monthly mean ZSD greatly increased to 0.41 and 0.58 m, respectively.ZSD remained relatively high from June to September, with an average of 0.70 m, and it began to decrease in October and November, with values of 0.56 and 0.41 m, respectively.Possible explanations for these patterns are discussed below.
In terms of the spatial distribution, the ZSD values of open areas were lower, while the ZSD values of bays to the east were relatively higher.Monthly changes in ZSD values from 2003 to 2011 were determined for six regions of Lake Taihu (Zhushan Bay, Meiliang Bay, Gonghu Bay, Xukou Bay, East Lake Taihu, and the Open Area).The distribution pattern of ZSD in January, February, March, and December were similar (Figures 11 and 12): the ZSD values were low in each region.From April to June, the ZSD values in Xukou Bay and East Lake Taihu significantly increased from 0.35 to 0.82 m and from 0.52 to 0.98 m, respectively, followed by Meiliang Bay (from 0.39 to 0.59 m), Gonghu Bay (from 0.44 to 0.57 m), the open area (0.25 to 0.38 m), and Zhushan Bay (0.54 to 0.60 m) (Figures 11 and  12).Moreover, the area with high ZSD in East Lake Taihu rapidly expanded during this period (Figure 11).The ZSD values maintained their high levels in each region from July to October.From October to December, the ZSD values in East Lake Taihu and Xukou Bay remarkably decreased from 1.00 m to 0.37 m and from 0.84 m to 0.28 m, respectively.Compared to East Lake Taihu and Xukou Bay, the decreases in ZSD were lower in the other four regions from October to December.

Factors Affecting ZSD in Lake Taihu
Light attenuation in water is controlled in general by four optically important substances: pure water, CDOM, tripton, and phytoplankton [40].Among these factors, tripton and phytoplankton form SPM. The concentrations of CDOM, tripton, and phytoplankton may vary in different waters or within the same water during different seasons [30,41,42].
Numerous studies have been conducted concerning the controlling factors of the temporalspatial variations in light attenuation in Lake Taihu.Previous studies revealed that high SPM concentrations dominate the light attenuation in Lake Taihu [43,44].The temporal-spatial distribution pattern of ZSD in our study is roughly the same as that from SPM in other studies [44].In Lake Taihu, the SPM dynamics were found to be controlled by three factors: the lake topography, wind-induced sediment resuspension, and the distribution of macrophytes [44].
The water depth is relatively shallower in the bays and deeper in the open areas.The water depth in each region increases in the following order: East Lake Taihu (1.76 ± 0.37 m, average ± standard error), Xukou Bay (1.82 ± 0.50 m), Zhushan Bay (1.97 ± 0.54 m), Gonghu Bay (1.98 ± 0.44 m), Meiliang Bay (2.11 ± 0.63 m), and the open area (2.63 ± 0.38 m).ZSD remained low in each month in regions with deeper depths, such as Meiliang Bay and the open area (Figures 11 and 12).In contrast, ZSD was significantly higher in shallower regions such as Xukou Bay and East Lake Taihu compared to the other regions throughout the year, except for January, February, and March (Figures 11 and  12).Unlike patterns in the ocean where lower ZSD are in the shallow coastal waters because of sediment resuspension and river input, ZSD is found to decrease with the increase of water depth in Lake Taihu.There are two mechanisms to explain this pattern.First, in the littoral zones and lake bays (shallow waters), the wind fetch was shorter than in the open area (deep waters), which decreased the intensity of wind-driven sediment resuspension and resulted in a higher ZSD than in the open area.Second, there are macrophytes distributed in the shallow area of Lake Taihu, which inhabit the resuspension of sediment, and, thereby, increase water transparency.We analyzed the correlations between the monthly average ZSD and mean water depth in each region and found negative correlations in each month, among which the correlations for November and December were most significant (p < 0.05, Supplementary Figure 1).Zhang, et al. [44] found a positive correlation between the daily wind speed and corresponding SPM concentration in 2003 and 2011 (r 2 = 0.69, p < 0.001).Shi et al. found a positive correlation between the diffuse attenuation coefficient of photosynthetically active radiation (Kd(PAR), m −1 ) and the wind speed [45] (r 2 = 0.85, n = 37, p < 0.005).Therefore, we further investigated the correlation between the monthly average ZSD and wind speed in Lake Taihu.All the valid pixels (not covered by cloud and

Factors Affecting Z SD in Lake Taihu
Light attenuation in water is controlled in general by four optically important substances: pure water, CDOM, tripton, and phytoplankton [40].Among these factors, tripton and phytoplankton form SPM. The concentrations of CDOM, tripton, and phytoplankton may vary in different waters or within the same water during different seasons [30,41,42].
Numerous studies have been conducted concerning the controlling factors of the temporal-spatial variations in light attenuation in Lake Taihu.Previous studies revealed that high SPM concentrations dominate the light attenuation in Lake Taihu [43,44].The temporal-spatial distribution pattern of Z SD in our study is roughly the same as that from SPM in other studies [44].In Lake Taihu, the SPM dynamics were found to be controlled by three factors: the lake topography, wind-induced sediment resuspension, and the distribution of macrophytes [44].
The water depth is relatively shallower in the bays and deeper in the open areas.The water depth in each region increases in the following order: East Lake Taihu ( 11 and 12).In contrast, Z SD was significantly higher in shallower regions such as Xukou Bay and East Lake Taihu compared to the other regions throughout the year, except for January, February, and March (Figures 11 and 12).Unlike patterns in the ocean where lower Z SD are in the shallow coastal waters because of sediment resuspension and river input, Z SD is found to decrease with the increase of water depth in Lake Taihu.There are two mechanisms to explain this pattern.First, in the littoral zones and lake bays (shallow waters), the wind fetch was shorter than in the open area (deep waters), which decreased the intensity of wind-driven sediment resuspension and resulted in a higher Z SD than in the open area.Second, there are macrophytes distributed in the shallow area of Lake Taihu, which inhabit the resuspension of sediment, and, thereby, increase water transparency.We analyzed the correlations between the monthly average Z SD and mean water depth in each region and found negative correlations in each month, among which the correlations for November and December were most significant (p < 0.05, Supplementary Figure S1).
Zhang et al. [44] found a positive correlation between the daily wind speed and corresponding SPM concentration in 2003 and 2011 (r 2 = 0.69, p < 0.001).Shi et al. found a positive correlation between the diffuse attenuation coefficient of photosynthetically active radiation (K d (PAR), m −1 ) and the wind speed [45] (r 2 = 0.85, n = 37, p < 0.005).Therefore, we further investigated the correlation between the monthly average Z SD and wind speed in Lake Taihu.All the valid pixels (not covered by cloud and surface scum) over the entire lake within each month were spatially averaged to generate the corresponding Z SD value.The correlation between the monthly average Z SD and wind speed was statistically significant (p < 0.05), which suggests that the wind speed is a major factor for the phenological variations in Z SD in Lake Taihu (Figure 13), as shown by Shang et al. [46] where wind is a major factor for lower clarity in Bohai Sea.
Based on field and remote-sensing observations, macrophytes are abundantly distributed in the littoral regions of East Lake Taihu, Xukou Bay, and Gonghu Bay [43,45].The growth habits of aquatic vegetation are largely responsible for the monthly variations in Z SD in these regions.From December to the following March, the death and degradation of aquatic vegetation produce lower Z SD values (Figure 11).The growth of aquatic vegetation can improve the underwater light climate by allelopathy and suppress sediment resuspension [44].Therefore, the area with high Z SD increased in April and May as aquatic vegetation grew.From June to October, the aquatic vegetation reached its peak density and extent [43].Correspondingly, the highest Z SD values of the year were found in the regions that were dominated by aquatic vegetation during this period (Figure 11).Based on field and remote-sensing observations, macrophytes are abundantly distributed in the littoral regions of East Lake Taihu, Xukou Bay, and Gonghu Bay [43,45].The growth habits of aquatic vegetation are largely responsible for the monthly variations in ZSD in these regions.From December to the following March, the death and degradation of aquatic vegetation produce lower ZSD values (Figure 11).The growth of aquatic vegetation can improve the underwater light climate by allelopathy and suppress sediment resuspension [44].Therefore, the area with high ZSD increased in April and May as aquatic vegetation grew.From June to October, the aquatic vegetation reached its peak density and extent [43].Correspondingly, the highest ZSD values of the year were found in the regions that were dominated by aquatic vegetation during this period (Figure 11).

Algorithm Uncertainties
Even if we assume that the in-situ Rrs measurements were completely error-free, multiple sources of uncertainty still affect the accuracy of the final ZSD value.In the new ZSD scheme, the uncertainty of derived ZSD mainly come from a(λ0) and η.According to a systematic study of the sources of uncertainties in IOP products from the QAA, the a(λ0) modeling errors for waters with high absorption coefficients have a larger effect than the η modeling errors on the a(λ) estimation [47].In our study, we attempted to minimize the uncertainties that are introduced from these two parameters by refining the corresponding steps in QAA.In this step (Step 4), the values of η were derived from a combination of Rrs and a, where a was calculated from ap and aCDOM.Thus, the uncertainties of measured ap and aCDOM as well as those from measured Rrs could propagate to the estimated η.The ap in this study was quantified by utilizing the standard transmittance/reflectance (T-R) method by using the integrating sphere, and the uncertainty for this measurement procedure was found to be ~10% [48].The detection limit of Ultraviolet (UV) spectrophotometers for the measurement of aCDOM was ~0.0001 absorbance units (equal to ~0.005 m −1 in absorption coefficients) [30].The uncertainty of this measurement could be relatively low, considering the high aCDOM in eutrophic inland waters (generally greater than ~0.2 m −1 in violet and blue bands).
The power coefficient of bbp, η, is sensitive to the shape of the particle size distribution (PSD), according to theoretical calculations and field measurements [49][50][51].It was found that small particles generally have enhanced backscattering toward shorter wavelengths, whereas larger particles tend to have a flatter backscattering spectrum.In oceanic water, the water optical properties are determined primarily by phytoplankton and related CDOM and degradation products.Previous

Algorithm Uncertainties
Even if we assume that the in-situ R rs measurements were completely error-free, multiple sources of uncertainty still affect the accuracy of the final Z SD value.In the new Z SD scheme, the uncertainty of derived Z SD mainly come from a(λ 0 ) and η.According to a systematic study of the sources of uncertainties in IOP products from the QAA, the a(λ 0 ) modeling errors for waters with high absorption coefficients have a larger effect than the η modeling errors on the a(λ) estimation [47].In our study, we attempted to minimize the uncertainties that are introduced from these two parameters by refining the corresponding steps in QAA.In this step (Step 4), the values of η were derived from a combination of R rs and a, where a was calculated from a p and a CDOM .Thus, the uncertainties of measured a p and a CDOM as well as those from measured R rs could propagate to the estimated η.The a p in this study was quantified by utilizing the standard transmittance/reflectance (T-R) method by using the integrating sphere, and the uncertainty for this measurement procedure was found to be ~10% [48].The detection limit of Ultraviolet (UV) spectrophotometers for the measurement of a CDOM was ~0.0001 absorbance units (equal to ~0.005 m −1 in absorption coefficients) [30].The uncertainty of this measurement could be relatively low, considering the high a CDOM in eutrophic inland waters (generally greater than ~0.2 m −1 in violet and blue bands).
The power coefficient of b bp , η, is sensitive to the shape of the particle size distribution (PSD), according to theoretical calculations and field measurements [49][50][51].It was found that small particles generally have enhanced backscattering toward shorter wavelengths, whereas larger particles tend to have a flatter backscattering spectrum.In oceanic water, the water optical properties are determined primarily by phytoplankton and related CDOM and degradation products.Previous studies showed a general decrease in η from oligotrophic to eutrophic regimes in oceanic waters, with higher values of η (>2) in oligotrophic waters where smaller phytoplankton dominate and lower values of η (<1) in more eutrophic waters where larger phytoplankton dominate [52].In oceanic waters, small-sized phytoplankton (e.g., cyanobacteria, prochlorophytes, and chlorophytes) are incapable of growing beyond a particular chlorophyll concentration, with an upper limit imposed possibly from a combination of a bottom-up (nutrient control) and a top-down (grazing) process [53].Therefore, an increasing trend of phytoplankton cell size with Chla can be observed in oceanic water and could explain the lower η in eutrophic regimes than in oligotrophic regions.However, in the eutrophic freshwaters, cyanobacteria are the most wide-spread algal taxa [54].Typically, this species is classified as pico-phytoplankton (0.2-2 µm) and is usually characterized by a cell size of <2 µm [55].In our study, the derived η value (of a highly eutrophic lake) ranged from −1.6 to 2.8 with a mean value of 1.7, because a large number of samples (~34%) used in this study were collected in regions with highly abundant phytoplankton (Chla > 30 µg/L), where the dominant species was cyanobacteria [56].Thus, these results could explain the higher η values for some data points in this study.Our results are consistent with reports in the literature regarding particle backscattering in turbid eutrophic waters with high Chla concentrations.Ma et al. [57] showed that the η values calculated from the Hydroscat-6 (HOBI Labs, Watsonville, California) measured backscattering coefficients in Lake Taihu and have a mean value of 2.19 (SD = 0.39).A study in the coastal region of Oslo Fjord also showed that algal bloom (when Chla is greater than 2 µg/L) always associates with steep spectral curves [58].This finding was further supported by a study in tropical water reservoirs in Queensland, Australia [59].In addition, in a modelling study of five different cultured phytoplankton species, decreasing backscattering toward longer wavelengths was most significant in the case of cyanobacteria [60], and the η showed an increasing trend with the increase of Chla [61].
For the Z SD maps derived from satellite measurements, another uncertainty is potentially associated with floating algae on the surface of the water.The floating algae is prevalent during the summer months in Lake Taihu when the cyanobacterial blooms occur [30,62].These masses of floating algae inhibit the reflectance signal of the water column from satellite reception, which results in questionable retrieval of Z SD in these areas during periods with heavy algae.However, because cyanobacteria (the dominant bloom species in Lake Taihu) contain phycocyanin (an accessory pigment with a diagnostic absorption peak at ~620 nm), cyanobacterial blooms can be detected with an algorithm that was developed for MERIS called PCI [35].By applying this simple algorithm, we detected and removed pixels that were greatly influenced by floating algae.This step can avoid error that is introduced by the influence of thick floating algae.However, the Z SD values in bloom areas are typically low (often <0.2 m in observations).Thus, removing these pixels would cause Z SD to be overestimated in such regions during long-term analysis.The influence of cyanobacterial blooms on the retrieval of water color parameters has not been considered in many previous studies of Lake Taihu, such as those that investigated SPM and K d (PAR) [43,44,63].Therefore, this step is expected to be helpful for remote sensing in other very eutrophic waters.

Implications of the Proposed Model
Many studies of the optical properties of coastal and estuarine regions have revealed that significant correlations exist between SPM and Z SD [64,65].For example, Lake Frisian in the Netherlands is a shallow water lake with a water depth of 1 to 2 m, and its Z SD is mainly influenced by SPM [66].Shi et al. [43] evaluated the correlation between Z SD and SPM in Lake Taihu based on 2432 in-situ observations and found a relationship of Z SD = 3.07•SPM −0.60 (r 2 = 0.73, p < 0.001).The significant correlation between SPM and Z SD demonstrated that the SPM concentration can be obtained from a semi-analytical model in waters where SPM dominates if the quantitative relationship between SPM and Z SD is determined.
The spatial distributions of Z SD in Lake Taihu were different from that of SPM presented in some previous studies, especially during the summer months.In those studies, the SPM exhibited low values in both phytoplankton-dominated turbid regions (such as Meiliang Bay and Zhushan Bay) and in clear regions that were dominated by submerged macrophytes (such as Xukou Bay and East Lake Taihu) when algal blooms occurred during the summer [43,44,67].Note that there could be misinterpretations of SPM in phytoplankton-dominated regions, where the absorption and scattering signals in the water column are inhibited when the waters are covered by massive floating algal.Therefore, SPM exhibited similar distribution trends in these two different regions.In our study, the regions covered by serious floating algae were removed by applying a PCI index [35].Thus, we avoided the influence of floating blooms when retrieving Z SD .

Figure 1 .
Figure 1.Station maps of datasets Lake Taihu (a) and Lake Erie (b).

Figure 1 .
Figure 1.Station maps of datasets Lake Taihu (a) and Lake Erie (b).

Figure 4 .
Figure 4. Examples of (a) a(560) and corresponding (b) Rrs spectrum at different locations.

Figure 4 .
Figure 4. Examples of (a) a(560) and corresponding (b) R rs spectrum at different locations.

Figure 6 .
Figure 6.Comparing the modeled and the measured a(λ) at MERIS bands using the validation (LT-B) dataset.

Figure 7 .
Figure 7. Comparing the estimated and the measured a(λ) at MERIS bands using the Lake Erie dataset.

Figure 6 .
Figure 6.Comparing the modeled and the measured a(λ) at MERIS bands using the validation (LT-B) dataset.

Figure 6 .
Figure 6.Comparing the modeled and the measured a(λ) at MERIS bands using the validation (LT-B) dataset.

Figure 7 .
Figure 7. Comparing the estimated and the measured a(λ) at MERIS bands using the Lake Erie dataset.

Figure 7 .
Figure 7. Comparing the estimated and the measured a(λ) at MERIS bands using the Lake Erie dataset.

Figure 8 .
Figure 8.Comparison of estimated ZSD with in-situ measurements in LT-A and LT-B datasets.

Figure 9 .
Figure 9.Comparison of estimated and in-situ ZSD with the Erie dataset.

Figure 8 .
Figure 8.Comparison of estimated Z SD with in-situ measurements in LT-A and LT-B datasets.

Figure 8 .
Figure 8.Comparison of estimated ZSD with in-situ measurements in LT-A and LT-B datasets.

Figure 9 .
Figure 9.Comparison of estimated and in-situ ZSD with the Erie dataset.Figure 9. Comparison of estimated and in-situ Z SD with the Erie dataset.

Figure 9 .
Figure 9.Comparison of estimated and in-situ ZSD with the Erie dataset.Figure 9. Comparison of estimated and in-situ Z SD with the Erie dataset.

Figure 10 .
Figure 10.Comparison between in-situ and MERIS-derived ZSD, where the color denotes the wavelength of minimum Kd.

3. 3 .
Temporal Variation of Water Clarity from MERIS Measurements: Lake Taihu as an Example 3.3.1.Phenological Distribution of ZSD from the MERIS Images The monthly climatology of ZSD was produced from the MERIS-estimated ZSD data from 2003 to 2011, as shown in Figures 11 and 12 .

Figure 10 .
Figure 10.Comparison between in-situ and MERIS-derived Z SD , where the color denotes the wavelength of minimum K d .

Figure 11 .
Figure 11.Phenological distribution of ZSD in Lake Taihu, calculated from all valid pixels of the 74 selected MERIS images using QAATI.

Figure 11 .
Figure 11.Phenological distribution of Z SD in Lake Taihu, calculated from all valid pixels of the 74 selected MERIS images using QAA TI .

Figure 12 .
Figure 12.Phenological variation of ZSD in each region of Lake Taihu.

Figure 12 .
Figure 12.Phenological variation of Z SD in each region of Lake Taihu.

Figure 13 .
Figure 13.Linear relationship between monthly averaged wind speed and MERIS-derived ZSD of the entire lake, from 2003 to 2011.

Figure 13 .
Figure 13.Linear relationship between monthly averaged wind speed and MERIS-derived Z SD of the entire lake, from 2003 to 2011.

Table 1 .
List of datasets used in this study.

Table 2 .
Sampling date, number of samples, and distribution of concurrent sampling sites with MERIS for six cruises in the SAT dataset.

Table 2 .
Sampling date, number of samples, and distribution of concurrent sampling sites with MERIS for six cruises in the SAT dataset.