1 INTRODUCTION

Water clarity (or water transparency) is an important parameter that indicates the degree of turbidity, a key ecological index that is indicative of the light transmission capacity of water (Lee et al., 2018; Shi et al., 2018). It also serves as a basic parameter in water quality investigations. Monitoring temporal and spatial changes in water clarity is of great importance to studies on water optical parameters, physical and chemical properties, water quality changes, and water eutrophication (Gong, 2012; Bukata, 2013; Yu et al., 2014). Water clarity is traditionally monitored using the Secchi disk method, which is also known as the Secchi disk depth (Z_sd, unit: m). The black-white Secchi disk with a diameter generally about 30 cm is gradually lowered into the water at the surface until it disappears out of the observer's sight. Here, the vertical distance between the observer's sight and water surface is the Z_sd of the water body. Although this method is simple to perform, the obtained data are discontinuous, the experiments are time-consuming and expensive, and real-time synchronization suffers from poor quality (Gong, 2012; Lee et al., 2018). However, remote sensing technology can acquire macroscopic, timely, accurate, and periodic data for a wide range of Z_sd distribution characteristics, especially in areas where access is diffcult to monitor. Therefore, remote sensing technologies have gradually become an important and effective means to monitor water clarity (Chen et al., 2007; Gong, 2012; Bukata, 2013; Lee et al., 2018).

Over the past 30 years, numerous studies have shown that the Z_sd distribution in both open ocean and inland waters can be obtained from different types of remote sensing satellite data, including SeaWiFS, MODIS, MERIS, and Landsat TM/ETM+/OLI. At present, these studies have used remote sensing satellite data to retrieve Z_sd characteristics and proposed various remote sensing methods. One of these methods is a semi-analytical model based on physical theories, such as a model based on the classic underwater visibility theory (Duntley, 1960), a second based on the diffuse attenuation coeffcient (K_d, /m) (Bukata et al., 1988; Kratzer et al., 2003; Aas et al., 2014; Alikas et al., 2015), and a third based on the remote sensing reflectance (R_rs(λ)), which corresponds to the spectral wavelengths from the visible domain (410-665 nm) (Lee et al., 2015, 2016). Finally, there is an empirical or semi-empirical model based on the use of band or band combinations, the band/band combination, or spectral coeffcients, which include green band algorithms at approximately 550 nm (Binding et al., 2007), red band algorithms at approximately 645 nm (Shi et al., 2018), blue-red band ratio algorithms (Kratzer et al., 2008; Olmanson et al., 2008), spectral coeffcient algorithms based on the surrounding baseline area and spectral curves in the 430- and 750-nm bands (Thiemann and Kaufmann, 2002), and algorithms based on the Forel-Ule index (FUI) and hue angle (α) (Van der Woerd and Wernand, 2015; Wang et al., 2015). Some Z_sd retrieval algorithms have been applied to analyze the temporal and spatial variations of water clarity in many areas (Barnes et al., 2013; Shang et al., 2016; Hou et al., 2017). Remote sensing technology plays an important role in the analysis of temporal and spatial variations of Z_sd, and it can also explore the causes for such changes.

Jiaozhou Bay is a semi-enclosed bay located in the southern coast of the Shandong Peninsula, China. It is the largest semi-enclosed bay of China, in Qingdao City, one of the biggest industrial cities in the northern China. With rapid coastal economic development, a large number of artificial buildings, such as aquaculture ponds, salt ponds, ports, wharves, and bridges, have gradually been built in the Jiaozhou Bay area, which has resulted in the gradual transformation of most of the coastline from natural to artificial areas (Zhang et al., 2015). In addition, Jiaozhou Bay has been eutrophicated over the last 30 years. (Jiao, 2001; Shen, 2001; Han et al., 2011). The Jiaozhou Bay Bridge, which was built in 2007, has also had a major impact on the hydrodynamic environment of Jiaozhou Bay. Collectively, Jiaozhou Bay is an ideal study site to investigate the natural and anthropogenic influences on the variation of water quality. Z_sd distribution characteristics will provide additional information regarding water mass analysis, stream identification, marine primary productivity, and regional oceanography for the waters that surrounding Jiaozhou Bay, as well as for environmental assessments and ecological protection planning. However, there is a lack of information on the temporal and spatial distribution of Z_sd in Jiaozhou Bay, and its natural and anthropogenic driving forces.

Based on MODIS data, we monitored the temporal and spatial changes in Z_sd in Jiaozhou Bay from 2000- 2018 and analyzed its main influencing factors. First, we tested and optimized the Z_sd retrieval model suitable for Jiaozhou Bay and then produced the Z_sd product for Jiaozhou Bay from 2000-2018. In addition, we analyzed the influence that the environment (water depth), climate (wind speed, precipitation, and river discharge), and human activities have on the annual and spatial variations of Z_sd in Jiaozhou Bay.

2 DATA AND METHOD 2.1 Data 2.1.1 Study area

Jiaozhou Bay (36°06'N-36°25'N; 120°10'E- 120°37'E; Fig. 1) is located in Qingdao, Shandong Peninsula, China. The coastal administrative regions include districts of Shinan, Shibei, Licang, Chengyang, Laoshan, Jimo, Huangdao, and Jiaozhou. Jiaozhou Bay covers an area of nearly 500 km² and is characterized by a complex underwater topography. The average water depth is approximately 7 m, with a maximum water depth of approximately 64 m (Song et al., 2015; Zhao et al., 2015). The main rivers flowing to Jiaozhou Bay include the Dagu, Moshui, Baisha, and Yanghe rivers. Due to the advancement of economic activities, and Qingdao City's recent population increase, these rivers have become discharge trenches for industrial and domestic waste waters (Shen, 2001; Liu et al., 2011).

2.1.2 In-situ data

The observed Z_sd data in Jiaozhou Bay were obtained from the Jiaozhou Bay Marine Ecosystem Research Station (http://jzw.qdio.cas.cn/). A total of 29 water samples from stations S1-S12 were collected on the following dates: December 11, 2008 (7 Z_sd measurements), July 15, 2009 (7 Z_sd measurements), December 17, 2010 (7 Z_sd measurements), August 19, 2013 (3 Z_sd measurements), and September 10, 2015 (5 Z_sd measurements) (Fig. 1). The Z_sd values were measured using a standard black-white disk of 30 cm.

2.1.3 Remote sensing data

The remote sensing data selected in this study was derived from MODIS-Terra data, which was launched in 2000. Therefore, we can perform long-term time series analysis using this data, which is the case for previous studies that have successfully applied this data to analyze turbid offshore and inland waters (Feng et al., 2012; Li et al., 2017; Wang et al., 2018). The MODIS surface reflectance product (MOD09) provided 500-m spatial resolution data at seven bands across the visible, near-infrared, and short-wave infrared wavelengths. Five scenes of MOD09 images, which were concurrent with the in-situ data, were used to calibrate and validate the Z_sd estimation model. The calibrated and validated Z_sd estimation model was then applied to the 8-day composite MOD09 products, which are called MOD09A1 products. The reason why MOD09A1 was used in this paper is that it can provide data over a longer time range than MYD09A1. Also, it has been used for long-term and large-area water quality monitoring research as this dataset is well georeferenced, synthesized, and cloud marked (Feng et al., 2012; Wu et al., 2013; Li et al., 2016; Hou et al., 2017; Klein et al., 2017; Wang et al., 2018). The MOD09A1 data, which covered Jiaozhou Bay from 2000-2018, was obtained from the U.S. NASA Goddard Space Flight Center (GSFC, http://oceancolor.gsfc.nasa.gov/). The MOD09A1 data contains 868 good-quality scenes (spring: 209 scenes, summer: 209 scenes, autumn: 228 scenes, winter: 222 scenes) that were selected after visual examination and exclusion of scenes that were affected by clouds, sun glint, and thick aerosols.

2.1.4 Meteorological and hydrological data

We obtained daily average wind speed and precipitation data from 2000-2017 from the China Meteorological Data Network (http://data.cma.cn), which was collected from the nearest meteorological station (Qingdao Station) in Jiaozhou Bay. This was done in order to analyze the factors that affect the temporal and spatial variations in Z_sd. The daily average wind speed during the day of in-situ measurements was less than 4 m/s, and the water surface was relatively calm. Therefore, it can be considered that the value of Z_sd did not change considerably within hours of the difference between the in-situ measured data and the MODIS retrieved data. We also obtained monthly average river discharge data from the Nancun Hydrological Station in Dagu River from 2000-2013.

2.2 Method

To determine an appropriate Z_sd retrieval algorithm for Jiaozhou Bay, the MOD09A1 images first required preprocessing, followed by calibration and the testing of various retrieval methods based on satellite-ground synchronous data.

2.2.1 Data pre-processing

The MOD09A1 data required a series of pre-processing, including water-leaving reflectance correction and water body identification and extraction.

MODIS surface reflectance products (MOD09) were first atmospherically corrected for effects from atmospheric gases, aerosols, thin cirrus clouds, and for land adjacency effects (Vermote and Vermeulen, 1999). But the MOD09 data was still inadequate to monitor water quality parameters because of the effects of residual aerosol scattering, solar flares, and skylight reflection. In this study, a method was adopted to reduce noise in the MOD09 data and convert it into R_rs(λ) by subtracting the near-infrared (NIR) minimum value from the shortwave infrared (SWIR) band and dividing by π (Wang et al., 2016). Previous studies have shown that this method easily and steadily yields good performance when calculating R_rs(λ) for different inland water bodies with exposure to different conditions (Wang et al., 2016, 2018).

Due to influences from changes in water volume, the water boundary, to a certain extent, is characterized by dynamic changes. Therefore, we used the modified histogram bimodal method (MHBM) to segment the water body (Zhang et al., 2018) based on the reflectance of the short wave infrared remote sensing data in the 1 640 nm band. In order to reduce the effect of adjacent land pixels (which have a higher reflectance) the coastline was extended seawards by two pixels (Feng and Hu, 2017).

2.2.2 Calibration and validation of Z_sd models

Out of the 29 measured Z_sd data samples, 20 were selected as the training samples to construct the model. The nine remaining samples were used as test samples to verify the model. Satellite products and field measurements have different temporal and spatial resolutions. It is necessary to determine a reasonable spatiotemporal window according to the spatial resolution of satellite products, as well as the spatiotemporal variation and uniformity of water bodies. Means of effective pixels in the spatial window and effective field measurements in the time window were taken as matching data pairs and incorporated into the validation data set. Several criteria were followed in obtaining matching pairs (Le et al., 2013; Feng et al., 2015). First, the time difference between field measurements and remote sensing images should be within ±3 h, and the wind speed during the day should be less than 4 m/s. Second, the water around the observation station should be uniform. Thus, a 3×3 pixel window centered on the measured site was selected. If the number of effective pixels, not including cloud and sun-glint, in the window exceeds five and the coeffcient of variation (CV) of the effective pixels is less than 0.4, the average pixel value of the window is calculated as the final value.

To obtain the optimal model to retrieve the Z_sd values in Jiaozhou Bay based on the MODIS data, we used the 20 training samples to calibrate several commonly used Z_sd estimation models, including the red band algorithms (Shi et al., 2018), blue-red band ratio algorithms (Kratzer et al., 2008; Olmanson et al., 2008), green band algorithms (Binding et al., 2015), red-green band ratio algorithms (Duan et al., 2009; Ren et al., 2018), the Forel-Ule index (FUI), and the hue angle (α) algorithms (Van der Woerd and Wernand, 2015; Wang et al., 2015). The R² value for each model was calculated. Using the test samples, we then calculated the mean relative error (MRE) and correlation coeffcient (r) of the measured and retrieved data (Table 1). MRE was defined as follow:

where Y_i and X_i are the derived Z_sd and field measurements Z_sd, respectively. N is the number of matchups.

Out of these algorithms, the methods based on the FUI and α were the most complex. To obtain these two parameters, we used the Red, Green, Blue (RGB) conversion method to calculate X, Y, Z in the International Commission on Illumination (Commission Internationale de l'Eclairage, CIE) color system using the R_rs(λ) of three visible bands (469, 555, and 645 nm) (CIE, 1931; Wang et al., 2015; Li et al., 2016).

The CIE chromaticity coordinates (x, y) were obtained by normalizing X, Y, and Z to 0 and 1 as follows (CIE, 1931; Wang et al., 2015).

The hue angle α were calculated by the new definition in Wang et al. (2018), which has positive range from 0 to 2π from the downward x axis clockwise:

To eliminate systematic deviation caused by hyperspectral R_rs(λ) and the equivalent MODIS bands, the corrected hue angle α was obtained using the method reported in Wang et al. (2018), which was based on the methods proposed in Van der Woerd and Wernand (2015), but with the new definition of α start axis and direction. Finally, the FUI was calculated using the lookup table (Table 2) based on the corrected hue angle α, which is transformed from the look-up-table in Novoa et al. (2013) with the new definition of α start axis and direction.

By comparison, we observed that the Z_sd retrieval model based on α had the highest accuracy, with a model MRE of 36%, as well as elevated model and test R² values compared with the other models. This indicates that, compared with other commonly used Z_sd retrieval models, the model based on α had higher retrieval accuracy and was more suitable for Z_sd retrieval in the Jiaozhou Bay area. Z_sd showed a linear decrease with increasing hue angle α (Fig. 2a). There was a negative correlation between Z_sd and α, with a model R² of 0.61. Equation 1 defines the developed model:

The paired data were characterized by a significant correlation that was predominantly located near the 1:1 line (Fig. 2b; R²=0.518), which demonstrated that the model was effective and can be used to derive the Z_sd in Jiaozhou Bay. Therefore, the model based on α was chosen for the derivation and the further analysis of the Z_sd in Jiaozhou Bay.

3 RESULT

Based on the MOD09A1 products that cover Jiaozhou Bay from 2000-2018, the long-term Z_sd product for Jiaozhou Bay was derived using the established retrieval model based on α. Based on this product, we analyzed the Z_sd temporal and spatial distribution and variation characteristics in Jiaozhou Bay.

3.1 Temporal variation in Z_sd 3.1.1 Seasonal variation

The Z_sd in Jiaozhou Bay varied significantly, ranging from a minimum value of 1.8 m in the spring of 2000 to a maximum value of 2.5 m in the winter of 2012, with a mean long-term Z_sd mean value of 2.1 m (Fig. 3). Over the past 19 years, the Z_sd of Jiaozhou Bay has been slowly increasing, but there were no obvious seasonal variations. The Z_sd in summer (June- August) and winter (December-February) were higher than that in spring (March-May) and autumn (September-November), but the seasonal variations were not stable from year to year.

3.1.2 Inter-annual variation

The mean annual Z_sd was lowest in 2000, which was 1.8 m, and the highest in 2016, which was 2.2 m (Figs. 4 & 5). The mean annual Z_sd in Jiaozhou Bay from 2000-2018 showed a significant upward trend (R²=0.78, P < 0.01), with a mean annual increase of 0.02 m. However, there was a sharp decline in the Z_sd in 2008 that deviated from the inter-annual trend line (Fig. 5). Annual Z_sd in Jiaozhou Bay was inversely correlated with wind speed (R²=0.65, P < 0.01).

3.2 Spatial distribution characterization of Z_sd

The MODIS derived Z_sd data from 2000-2018 were averaged to calculate the Z_sd distribution in Jiaozhou Bay (Fig. 6a). Z_sd in Jiaozhou Bay gradually increased from the northwest to the southeast. The Z_sd was highest in the southeastern region near the open sea area while the lowest value was in the northwest coastal estuary. The Z_sd spatial distribution in Jiaozhou Bay was consistent for each year (Fig. 4).

The spatial distribution of CV was also calculated from the MODIS derived Z_sd values from 2000-2018 (Fig. 6b). The CV increased from the inside to the outside of Jiaozhou Bay, ranging from 0 to 82%, with an average of 27%. This demonstrated that there was a certain Z_sd spatial variation. The lowest Z_sd values and the largest CV values were found in the coastal areas of Jiaozhou Bay that are connected to inland regions by several rivers and streams, while the lowest values were observed in the middle of Jiaozhou Bay. This indicates that Z_sd variations in the coastal areas of Jiaozhou Bay which are strongly influenced by humans activities and coastal runoff.

4 DISCUSSION 4.1 Driving factor

The above results demonstrate that Z_sd values in Jiaozhou Bay had clear temporal and spatial variations. We attempted to analyze the influence that the topography (water depth), hydrological (wind speed, precipitation, and river discharge), and human activities have on the annual, seasonal, and spatial variations of Z_sd in Jiaozhou Bay.

4.1.1 Water depth

Xue et al. (2015) reported that water depth and runoff had the greatest impact on the Z_sd distribution in Chinese offshore waters. Based on the water depth distribution map of Jiaozhou Bay reported in Ye et al. (2009), we divided it into five levels from low to high (Fig. 6a) and superimposed it with the mean Z_sd from 2000-2018 to calculate the mean Z_sd of each water depth level area.

Z_sd in Jiaozhou Bay noticeably increased with increasing water depths (Fig. 6a). The mean Z_sd of the five regions in Fig. 6a was 1.6, 1.7, 1.9, 2.0, and 2.1 m, with low and high Z_sd values at nearshore and offshore locations, respectively, which indicates that water depth is an important factor that affects the Z_sd spatial distribution in Jiaozhou Bay. Shallow depths are usually found in near-coastal areas while deeper waters are found in the inner bay. Coastal runoff transports a large amount of suspended sediment to the sea. In shallow coastal areas bottom sediments are easily re-suspended when exposed to waves and currents, which increases the concentration of suspended solids in the water column. The water depth mainly affects the spatial distribution of Z_sd but seems to have little effect on the temporal variations in Z_sd (Figs. 4 & 6).

4.1.2 Wind speed

The inter-annual variations in the Z_sd values had a negative correlation with wind speed (R²=0.649, P < 0.01; Fig. 7a) indicating that lower wind speeds corresponded to higher Z_sd values. This is possibly due to the fact that increased wind speeds in Jiaozhou Bay facilitate resuspension and transportation of sediments deposited on the seabed, thus leading lower Z_sd values (Booth et al., 2000; Gao et al., 2017; Feng et al., 2019). Therefore, the significant inter-annual Z_sd fluctuations in Jiaozhou Bay may also be due to inter-annual variations in wind or sand resuspension. Based on the annual Z_sd variations in Jiaozhou Bay (Fig. 5), Z_sd has been significantly higher since 2015 compared with the years before 2015, which is due to lower wind speeds between 2015 and 2018 compared with previous years. The results clearly suggest that wind speed is an important factor affecting the inter-annual Z_sd variations in Jiaozhou Bay.

4.1.3 Precipitation and river discharge

Precipitation data were averaged with respect to year and compared with the mean annual Z_sd (Fig. 7b). We observed that the correlation between these parameters was not significant for the annual time scales (R²=0.044, P=0.40). The analysis showed that precipitation has a relatively small effect on the Z_sd in Jiaozhou Bay and thus is not a key factor influencing inter-annual Z_sd variations in Jiaozhou Bay. This is normal, because precipitation mainly affects the short-term change in Z_sd, but not the long-term change such as inter-annual change (Song et al., 2015).

For the relationship between river discharge and Z_sd in Jiaozhou Bay, we must consider the 16 rivers that flow into Jiaozhou Bay, which include the Dagu, Moshui, Baisha, Yanghe, and Licun rivers. Out of these rivers, Dagu River has the largest catchment area, which occupies 82% of the total catchment area of the Jiaozhou Bay river system (Jiang and Wang, 2013). Average monthly river discharge from Nancun Hydrological Station, which exists for 2000-2013, was used to analyze the relationship. Nancun Hydrological Station is located in the middle and lower reaches of the main stem of the Dagu River. Nancun is the closest station to the lower reaches of the Dagu River, with long-term river discharge data, and is representative of the overall flow characteristics in the basin. Dagu River is a river that is mostly fed by precipitation, and thus is subject to strong fluctuations in water levels. The flood season for the Dagu River occurs in the summer.

Using the same method to analyze the relationship between the mean annual river discharge and Z_sd (Fig. 7c), we found that there was no clear correlation between these parameters (R²=0.018, P=0.65). This suggests that the influence that river discharge has on Z_sd in Jiaozhou Bay is not significant, as well as the fact that river discharge is not the main factor that affects annual variations in Jiaozhou Bay.

4.1.4 Human activity

Human activity has an important impact on the temporal and spatial Z_sd variations in Jiaozhou Bay. Z_sd can also reflect changes, to a certain extent, in ecological environments in Jiaozhou Bay, which is of great significance to urban development planning. Z_sd variations in the bay were mainly concentrated in coastal areas, whereas variations in the center of the bay were relatively small. This is likely due to the fact that coastal areas are more directly affected by anthropogenic factors, such as industrial pollution.

During the construction of the Jiaozhou Bay Bridge (2007-2010), Z_sd values in the Jiaozhou Bay decreased compared with those of previous years, reaching their lowest level in 2008, which was only 1.9 m. After the bridge was completed and opened to traffc, the Z_sd values displayed an upward trend. The mean Z_sd during bridge construction was 2.0 m, while that after bridge construction was 2.1 m, mainly due to the construction of Jiaozhou Bay Bridge, which offcially began in 2007 and was completed the following year. A large number of operations associated with bridge foundations churned the water, which resulted in a sharp decline in the mean annual Z_sd value in 2008. With the completion of the bridge, the Z_sd value in the Jiaozhou Bay area began to gradually rise beginning in 2009. Relevant studies have shown that the construction of the Jiaozhou Bay Bridge has affected the hydrodynamic conditions in Jiaozhou Bay (Zhang et al., 2015) and further affected the spatial and temporal distribution of sea ice in Jiaozhou Bay in the winter (Huang et al., 2019). We infer that the changes in the hydrodynamic conditions in Jiaozhou Bay caused by the construction of the sea-crossing bridge, as well as operations, such as bridge piling, will cause sediment resuspension, which also has an impact on the Z_sd distribution. Based on our analyses, it is clear that human activity rapidly changes the Z_sd, such that human activity has an important impact on short-term annual Z_sd variations in Jiaozhou Bay.

4.1.5 Summary of driving factor

The spatial Z_sd distribution in Jiaozhou Bay displayed a trend of elevated values in the southeast and lower values in the northwest. The water depth distribution had a trend of high values in the southeast and decreased values in the northwest. In contrast, the northwestern region of Jiaozhou Bay is contiguous with inland areas, and is significantly influenced by human activity, such as industrial pollution, while the southeastern areas connect with the sea, such that the water is relatively clean. In addition, rivers along the coast, such as the Dagu, Moshui, and Baisha rivers, are mainly located in the northwestern and northern parts of Jiaozhou Bay. Therefore, when it rains, rivers have a greater effect on the northwest region, where changes in the river discharge contribute to changes in Z_sd. Compared with wind, rainfall and river discharge, it was found that wind speed had a greater impact on the inter-annual Z_sd variations in Jiaozhou Bay. Meanwhile, we explored the factors governing seasonal Z_sd variations. We calculated the relationships between the seasonal average wind speed, precipitation, and river discharge, and the seasonal average Z_sd, and found that the correlations are not strong (R²=0.3, 0.018, and 0.004, respectively). Therefore, there are possibly other factors on the seasonal scale that result in unstable seasonal Z_sd variations. In addition, human activity has a significant influence on short-term annual variations in Z_sd.

4.2 Uncertainty and future application

The biggest challenge associated with the long-term-scale Z_sd retrieval model lies in its application to regional and seasonal scales, especially the stability of retrieval model in turbid waters. The commonly used quasi-analytical algorithm (QAA) based Z_sd retrieval model requires R_rs products (Lee et al., 2002), but because of the highly turbidity of Jiaozhou Bay, there is no available R_rs products. Therefore, only the surface reflectance product MOD09A1 can be used, but MOD09A1 only has 469, 555, 645, and 859 available bands, so QAA algorithm cannot be used. Therefore, other empirical or semi-empirical methods should be used. The theoretical basis for Z_sd retrieval using the hue angle α, is that water with different Z_sd values exhibiting different hues, where α can more accurately capture this change in information. Compared with other commonly used models, we observed that the model based on the hue angle had a higher accuracy with respect to Z_sd retrieval in Jiaozhou Bay. Therefore, a regional model suitable for the inversion of Z_sd in Jiaozhou Bay was proposed based on α.

Some studies have shown that there is good correlation between water color and Z_sd (Novoa et al., 2013; Van der Woerd and Wernand, 2015; Li et al., 2016; Wang et al., 2018). However, there are some limitations in the method of Z_sd based on FUI or α. The color of water is mainly determined by the total absorption spectrum of suspended matter, chlorophyll a, colored dissolved organic matter (CDOM), and pure water. When the optical properties of water are dominated by suspended matter, CDOM, or pure water, it is not easy to have the same color corresponding to different Z_sd values. However, when algae blooms occur, phytoplankton dominates the total absorption spectrum of the water, such that the color of the water does not change monotonously with the wavelength and the same color may correspond to different Z_sd values. Fortunately, for Jiaozhou Bay, which has no large-scale and long-term water bloom and whose absorption spectrum is dominated by suspended matter, the same water color does not correspond to different Z_sd. The atmospheric correction effect for remote sensing images is another challenge to water quality parameter retrieval. In this study, we adopted a correction method using the water-leaving reflectance based on the MODIS surface reflectance product, which essentially meets the retrieval requirements for the Jiaozhou Bay water-leaving reflectance. Furthermore, several previous and related studies have indicated that α was less sensitive to observation conditions and aerosol types (Wang et al., 2018; Pitarch et al., 2019). In addition, the α value used in this model was derived from the corrected synchronous MODIS surface reflectance products, rather than the measured water-leaving reflectance, which could further reduce the impact that the atmosphere has on the retrieval results. Therefore, the Z_sd retrieval model using the MODIS hue angle α, can largely avoid the uncertainty associated with atmospheric correction.

Note that the accuracy of the measured Z_sd data may affect the model, especially in water characterized by high Z_sd values. The Z_sd measurements are vulnerable to wind speed, wave, and even ship navigation, which may lead to associated uncertainties. Therefore, high-precision data measurements will reduce model uncertainty. More measured data can further improve the Z_sd retrieval model. At the same time, the satellite transit time may be different from the measurement time (e.g., more than 3 h), which will also lead to certain errors in the Z_sd retrieval results (Feng et al., 2012).

Due to the complex water environment in Jiaozhou Bay, it is diffcult to accurately obtain the characteristic changes in Z_sd using traditional field surveys. In this study, for the first time through the analysis of Z_sd MOD09A1 images in Jiaozhou Bay from the past 19 years, we were able to decipher the significant temporal and spatial changes in Z_sd in Jiaozhou Bay and its driving factors. Although the MODIS data is characterized by strong stability and the advantages of a global data product service, its spatial resolution is not suffcient to monitor small bodies of water. With the atmospheric correction for high spatial and spectral resolution remote sensing data, such as Landsat-8 OLI and the production of water products, the algorithm established in this study can be used with new sensor images to cover more water bodies, which can provide new possibilities and challenges when evaluating Z_sd changes in Jiaozhou Bay, as well as providing certain theoretical support to environmental assessments and ecological protection decision-making for the relevant governmental departments.

5 CONCLUSION

In this study, using an eight-day composite MODIS-Terra surface reflectance product (MOD09A1), we established a model to derive the Z_sd in Jiaozhou Bay using the CIE color space parameter of hue angle α. Using this model, we obtained the Z_sd products in Jiaozhou Bay from 2000 to 2018. We investigated the temporal and spatial distribution characteristics of the Z_sd in Jiaozhou Bay, as well as an analysis of its main influencing factors. The main conclusions of this study are as follows:

(1) The mean inter-annual Z_sd in Jiaozhou Bay from 2000-2018 ranged from 1.8 m to 2.5 m, with a mean value of 2.1 m. The variations in Z_sd with time showed an upward trend, with a mean annual increase of 0.02 m. Comparing wind, precipitation, river discharge, and other factors, we found that wind speed has a greater impact on Z_sd.

(2) The spatial distribution of Z_sd in Jiaozhou Bay was stable, showing high values in the southeast and low values in the northwest, where variations in Z_sd in coastal areas were much larger than those in the inner bay. Water depth was an important factor that caused the spatial distribution of Z_sd in Jiaozhou Bay. Deeper waters yield higher Secchi depths whilst more shallow areas are more influenced by resuspension of bottom sediments.

(3) Human activity has also caused changes in Z_sd to a certain extent, where human activity, such as the large-scale construction of artificial facilities, has an impact on the short-term annual variations in Z_sd.

This study shows that the Z_sd retrieval model based on hue angle α can be applied to turbid water. Although we only selected Jiaozhou Bay to test this model, we suggest that this method can be extended to other water bodies with similar optical characteristics and ranges in Secchi depth, such as Liaodong Bay, Bohai Bay, and Laizhou Bay to explore their Z_sd distributions. This study also shows that this method is a support tool for environmental managers to improve decision-making and governmental policies in Jiaozhou Bay.

6 DATA AVAILABILITY STATEMENT

The data of this study are available from the Jiaozhou Bay Marine Ecosystem Research Station but restrictions may apply due to the availability of some of the data that were used under license for this study, and thus, they are not publicly available. The data are however available from the authors upon reasonable request and with permission from the Jiaozhou Bay Marine Ecosystem Research Station.

7 ACKNOWLEDGMENT

The authors would like to thank NASA GSFC for providing MODIS data, and we are also grateful to Jiaozhou Bay Marine Ecosystem Research Station for providing in-situ Z_sd data.