Steady increase in water clarity in Jiaozhou Bay in the Yellow Sea from 2000 to 2018: Observations from MODIS*

Ziyao YIN , Junsheng LI ,**, Jue HUANG , Shenglei WANG, Fangfang ZHANG,Bing ZHANG

1 Key Laboratory of Digital Earth Science, Aerospace Information Research Institute, Chinese Academy of Sciences, Beijing 100094, China

2 University of Chinese Academy of Sciences, Beijing 100049, China

3 College of Geomatics, Shandong University of Science and Technology, Qingdao 266590, China

4 Institute of Remote Sensing and Geographic Information System, Peking University, Beijing 100871, China

Abstract The Moderate Resolution Imaging Spectroradiometer (MODIS) surface reflectance data were used to analyze the temporal and spatial distribution characteristics of water clarity ( Z sd) in the Jiaozhou Bay,Qingdao, China, in the Yellow Sea from 2000 to 2018. Z sd retrieval models were regionally optimized using in-situ data with coincident MODIS images, and then were used to retrieve the Z sd products in Jiaozhou Bay from 2000–2018. The analysis of the Z sd results suggests that the spatial distribution of relative Z sd spatial characteristics in Jiaozhou Bay was stable, being higher Z sd in the southeast and a lower Z sd in the northwest.The annual mean Z sd in Jiaozhou Bay showed a significant upward trend, with an annual increase of approximately 0.02 m. Water depth and wind speed were important factors aff ecting the spatial distribution and annual variation of Z sd in Jiaozhou Bay, respectively.

Keyword: water clarity; Jiaozhou Bay; Moderate Resolution Imaging Spectroradiometer (MODIS); spatial distribution

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 (Zsd,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 theZsdof 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 suff ers 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 ofZsddistribution characteristics, especially in areas where access is diffi cult to monitor. Therefore, remote sensing technologies have gradually become an important and eff ective 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 theZsddistribution in both open ocean and inland waters can be obtained from diff erent 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 retrieveZsdcharacteristics 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 diff use attenuation coeffi cient (Kd, /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 (Rrs(λ)), 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 coeffi cients, 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 coeffi cient 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). SomeZsdretrieval 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 ofZsd, 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.Zsddistribution 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 ofinformation on the temporal and spatial distribution ofZsdin Jiaozhou Bay, and its natural and anthropogenic driving forces.

Fig.1 Locations of sampling, meteorological, and hydrological stations in Jiaozhou Bay

Based on MODIS data, we monitored the temporal and spatial changes inZsdin Jiaozhou Bay from 2000–2018 and analyzed its main influencing factors. First,we tested and optimized theZsdretrieval model suitable for Jiaozhou Bay and then produced theZsdproduct 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 ofZsdin 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 km2and 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 observedZsddata 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 (7Zsdmeasurements), July 15, 2009 (7Zsdmeasurements),December 17, 2010 (7Zsdmeasurements), August 19,2013 (3Zsdmeasurements), and September 10, 2015(5Zsdmeasurements) (Fig.1). TheZsdvalues 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 off shore 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 theZsdestimation model.The calibrated and validatedZsdestimation 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 aff ected 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 aff ect the temporal and spatial variations inZsd. The daily average wind speed during the day ofin-situ measurements was less than 4 m/s, and the water surface was relatively calm. Therefore, it can be considered that the value ofZsddid not change considerably within hours of the diff erence 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 appropriateZsdretrieval 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 preprocessing, including water-leaving reflectance correction and water body identification and extraction.

MODIS surface reflectance products (MOD09)were first atmospherically corrected for eff ects from atmospheric gases, aerosols, thin cirrus clouds, and for land adjacency eff ects (Vermote and Vermeulen,1999). But the MOD09 data was still inadequate to monitor water quality parameters because of the eff ects 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 intoRrs(λ) 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 calculatingRrs(λ) for diff erent inland water bodies with exposure to diff erent 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 eff ect 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 ofZsdmodels

Out of the 29 measuredZsddata 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 diff erent 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 eff ective pixels in the spatial window and eff ective 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 diff erence 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 eff ective pixels, not including cloud and sun-glint, in the window exceeds five and the coeffi cient of variation (CV) of the eff ective 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 theZsdvalues in Jiaozhou Bay based on the MODIS data, we used the 20 training samples to calibrate several commonly usedZsdestimation 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). TheR2value for each model was calculated. Using the test samples,we then calculated the mean relative error (MRE) and correlation coeffi cient (r) of the measured and retrieved data (Table 1). MRE was defined as follow:

whereYiandXiare the derivedZsdand field measurementsZsd, respectively.Nis 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 calculateX,Y,Zin the International Commission on Illumination (Commission Internationale de l’Eclairage, CIE) color system using theRrs(λ) 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 normalizingX,Y, andZto 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 downwardxaxis clockwise:

To eliminate systematic deviation caused by hyperspectralRrs(λ) 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-uptable in Novoa et al. (2013) with the new definition ofαstart axis and direction.

By comparison, we observed that theZsdretrieval model based onαhad the highest accuracy, with a model MRE of 36%, as well as elevated model and testR2values compared with the other models. This indicates that, compared with other commonly usedZsdretrieval models, the model based onαhad higher retrieval accuracy and was more suitable forZsdretrieval in the Jiaozhou Bay area.Zsdshowed a linear decrease with increasing hue angleα(Fig.2a). There was a negative correlation betweenZsdandα, with a modelR2of 0.61. Equation 1 defines the developed model:

Table 1 Calibration and evaluation of the various tested Z sd retrieval models for Jiaozhou Bay

Table 2 21 levels of FUI and the corresponding hue angle α values in the Forel-Ule scales

Fig.2 Calibration (a) and validation (b) of the proposed model to estimate Z sd in Jiaozhou Bay

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

3 RESULT

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

Fig.3 Mean MODIS seasonal time series and Z sd linear regression line from 2000–2018

3.1 Temporal variation in Z sd

3.1.1 Seasonal variation

TheZsdin 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-termZsdmean value of 2.1 m(Fig.3). Over the past 19 years, theZsdof Jiaozhou Bay has been slowly increasing, but there were no obvious seasonal variations. TheZsdin 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 annualZsdwas lowest in 2000, which was 1.8 m, and the highest in 2016, which was 2.2 m(Figs.4 & 5). The mean annualZsdin Jiaozhou Bay from 2000–2018 showed a significant upward trend(R2=0.78,P＜0.01), with a mean annual increase of 0.02 m. However, there was a sharp decline in theZsdin 2008 that deviated from the inter-annual trend line(Fig.5). AnnualZsdin Jiaozhou Bay was inversely correlated with wind speed (R2=0.65,P＜0.01).

3.2 Spatial distribution characterization of Z sd

The MODIS derivedZsddata from 2000–2018 were averaged to calculate theZsddistribution in Jiaozhou Bay (Fig.6a).Zsdin Jiaozhou Bay gradually increased from the northwest to the southeast. TheZsdwas highest in the southeastern region near the open sea area while the lowest value was in the northwest coastal estuary. TheZsdspatial distribution in Jiaozhou Bay was consistent for each year (Fig.4).

The spatial distribution of CV was also calculated from the MODIS derivedZsdvalues 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 certainZsdspatial variation. The lowestZsdvalues 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 thatZsdvariations 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 thatZsdvalues 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 ofZsdin Jiaozhou Bay.

4.1.1 Water depth

Xue et al. (2015) reported that water depth and runoff had the greatest impact on theZsddistribution in Chinese off shore 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 meanZsdfrom 2000–2018 to calculate the meanZsdof each water depth level area.

Fig.4 MODIS-derived mean annual Z sd in Jiaozhou Bay from 2000 to 2018

Fig.5Mean annualMODIS-derivedZ sd from2000–2018 andmeanannual windspeeds measuredat Qingdao meteorological station from 2000–2017

Zsdin Jiaozhou Bay noticeably increased with increasing water depths (Fig.6a). The meanZsdof the five regions in Fig.6a was 1.6, 1.7, 1.9, 2.0, and 2.1 m,with low and highZsdvalues at nearshore and off shore locations, respectively, which indicates that water depth is an important factor that aff ects theZsdspatial 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 resuspended when exposed to waves and currents, which increases the concentration of suspended solids in the water column. The water depth mainly aff ects the spatial distribution ofZsdbut seems to have little eff ect on the temporal variations inZsd(Figs.4 & 6).

4.1.2 Wind speed

Fig.6 Spatial distribution of Z sd in Jiaozhou bay based on the average Z sd estimates from all MODIS data collected during 2000–2018 (a) and the spatial distribution of the mean CV in Jiaozhou bay based on all Z sd estimates from the MODIS data collected during 2000–2018 (b)

The inter-annual variations in theZsdvalues had a negative correlation with wind speed (R2=0.649,P＜0.01; Fig.7a) indicating that lower wind speeds corresponded to higherZsdvalues. 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 lowerZsdvalues (Booth et al., 2000; Gao et al., 2017; Feng et al., 2019). Therefore, the significant inter-annualZsdfluctuations in Jiaozhou Bay may also be due to inter-annual variations in wind or sand resuspension.Based on the annualZsdvariations in Jiaozhou Bay(Fig.5),Zsdhas 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 aff ecting the interannualZsdvariations in Jiaozhou Bay.

4.1.3 Precipitation and river discharge

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

For the relationship between river discharge andZsdin 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 andZsd(Fig.7c), we found that there was no clear correlation between these parameters (R2=0.018,P=0.65). This suggests that the influence that river discharge has onZsdin Jiaozhou Bay is not significant, as well as the fact that river discharge is not the main factor that aff ects annual variations in Jiaozhou Bay.

4.1.4 Human activity

Human activity has an important impact on the temporal and spatialZsdvariations in Jiaozhou Bay.Zsdcan also reflect changes, to a certain extent, in ecological environments in Jiaozhou Bay, which is of great significance to urban development planning.Zsdvariations 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 aff ected by anthropogenic factors, such as industrial pollution.

During the construction of the Jiaozhou Bay Bridge(2007–2010),Zsdvalues 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 traffi c, theZsdvalues displayed an upward trend. The meanZsdduring 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 offi cially 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 annualZsdvalue in 2008.With the completion of the bridge, theZsdvalue 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 aff ected the hydrodynamic conditions in Jiaozhou Bay (Zhang et al., 2015) and further aff ected 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 theZsddistribution. Based on our analyses, it is clear that human activity rapidly changes theZsd, such that human activity has an important impact on short-term annualZsdvariations in Jiaozhou Bay.

4.1.5 Summary of driving factor

The spatialZsddistribution 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 eff ect on the northwest region, where changes in the river discharge contribute to changes inZsd. Compared with wind, rainfall and river discharge, it was found that wind speed had a greater impact on the interannualZsdvariations in Jiaozhou Bay. Meanwhile,we explored the factors governing seasonalZsdvariations. We calculated the relationships between the seasonal average wind speed, precipitation, and river discharge, and the seasonal averageZsd, and found that the correlations are not strong (R2=0.3,0.018, and 0.004, respectively). Therefore, there are possibly other factors on the seasonal scale that result in unstable seasonalZsdvariations. In addition,human activity has a significant influence on shortterm annual variations inZsd.

4.2 Uncertainty and future application

The biggest challenge associated with the longterm-scaleZsdretrieval 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) basedZsdretrieval model requiresRrsproducts (Lee et al.,2002), but because of the highly turbidity of Jiaozhou Bay, there is no availableRrsproducts. 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 forZsdretrieval using the hue angleα, is that water with diff erentZsdvalues exhibiting diff erent 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 toZsdretrieval in Jiaozhou Bay. Therefore, a regional model suitable for the inversion ofZsdin Jiaozhou Bay was proposed based onα.

Some studies have shown that there is good correlation between water color andZsd(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 ofZsdbased on FUI orα.The color of water is mainly determined by the total absorption spectrum of suspended matter, chlorophylla, 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 diff erentZsdvalues. 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 diff erentZsdvalues. 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 diff erentZsd. The atmospheric correction eff ect 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, theZsdretrieval model using the MODIS hue angleα, can largely avoid the uncertainty associated with atmospheric correction.

Note that the accuracy of the measuredZsddata may aff ect the model, especially in water characterized by highZsdvalues. TheZsdmeasurements 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 theZsdretrieval model. At the same time, the satellite transit time may be diff erent from the measurement time (e.g., more than 3 h), which will also lead to certain errors in theZsdretrieval results (Feng et al., 2012).Due to the complex water environment in Jiaozhou Bay, it is diffi cult to accurately obtain the characteristic changes inZsdusing traditional field surveys. In this study, for the first time through the analysis ofZsdMOD09A1 images in Jiaozhou Bay from the past 19 years, we were able to decipher the significant temporal and spatial changes inZsdin 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 suffi cient 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 evaluatingZsdchanges 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 theZsdin Jiaozhou Bay using the CIE color space parameter of hue angleα. Using this model, we obtained theZsdproducts in Jiaozhou Bay from 2000 to 2018. We investigated the temporal and spatial distribution characteristics of theZsdin Jiaozhou Bay, as well as an analysis ofits main influencing factors. The main conclusions of this study are as follows:

(1) The mean inter-annualZsdin Jiaozhou Bay from 2000–2018 ranged from 1.8 m to 2.5 m, with a mean value of 2.1 m. The variations inZsdwith 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 onZsd.

(2) The spatial distribution ofZsdin Jiaozhou Bay was stable, showing high values in the southeast and low values in the northwest, where variations inZsdin coastal areas were much larger than those in the inner bay. Water depth was an important factor that caused the spatial distribution ofZsdin 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 inZsdto a certain extent, where human activity, such as the large-scale construction of artificial facilities, has an impact on the short-term annual variations inZsd.

This study shows that theZsdretrieval 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 theirZsddistributions.This study also shows that this method is a support tool for environmental managers to improve decisionmaking 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-situZsddata.