The tidal flats along the Zhejiang coast face directly toward the East China Sea, making them highly susceptible to the combined influence of macrotidal currents and significant wave and swell activity, which contribute to vigorous sediment dynamics. These flats provide an ideal setting for investigating the impact of tidal and wave forces on sediment transport processes. The field study was conducted on the Oufei intertidal mudflat, located between the Oujiang River estuary and the Feiyun River estuary. The surface sediment composition ranges from 3.6% to 5.0% sand, 76.2% to 77.7% silt, and 17.3% to 20.2% clay, with a median grain size (D₅₀) of 14 μm. The mudflat extends 4 to 6 kilometers seaward, exhibiting a gentle slope that varies from 1° in the upper intertidal region to 10° at the outer edge. The intertidal flat is devoid of vegetation and lacks a prominent tidal creek system or significant gullies.

Coastal hydrodynamics in the study area are primarily influenced by a semidiurnal macrotidal regime, with an average tidal range of approximately 4 meters (Wang and Ke, 1997; Xing et al., 2012). Due to the progressively shallowing water depths towards the coastline and the river discharge from the Oujiang River, the M₄ tidal component becomes prominent in the nearshore zone, indicating the development of a flood-dominated tidal wave. Swell and surface waves originating from the East China Sea can propagate onto the study mudflats, which are not sheltered by offshore islands. Observations from gauge stations off the Oujiang estuary indicate that the mean wave period ranges from 4 to 10 seconds, with significant wave heights varying between 0.5 meters and 5 meters under typical weather and cold wave conditions (Ye et al., 2022).

The field survey was conducted from 11 to 15 March 2015, spanning six semidiurnal tidal cycles. A variety of time-series data, including water depth, wave level, turbidity, and near-bed velocity, were collected using instruments affixed to a tripod ( Figure 1d ). The observation site was established in the southern region of a reclamation area, where the sheltered conditions eliminated the alongshore currents at the site ( Figures 1b, c ).

A 6 MHz Nortek Acoustic Doppler Velocimeter (ADV) was employed to measure high-frequency near-bed velocity and water level variations. The ADV probe was mounted on a tripod, positioned 0.25 m above the seabed, and directed downward at a standoff distance of 0.15 m from the probe head. Consequently, the sample volume of the ADV transducer was located 0.1 m above the bottom (mab). The ADV operated in burst mode at 10-minute intervals, with each burst recording 1000 samples at a frequency of 8 Hz. Velocity profiles were acquired using a downward-looking 2 MHz High-Resolution Acoustic Doppler Current Profiler (HR-ADCP). The HR-ADCP was positioned 0.8 m above the seabed, with velocity data collected in 3 cm bins at a sampling frequency of 1 Hz.

Water turbidity was measured at four vertical layers positioned at 0.1 m, 0.2 m, 0.3 m, and 0.6 m above the seabed using Optical Backscatter Sensors (OBS-3A, self-recording turbidity–temperature monitoring instruments), with a sampling period of 2.5 minutes ( Figure 1d ). Water samples were simultaneously collected using four custom-made samplers (600 mL volume) placed at the same heights as the OBS-3A sensors to calibrate turbidity measurements, thereby enabling the determination of suspended sediment concentrations.

In order to characterize the tidal currents in the study area, long-term ADCP data were gathered from four mooring stations (C1–C4) along the Zhejiang coast ( Figure 1a ). The ADCP measurements were conducted between December 2008 and March 2009, along a transect perpendicular to the isobath, covering water depths from approximately 20 m to 70 m. The velocity profile was sampled every 20 minutes, with bin sizes ranging from 0.5 m to 2 m.

The significant wave height, Hs, can be estimated from the high-frequency fluctuations in water level recorded by the near-bed ADV. The high-frequency fluctuations of the water level are the difference between the observed water level and the mean water level during a burst. Under the assumption that surface wave heights generally adhere to a Rayleigh distribution, as proposed by Longuet-Higgins (1952), the value of Hs is calculated following the methodology outlined by Wiberg and Sherwood (2008).

where S_η is the spectral density of surface elevation as a function of frequency, f (=1/T).

The mean wave period T is estimated from the spectra of surface elevation as follows (Equation 2):

where the symbols have the same meaning as in Equation 1.

Each measured instantaneous velocity (u, v, w) can be considered as the summation of mean flow velocity ( u ¯ , v ¯ , w ¯ ), wave oscillatory ( u ˜ , v ˜ , w ˜ ), and turbulent (u’, v’, w’) components, indicating (Equation 3)

In general, the turbulent energy spectra within the bottom boundary layer of intertidal mudflats exhibit characteristics consistent with Kolmogorov’s −5/3 law in the inertial subrange. However, wave-induced motion can significantly distort the original turbulent energy spectra in the energy-containing subrange. When wave oscillations are superimposed on turbulent components, the resulting turbulent parameters, such as Reynolds stress and turbulent kinetic energy, may be substantially overestimated.

Several methods have been proposed for separating turbulence and wave-induced fluctuations, including the paired differencing method (Trowbridge, 1998; Shaw and Trowbridge, 2001; Feddersen and Williams, 2007), the cohesive spectral method (Benilov and Filyushkin, 1970; Bendat and Piersol, 2000), and the empirical mode decomposition (EMD) method (Qiao et al., 2016; Bian et al., 2018, 2020). In this study, we employed the cohesive spectral method (Bricker and Monismith, 2007) to identify wave motions by examining the correlations between mutually independent velocity and pressure measurements. The coherence between ADV pressure and high-frequency velocity data (p’ and ui’) was calculated as follows (Equation 4):

where C p u i is the cross spectra of p ′ and u i ' , C p u i * is the complex conjugate of C p u i , and S p p and S u i u i are power spectra of p ′ and u i ' , respectively. Then, the power spectra of wave-removed turbulence S u i u i ' were obtained following Equation 5

To assess the effectiveness of the wave-removal process, we used velocity bursts from the eastward and vertical components collected during high tide in the second tidal cycle as a case study. The wave-induced peak appears at frequencies lower than 1 Hz (represented by the gray line in Figure 2A ), suggesting that the low-frequency band is particularly susceptible to wave fluctuations. For frequencies exceeding 1 Hz, the spectral slope generally adheres to Kolmogorov’s −5/3 law, indicating that the influence of waves on the high-frequency band is negligible. Following the wave-turbulence decomposition procedure, the peaks in the low-frequency band (< 1 Hz) are largely eliminated, with the spectral slope in this range becoming relatively horizontal. Furthermore, the wave-enhanced spectrum exhibits a significantly lower peak in the vertical velocity component compared to the horizontal velocity component ( Figure 2B ). This finding is consistent with prior research indicating that surface waves have a less pronounced effect on vertical velocity than on horizontal velocity (Stapleton and Huntley, 1995). In accordance with Bricker and Monismith (2007), we assume that waves primarily govern the phase of the Fourier coefficients, as waves dominate velocity fluctuations at the wave crest. The time series of wave oscillations ( u ˜ ) can be reconstructed by inverse Fourier transform with the modules ( | u j ˜ | ) of the Fourier coefficients of the wave part and the original phase ( e i ∠ u j ˜ ) by u j ˜ = | u j ˜ | e i ∠ u j ˜ . The corresponding turbulent velocity can be estimated by subtracting the mean flow and wave velocity from the original velocity ( u j ' = u j − u j ¯ − u ˜ j ) .

The bottom shear stress due to waves (τw) was estimated according to van Rijn et al. (1993) (Equation 6):

where ρ_w is seawater density (= 1030 kg m^-3). The wave orbital velocity ( u ˜ b ) at the edge of the wave boundary layer is given by (Equation 7)

where ω (= 2π/T) is the angular velocity and A ˜ b is the peak value of the orbital excursion, given by (Equation 8)

where Hs is the significant wave height, h is the water depth, L is the wavelength (= ( g T 2 / 2 π ) t a n h ( 2 π h / L ) )), g is the gravitational acceleration (9.8 m s^-2), and T is the mean period.

Based on Soulsby (1995), the wave friction coefficient fw depends on the hydraulic regime and can be expressed as (Equation 9)

where R e w (= u ˜ b A ˜ b / v ) is the wave Reynolds number and r (= A ˜ b / k s ) is the relative roughness. For the hydrodynamic condition of the observed intertidal mudflat, R e w has values smaller than ~ 4×10⁴, which indicates a laminar boundary layer induced by waves over the intertidal flat (Jensen et al., 1989).

Since vertical velocity is less affected by surface waves, the bottom shear stress induced by tidal currents (τc, Pa) can be more accurately estimated using near-bed vertical turbulent velocities (Huntley and Hazen, 1988; Kim et al., 2000), as follows (Equation 10):

where C₀ is a constant of 0.9, as recommended by Kim et al. (2000). ρ is the seawater density (kg m^-3), and w " 2 ¯ is the variance in the vertical turbulent velocity (m² s^-2). w’ is the wave-removed vertical turbulent velocity.

The bed shear stress due to combined current–wave action τcw can be estimated through the van Rijn et al. (1993) model, Soulsby (1995) model, and Grant and Madsen (1979) model. All three models yield very similar results for bed shear stress. Here, we use the Grant and Madsen (1979) model to estimate the current-wave shear stress (Equation 11).

where τ_w and τ_c are the bed shear stresses due to waves and currents (N m^-2), respectively, and φ_cw is the angle between the current direction and wave direction. To determine the wave direction, we combine horizontal velocities and pressure data from the ADV using a standard PUV method (available at http://www.nortekusa.com/usa/knowledge-center/table-of-contents/waves) and some MATLAB tools in the Nortek web page that can compute wave directional spectra from data using Nortek Vector. These models have been widely employed in estimating bed shear stress due to current–wave interaction (e.g., Feddersen and Guza, 2003; Keen and Glenn, 2002; Shi et al., 2012; Styles and Glenn, 2002).

The critical shear stress ( τ c r ) is essential for evaluating bottom sediment erosion. For cohesive sediments, τ c r is primarily governed by cohesive forces resulting from electrochemical interactions between sediment particles, making it challenging to estimate precisely. In this study, we employ two methods to estimate the lower and upper bounds of τ c r . The lower band of τ c r is determined using the extended Shields curve, as proposed by Mantz (1977) and Govers (1987).

where D * ( = d 50 [ ( s − 1 ) g / v 2 ] 1 / 3 ) is the dimensionless grain size and θ_cr and τ_cr are the dimensionless critical shear stress and critical shear stress, respectively. ρ and ρ_s represent the fluid density (1030 kg m^-3) and sediment density (2650 kg m^-3), respectively. v is the kinematic viscosity coefficient (=1×10–⁶ m² s^-1). g (= 9.8 m s^-2) is the gravity acceleration. d₅₀ is the median diameter of the surface sediment and has a value of 12 μm. The upper band of τ c r is estimated based on a formula derived from the experimental results of Taki (2000) that is applicable to fine sediments (typically less than several tens of microns) and relatively high-water contents.

The average inundation time at the observation site is approximately 4 hours for each semidiurnal tidal cycle. The M₂ semidiurnal tide dominates the variance in water level at the intertidal mudflat, with a maximum water depth of 1.8 m ( Figure 4A ). The time series of water level during the first two tidal periods clearly exhibit diurnal variations, while the subsequent tidal cycles show nearly equal amplitude (~1.5 m). The rate of water-level change (dη/dt) decreases from the beginning of flood tides to the end of ebb tides. Extremum values are observed at the start and end of each tidal cycle.

Tidal currents on the intertidal mudflat are nearly rectilinear ( Figure 4C ). The direction of the streamwise velocity axis is 35.4° west of north, almost perpendicular to the shoreline. The velocity along the mudflat is considerably smaller than the velocity across it. Two speed extremes are observed in each semidiurnal tidal cycle when tidal currents flow onto and retreat from the observation station ( Figure 4B ), indicating the occurrence of strong surges during the very shallow water stages. Surges during such periods have been shown to play a significant role in mudflat sediment transport and bottom erosion (Wang et al., 2013; Shi et al., 2019; Zhang et al., 2021). The average duration of flood currents is 1.4 hours (red segments in Figure 4B ), whereas ebb currents have an average duration of 2.7 hours (blue segments in Figure 4B ), highlighting tidal asymmetry in flood-ebb duration.

The maximum significant wave height and the mean wave period during the six semidiurnal tidal cycles are approximately 0.18 m and 6 s, respectively ( Figures 5B, C ). These two wave parameters exhibit a clear positive correlation with local water depth, gradually increasing from the beginning of the flood tide and decreasing after high tide ( Figure 5A ). Both wave periods, significant wave height, and maximum wave orbital velocity (U_w ​) are relatively small during flood and ebb surges. The value of U_w ​​ is approximately 0.14 m/s outside of flood and ebb surges and does not exhibit significant peaks during high tide, as H_s and water depth exert opposing effects on U_w ​​ ( Figure 5D ). The intratidal variation in wave parameters suggests that as the tide level rises, the deeper water facilitates the propagation of larger amplitude waves with longer wavelengths onto the observed mudflat. The mudflat of our observation site has flat seabed, which induces the clear relationship between wave parameter and water depth. For complex topography, the wave climates may be complicated due to wave refraction and reflection (Zhu et al., 2024; Gao et al., 2024)

The temporal variation in suspended sediment concentration (SSC) in the near-bed layer at 0.1 mab reveals three distinct SSC peaks, occurring at the beginning of the flood tide, the end of the ebb tide, and during high tide ( Figures 6A, B ). The W-shaped variation in SSC is particularly pronounced in tidal cycles 3, 4, and 5. SSC variation is similar in the 0.2 mab layer, although the concentration during high tide is lower compared to the 0.1 mab layer ( Figure 6C ). At the 0.3 mab layer, SSC only shows a slight increase during high tide, with extreme SSC values primarily occurring at the beginning of the flood tide and the end of the ebb tide ( Figure 6D ). The highest SSC at the 0.6 mab layer is observed as the flood tide approaches the observation site. SSC gradually decreases until the end of the ebb tide, without significant high concentrations during high tide ( Figure 6E ). While Tides 1 and 2 capture the increase in SSC at the end of the ebb tide, this elevated SSC does not develop in subsequent tidal cycles.

The elevated SSC at the beginning of the flood tide and the end of the ebb tide may be closely associated with tidal surges, which not only facilitate the horizontal transport of SSC but also significantly resuspend sediment from the seabed. However, no direct evidence exists to confirm the mechanism responsible for the high SSC observed in the near-bed layer during high tide. Given that the mean flow is relatively weak during high tide, horizontal transport is unlikely to account for the observed increase in SSC. The significant SSC during high tide is more likely attributed to either sediment settling from the upper layers or vertical diffusion and resuspension from the lower layers. The occurrence of high SSC during high tide, as observed in this study, is not unique; similar phenomena have been reported in previous studies, such as those conducted on intertidal flats along the Jiangsu coast, China (Wang et al., 2013; Shi et al., 2019), on intertidal mudflats in Brouage, France (Bassoullet et al., 2000), and on intertidal mudflats in the Humber Estuary, UK (Christie et al., 1999).

The instantaneous suspended sediment flux per unit width (ISF, kg m−1 s−1) on intertidal flats is the product of the mean velocity and the corresponding instantaneous sediment concentration. The layers at 0.1 mab and 0.6 mab were selected to represent sediment transport in the near-bed and middle layers, respectively. The subtidal sediment flux (SSF) can be estimated as (Equation 20):

where the overbar denotes the tidal-averaged process.

In the middle layer of the water column (0.6 mab), the ISF is most pronounced during flood surges, as maximum tidal currents and elevated SSC occur simultaneously during these periods ( Figure 7B ). The ISF then decreases substantially and remains at a relatively low level throughout the rest of the tidal cycles. Since high SSC at the end of the ebb tide is not consistently observed in the middle water column, ebb surges do not always lead to significant ISF. Specifically, landward sediment flux is consistently induced by flood surges, while ebb surges may not consistently generate substantial seaward sediment flux to counterbalance the landward transport ( Figure 7C ). As a result, sediment transport tends to be landward in the 0.6 mab layer due to the tidal asymmetry of SSC on a subtidal time scale.

In the near-bed layer (0.1 mab), the instantaneous suspended sediment flux (ISF) is significant during both flood and ebb surges. Unlike the 0.6 mab layer, the ISF during flood and ebb surges are nearly balanced. The cumulative effect of ISF during high tide is substantial, despite the weak slack-tide currents. The peak SSC during high tide largely corresponds with ebb currents, leading to seaward sediment flux. To highlight the significance of sediment transport during high tide, the periods of high-water peaks of SSC are indicated by vertical dashed lines in Figures 7A, B . As shown in Figure 7D , sediment flux during high tide significantly contributes to the seaward sediment transport, accounting for approximately 30% of the total seaward flux. In summary, the transport of high SSC with weak ebb currents during high tide promotes seaward subtidal sediment flux. The subtidal sediment flux exhibits a two-layered structure, with the middle and upper layers directed landward, while the near-bed flux is directed seaward.

Sediment resuspension is closely related to bottom shear stress induced by currents and waves. According to Equation 12 and Equation 13, the lower bound of critical shear stress for sediment resuspension is estimated to be 0.0452 Pa, while the upper bound is 0.1 Pa. Significant current-induced bottom shear stress occurs during flood and ebb surges ( Figure 8 ). After these surges, the current-induced bottom shear stress rapidly drops to the lower bound of the critical shear stress for the remainder of the tidal cycles ( Figure 8 ). Due to the low orbital velocity, wave-induced bottom shear stress is minimal during flood and ebb surges ( Figure 8 ). However, during non-surge periods, wave-induced shear stress exceeds the critical shear stress, suggesting that wave-induced shear stress may prevent suspended sediment from settling and can resuspend sediment from the seabed ( Figure 8 ). In general, tidal currents generate high bottom shear stress during flood and ebb surges, while waves play a more significant role in the middle of tidal cycles.

Tidal duration asymmetry has been shown to enhance seaward sediment transport, as the near-bed high SSC during high tide closely aligns with ebb currents. This asymmetry can arise from large-scale influences, such as tide propagation in a shallowing coastal ocean, or from local asymmetries generated by the topography of flat regions (Le Hir et al., 2000). In coastal oceans primarily governed by the M₂ semidiurnal tidal cycle, interactions between the M₂ and M₄ harmonics (2ω_M2 = ω_M4) are widely recognized as the primary drivers of tidal wave deformation and associated tidal asymmetry (Friedrichs and Aubrey, 1988; Aubrey and Speer, 1985). A phase difference of 2θ_M2 − θ_M4 in the range of 0° to 180° results in flood tide dominance, with the flood tide being shorter than the ebb tide, while a phase difference in the range of 180° to 360° leads to ebb tide dominance, with the ebb tide being shorter. A phase difference of exactly 0° or 180° between 2θ_M2‐θ_M4 results in equal durations for both flood and ebb tides, thereby eliminating tidal asymmetry, though the wave shape remains statistically skewed. The amplitude ratio A_M4/A_M2 (where A is tidal amplitude) is used to quantify the magnitude of tidal asymmetry under a given phase difference.

Harmonic analysis was performed on depth-mean velocity series measured by moored Acoustic Doppler Current Profilers (ADCPs) over a 30-day period in the southern Zhejiang coastal ocean, near the observed mudflats ( Figure 1 ). Two indicators were considered: the phase differences and amplitude ratios between the M4 and M2 tidal constituents. The amplitudes and phase angles were calculated using the MATLAB program T_TIDE (Pawlowicz et al., 2002).

The harmonic analysis reveals that the phase difference of 2θ_M2‐θ_M4 ranges from approximately 70° to 110°, suggesting a shorter flood duration in the southern Zhejiang coastal ocean. The amplitude ratio between the M4 and M2 tidal constituents increases significantly from offshore to nearshore ( Table 1 ), indicating that the nonlinear effects of tidal currents intensify as water depth decreases. Tidal asymmetries established offshore of the Zhejiang coast directly contribute to tidal duration asymmetries in the intertidal mudflat.