The historical bathymetrics of the depth profile perpendicular to the shoreline during 2000–2018 were derived from the archived dataset of the Yellow River Estuary Hydrology and Water Resources Survey Bureau. The 2020 water depth was measured using a single-beam echo-sounder on January 13–14, 2020, with tidal corrections made based on tidal elevation data from the Gudong tidal observation station (Xu et al., 2022a). The topographic map based on water depth data in 2018 is shown in Figure 1C.

The interannual variation in the water depth profile of the survey line is shown in Figure 1D . From 2000 to 2020, the water depth at the root of the seawall increased from 2 m to 5.7 m, and the depth of the profile showed a deepening trend. To analyze scour dynamics in front of the seawall, a GD bottom observation ( Figure 1B ) system was deployed near the survey line (37.9363° N, 119.0576° E). Additional information on the distribution and interannual variation of water depth near the observation station, as well as the layout of bathymetric survey lines is provided in Xu et al. (2022a). The observation station was located approximately 100 m from the shore, and all the instruments were mounted on a stainless-steel platform with a height of ~0.5 m. Observations lasted from November 4 to December 4, 2021 (31 days). The stainless-steel platform was equipped with a Nortek 6 MHz Acoustic Doppler Velocimetry (ADV) to measure near-bed (0.8 m) three-dimensional flow velocity at a sampling frequency of 16 Hz. A Nortek 600 kHz acoustic wave and current (AWAC) was used to monitor surface waves and currents in the upper water column. Optical backscatter sensors (OBS, Seapoint Sensors, Inc), powered by AWAC were installed 0.4 m above the seabed. The current and OBS measurements were configured to simultaneously sample every 10 min, while wave measurements were set to sample hourly.

From January 14 to April 19, 2020, we conducted in-situ observations at GD (37.9359° N, 119.0575° E) and KD (37.9200° N, 119.1379° E) simultaneously to compare the discrepancies of hydrodynamics in front of the seawall and the far field ( Figure 1B ). An AWAC were deployed at each station with the same parameter settings as the field observations from 2021, and no additional instruments were deployed.

Wind data at 10 m above the sea surface (W ₁₀, unit: m/s) near the observation station were obtained from National Centers for Environmental Prediction climate forecast system version 2 (NCEP CFSv2) at 1 h intervals (URI: https://rda.ucar.edu/datasets/ds094.1/). Owing to the difficulty of collecting water samples during high winds, we collected bottom sediment near the observation station. Meanwhile, OBS sensor calibration was performed in the laboratory to establish the relationship between OBS measured turbidity values and SSC. Table 1 summarizes the instruments deployed on the stainless-steel platform.

The type of bottom sediment type near the site is silty (d ₅₀ = 0.068 mm, clay content CC = 2.85%) and the grain size accumulation curve for surface sediment near the observation site is shown in Figure 2 .

The accuracy of AWAC wave measurement benefits from a combination of acoustic surface tracking, velocity data, and pressure data. Acoustic surface tracking uses a vertical acoustic beam to measure water surface height. This has the advantage of reducing the impact of velocity and water depth compared to pressure sensors (Pedersen and Lohrmann, 2004). The surface wave pressure spectrum Sp(f) and directional wave spectrum were calculated using Storm software—a Nortek software for analyzing wave data measured by AWAC. Wave parameters, namely, significant wave height (Hs) and bottom wave orbital velocity (U_w ) were expressed via Equations (1) and (2), respectively (Wiberg and Sherwood, 2008)

Where w represents the angular wave frequency, Δf is the frequency band of wave pressure spectrum, i is the number of bands, and k represents the wave number. The peak period (T_p ) is defined as the period corresponding to the highest energy of Sp(f).

Given that near-bed velocity measured by ADV is susceptible to environmental interference (Fugate and Friedrichs, 2002), strict post-processing was performed before analysis. The quality of ADV data was controlled by removing ADV-measured data points with a correlation< 70% and Signal-to-Noise ratio (SNR, unit: dB)< 5 dB. The “phase space method” was then used to detect and replace spikes in the data (Goring and Nikora, 2002).

High-frequency flow velocity measured via ADV captured changes in mean flow velocity, wave motion, and turbulence. Accurate extraction of turbulence and calculation of turbulent parameters, such as turbulent kinetic energy (TKE) and turbulence-induced sediment vertical diffusion flux (TSF) requires decomposition of these motions from the original velocity signal. Bian et al. (2018) evaluated the existing wave–turbulence decomposition methods and determined the advantages of the SWT method (Daubechies et al., 2011; Thakur et al., 2013; Bian et al., 2018). Taking the eastward velocity (u) measured via ADV as an example, it can be broken down into mean velocity ( u ¯ ) and velocity fluctuation (u ^′) . In the wave environment, (u ^′) is the linear superposition of further wave orbital velocity ( u w ′ ) and turbulent fluctuation ( u t ′ ) as follows:

Figure 3 shows the results of wave–turbulence decomposition (Burst = 28, Hs = 2.4 m). The SWT of u ^′ ( Figure 3A ) and power spectrum analysis ( Figure 3C ) shows that u ^′ is significantly affected by waves. At frequencies > 0.5 Hz, the power spectrum of u ^′ shows a distinct inertial subrange characterized by the Kolmogorov’s “−5/3” theoretical spectrum (Bian et al., 2018; Fan et al., 2019). SWT was used to reconstruct the velocity fluctuation component at wave-dominant frequencies ( Figure 3A ) to obtain u w ′ ( Figure 3B ). Meanwhile, the remainder was assigned to the turbulent component u t ′ ( Figure 3B ). Therefore, the power spectrum of each component after wave–turbulence decomposition indicates that the components of the wave in the original velocity fluctuation were effectively removed. Meanwhile, the turbulence spectrum exhibited a pronounced “energy groove” in the wave-dominant frequencies ( Figure 3C ). However, the frequency range affected by the wave is relatively narrow, meaning that turbulence energy loss at the wave-dominant frequencies was negligible (Bian et al., 2018).

Bottom bed shear stress is a critical parameter controlling erosion, deposition, and resuspension of seabed sediments. We calculated the bottom bed shear stress induced by currents, waves, and TKE, expressed as τ_c , τ_w , and τ_TKE , respectively.

Current-induced bottom shear stress (τ_c ,unit: N/m²) was calculated using the formula reported by Soulsby (1997), which relies on a specific log profile and estimates the shear stress from the first moment statistics. The bed shear stress τ_c is expressed as follows:

where U (unit: m/s) is the burst-averaged velocity at height z (z = 0.8 m) collected by ADV, and z ₀ = d ₅₀/12 (unit: μm) indicates the roughness of the seabed (d ₅₀is the median particle size, 68 μm). ρ = 1025 kg/m³ and κ = 0.4 represent the seawater density and von Karman constant, respectively.

Wave-induced bottom shear stress (τ_w , unit: N/m²) can be expressed as a function of the wave friction coefficient (f_w ) and wave orbital velocity (U_w ) (Grant and Madsen, 1979; Soulsby, 1997). For a wave with period T, and orbital velocity U_w , the bed shear stress τ_w is expressed as follows:

Combined shear stress of currents and waves play an important role in seabed erosion in estuaries and coastal areas (Zhu et al., 2016). Bed shear stress caused by combined wave–current action (τ_cw , unit: N/m²) is described by (Grant and Madsen, 1979) as follows:

where φ is the angle between the current direction and wave propagation direction.

Along with currents and waves, turbulence is also critical for resuspended sediment to enter the water column. Turbulence strength determines the vertical diffusion of suspended sediment. The turbulence-induced bottom shear stress τ_TKE (unit: N/m²) is expressed as follows:

where C is a constant (0.19) (Stapleton and Huntley, 1995; Kim et al., 2000), and TKE is calculated from the turbulence components described in Eq (3).

Spectral and wavelet analysis are widely used in signal processing and analysis (Elsayed, 2008; Pomeroy et al., 2015). They were used in the analysis to process flow velocity and SSC data. Spectral analysis can be used to determine the frequency range of fluid motion. Wavelet analysis allows the extraction of time and frequency information to determine the intensity of each frequency component at each time of measurement. That is, the instantaneous distribution of each signal in the frequency domain (Liu and Babanin, 2004). In the present study, wavelet power spectrum was applied to turbulence-induced sediment vertical diffusion flux and momentum flux to visualize the coherent structures, which is shown in Section 5.2.

Variation of SSC in the boundary layer can be obtained by combining turbidity measured using optical instruments, acoustic reflection signals from acoustic instruments, and field water samples (Yuan et al., 2008; Yuan et al., 2009; Zhao et al., 2016; Fan et al., 2019). A similar approach has been used in this study. The in-situ near-bottom water samples obtained were dried in an oven at a constant temperature of 60°C. The OBS turbidimeter, and a certain volume of distilled water, were placed in the calibration tank with a volume of 30 L, where turbidity and SSC were near 0. The dried mud samples were then gradually added to the tank and stirred for 5 min to ensure complete mixing of the suspended sediment, and the OBS recorded the turbidity simultaneously. A total of 29 sets of corresponding SSC and turbidity values were obtained. Correlation was determined using the regression method and the coefficient of determination (R²) reached 0.9928 ( Figure 4 ). According to the calibration formula obtained, the turbidity value recorded by OBS was converted to SSC (denoted as SSC_OBS, unit: mg/L).

The SNR recorded by ADV can reflect SSC variation (Fugate and Friedrichs, 2002; Voulgaris and Meyers, 2004). We established linear regression between the logarithm of SSC_OBS and SNR value, as shown in Figure 4B . Results showed that log₁₀(SSC_OBS) was significantly correlated with SNR with an R² of 0.8689. The SNR was then transformed into high frequency SSC (denoted as c, unit: mg/L). The burst averaged value of c, was denoted as SSC_ADV (unit: mg/L).

During the observation periods, the maximum and average values of W ₁₀ were 18.0 and 7.3 m/s, respectively, with a dominant northwest wind ( Figure 5A ). Three weather events with relatively strong winds (W ₁₀ > 10.0 m/s) occurred on November 6–11, November 21–22, and November 30–December 1, 2021 ( Figure 5 ).

At the beginning of the first strong wind weather event, the northeast wind gradually strengthened, with a W ₁₀ range of 2.2–18.0 m/s. Hs rapidly increased to 3.19 m, reaching the maximum value during the observation period. From November 7 to 8, the wind direction changed to the northwest, while the W ₁₀ continued to range between 10.0 and 11.0 m/s. The Hs value rapidly decreased to approximately 0.5 m, and the Tp was stable at 4–6 s. At the beginning of the second weather event, a moderate intensity northeast wind was observed (W ₁₀ range: 2.4–11.9 m/s) with a rapid rise in Hs (2.27 m). On November 21, the wind direction shifted to northwest and wind speed increased to 16.4 m/s. However, the Hs was lower than that at the beginning of the second weather event. Wave growth was then limited by the transient northwest wind. In the third weather event, W ₁₀ had a range of 2.4–13.5 m/s, Hs increased from 0.26 m to 1.8 m, and Tp was stable at 5–7 s. The changes in Hs and Tp maintained a relatively high level of consistency, and increased rapidly under the influence of the northeast wind and decreased slowly under the northwest wind. Under the south or northwest wind, Hs remained low, while the high wave height condition was caused by northeast winds.

The time series for water level and flow velocity during the observation period are shown in Figure 6 . Hs was used to reflect wave intensity and is shown in Figure 6A . The measured water level (η, unit: m) ranged from −1.3 to 1.4 m ( Figure 6B ). Its variation notably increased during high waves. Harmonic analysis on the time series of η using the T_tide package (Pawlowicz et al., 2002) showed that the tidal-induced water level (η_tidal , unit: m, Figure 6B ) was consistent with the variation in the water level measured under normal weather conditions. During high waves, variation in the residual water level (η–η _tidal) caused by wind stress is notably greater than η_tidal (−0.5 to 0.6 m), while the maximum residual water level reaches 1.28 m. A comparison of residual water level and wind conditions shows that the rise in water level was caused by offshore winds (NE direction), while the fall in water level was associated with onshore winds (NW direction).

In calm weather conditions, the measured flow velocity (U) was< 0.3 m/s, which is nearly equivalent to the tidal current (U_tidal )value, the maximum of which was 0.28 m/s ( Figure 6C ). During periods of high wave activity, the flow velocity was significantly enhanced, reaching a maximum of 1.11 m/s. The enhanced velocity was primarily attributed to residual current generated by high waves. The flow direction also changed notably, indicating that strong longshore currents exceeded the reciprocating tidal currents. The average azimuth of the flow direction observed during high waves was 133°. The azimuth is the northward reference angle, and clockwise direction is denoted as positive. This is predominantly related to the orientation of the Gudong seawall located in the NW–SE direction with an average azimuth of 120° ( Figure 6D ).

The SSC_OBS and SSC_ADV obtained using optical and acoustic instruments after calibration was shown in Figure 7 . The SSC exhibited robust coupling with Hs. In calm weather conditions, the average SSC_OBS and SSC_ADV values were approximately 300 mg/L. When wave conditions strengthened, the trends in SSC_OBS and SSC_ADV were relatively consistent, and rapidly increased to 7,137 mg/L and 12,222 mg/L, respectively. During the gale from November 7–10, the SSC_OBS value was out of range during high waves with Hs > 2.4 m, and the SSC_OBS and SSC_ADV values differed due to distinct behavior of coarse and fine particles during the decay stage of the gale. Coarse particles settled while the fine particles may have remained in suspension when Hs dropped below 0.5 m. Considering that the SSC obtained by ADV acoustic inversion is not sensitive to the dependence on particle size variation (Fugate and Friedrichs, 2002), the optical signal was more susceptible to the influence of unsettled fine particles. Therefore, SSC_OBS may be higher than SSC_ADV during the decay stage of the gale. The intra-burst standard deviation of high frequency SSC was typically< 1,000 mg/L, while a previous study reported that the quality of acoustic backscattering was superior to optical backscattering (MacVean and Lacy, 2014). Therefore, SSC_ADV was used to characterize the SSC of the bottom boundary layer.

Formation of high SSC is predominantly attributed to horizontal convective transport or the local resuspension process (Zhang et al., 2021). Horizontal convective transport is induced by the horizontal SSC gradient, while the local resuspension process is primarily induced by turbulence in the form of vertical sediment diffusion (Yuan et al., 2009). Band-pass filtering and SWT have been applied to extract the wave ( c w ′ ) and turbulent ( c t ′ ) components of high-frequency SSC fluctuation (c ^′ , c ′ = c − c ¯ , the overbars (^—) represents the average of each burst) (Fan et al., 2019; Li et al., 2022). The power spectrum of c ^′ exhibited a clear peak near the wave-dominant frequency, which was in line with velocity fluctuations. However, at frequencies > 3 Hz the power spectrum appeared to be contaminated by white noise ( Figure 8 ). Therefore, in this study, the SWT method was applied for wave-turbulence decomposition of c ^′ TSF is a more direct indicator of resuspension dynamics than SSC (Brand et al., 2010). It can be expressed as c t ′ w t ′ ¯ , where w t ′ represents the vertical velocity fluctuation extracted using the SWT method. SSC_ADV and TSF exhibited the same variation trend. High SSC_ADV is typically accompanied by high TSF values ( Figure 9 ). This indicates that high SSC_ADV is mainly generated by local resuspension during gale events. Therefore, SSC_ADV was applied to characterize scour intensity in the study area with a focus on sediment resuspension.

To analyze the principal factors controlling sediment resuspension, we calculated the bottom shear stresses generated by currents, waves, turbulence, and the combined action of waves and currents. In calm weather conditions, τ_c , τ_w , τ_cw and τ_TKE were< 0.1 N/m². Following the increase in wave energy, τ_c , τ_w , τ_cw and τ_TKE increased by 1–2 orders of magnitude reaching maximum values of 1.41 N/m², 5.48 N/m², 5.68 N/m² and 4.44 N/m² on November 7, 2021, respectively ( Figure 9 ). The maximum τ_w value was 3.89 times that of τ_c , and τ_w dominated the change in τ_cw . Therefore, wave and current acted together in the study area, among which the effect of wave-induced bottom shear stress on the seabed exceeded that of the current.

To discuss the effects of seawalls on hydrodynamic conditions, we analyzed the waves and currents at GD and KD stations that were simultaneously observed in 2020. The time series of Hs showed similar variation patterns at the two observation sites ( Figure 10A ), although the statistical wave condition of KD was slightly stronger than that of GD, with the difference generally fluctuating around 0.4 m. The strong waves at KD were mainly from the NE ( Figure 10B ), with a mode azimuth of 32˚. However, the mode azimuth for strong waves at GD was 12˚. This was mainly due to wave refraction in front of the Gudong seawall. Furthermore, the flow velocity of GD was stronger and more sensitive to Hs than KD ( Figure 10C ). During periods with high waves, the flow direction at GD changed substantially with an average azimuth angle of 71° ( Figure 10D ). The flow direction were different from the observation results from 2021, which may have been related to the stations set up in 2020 being closer to the shore and the stronger influence of wave reflection on the hydrodynamic conditions.

To elucidate net transportation of sediment in the study area, the horizontal net flux (HNF, unit: t/m²) of sediment in the bottom boundary layer was calculated. The horizontal velocity component of the bottom boundary layer observed using ADV was projected to the coastal direction. The southeast coastal transportation was regarded as positive, while the northwest coastal transportation was negative. The horizontal net flux was then calculated, with integration of SSC_ADV and the longshore currents vector being performed in the diurnal tidal cycle.

In the small wave stages, local resuspension of sediment was relatively weak and the horizontal net flux was low (maximum< 0.3 t/m²) ( Figure 11B ). At this time, the net sediment transport level was in relatively dynamic balance. In contrast, during significant wave effects, most suspended sediments were transported to the southeast ( Figure 11A ), and the transportation flux of resuspending sediment was high, with horizontal net flux reaching a maximum of 769 t/m² per day ( Figure 11B ). These results highlight the significant influence of wave-induced longshore currents on transporting suspended sediment. Specifically, copious amounts of local resuspended sediments were transported southeast from the coastline of the study area by longshore currents, and local seabed erosion persisted during gales.

Sediment resuspension is a dynamic process frequently occurring in the bottom boundary layer that drives material exchange between the seabed and seawater. Sediments are mainly resuspended by waves and currents and the bottom shear stress is an key parameter to control the sediment resuspension (Zhu et al., 2016; Niu et al., 2020; Li et al., 2022). The relationship between SSC and bottom shear stress is shown in Figure 12 . For all calm weather periods (Hs< 0.2 m), a linear regression was observed between SSC and τ_c , with R²< 0.1 ( Figure 12A ), indicating that tidal shear stress is not sufficient to resuspend the bottom sediment. For the wave-influenced episodes with Hs > 0.2 m, the R² between SSC and wave-induced bed shear stress τ_w was 0.68 ( Figure 12B ). During wave periods, τ_w has a larger contribution in τ_cw than in τ_c ( Figure 9A , Figure 12 ) and SSC has a similar R² to τ_cw and τ_w . This demonstrates that wave-induced shear stress is the main contributor to total bed shear stress during large gales. This means that waves generate sufficient shear stress to resuspend sediment, which is consistent with previous studies with remote sensing (Chu et al., 2006; Zhang et al., 2018; Li et al., 2021), field investigations (Yang et al., 2011; Liu et al., 2020; Niu et al., 2020; Zhang et al., 2021), and numerical simulations (Jiang et al., 2004; Fan et al., 2020).

In the present study, τ_TKE was found to be notably enhanced during high waves. Therefore, to further characterize the effect of waves on SSC, the influence of turbulence on SSC changes was further examined. During high shear stress induced by wind and waves, SSC notably increased, with a maximum value of 12,222 mg/L. Enhanced SSC and seabed erodibility during wave periods was reflected in the response of TSF to Reynolds stress ( U t ′ w t ′ ¯ , where   U t ′ and w t ′ represent horizontal and vertical turbulence fluctuations, respectively) (MacVean and Lacy, 2014). We regarded an Hs of 0.2 m as being the boundary between calm weather and periods with high waves and grouped all the data according to Reynolds stress values with step sizes of 0.0004 N/m². TSF and Reynolds stress generally show strong coupling, with low values in calm weather periods ( Figure 13 ). However, TSF decreases when Reynolds stress is > 0.004 N/m². There are two likely explanations for this phenomenon. (1) Turbulence can cause vertical distribution of suspended sediments to become more uniform by enhancing the vertical diffusion, thereby decreasing the bottom sediment concentration and TSF in windy conditions. (2) Wave-induced pore pressure causes liquefaction in the seabed (Xu et al., 2022b), thereby reducing turbulence intensity. This is in line with that reported by Li et al. (2022), who indicated that SSC increased first and subsequently decreased in a typhoon event.

CWT has been used to visualize the intermittency of instantaneous Reynolds stress and TSF effectively (Torrence and Compo, 1998; Salmond, 2005; Keylock, 2007; Yuan et al., 2009; Li et al., 2022). Wavelet analysis results of instantaneous Reynolds stress and TSF signals during a gale (Burst = 28, Hs = 2.4 m) are shown in Figure 14 , where the wavelet analysis passed the confidence test (95% confidence interval). The energy of Reynolds stress and TSF in the wave frequency range are relatively low. This suggests that the wave energy in the original signal was filtered out via the SWT method. Distinct plume strip structures with distribution regularity were observed between 0.2 and 8 Hz. These plume structures contained the strongest TKE, which contributed most to Reynolds stress and TSF. The high-energy bands in the Reynolds stress and TSF wavelet coefficients are consistent over time, which highlights the significant impact of turbulence on sediment resuspension.