Seasonal Variability of Di ﬀ use Attenuation Coe ﬃ cient in the Pearl River Estuary from Long-Term Remote Sensing Imagery

: We evaluated six empirical and semianalytical models of the di ﬀ use attenuation coe ﬃ cient at 490 nm ( K d (490)) using an in situ dataset collected in the Pearl River estuary (PRE). A combined model with the most accurate performance (correlation coe ﬃ cient, R 2 = 0.92) was selected and applied for long-term estimation from 2003 to 2017. Physical and biological processes in the PRE over the 14-year period were investigated by applying satellite observations (MODIS / Aqua data) and season-reliant empirical orthogonal function analysis (S-EOF). In winter, the average K d (490) was signiﬁcantly higher than in the other three seasons. A slight increasing trend was observed in spring and summer, whereas a decreasing trend was observed in winter. In summer, a tongue with a relatively high K d (490) was found in southeastern Lingdingyang Bay. In Eastern Guangdong province (GDP), the relatively higher K d (490) value was found in autumn and winter. Based on the second mode of S-EOF, we found that the higher values in the eastern GDP extended westward and formed a distinguishable tongue in winter. The grey relational analysis revealed that chlorophyll-a concentration ( C chla ) and total suspended sediment concentration ( C tsm ) were two dominant contributors determining the magnitude of K d (490) values. The C tsm -dominated waters were generally located in coastal and estuarine turbid waters; the C chla -dominated waters were observed in open clear ocean. The distribution of constituents-dominated area was di ﬀ erent in the four seasons, which was a ﬀ ected by physical forces, including wind ﬁeld, river runo ﬀ , and sea surface temperature. of the Several approaches, including and semianalytical models, were applied to retrieve the K d (490) in PRE water. The results showed that Wang’s model was more accurate and is most suitable for PRE water, which uses a combination of di ﬀ erent algorithms for clear and turbid waters. Hence, Wang’s model was selected for deriving K d (490) products from long-term MODIS / Aqua imagery. Aqua the variability and


Introduction
The light diffuse attenuation coefficient (K d (λ)) in aquatic systems is defined by the exponential decrease in the irradiance with depth [1,2]. K d (λ) is an ecologically important water property that provides an estimate of the availability of light to underwater communities, which influences ecological processes and biogeochemical cycles in natural waters [3,4]. The estimation of K d (λ) is also critical for understanding physical processes such as sediment resuspension and heat transfer in the upper layer of the ocean [5][6][7].
The in situ K d (λ) is traditionally measured by the ocean color scientific community at 490 nm, K d (490), following the primary studies in the 1970s [8]. Traditional field measurement of K d (λ) is costly and time consuming, but recent advances in satellite sensors have provided synoptic and frequent A cruise was conducted on 5 June 2012 to collect water samples and the water spectrum. Positions for all sampling stations are plotted in Figure 1. The field spectral measurements were composed of two parts: the above-water remote sensing reflectance (R rs ) and the downwelling irradiance within the water column. To obtain the background water column conditions, water samples from the 15 sampling stations were used for measurement of chlorophyll-a (C chla ), total suspended sediment (C tsm ), absorption coefficient for phytoplankton (a p (λ)), and colored dissolved organic matter (CDOM, a g (λ)) ( Table 1).
The above-water R rs was measured using a spectroradiometer (USB4000, Ocean Optics, Inc., Dunedin, FL, USA) following the National Aeronautics and Space Administration (NASA) ocean-optics standard protocol [21]. The upward radiance (L u ), downward sky radiance (L sky ), and radiance from standard spectra on a reference plaque (L pla ) were measured, and R rs was calculated using the following equation: R rs (λ) = ρ pla (λ) L u (λ) − ρ f (λ)L sky (λ) / πL pla (λ) (1) where λ is the wavelength, ρ pla is the reflectance of the plaque provided by the manufacturer (Ocean Optics, Inc., Dunedin, FL, USA), ρ f is the water surface Fresnel reflectance, where a value of 0.028 was taken for wind speeds of less than 5 m·s −1 .
Remote Sens. 2020, 12, 2269 3 of 18 To evaluate the MODIS-based K d (490) retrieval models, in situ R rs was aggregated to simulate MODIS/Aqua R rs according to the following equation [22][23][24]: R rs (Bi) = λ n λ m RSR(λ) * R rs_meas (λ)dλ λ n λ m RSR(λ)dλ (2) where R rs (Bi) denotes the simulated R rs for the ith band of MODIS/Aqua, with integration from λ m to λ n ; R rs_meas (λ) denotes the field-measured R rs (λ); and RSR(λ) denotes the MODIS/Aqua spectral response function. Remote Sens. 2020, 12 The above-water Rrs was measured using a spectroradiometer (USB4000, Ocean Optics, Inc., Dunedin, FL, USA) following the National Aeronautics and Space Administration (NASA) oceanoptics standard protocol [21]. The upward radiance (Lu), downward sky radiance (Lsky), and radiance from standard spectra on a reference plaque (Lpla) were measured, and Rrs was calculated using the following equation: where λ is the wavelength, ρpla is the reflectance of the plaque provided by the manufacturer (Ocean Optics, Inc., Dunedin, FL, USA), ρf is the water surface Fresnel reflectance, where a value of 0.028 was taken for wind speeds of less than 5 m·s −1 .
To evaluate the MODIS-based Kd(490) retrieval models, in situ Rrs was aggregated to simulate MODIS/Aqua Rrs according to the following equation [22][23][24]:   Downwelling irradiance within the water column was measured with a TriOS-RAMES hyperspectral spectroradiometer (TriOS GmbH, Oldenburg, Germany). The spectroradiometer recorded irradiance signal in the range of 320 to 950 nm with a wavelength resolution of 3.3 nm. The TriOS-RAMES instrument was slowly hand-lowered at a stable speed from the surface to a water depth of about 5 m and set to a sampling rate of one sample every five seconds. Meanwhile, a pressure sensor recorded the corresponding depth of water. By releasing the TriOS-RAMES instrument (TriOS GmbH, Oldenburg, Germany) into water twice, two profiles of the downwelling irradiance were collected. The two profiles were averaged to minimize the effect of near-surface wave focusing. The natural logarithm of the measured irradiance was plotted against depth, and an estimate of K d (λ) was acquired from the resulting slope [25]:

of 18
where λ is the wavelength, E d (z) is the downwelling irradiance at depth z, and ∆z is the infinitesimal thickness at depth z.

MODIS/Aqua Imagery
The Level-1B MODIS/Aqua ocean color dataset and the geolocation dataset from 2003 to 2017 were obtained from the Level-1 and Atmosphere Archive and Distribution System (LAADS) Distributed Active Archive Center (DAAC). Imagery was preprocessed using the SeaWiFS data analysis system (SeaDAS, version 7.5.1). The Management Unit of the North Seas Mathematical Models (MUMM)-based atmospheric correction [26] and an iterative f /Q Bidirectional Reflectance Distribution Function (BRDF) correction [27][28][29][30] were used to acquire accurate R rs values. Flags were used to mask contamination from land, clouds, sun glint, and other potential disturbances to the radiance signal.

Ancillary Data
The wind field dataset was obtained from the National Centers for Environmental Prediction (NCEP) Climate Forecast System Version 2 (CFSv2). The model is fully coupled, representing the Earth's atmosphere, oceans, land, and sea ice [31]. The mixed layer depth (MLD), defined as the depth where the density is equal to the sea surface density plus an increase in density equivalent to 0.8 • C, was acquired from the global ocean Argo gridded dataset (BOA_Argo, provided by the China Argo Real-time Data Center, ftp://data.argo.org.cn/pub/ARGO/BOA_Argo/) [32]. The monthly river runoff was acquired from the Chinese River Sediment Bulletin. The Level-3 MODIS/Aqua sea surface temperature (SST) dataset was obtained from the Ocean Color Website (https://oceancolor.gsfc.nasa.gov/l3/), a website that provides the derived geophysical variables that have been aggregated/projected onto a well-defined spatial grid during a well-defined time period.

Models for K d (490) Retrieval
At present, the standard methods for K d (490) estimation are roughly classified into three types: (1) empirical relationship between K d (490) and apparent optical properties (AOP), including water-leaving radiance or reflectance [11,33,34]; (2) empirical relationship between K d (490) and chlorophyll-a based on regression analyses [35]; and (3) semianalytical approaches based on radiative transfer models [1,36]. These three types of models, six models in total (Table 2), were evaluated in the PRE waters using the in situ dataset. analysis are as follows: Firstly, the time series of seasonal K d (490) anomaly was calculated. Secondly, the EOF was analyzed based on the matrix composed of the four seasons. Finally, each S-EOF mode containing four spatial modes, which represented the spatial patterns of K d (490) in the four seasons, and a corresponding principal component time series were obtained.
GRA is an important part of grey system theory, which is used to determine the relational degree among factors according to the similarities in their geometry [38]. The GRA was applied here to identify the dominant water constituents (total suspended matter, phytoplankton, and dissolved matter) affecting the spatial distribution and temporal variation of K d (490). In GRA, the reference series of K d (490) and comparison sequences (water constituents, including C tsm , C chla , and a dg (443)) were constructed in advance to calculate the grey relational grade (GRG), which is a measure of similarity between the reference sequence and comparison sequences. Details about the calculation of GRG were described by Liu and Lin (2005)

Performance Assessment
To compare the performance of different K d (490) retrieval models, several statistical parameters were used: the determination coefficient (R 2 ), root mean square error (RMSE), mean absolute difference (MAD), and mean absolute percentage difference (MAPD), which are calculated as: where x m and x p denote the measured and predicted samples, respectively; x m denotes the mean value of the measured samples; and N is the number of samples.

Model Performance
We evaluated the six different models with MODIS/Aqua spectral bands or C chla . The evaluation was based on the comparison of the model-derived K d (490) with in situ measured K d (490) collected from the PRE on 5 June 2012. Figure 2 shows scatterplots between the in situ measured and different models' K d (490) retrievals, and Table 3 lists the statistical parameters. The results provided by both Mueller's and Morel's models constantly underestimated the K d (490) compared with the in situ dataset for the PRE, with RMSEs higher than 1.1 m −1 , MADs close to 1.0 m −1 , and MAPDs up to 70%. The Morel (empirical model with C chla ) and Mueller (empirical model with water-leaving radiance) models not only underestimated the in situ values of the PRE, they had little to no sensitivity along a broad gradient of in situ values.
By comparison, the other four models appeared to be more effective when applied in the PRE waters. These four models performed well with R 2 values higher than 0.9, RMSEs ranging from 0.31 to 0.70 m −1 , MADs ranging from 0.27 to 0.54 m −1 , and MAPDs ranging between 25.51% and 37.10%. We found that Wang's model, combining Lee's algorithm for turbid waters and Mueller's algorithm for clear waters, was a better choice for K d (490) retrieval in the PRE waters. Comparison of Wang's model to the other models showed that Wang's model had considerably lower RMSE and MAD values and outperformed the other models, especially at relatively higher K d (490) levels. and different models' Kd(490) retrievals, and Table 3 lists the statistical parameters. The results provided by both Mueller's and Morel's models constantly underestimated the Kd(490) compared with the in situ dataset for the PRE, with RMSEs higher than 1.1 m −1 , MADs close to 1.0 m −1 , and MAPDs up to 70%. The Morel (empirical model with Cchla) and Mueller (empirical model with water-leaving radiance) models not only underestimated the in situ values of the PRE, they had little to no sensitivity along a broad gradient of in situ values. By comparison, the other four models appeared to be more effective when applied in the PRE waters. These four models performed well with R 2 values higher than 0.9, RMSEs ranging from 0.31 to 0.70 m −1 , MADs ranging from 0.27 to 0.54 m −1 , and MAPDs ranging between 25.51% and 37.10%. We found that Wang's model, combining Lee's algorithm for turbid waters and Mueller's algorithm for clear waters, was a better choice for Kd(490) retrieval in the PRE waters. Comparison of Wang's

Spatial Distribution and Temporal Variation
Given its superior performance of the six considered models, the long-term MODIS K d (490) products were derived based on Wang's model. Significant seasonal variation was identified over the entire study area from 2003 to 2017 ( Figure 3). The mean values for the entire study area were 0.13 m −1 in spring, 0.12 m −1 in summer, 0.14 m −1 in autumn, and 0.21 m −1 in winter. In the coastal area, the relatively high K d (490) was observed in Lingdingyang Bay (LB) and western Guangdong Province (GDP), where the highest value exceeded 4.0 m −1 . In summer, the river plume extends from LB southeastward into the coastal region, resulting in a wider distribution of high-value K d (490). The plume waters formed a tongue along the eastern GDP (located 113-115 • E, 22-22.5 • N) in some specific years, though this feature was not so remarkable in the seasonal climatological imagery due to long-term average smoothing. The influence of the terrestrial input of nutrients from the Pearl River is highest in summer. In addition, a southeast wind prevails and rainfall mainly occurs in summer. The winds blow from PRE to the middle shelf. The winds control the spatial pattern of K d (490) distribution in the PRE. In the eastern GDP, relatively higher K d (490) values were found in autumn and winter, whereas lower values were observed in spring and summer. This phenomenon is closely correlated with the coastal upwelling along the eastern GDP coast [41,42]. Chen et al. (1982) [43] reported that a radiating current could generate an upwelling in winter near the Jieshi Bay in the eastern GDP. In open ocean areas, the average K d (490) values in winter and spring were higher than that in summer and autumn. The prevailing northeasterly monsoon is stronger in the northern South China Sea in winter, so the MLD was deeper. The mixing effects are relatively stronger in winter.

GRG of Water Constituents
The optical properties were determined using the absorption or backscattering of different water constituents. Here, three types of water constituents were considered; the monthly average C tsm , C chla , and a dg (443) were obtained based on the same atmospheric-corrected MODIS/Aqua R rs dataset that was used for K d (490) retrieval. A band ratio algorithm was adopted for C tsm retrieval [17]. The OC3M algorithm was used for C chla retrieval [44]. The generalized inherent optical property (GIOP) model was applied for a dg (443) retrieval [45,46]. The GRGs, which can be used to measure the relationships between the K d (490) and the three water constituents, were calculated pixel by pixel for the four seasons ( Figure 5).

GRG of Water Constituents
The optical properties were determined using the absorption or backscattering of different water constituents. Here, three types of water constituents were considered; the monthly average Ctsm, Cchla, and adg(443) were obtained based on the same atmospheric-corrected MODIS/Aqua Rrs dataset that was used for Kd(490) retrieval. A band ratio algorithm was adopted for Ctsm retrieval [17]. The OC3M algorithm was used for Cchla retrieval [44]. The generalized inherent optical property (GIOP) model was applied for adg(443) retrieval [45,46]. The GRGs, which can be used to measure the relationships between the Kd(490) and the three water constituents, were calculated pixel by pixel for the four seasons ( Figure 5). Spatially, GRGs between Kd(490) and Cchla or adg(443) were higher in the clear open ocean region than in the coastal region, but the values contrasted between Kd(490) and Ctsm. The GRGs gradually decreased from nearshore to offshore, similar to the distribution of Ctsm. Seasonally, the GRG was higher in summer and autumn than in spring and winter between Kd(490) and Cchla, with most of the pixels' values being above 0.8. Similar phenomena were observed in the GRGs between Kd(490) and adg(443), although the average value was lower than between Kd(490) and Cchla. Spatially, GRGs between K d (490) and C chla or a dg (443) were higher in the clear open ocean region than in the coastal region, but the values contrasted between K d (490) and C tsm . The GRGs gradually decreased from nearshore to offshore, similar to the distribution of C tsm . Seasonally, the GRG was higher in summer and autumn than in spring and winter between K d (490) and C chla , with most of the pixels' values being above 0.8. Similar phenomena were observed in the GRGs between K d (490) and a dg (443), although the average value was lower than between K d (490) and C chla .

S-EOF Analysis
An S-EOF analysis was performed after subtracting the long-term monthly climatological average K d (490). The first two modes and the corresponding principal components (PC) were separated, which accounted for approximately 81.16% of the total variance (Table 4).

S-EOF Analysis
An S-EOF analysis was performed after subtracting the long-term monthly climatological average Kd(490). The first two modes and the corresponding principal components (PC) were separated, which accounted for approximately 81.16% of the total variance (Table 4).    Figure 7 shows the spatial distribution of the first mode of S-EOF, which explained approximately 56.7% of the total variance. Relatively high values were observed in LB and the western coast of GDP in the four seasons, and the average values for spring were significantly lower than in the other seasons. In summer, the high value area tended to expand to the southeastern LB. From autumn to winter, we observed a high value area along the east coast of the GDP, whereas this high value area disappeared in spring and summer.  Figure 7 shows the spatial distribution of the first mode of S-EOF, which explained approximately 56.7% of the total variance. Relatively high values were observed in LB and the western coast of GDP in the four seasons, and the average values for spring were significantly lower than in the other seasons. In summer, the high value area tended to expand to the southeastern LB. From autumn to winter, we observed a high value area along the east coast of the GDP, whereas this high value area disappeared in spring and summer.
The second mode of S-EOF explained 24.5% of the total variance. A relatively higher value area was observed in LB during the four seasons ( Figure 8). In summer, the high value area extended eastward and formed a distinguishable tongue. Compared with spring and summer, higher values were observed along the whole coastal zone of the PRE. Based on the second mode of S-EOF, we found that the higher values in the eastern GDP extended westward and formed a distinguishable tongue in winter. Remote Sens. 2020, 12, x FOR PEER REVIEW 12 of 20 The second mode of S-EOF explained 24.5% of the total variance. A relatively higher value area was observed in LB during the four seasons ( Figure 8). In summer, the high value area extended eastward and formed a distinguishable tongue. Compared with spring and summer, higher values were observed along the whole coastal zone of the PRE. Based on the second mode of S-EOF, we found that the higher values in the eastern GDP extended westward and formed a distinguishable tongue in winter.

Evaluation of Kd(490) Models
Both Mueller and Morel's models were unsuitable for the PRE waters because these two models were established for clear waters and only use the spectral information from the blue and green bands. For clear waters where the downwelling attenuation is mainly determined by phytoplankton, the blue-green band ratio is sensitive to the variability of the Cchla, resulting in a high accuracy for Kd(490) retrieval. However, for turbid waters where the optical properties are

Evaluation of K d (490) Models
Both Mueller and Morel's models were unsuitable for the PRE waters because these two models were established for clear waters and only use the spectral information from the blue and green bands. For clear waters where the downwelling attenuation is mainly determined by phytoplankton, the blue-green band ratio is sensitive to the variability of the C chla , resulting in a high accuracy for K d (490) retrieval. However, for turbid waters where the optical properties are more complex, the blue-green band ratio demonstrates a lower sensitivity to the variability in K d (490). The strong absorption of phytoplankton and CDOM could lead to relatively smaller R rs in the blue and green bands [47,48]. In the PRE, the water constituents are from river inputs and coastal erosion. The C chla , CDOM, and C tsm are very high, which may result in both Mueller and Morel's models being inapplicable in the PRE. Tiwari's model uses the reflectance ratio at 490 and 670 nm, R rs (490)/R rs (670), to derive K d (490). Zhang's model is composed of two independent algorithms: one based on the ratio R rs (490)/R rs (555) for clear waters and another based on the ratio R rs (490)/R rs (665) for turbid waters. When tested with the independent PRE dataset, the predictions of these two models were statistically better compared to both Mueller's and Morel's models. However, the two models also showed a pronounced underestimation for higher K d (490) values (>1.0 m −1 ), which might be due to the strong backscattering of suspended sediments in the more turbid waters in the PRE. Lee's model produced a suitable estimation, which uses a relationship relating the backscattering coefficient at 490 nm to the irradiance reflectance just beneath the surface within the red band. The performance of Wang's model was the same as Lee's, which was attributed to the same approach used in both models for highly turbid waters. In Wang's model, the retrieval method switched to Mueller's in clear waters, and the bridging of the two types of models is based on a certain weighting function (W). The spectral information within the red band cannot be ignored when retrieving K d (490) for turbid waters. However the values of R rs (670)/R rs (490) tended to be very low and, therefore, values of K d (490) for clear waters were inaccurate. Therefore, Wang's model, which uses a combination of different algorithms for clear and turbid waters, is a better choice for K d (490) retrieval in the PRE waters.

Dominant Contributor to K d (490) of Water Constituents
Attenuation of light in water depends on concentrations of particulate matter and dissolved matter, which can be expressed by C tsm , C chla , and the absorption coefficient of CDOM [7,49]. The contribution of these constituents varies for different types of water and within the same water body in different seasons [50][51][52]. Since the calculated GRGs between K d (490) and a dg (443) were significantly lower than the other two water constituents, only the GRGs of C tsm and C chla were considered. Figure 9 depicts the subtraction of both GRGs. Positive values indicate the GRGs of C tsm were higher than those of C chla , which means that C tsm played a dominant role in K d (490) variability. In contrast, negative values indicate that the C chla had a greater influence. The C tsm -dominated were waters generally located in coastal and estuarine turbid areas, whereas the C chla -dominated waters were observed in open clear ocean. Notably, waters dominated by a dg (443) were rare. The strong absorption of CDOM in the blue bands influenced the variability of K d (490), particularly in waters with high CDOM concentrations. The major sources of CDOM in the PRE were the river water and the human and industrial sewage [53,54]. However, in coastal or estuarine areas with highly turbid waters, C tsm can reach over 100 g·m −3 . During the survey conducted on 5 June 2012, the range of measured a g (443) was 0.12 to 0.58 and the range of measured a p (443) was 0.31 to 1.61. The latter was approximately three orders higher than the former, indicating that the influence of total suspended sediments on K d (490) was far greater than that of CDOM.
with high CDOM concentrations. The major sources of CDOM in the PRE were the river water and the human and industrial sewage [53,54]. However, in coastal or estuarine areas with highly turbid waters, Ctsm can reach over 100 g˖m −3 . During the survey conducted on 5 June 2012, the range of measured ag(443) was 0.12 to 0.58 and the range of measured ap(443) was 0.31 to 1.61. The latter was approximately three orders higher than the former, indicating that the influence of total suspended sediments on Kd(490) was far greater than that of CDOM. The distribution of dominant constituents showed some seasonality. In spring and summer, the Ctsm-dominated waters were mainly distributed in LB and the western GDP. The Ctsm-dominated area was confined close to the nearshore areas in the eastern GDP, indicating the impact of Cchla can extend from offshore to nearshore regions. Compared with other seasons, the most significant feature in summer was the southward extension of Ctsm from LB to the open ocean, which can The distribution of dominant constituents showed some seasonality. In spring and summer, the C tsm -dominated waters were mainly distributed in LB and the western GDP. The C tsm -dominated area was confined close to the nearshore areas in the eastern GDP, indicating the impact of C chla can extend from offshore to nearshore regions. Compared with other seasons, the most significant feature in summer was the southward extension of C tsm from LB to the open ocean, which can probably be attributed to the increase in river runoff. In autumn and winter, the C tsm -dominated area was wider in the eastern GDP than in spring and summer. The underlying reason for the change in area still requires future research. Currently, the change in area in the eastern GDP during autumn and winter might be indirectly caused by the decrease in C chla rather than the variability of C tsm . In autumn and winter, the entire eastern GDP is influenced by monsoons. The northeasterly wind-induced downwelling appears to decrease the amount of resuspension, resulting in the slight decrease in surface C tsm , which seems to contradict the expansion of the C tsm -dominated area. However, the downwelling also inhibits the growth of phytoplankton. The decrease in surface C chla may prevent it from becoming the primary factor affecting the variability of K d (490).  Figure 9 shows that the spatial variations can be attributed to the changes in C chla and C tsm .

Influence of Physical Factors on K d (490) Variability
To understand the mechanism through which the seasonal K d (490) varies, correlation analysis was performed between several types of physical factors, including wind field, river runoff, MLD, SST, and seasonal average K d (490). The results showed that the average K d (490) was highly correlated with the wind speed (u-component) in summer, with an R 2 of about 0.69. During winter, we found a significant negative correlation between K d (490) and SST, with an R 2 of -0.66 ( Figure 10). Seasonal anomalies were also obtained by subtracting the seasonal climatological average. In 2007 and 2015, when the wind speed anomaly (u-component) reached its peak, a distinguishing tongue of K d (490) anomaly was observed near the southeastern LB ( Figure 11). Inside this tongue region, the K d (490) anomaly in the west was higher than in the east, indicating that the variability can be attributed to the high turbid river plume waters in the surface layer, which are driven by the intense eastward wind.
SST, and seasonal average Kd(490). The results showed that the average Kd(490) was highly correlated with the wind speed (u-component) in summer, with an R 2 of about 0.69. During winter, we found a significant negative correlation between Kd(490) and SST, with an R 2 of -0.66 ( Figure 10). Seasonal anomalies were also obtained by subtracting the seasonal climatological average. In 2007 and 2015, when the wind speed anomaly (u-component) reached its peak, a distinguishing tongue of Kd(490) anomaly was observed near the southeastern LB ( Figure 11). Inside this tongue region, the Kd(490) anomaly in the west was higher than in the east, indicating that the variability can be attributed to the high turbid river plume waters in the surface layer, which are driven by the intense eastward wind.   The winter SST cooling in 2004 was the most significant during the whole study period, and was located in the southeastern PRE, which was about 0.4 °C cooler than the winter climatological average. Within these cooling regions, Kd(490) values higher than the average values were observed, with anomalies ranging approximately from 0.1 to 0.35 m −1 . The observed winter variations in Kd(490) in the southeastern PRE were strongly consistent with the changes in SST anomalies, and higher values coincided with lower SST (Figure 12). The variability of Kd(490) in the southeastern PRE was mainly determined by Cchla. The average values of Kd(490) were higher in winter than in summer. This seasonal variability might be attributed to the deepening of MLD in winter. In marine systems, MLD is generally deeper in winter than in summer [55]. Nutrients are brought from the bottom of the ocean to the surface or subsurface, which may enhance phytoplankton growth. The strong mixing in winter was The winter SST cooling in 2004 was the most significant during the whole study period, and was located in the southeastern PRE, which was about 0.4 • C cooler than the winter climatological average. Within these cooling regions, K d (490) values higher than the average values were observed, with anomalies ranging approximately from 0.1 to 0.35 m −1 . The observed winter variations in K d (490) in the southeastern PRE were strongly consistent with the changes in SST anomalies, and higher values coincided with lower SST (Figure 12).  The winter SST cooling in 2004 was the most significant during the whole study period, and was located in the southeastern PRE, which was about 0.4 °C cooler than the winter climatological average. Within these cooling regions, Kd(490) values higher than the average values were observed, with anomalies ranging approximately from 0.1 to 0.35 m −1 . The observed winter variations in Kd(490) in the southeastern PRE were strongly consistent with the changes in SST anomalies, and higher values coincided with lower SST (Figure 12). The variability of Kd(490) in the southeastern PRE was mainly determined by Cchla. The average values of Kd(490) were higher in winter than in summer. This seasonal variability might be attributed to the deepening of MLD in winter. In marine systems, MLD is generally deeper in winter than in summer [55]. Nutrients are brought from the bottom of the ocean to the surface or subsurface, which may enhance phytoplankton growth. The strong mixing in winter was demonstrated by the deepening of MLD, and a significant relationship between SST and MLD The variability of K d (490) in the southeastern PRE was mainly determined by C chla . The average values of K d (490) were higher in winter than in summer. This seasonal variability might be attributed to the deepening of MLD in winter. In marine systems, MLD is generally deeper in winter than in summer [55]. Nutrients are brought from the bottom of the ocean to the surface or subsurface, which may enhance phytoplankton growth. The strong mixing in winter was demonstrated by the deepening of MLD, and a significant relationship between SST and MLD provided evidence that nutrients were supplied from the bottom waters ( Figure 13). The variability of Kd(490) in the southeastern PRE was mainly determined by Cchla. The average values of Kd(490) were higher in winter than in summer. This seasonal variability might be attributed to the deepening of MLD in winter. In marine systems, MLD is generally deeper in winter than in summer [55]. Nutrients are brought from the bottom of the ocean to the surface or subsurface, which may enhance phytoplankton growth. The strong mixing in winter was demonstrated by the deepening of MLD, and a significant relationship between SST and MLD provided evidence that nutrients were supplied from the bottom waters ( Figure 13).

Conclusions
Accurate estimation of K d (490) using ocean color remote sensing imagery is challenging in turbid coastal waters due to the optical complexity of the water. Several approaches, including empirical and semianalytical models, were applied to retrieve the K d (490) in PRE water. The results showed that Wang's model was more accurate and is most suitable for PRE water, which uses a combination of different algorithms for clear and turbid waters. Hence, Wang's model was selected for deriving K d (490) products from long-term MODIS/Aqua imagery.
Derived from long-term MODIS/Aqua imagery, the temporal variability and spatial distribution of K d (490) were tracked using S-EOF analysis. The results of GRA showed that both phytoplankton and suspended sediments were the two dominant contributors to the variability in K d (490). The C tsm -dominated waters were generally located in coastal and estuarine turbid area, whereas the C chla -dominated waters were observed in clear open ocean. The influence of wind field on the variability of K d (490) was significant near the coastal and estuarine regions in summer. With the strengthening of the eastward wind, a water tongue of relatively higher K d (490) values formed in the southeastern PRE. In winter, the location of the negative SST anomaly and positive K d (490) anomaly was strongly consistent, indicating that the sea surface cooling was related to the positive K d (490) anomaly. The winter variability might be attributed to the strong mixing, which brought nutrients from the bottom layer to the surface to enhance phytoplankton growth.
Estuarine and coastal regions are complex ecosystems. To better examine the biogeochemical responses to physical events, a combination of remote sensing and coupled hydrodynamic-biological models should be applied in future research.