Mapping Dynamic Turbidity Maximum Zone of the Yellow River Estuary from 38 Years of Landsat Imagery

Abstract

Turbidity maximum zone (TMZ) plays a crucial role in estuarine ecosystems, exerting effects on erosion, environment evolution and socioeconomic activities in the coastal area. However, the long-term understanding of the TMZ in large river estuary such as the Yellow River estuary is still lacking. In this study, we focus on the TMZ distribution, variation and regulation mechanisms in the Yellow River estuary from different time scales. Based on time series Landsat images during the period 1984 to 2021 and Google Earth Engine (GEE), we proposed a TMZ extracting method in the Yellow River estuary to generate 322 TMZ maps. The overall accuracy of our algorithm reached 97.4%. The results show that there are clear decadal and seasonal TMZ variations during the 38-year period in the Yellow River estuary. Morphology, currents and wind speeds combined with seawater stratification have direct effects on TMZ at different time scales, while the direct impacts of tides and fluvial output of the Yellow River on TMZ are limited. In this article, the highly robust method provides a cost-effective alternative to accurately map the TMZ in global large river estuaries and systematically reveals the spatiotemporal evolution of TMZ, shedding light on the response mechanism of coastal geomorphology, marine ecological environment and biogeochemical cycle.

1. Introduction

As the feedback of complex interactions among hydrodynamics, sediment transport and morphodynamics [1], turbidity maximum zone (TMZ) is the region within an estuary where the suspended sediment concentration (SSC) is much higher than that in surrounding waters [2,3]. The suspended sediment dynamics in the TMZ plays an important role in estuarine ecosystems and biogeochemical cycles [4], affecting the light propagation in water [5] and transportation of nutrients and pollutants [6], exerting influences on erosion and biotic community in the coastal fringe. It also has socioeconomic implications, such as the impact on coastal engineering stability and production activities [7,8]. Thus, to sustain its ecological and social services, it is critical to identify the extent, formation and variation of TMZ.

The Yellow River has been ranked the second largest sediment-laden river in the world [9] and is well-known for its frequent channel migrations in history. The large amount river outflow into sea causes huge carbon flux and strongly affects the marine environment in the river estuary [10]. Meanwhile, as the generalization of highly turbid water, TMZ could reflect the status of sediment dynamics. Studies on the TMZ in global large river estuaries, such as the Yellow River estuary, have significant scientific implications in promoting cycles of carbon, oxygen and nutrients in coastal and surrounding sea areas [11].

The spatial and temporal variation of estuarine SSC is drastic. This heterogeneity makes it difficult for in situ monitoring to obtain a synoptic view of estuarine SSC and TMZ [12]. Remote sensing methods based on satellite imagery have the advantages of large monitoring extent and regular revisit period. Ocean color satellite sensors, such as Moderate Resolution Imaging Spectroradiometer (MODIS), Medium Resolution Imaging Spectrometer (MERIS) and Geostationary Ocean Color Imager (GOCI), have been widely utilized in estuarine turbidity monitoring [12,13,14,15]. Nevertheless, the coarse spatial resolution (250 m to 1 km) of these sensors is incompetent in capture the detailed information of SSC in relatively small areas such as estuaries [12]. Thus, the higher spatial resolution of 30 m for Landsat sensors makes them more suitable in estuary application [16].

Both SSC and suspended particulate matter (SPM) are important aspects of water turbidity. SSC or SPM assessment based on satellite imagery have been carried out in the Yellow River estuary in recent years [15,17,18,19,20,21,22], identifying their pattern and influencing factors. By means of developing model from in situ data, the SSC in the Yellow River estuary was retrieved [17,19,21]. Interannual and seasonal variations plus spatial distributions of SSC in the Yellow River estuary from 2000 to 2010 were investigated, showing the strong influence of the Yellow River sediment output on the estuary [17]. Correlation analysis revealed that wind speed reduction was critical in the decreasing trend of SSC from 2003 to 2016 in the Yellow River estuary [21]. In addition, SPM in the Yellow River estuary was also estimated from tuned model [18,20,22]. The SPM variance in dry and wet seasons has been observed during 2013 to 2016 [18] and hourly river mouth SPM dynamics has been explored [20]. Using a clustering method, Wang et al. [15] further identified six turbidity patterns. Using time series satellite images, Li et al. [22] studied the spatiotemporal dynamics of SPM in the Yellow River estuary from 1997 to 2018. However, all these studies did not extract the extent of TMZ and did not treat TMZ as an integral phenomenon.

The retrieval model requires cruise data, which needs much financial support and labor power. Theoretically, in consideration of the high variability of sediment dynamics, retrieval model demands concurrent field samplings and satellite overpass, while this is almost impossible [13]. Moreover, signals of satellite sensor bands could saturate in high turbidity waters [23], resulting unreliable retrieval model results. In addition, current remote sensing sediment dynamics studies in the Yellow River estuary rarely utilize images for more than 30 years, prohibiting the understanding from a larger picture. Finally, the regulation mechanism of sediment dynamics has not been thoroughly understood and is not unanimous [21].

TMZ has been systematically studied in other estuaries around the world [8]. However, long time series patterns and mechanisms of TMZ in the Yellow River estuary are not clear at present. It needs to be noted that the relatively broad TMZ range in the Yellow River estuary makes it closer to the coastal turbidity maximum. To solve this problem, there are two clear obstacles: a high computing power demand for efficient large amounts of images procession and a robust TMZ extraction algorithm. Google Earth Engine (GEE) is a cloud-based platform that possesses massive computational capabilities and archives multi-petabyte analysis-ready data with an intrinsically parallel computation service [24]. The Landsat images and computing power of GEE could greatly facilitate the long time series TMZ research.

Considering these challenges, we aim to develop a long time series TMZ extracting algorithm based on Landsat images using GEE. This algorithm can avoid the dependence on field data and model development with high efficiency and reproductivity. Then, we will investigate the time series variation, decadal and seasonal pattern of TMZ from 1984 to 2021 in the Yellow River estuary. In addition, we will also discuss the various interactive factors affecting TMZ distribution at different time scales.

2. Materials and Methods

2.1. Study Area

With a length of 5464 km, the Yellow River originates from the northern Qinghai-Tibetan plateau and meanders eastward before entering into the Bohai Sea [25] (Figure 1). In history, about 1.64 × 10⁹ tons of sediment and 5.74 × 10¹⁰ m³ of freshwater were discharged into the Bohai Sea annually, accounting for approximately 90% of sediment and 50–60% freshwater the Bohai Sea received [18,26]. The currently active Yellow River Delta (YRD) has been quickly built since the diversion of the Yellow River from the Diaokou course to the Qingshuigou course in 1976. In addition, to facilitate the construction of local oilfields, the river channel was slightly redirected to the current Qing 8 course in 1996 [25]. Sediments discharged from the Yellow River mouth tend to migrate southward [27].

The Yellow River estuary belongs to the temperate monsoon climate zone. Weak southerly winds prevail in summer, while frequent strong northly winds dominate in winter [20]. Furthermore, most local waves belong to wind waves and thus have clear seasonal variability. This microtidal estuary is dominated by an irregular semi-diurnal tide, with an average tidal range of 0.6 to 0.8 m at the river mouth and 1.5 to 2.0 m in the Laizhou Bay [28]. The water depth of Laizhou Bay is less than 15 m. More than half of it is shallower than 10 m, especially along the coast [29]. In general, currents off the Yellow River estuary head south [18]. Most suspended sediments of the Yellow River estuary are clay, silt and fine sand. Since the channel diversion in 1976, most of these riverine sediments have been deposited in Laizhou Bay [19].

2.2. Datasets

2.2.1. Landsat Images

Currently, the Landsat surface reflectance dataset on GEE comes from five Landsat satellites: Landsat 4, 5, 7, 8 and 9. They have been successively launched from 1982 to 2021, carrying Thematic Mapper (TM), TM, Enhanced Thematic Mapper Plus (ETM+), Operational Land Imager (OLI) and OLI-2 sensors, respectively. The Landsat series run around the same sun-synchronous track, which means two satellites in orbit at a given time could achieve an 8-day revisit frequency [30]. The Landsat 4, 5, 7, 8 and 9 share the same global notation system (Worldwide Reference System-2) and take advantage of the similar Path/Row system. All the spectral mappers or imagers loaded on these satellites could achieve a 30 m spatial resolution in red, green, blue and NIR bands.

GEE archives the whole Landsat time series surface reflectance images. In addition, the surface reflectance data of Landsat 4, 5, 7, 8 and 9 has been thoroughly pre-processed on GEE [24], which makes them analysis-ready. In this study, we used the tile of Path 121 and Row 34, since this tile could entirely cover the study area. In total, 761 Landsat images acquired during the period 1984–2021 were processed for water–land boundary generating, including Landsat 4 (2 images), 5 (325 images), 7 (299 images), 8 (133 images) and 9 (2 images). After screening, we obtained 322 dense cloud and cloud-shadow-free Collection-1 Landsat images from 1984 to 2021 for TMZ mapping, including Landsat 4, 5, 7 and 8 (Figure 2). The whole procedures of Landsat data processing were implemented on GEE.

2.2.2. Tide Data

As shown in Figure 1, the tide height data at two tidal stations, namely Yangjiaogou (abandoned in 2011) and Weifanggang (2011 to present), were obtained from tide tables published by National Marine Data and Information Service (http://global-tide.nmdis.org.cn/, accessed on 1 June 2022). The tidal datum has been united to 107 cm below mean sea level.

2.2.3. Water Discharge and Sediment Load

The runoff and sediment load data of the Yellow River to the sea are recorded at the Lijin gauging station (Figure 1), which is the closest main hydrometric station to the Yellow River estuary (about 100 km upstream the estuary). These data were collected from the Yellow River Sediment Bulletin published by the Yellow River Conservancy Commission of the Ministry of Water Resources (http://www.yrcc.gov.cn/, accessed on 31 May 2022).

2.2.4. Wind Speed Data

The wind speed data were obtained from the fifth generation European Centre for Medium-Range Weather Forecasts reanalysis for the global climate and weather (ERA5), with temporal resolution of 1 h and spatial resolution of 0.25° × 0.25° (https://cds.climate.copernicus.eu/, accessed on 3 June 2022).

Considering that the acquisition time of all Landsat images of Path 121 and Row 34 was centered at around 2:00 in UTC, we used the average wind speed during 1:00 to 2:59 in UTC in this study.

2.3. Methods

Our TMZ mapping algorithm was mainly composed of three parts: (1) high-quality images collection, (2) land–water separation and (3) TMZ extraction through red band thresholding (Figure 3).

2.3.1. High-Quality Imagery Selection

In this study, high-quality images mean images free of dense cloud and cloud shadows that could lead to severe ground information loss [31]. First, based on the ‘CLOUD_COVER’ property of Landsat images on GEE, the original image archive was filtered with the restriction of cloud cover less than 60% in the tile Path 121 Row 34 (Figure 3a). Then, these filtered images were clipped to the study area. These images were used for two purposes, that is, water–land boundary generation and TMZ delineating (Figure 3b).

For the water–land boundary generation, the bad-quality pixels with cloud and cloud shadows were masked out in images based on the quality assessment bands (QA bands). The QA bands were integrated in the Landsat surface reflectance dataset packaged by GEE, flagging bad-quality pixels. The derived cloud-masked images were ready for Normalized Difference Water Index (NDWI) compositing (Section 2.3.2).

For single-image TMZ delineating, images are required to be generally intact and free of cloud coverage. Within the extent of the study area, we managed to calculate the cloud and cloud shadow percentage through QA bands and limited this percentage to less than 10% through programming on GEE (Figure 3b). In addition, small cloud and cloud shadow patches were still masked out. The obtained images could be regarded as intact and clear images.

2.3.2. Land–Water Boundary Generation

The accuracy of TMZ mapping could be affected by inland features such as tidal flats, inland turbid water bodies and aquacultures, as they may share spectral characteristics with TMZ. Previous related studies often ignored the lowest tide line during long time periods in the Yellow River estuary [15,17,18,19], causing inland features (especially tidal flats) to be mixed into the sea water turbidity mapping results. Thus, to avoid misleading outcomes, it is crucial to mask out inland features before TMZ mapping.

Image compositing methods select suitable pixels for each pixel location from an image stack and can robustly generate composite images representing different tidal inundation statuses [32]. Tidal flats are non-vegetated areas between the highest and lowest tidal inundations [33]. The lowest tide composite images represent the largest land extent, which includes the whole tidal flats. Therefore, we used the lowest tide composite images to generate the land–water boundary, which could maximize the removal of inland features from sea waters.

Considering that the NDWI extremum composites have been successfully applied in representing the tidal inundation status, that is, maximum composite denotes the highest tide condition and the minimum composite denotes the lowest tide condition [32,33,34], we selected the NDWI minimum composites to derive the land–water boundaries. In light of the relatively low revisit frequency of Landsat satellites (16 days when alone in orbit), we referred to the scheme adopted by Murray et al. [34] in which the composites were derived from Landsat image stacks of every three years. From 1 January 1983–1 January 1986 to 1 January 2019–1 January 2022, 13 sets of Landsat image stack were built. To obtain 13 NDWI minimum composites, it is theoretically reasonable with all images in every image stack to composite the NDWI minimums into one image.

However, this method could accumulate lots of noise in the composite images [32,33]. Previous researchers have found positive relationships between the NDWI of tidal flat and tidal height [32]. Taking advantage of the time series tidal flat datasets beginning in 1984 (http://intertidal.app, accessed on 1 April 2022) [34], we flagged each Landsat image with the sum of NDWI value within the tidal flat regions and sorted them in ascending order. The top 20% ranked scenes were regarded as the low tide images [32]. The NDWI minimum composite image was generated through minimum NDWI compositing on these 20% images. In total, 13 composite images were derived. In this way, the occurrence of noise was greatly reduced in the results.

In order to separate land and water in the NDWI composites, the adaptive binary classification algorithm, namely the Otsu method [35], was applied. The Otsu method could automatically select the optimal dividing value of two categories and has shown good performance in separating land and water [33,36,37]. However, in a small portion of the study area, the Otsu-derived boundary could be ambiguous due to inland water and coastal aquacultures, so some minor adjustments were made to delineate the final land–water line (Figure 3b).

2.3.3. Red Band Thresholding

After clipping the imagery using the land–water boundary, only the sea water area was retained for TMZ mapping (Figure 3c). We tested random forest (RF) classifiers [38] to carry out a supervised classification of TMZ owning to its robustness in ground feature mapping [39].

The result showed that band 4 (red) has the highest importance in the process of RF classification (Figure S1b), meaning that it could act as a good proxy in TMZ delineating. The red band image also shares a similar pattern (especially the extent) with the TMZ within the sea area (Figure 3c). In addition, the red band has shown strong reliability and sensitivity in Case 2 water-turbidity-related mapping in previous studies [23,40,41,42]. Therefore, the red band was used in this study. Furthermore, to delineate the TMZ, proper thresholds were needed for different Landsat sensors.

Based on 15 random-selected and land-masked clear images (5 TM, 5 ETM+ and 5 OLI), we manually drew multiple polygon samples which contained only TMZ or other water area on the red band images. As a result, a total of nearly 4 × 10⁷ pixels were selected to build the histogram in Figure 4. We found that the threshold 1400 could effectively identify TMZ from other water area on TM and ETM+, and the threshold 1300 could be applied on OLI (Figure 4). Thus, the TMZ was defined as red band >1400 on TM and ETM+ or red band >1300 on OLI.

$$\mathrm{TMZ}=\{\begin{array}{c}\mathrm{red} \mathrm{band}>1400; \mathrm{TM} \mathrm{and} \mathrm{ETM}+\\ \mathrm{red} \mathrm{band}>1300; \mathrm{OLI} \end{array}$$

(1)

The difference in thresholds was attributed to the different wavelength range of red bands on different sensors (0.63–0.69 μm for TM and ETM+, and 0.636–0.673 μm for OLI).

2.3.4. Accuracy Assessment

Firstly, all the derived 322 TMZ maps from the period 1984–2021 were checked using visual inspection. In general, the TMZ maps could reflect the real condition of TMZ in the Yellow River estuary very well despite its drastic variations during such a long time period.

We randomly selected 30 images to be applied to quantitatively accuracy assessments. These images involved various TMZ patterns, covered four seasons and roughly the whole study time period and included sensors of TM, ETM+ and OLI. The 30 images provided wide representation and could be used to reflect the accuracy of the whole. In total, 3552 validation samples were selected based on these images and divided into two categories, that is, ‘TMZ’ and ‘other’, with the number of 1746 and 1806, respectively (Figure 5). The confusion matrix of each image and entire samples was calculated (

Table 1. Confusion matrix of the whole validation samples. OA, UA and PA represent overall accuracy, user accuracy and producer accuracy, respectively.

		Reference
		TMZ	Other	UA
Classified	TMZ	1720	69	96.1%
	Other	26	1739	98.5%
PA		98.5%	96.2%
OA				97.4%

and Table S1).

3. Results

3.1. Outcome of Evaluation

Through visual inspection of the full time series TMZ maps, our algorithm has high practicability and robustness in capturing the extent of TMZ (30 examples were displayed in Figure 5). Furthermore, all the overall accuracy (OA) of the random selected images have achieved more than 91%, most of which were higher than 95% (Table S1). In a holistic view, the integral OA of the complete 3552 samples was 97.4%, accompanied with high user accuracy (UA) and producer accuracy (PA) of ‘TMZ’ and ‘other’ categories (Table 1).

3.2. Time Series TMZ Maps

Time series TMZ extent results were analyzed in two aspects: area (Figure 6) and spatial pattern (Figure 7). For the spatial pattern, all derived TMZ maps in this study were composited based on seasons and years on GEE, as shown in Figure 7.

Notably, we did not stretch the 25 legends (the color bar) in Figure 7 to the same range. Theoretically, the range of TMZ occurrences during the same time spans should be approximately equivalent to each other. However, higher cloud cover in the wet season and the dynamic running status of Landsat constellation make it difficult to acquire a similar number of clear images during the same time spans (e.g., Figure 7u,v). Furthermore, the true TMZ occurrence pattern in one period does not necessarily depend on the absolute number of TMZ maps. In other words, due to the randomness of clear images, thin TMZ maps stacks could still generate a composite that reveals the true relative pattern of TMZ occurrences during that corresponding period. In addition, if all the ranges of legends in Figure 7 were adjusted to the same scale, then small-dataset composites (e.g., Figure 7a–e) may only exhibit one color, and their relative TMZ occurrence patterns would be completely phased out. Therefore, to better highlight the relative TMZ occurrence pattern during every selected period, the range of legends in Figure 7 was not scaled to the same but suited the actual ranges of every composite image instead.

3.2.1. Variation of Area

The area of TMZ in the Yellow River estuary has varied not only abruptly but cyclically from 1984 to 2021. The largest area of TMZ could reach 6360 km² on 1 March 1992, while the smallest only covered an area of 1 km² on Jul 16, 2001 (Figure 6). The amplitude of TMZ variation has also experienced evolution. Based on the amplitude, we divided the study period into four decadal periods: 1984–1989, 1990–1999, 2000–2009 and 2010–2021.

In general, the decadal average TMZ area presents an increasing momentum during 1984–2009 but declines sharply from 2010. From 2000 to 2009, the average TMZ area is 1785 km², which ranked the highest among the four periods. The period 1990–1999 has an average TMZ area of 1682 km², which is slightly lower than that of the period 2000–2009 and significantly higher than that of the period 1984–1989 (1011 km²) and 2010–2021 (1117 km²).

3.2.2. Decadal Variation

Figure 7y is a composite of all 322 TMZ maps, representing the overall TMZ spatial distribution. It reveals that TMZ occurred most frequently in an ovate area south of the Qing 8 course from 1984 to 2021. Furthermore, based on the four decadal periods mentioned above, another four TMZ occurrences images were composited (Figure 7e,j,o,t). The number of TMZ maps used in the TMZ map stacks to generate Figure 7e,j,o,t was 28, 73, 115 and 106, respectively. Compared with the other three periods, 28 maps of the period 1984–1989 were relatively sparse. This is caused by two factors. First, the period 1984–1989 is shorter than other periods. Additionally, during the period 1984 to 1989, Landsat 4 and Landsat 5 were in track at the same time. However, Landsat 4 TM just provided two images for this study (Figure 2) owing to its constant device failures.

In the period 1984–1989, TMZ mainly occurred along the coast of the YRD, and its extent was mostly limited to a narrow belt near the coast (Figure 7e). For the following 10 years, the TMZ significantly expanded its extent south of the Qingshuigou course (Figure 7j). Compared with the period 1990–1999, the TMZ in the period 2000–2009 increased in the southern part of the Yellow River mouth except for the landward intrusion, while it also extended further seaward following the Qing 8 course (Figure 7o). In the period 2010–2021, the TMZ appeared more frequently around the current Yellow River mouth, assembling in an ovate area surrounding the mouth of the Yellow River (Figure 7t).

3.2.3. Seasonal Variation

TMZ in the Yellow River estuary indicated clear seasonal variations. In both four decadal and overall periods, the maximum range and high occurrence area of TMZ in summer were the smallest (Figure 7b,g,l,q,v), with an average TMZ area of 364 km² (Figure 8a). These two aspects of TMZ grew larger in spring (Figure 7a,f,k,p,u) and autumn (Figure 7c,h,m,r,w), with spring TMZ showing higher coverage than autumn TMZ (1357 km² to 814 km² in Figure 8a). Meanwhile, winter TMZ showed the largest pattern among four seasons in all time spans (Figure 7d,i,n,s,x), with an average TMZ area of up to 2730 km² (Figure 8a). In general, the TMZ in the dry season was more active and stronger than that in the wet season. In addition, the spectrum analysis revealed the prominent one-year period (370 days) of time series TMZ area in Figure 8b, which confirmed the statistical significance of the TMZ seasonal dynamics.

4. Discussion

4.1. Water Discharge and Sediment Load of Yellow River

The riverine output of the Yellow River is an important source of sediment in the Laizhou Bay. Water discharge and sediment load of the Yellow River has shown decreasing trends during the period 1984 to 2002 (Figure 9a). In 2002, the water-sediment regulation scheme (WSRS) was initiated by the Yellow River Conservancy Committee (YWCC). After the WSRS, the variation tendencies of water and sediment began to depart from each other [43]. From 2002 to 2021, the water discharge remained relatively stable in general, while the sediment load had a downward trend (Figure 9a).

However, compared with water discharge and sediment load, the variation of yearly averaged TMZ area displayed quite a different pattern. Basically, the yearly TMZ area has experienced an increasing–stable–decreasing trend (Figure 6 and Figure 9b). The correlation coefficient for TMZ and water and TMZ and sediment was −0.27 and −0.15, respectively. It means that no robust correlation for TMZ and riverine output existed on the yearly and decadal time scale.

In addition, most of the Yellow River sediment deposited quickly in the region less than 15 km from the river mouth, limiting its direct impact on the estuary [21,28]. Furthermore, the seasonal trends of the riverine output of the Yellow River and the TMZ were not consistent. The water discharge and sediment load from the Yellow River to the sea mostly centered in summer and autumn both before and after the start of WSRS in 2002, yet the TMZ composite results indicated a much smaller TMZ in summer and autumn compared with that in spring and winter. In winter, the water discharge and sediment load of the Yellow River into the sea have not changed distinctly during recent decades [21], yet the area of TMZ has markedly reduced from the period 1990–1999 to the period 2010–2021 (Figure 6 and Figure 7i,n,s).

Before 2000, the huge water extraction and low precipitation resulted in the drying up of the Yellow River. Specifically, starting Feb 7, the no-flow days of the Yellow River even reached 226 days in 1997, causing extremely low water and sediment output to the sea (Figure 9). However, the TMZ still showed high vitality with relatively large extent even during the no-flow period in 1997 (Figure 6 and Figure 9).

The impact of the fluvial discharge of the Yellow River on the YRD and estuary is complex and systemic. The fluvial discharge of the Yellow River into the sea is the fundamental source of the suspended sediment in the Yellow River estuary, offering the prerequisite of the formation of TMZ. However, based on all the above-mentioned evidence, it can be concluded that the TMZ variation in the Yellow River estuary is not directedly related to the water discharge and sediment load of the Yellow River.

4.2. Morphology Evolution of Yellow River Mouth

In 1984, the river mouth was heading east–southeast (Figure 10a). Since 1984 to 1996, this lobe has prograded seaward, and the river mouth direction has turned to the southeast (Figure 10b). After 1996, the Qing 8 delta lobe began to develop, and the river mouth started to head northeast because of the artificial channel diversion (Figure 10c). In the meantime, the abandoned Qingshuigou lobe was beginning to retreat mainly because of sediment starvation [25]. Around 2007, bank bursting started to develop a new north heading branch despite the fact that the river mouth was still northeastward (Figure 10d). Those small branches have finally evolved into a northward river mouth by 2010 (Figure 10e). From 2010 to 2020, the northward river mouth has remained relatively stable except for the formation of a sub channel heading east (Figure 10f).

From the period 1984–1989 to the period 1990–1999, the intrusion of the Qingshuigou lobe was accompanied by the seaward retreat of TMZ (Figure 7e,j). From 2000 to 2021, the development of Qing 8 and the gradually north-turning of the river mouth resulted in the northward shifting of TMZ (Figure 7o,t). The TMZ extent south of the Qingshuigou lobe also decreased in general. TMZ was more inclined to appear in the shape of an oval (Figure 7t). Therefore, a morphology evolution of the Yellow River mouth has a clear impact on the TMZ pattern to a large extent.

4.3. Currents

In the long term, currents are the major driver of sediment transport in Laizhou Bay [18]. In recent years, the coastal dynamics have barely changed near the Yellow River mouth [44]. Model results and previous research have shown that coastal circulation off the estuary were normally directed south and kept a clockwise circulation south of the Yellow River mouth [15,18,45]. As for the tidal currents, they are low-velocity reciprocating flows transported northward during ebb tides and southward during flood tides, forming a swirling water body south of the Yellow River mouth [28].

However, the impact of tidal currents has been, on average, minor in the long term [18]. In addition, the tidal shear front, located near the river mouth, acted as a barrier to offshore transport of fluvial sediment [46]. The highest velocity of currents was found near the Yellow River mouth, and the velocity decreased southward and northward [18]. As a result, a large amount of riverine sediment was carried south and gathered around the south of the river mouth, while a minority has been transported northwest [15,18,19]. Therefore, it explained the fact that, in the long term, TMZ in the Yellow River estuary tends to gather around south of the river mouth, and only a small part of TMZ moves northwest along the coast (Figure 7e,j,o,t,y).

4.4. Wind Speeds

Wind is one of the main driving forces of sediment resuspension in shallow seas such as Laizhou Bay and has the ability to reach the seabed of Laizhou Bay [29]. Previous studies have shown that decadal turbidity changes in the Laizhou Bay were related to the long-term wind speed variation trend, especially the variation in winter. The seasonality of turbidity patterns was attributed to the seasonal variation of wind speeds. Weak wind in summer was related to low SSC, while strong wind in winter was related to high SSC [18,21]. In addition, the impact of wind direction on the long-term turbidity pattern variation was relatively small [21].

To further investigate the effect of wind speeds on TMZ distribution, we explored the relationship between every single TMZ area and its synchronous wind speed during 1984 to 2021. In general, the TMZ area tended to increase as the wind speed grew. In summer, wind was weaker, corresponding to a smaller TMZ. Larger TMZ corresponded to relatively stronger wind in winter (Figure 11). However, it should be noted that the TMZ area could still vary greatly at similar wind speed, especially in the range of 3 to 7 m/s (Figure 11a). A higher wind speed did not fully guarantee a larger TMZ area on a singular basis.

We divided the results into four groups based on four seasons (Figure 11b–e). In a holistic view, wind speeds in spring and winter are higher than those in summer and autumn. Therefore, the stronger wind in spring and winter caused a more violent resuspension of sediment, leading to a larger TMZ distribution (Figure 7u–x). Previous studies also confirmed the different impacts of wind strength on the sediment resuspension during wet and dry seasons [18,21,22]. However, these studies did not pay much attention to the different response behaviors of sediment suspension to wind in four seasons. In fact, the response sensitivity of TMZ area to wind speeds varied with seasons (Figure 11b–e). Most points were circled in ellipses that excluded outliers, and the slope of major axis of the ellipse could roughly reflect the response sensitivity of TMZ area to wind speeds. For example, the TMZ area in summer barely changed even though the wind speed was increasing and the slope of the major axis of that ellipse was the smallest among four seasons (Figure 11c). On the other hand, the TMZ area expanded rapidly with the increase in wind speeds in winter, and its major axis slope was almost the largest among four seasons (Figure 11e). The major axis slope of ellipses in spring and autumn followed the similar pattern with winter and summer, respectively (Figure 11b,d).

Seawater stratification is an important oceanographical phenomena, which impedes the vertical matter exchange in the sea. In summer, the water column was strongly stratified [47], hindering the wind resuspension of seabed sediment. In winter, the water column was vertically uniform [29], allowing the effect of wind to easily penetrate through seawater to stir the underlying sediments. Spring and autumn were the transition period between summer and winter. Therefore, it indicated that the stratification condition in different seasons caused different seasonal responses of TMZ to winds (Figure 11b–e). In addition, water column stratification could partly explain the weak correlation between TMZ area and wind speed, with an overall correlation of 0.4.

4.5. Tides

We find that TMZ area can still vary largely at similar wind speeds of 4–5 m/s within the same or similar seasons (Figure 11a). We further explored the relationship between tide and TMZ area under the circumstances that they were in similar seasons and wind speeds (Figure 12).

Under the influence of similar wind speeds and seasons, TMZ area did not display a clear correlation with tide, which has sharp contrast with the condition in the Yangtze estuary where tide could exert significant influences on TMZ extent [3,48]. Regardless of whether it was in flood season (summer and autumn) or in dry season (spring and winter), the TMZ area did not change linearly or quasi-linearly with the absolute value of tide range (Figure 12a,b). At the same time, there were no significant difference between the impact of ebb and flood tide on TMZ area (Figure 12a,b). This evidence suggests that tides were not the essential factor influencing the area of TMZ on a singular basis in the Yellow River estuary.

4.6. Limitations and Uncertainties

It should be noted that the scan line corrector (SLC)-off issue of Landsat 7 and haze contamination images could lead to inaccurate TMZ mapping. In 2003, the SLC on board Landsat 7 failed, resulting in intermittent pixels missing. Although gap filling methods existed [49], these artificial filling values may have created false detections considering the dynamic spatiotemporal variability of water. These techniques should be used with caution [50]. Haze is similar to thin cloud and mostly arises from water droplets, smog, dust and other particles, diming the scene clarity. Haze images accounted for a small portion of our ready-to-map imagery stacks. It is difficult to remove haze in satellite imagery without information distortion. Therefore, to avoid missing TMZ area count, we manually delineated the outline according to TMZ extraction results for SLC-off images. For haze contamination, its similar spectral property with TMZ in the red band tended to enlarge the TMZ extent. Thus, we modified the TMZ extraction results for haze images.

Previous studies utilized cruise data to validate the time series water turbidity results from remote sensing retrieval models [12,17,18]. Cruise data are often acquired at specific locations and during short time periods, meaning that cruise data could only validate limited derived results. Like these studies, we did not display and validate all the 322 TMZ maps, which may have limited the quality of accuracy assessment. However, 30 images used in the accuracy assessment were random selected, covered a variety of conditions and achieved a high overall accuracy, meaning that these images were representative of all the results.

The sediment discharge of the Yellow River into the sea was the foundation of TMZ formation, while the TMZ was not directly related to the fluvial output of the Yellow River. Morphology, currents, wind speeds and water column stratification were the main factors regulating the behavior of TMZ in the Yellow River estuary. However, these factors still could not fully explain the behavior of TMZ. The coexistence of various interactive factors effecting sediment dynamics, e.g., biogeochemical processes, subaqueous topography, sediment grain size especially after the WSRS, and the combined effect of direction of wind and currents, makes it challenging for mechanism research of TMZ at different time scales. In addition, uncertainties intensified when trying to fully understand sediment dynamics only through the TMZ extent. Besides the extent of the TMZ, the inner variability could still bear information about sediment dynamics to some extent, especially when the TMZ was highly extensive. Thus, to better understand the TMZ in the large river estuary, details of the relationship between different factors, such as the hysteresis between discharge and SSC [51], and more statistical tools, such as the wavelet methodology [52], should be considered.

5. Conclusions

In this article, we proposed a new method to derive the long-term time series TMZ area based on time series Landsat images and GEE. This method has been successfully applied to the Yellow River estuary to generate a 38-year TMZ product. Visual and quantitative accuracy assessment show that our algorithm is robust and has great potential in large-scale TMZ mapping of global estuaries.

TMZ varies drastically in different time scales. In decadal time scale, the TMZ area shows an increasing trend from 1984 to 2009 but has sharply decreased since 2010. In general, as for seasonal variation, TMZ extent is largest in winter and smallest in summer.

Possible driving forces of TMZ variation were investigated. The morphology of the Yellow River mouth shifts the position of the TMZ. Currents aggregate the TMZ in an ovate area south of the Yellow River estuary. In terms of the Yellow River fluvial output, the water discharge and sediment load show negligible direct influences on the TMZ. The seasonal pattern of TMZ is caused by the different strength of wind in different seasons. However, on a singular basis, higher wind speeds do not necessarily lead to a larger TMZ area. The seasonality of stratification condition causes the different response sensitivity of TMZ area to wind speeds. In addition, we find that tides could barely affect the area of TMZ even in the case of similar wind speeds and seasons.