The Effect of Wave-Induced Current and Coastal Structure on Sediment Transport at the Zengwen River Mouth ^†

Abstract

The mechanisms that control estuarine sediment transport are complicated due to the interaction between riverine flows, tidal currents, waves, and wave-driven currents. In the past decade, severe seabed erosion and shoreline retreat along the sandy coast of western Taiwan have raised concerns regarding the sustainability of coastal structures. In this study, ADCPs (Acoustic Doppler Current Profiler) and turbidity meters were deployed at the mouth of the Zengwen river to obtain the time series and the spatial distribution of flow velocities and turbidity during the base flow and flood conditions. A nearshore circulation model, SHORECIRC, has been adapted into a hybrid finite-difference/finite-volume, TVD (Total Variation Diminishing)-type scheme and coupled with the wave-spectrum model Simulating Waves Nearshore (SWAN). Conventional finite-difference schemes often produce unphysical oscillations when modeling coastal processes with abrupt bathymetric changes at river mouths. In contrast, the TVD-type finite volume scheme allows for robust treatment of discontinuities through the shock-capturing mechanism. The model reproduces water levels, waves, currents observed at the mouth of the Zengwen River reasonably well. The simulated residual sediment transport patterns demonstrate that the transport process at the river mouth is dominated by the interaction of the bathymetry and wave-induced currents when the riverine discharge was kept in reservoirs. The offshore residual transport causes erosion at the northern part of the river mouth, and the onshore residual transport causes accretion in the ebb tidal shoals around the center of the river mouth. The simulated morphological evolution displays significant changes on shallower deltas. The location with significant sea bed changes is consistent with the spot in which severe erosion occurred in recent years. Further analysis of morphological evolution is also discussed to identify the role of coastal structures, for example, the extension of the newly constructed groins near the river mouth.

1. Introduction

Understanding hydrodynamics and sediment transport in the coastal zone is crucial for coastal protection, sustainable maintenance, and decision support, which are in accordance with social, economic, and ecological requirements. Numerical models are widely used by engineers to simulate bathymetric change due to the interactions of coastal processes and ocean structures [1,2,3,4]. Comprehensive predictions based on numerical simulations provide the information for engineering and ecosystem applications. Therefore, there is a need for the development of the modeling approach based on comprehensive field measurements. This study tends to combine field and numerical approaches that directly resolve the flood-ebb tidal cycles and energetic waves to simulate the transport processes due to the interaction of waves and currents over complex bathymetry at the river mouth.

Near mouths of estuaries and coastal zones, intense waves, and wave–current interactions play an important role in the exchange of sediment between the estuary and the adjacent coast. Prior numerical studies indicate that the riverine flow brings sediment at both sides of the river entrance and the adjacent coastal zone during the extreme riverine flow condition. For example, field evidence in the Choshui River and the Gaoping River in Taiwan shows that the sediment-laden river flow was first deposited on the ebb tidal shoal adjacent to the river mouth, subsequently re-suspended and transported via the oceanic circulation system [5,6]. During moderate to low riverine discharge conditions, shallower shoals at river mouth due to sediment deposition plays an important role on the evolution of wave energy, flow patterns, and the resulting sediment transport [7]. The prior study in the Delaware Estuary, USA, found that modeled incident wave energy fluxes at the shoreline have similar spatial distribution as the observed rates of shoreline retreat averaged over the past hundred years. Regions with some of the most severe shoreline retreat correspond with locally intensified wave energy fluxes [8]. Previous modeling studies in New River Inlet, NC, further demonstrate that the flows and water levels across the ebb shoal are modulated by breaking waves at the offshore edge of the ebb shoal. Tidally averaged flows on the ebb shoal were seaward directed owing to this asymmetry in flood versus ebb [9,10]. In the past decade, the evidence of seabed erosion near the Zengwen river mouth has raised concerns regarding coastal resilience, for example, coastal erosion occurred near the northern ebb tidal delta of the Zengwen River mouth (see the location of the newly constructed groins near the river mouth in Figure 1). Hence, the mouth of the Zengwen River provides an excellent environment for our understanding of sediment transport processes due to waves and wave–current interactions.

The Zengwen River has a total length of 138 km and a drainage area of 1176 km². The average slope of the Zengwen River in its lower reaches is only 1:3500, lower than most of the rivers in Taiwan. The natural environment of the Taiwan Strait has been known for its high current speed and the tidal wave surges into the strait and intersects at the center of the strait [11,12]. The averaged tidal range near the Zengwen River is 1 m, and up to 1.37 m during spring tides. Because of the cross-sectional contraction of river dykes, the mean tidal range may reach a maximum of 2.42 m within the system. The phase relationship between tidal elevation and tidal flow is close to standing wave characteristics [13]. Seawater intrudes upriver as a result of tidal dispersion and residual circulation [14,15,16]. Four reservoirs were built in the catchment of the Zengwen River, where occupied approximately 55% of the basin. The reduction of riverine flow and sediment supply due to the construction of the Zengwen reservoir may be one of the reasons which causes significant shoreline erosion at the study site. When the riverine discharge was kept in reservoirs, the coastal processes, such as waves and tidal currents may play an important role in the redistribution of sediment and the morphological evolution near the river mouth.

In this study, a numerical investigation on wave–current–bathymetry interactions at the mouth of the Zengwen River is carried out in conjunction with a field experiment. This study aims to understand coastal processes and their effects on sediment transport and the resulting morphological evolution. The alternative goal of the proposed research is to combine advanced numerical models with observations to facilitate the prediction and coastal management under the effect of waves and local anthropogenic change. It is conjectured that the reduction of riverine flow and sediment supply due to the construction of the Zengwen reservoir may reduce the riverine deposition and enhance coastal erosion processes. The river or sea dike, which is used to regulate water from a river and the ocean, may restrict the flood stream, speed up the flow velocity, and cause the sediment transport to the deep-sea instead of depositing in the nearshore region. A diagnostic study combing field and numerical methods is carried out to understand the characteristics of waves, wave-driven currents, and the resulting sediment transport in the study area. An array of co-located wave gauges, current meters (ADCPs), and turbidity meters (Lisst 100X) were deployed throughout the channels and ebb tidal shoals (Figure 2)to obtain the time series and spatial distribution of hydrodynamic and sediment transport conditions at the river mouth in the winter of 2018 (21–28 November) and the summer of 2019 (23 August–4 September). The river discharge is around 4 m³/s during 21–28 November and up to 200 m³/s during 23 August–25 August due to reservoir water release. A nearshore circulation model, SHORECIRC [17], has been adapted into a hybrid finite-difference/finite-volume, TVD (Total Variation Diminishing)-type scheme [18] and coupled with the wave-spectrum model Simulating Waves Nearshore (SWAN) [19]. The quasi-3D circulation model includes the effects of waves on the vertical structure of the currents based on theory [20]. The numerical model is fully parallelized with MPI and is suitable for long-term (months), large-scale (O (10–100) km) simulations. NearCoM-TVD was verified and validated with several coastal applications [21]: the modeled flow velocity through an idealized tidal inlet was compared with an analytical solution [22]; the modeled wave heights and current velocities over a rip channel system were compared with the measured data during a rip-current field experiment to evaluate the model’s capability for applications with wave–current interactions [23]. NearCoM-TVD was applied to studies in New River Inlet, USA and modeling results was validated with co-located wave gauges, current meters (ADCPs), and observations from autonomous underwater vehicle (AUV) [8,9,10]. Here, NearCoM-TVD is applied to study sediment transport and morphological evolution at the mouth of the Zengwen River under the interaction of tides and waves. The paper is organized as follows: the theory and formulations of the quasi-3D nearshore community model are discussed in Section 2. Section 3 andSection 4 describe the model setup and model calibration between simulation and observation data. The results of model simulation and the effect of waves and wave-induced currents on sediment transport are discussed in Section 5. The study is summarized in Section 6.

2. Numerical Model

NearCoM-TVD (Figure 3) is one of the coastal modeling systems, which integrates the wave model SWAN and a modified version of the circulation model SHORECIRC [24,25] with sediment transport formula [26]. The wave model Simulating Waves Nearshore (SWAN) has been developed [19] to estimate wave conditions in small-scale, coastal regions with shallow water, (barrier) islands, tidal flats, local wind, and ambient currents [27]. Wave can dominate a variety of processes associated with wave-breaking and wave-current interaction on shallower shoals [28]. In addition, waves can cause enhancement of bottom boundary-layer turbulence, wave orbital motions, apparent bed roughness, the increase of the vertical turbulent mixing, and suspended-sediment carrying capacity [29].

SWAN is a spectral wave model which solves the wave action balance equation. The governing equation can be written as

$$\frac{\partial N}{\partial t}+\frac{\partial \left(C_g^{\alpha }N\right)}{\partial x_{\alpha }}+\frac{\partial \left(C_g^{\sigma }N\right)}{\partial \sigma }+\frac{\partial \left(C_g^{\theta }N\right)}{\partial \theta }=\frac{S}{\sigma }$$

(1)

where N is the wave action defined by $E/\sigma$; E is the wave energy density; $\sigma$ is the relative angular frequency; $\theta$ is propagation direction of each wave component; $C_g^{\alpha }$ represents the propagation speed of wave energy, and S represents source and sink terms of energy density.

SHORECIRC is a two-dimensional horizontal (2DH) model that incorporates the mixing effect induced by the vertical variation of wave-induced horizontal circulations. The instantaneous total fluid velocity ${\mathrm{u}}_{\mathrm{ins}}^{\mathsf{\alpha }}$ in coordinates (x₁, x₂) is split into four components:

$${\mathrm{u}}_{\mathrm{ins}}^{\mathsf{\alpha }}=u_{\alpha }^{’}+{\mathrm{u}}_{\mathsf{\omega }}^{\mathsf{\alpha }}+{\mathrm{u}}^{\mathsf{\alpha }}+{\mathrm{u}}_1^{\mathsf{\alpha }}$$

(2)

where $u_{\alpha }^{’}$ is the turbulence component, ${\mathrm{u}}_{\mathsf{\omega }}^{\mathsf{\alpha }}$ is the wave component whose short-wave-averaged value (ū_w) is zero below trough level; ${\mathrm{u}}^{\mathsf{\alpha }}$ and ${\mathrm{u}}_1^{\mathsf{\alpha }}$ are two components of the short-wave-averaged velocity. In the direction terms, α = 1 represents the east–west and α = 2 represents the north–south. To get the third component (u^α)$,$ the depth-averaged short-wave-averaged velocity, “Lagrangian averaging” is used as

$$u^{\alpha }=\frac{1}{H}\overline{\underset{-h}{\overset{\xi }{{\displaystyle \int }}}{u}_{\alpha }^{ins}dz}$$

(3)

$$H=\mathsf{\eta }+\mathrm{h}$$

(4)

where η is the wave-averaged surface elevation, h is the still water level, and ξ is the instantaneous surface elevation.

Based on the depth-integrated, the formulations of SHORECIRC model [17] in Cartesian coordinates are given as

$$\frac{\partial \eta }{\partial t}+\frac{\partial H{u}_{\alpha }}{\partial x_{\alpha }}=0$$

(5)

$$\frac{\partial H{u}_{\alpha }}{\partial t}+\frac{H{u}_{\alpha }{u}_{\beta }}{\partial x_{\beta }}+f_{\alpha }+gH\frac{\partial \eta }{\partial x_{\alpha }}+\frac{1}{\rho }\frac{\partial T_{\alpha \beta }}{\partial x_{\beta }}+\frac{1}{\rho }\frac{\partial S_{\alpha \beta }}{\partial x_{\beta }}+\frac{{\tau }_{\alpha }^b}{\rho }-\frac{{\tau }_{\alpha }^s}{\rho }+ROT=0$$

(6)

Equation (6) shows the added Coriolis force term f_α caused a deflection of the path of an object that moves within a rotating coordinate system. The important terms on this equation T_αβ,S_αβ,${\tau }_{\alpha }^b, {\tau }_{\alpha }^s$, ROT are the depth-integrated Reynold stress, the wave-induced radiation stress [30], the bottom shear stress, the surface shear stress, and the rest associated with 3D dispersion, respectively.

Wind-induced current is one of the critical factors for hydrodynamics in the ocean especially for the transport of mass, energy, and momentum between air and sea. The SHORECIRC model also defines the wind-induced surface stress as

$${\tau }_{\alpha }^s=C_{dw}{\rho }_a\left|\mathit{W}\right|W_a$$

(7)

where C_dw, ρ_a, W, and W_α are the drag coefficient which follows the formula of The WAMDI Group [31], air density, the wind speed measured at 10 m elevation of the water surface, and stress (wind velocity vector), respectively.

The wave-averaged bottom stress equation in SHORECIRC model is given as

$${\tau }_{\alpha }^b=\frac{1}{2}f\rho u_0\left({\beta }_1{u}_{b\alpha }+{\beta }_2{U}_{w\alpha }\right)$$

(8)

where f, u_bα, are a constant friction factor [32] and the current velocity at the bottom. $u_0$ is the magnitude of current velocity ($u_{b\alpha }$) or the magnitude of wave velocity U_wα when U_wα > u_bα. The current velocity, u_bα, with weighting factor β₁ and the wave velocity, U_wα, with weighting factor β₂ are the other contributors to the shear stress ${\tau }_{\alpha }^b$. The weighting factors (β₁ and β₂) are the function of U_wα/u_bα and the angle between wave and current vectors taken from the laboratory experiment data under the monochromatic wave condition [33]. We set β₁ = 1.0 and β₂ = 0.5 in the present model. To integrate the effect of water depth H with a given Manning coefficient M, the friction factor f is calculated using the Manning formula.

$$f=2g{M}^2/H^{1/3}$$

(9)

The total load sediment transport model is closely connected to the current-dominated condition in which tidal currents and wave-induced currents determine the direction of net transport but waves can enhance the magnitude of transport through $u_{rms}^{}$ [26]. If sediment transport and the evolution of seabed in the coastal zone are driven by waves and wave induced-currents, the total sediment transport rate combining the effect of currents and waves on the bed shear stress is given as

$$q_{ɑ}=A_s{u}_{\alpha }{\left[\sqrt{{\left|u\right|}^2+\frac{0.018}{C_d}{u}_{rms}^2}-u^{cr}\right]}^{2.4}\left(1-1.6tan\beta \right)$$

(10)

where $u_{\alpha }$ is the short-wave-averaged velocity in coordinates (x₁, x₂), |u| is the magnitude of current velocity, u_cr is the critical velocity for erosion, which is determined by the water depth and the mean grain size of the sediment (D₅₀), u_rms is the root mean square of wave orbital velocity, and β is the bed slope. Among them, the sediment load (A_s = A_sb + A_ss) is calculated based on the characteristics of the bed load (A_sb) and suspended load (A_ss)$:$

$${\mathrm{A}}_{\mathrm{sb}}=\frac{0.005\mathrm{h}{\left({\mathrm{D}}_{50}/h\right)}^{1.2}}{{\left[\left(\mathrm{s}-1\right){\mathrm{gD}}_{50}\right]}^{1.2}}$$

(11)

$${\mathrm{A}}_{\mathrm{ss}}=\frac{0.012{\mathrm{D}}_{50}{\mathrm{D}}_{*}^{-0.6}}{{\left[\left(\mathrm{s}-1\right){\mathrm{gD}}_{50}\right]}^{1.2}}$$

(12)

where ${\mathrm{D}}_{*}={\left[\frac{\mathrm{g}\left(\mathrm{s}-1\right)}{v^2}\right]}^{1/3}{\mathrm{D}}_{50}$, s is relative density of sediment, and $v$ is kinematic viscosity of water. The drag coefficient (C_d) plays an important role for calculating the total sediment transport rate which is inversely parameterized to the wave stirring and can be defined as the equation below:

$$C_d={\left[\frac{0.40}{\ln \left(\frac{h}{z_0}\right)-1}\right]}^2$$

(13)

The magnitude of wave stirring is also influenced by the bed roughness (z_o = 0.001 m).

The seabed evolution can be described by using sediment transport flux in Cartesian coordinates:

$$\left(1-p\right)\frac{\partial h_1}{\partial t}+\frac{\partial f_{mor}{q}_{ɑ}}{\partial x_{\alpha }}=0$$

(14)

where $h_1$ is the bed evolution, $q_{ɑ}$ is the sediment transport rate based on the given hydrodynamic condition, p is the bed porosity, and $f_{mor}$ is the morphological factor based on the new method developed for modeling morphological change [34]. In this paper, we used $f_{mor}$= 12. Equation (14) is solved by using upwinding first order finite difference scheme. The dissipation is minimal due to the small-time step used in the circulation model.

3. Model Configuration

The Zengwen river is located on the southwest coast of Taiwan. The system includes a channel (3 < depth < 5 m) and shallower shoals (depth ≤ 3 m) on the northwestern side (Figure 4). The model domain, which includes the river channels and the edge of the continental shelf (depth ≤ 20 m), has a length of 20 km and a width of 9 km. Three DEMs are integrated to provide the bathymetry used in numerical simulations: the digital elevation model (DEM) provided by the Ministry of Science and Technology (MOST), Taiwan, has a 200-m resolution in coastal regions. The DEM surveyed in the year 2017 has a 10-m resolution in the surf zone, shallower deltas, and the adjacent lagoon. The DEM of the river channel is interpolated by surveyed cross-sections in the year 2017 using CCHE TOOL. The tidal boundary conditions are provided by the global tidal level forecast of Oregon State University (OSU) with eight tidal constituents (K1, O1, P1, Q1, K2, M2, N2, and S2). The dominant tidal constituent in the northern boundary is M2 (tidal amplitude ~ 0.22 m), with smaller amplitudes at other tidal constituents, including K1 (amplitude ~ 0.17 m), S2 (amplitude ~ 0.07 m), and O1 (amplitude ~ 0.17 m). The M2 tidal constituent in the southern boundary is greater and the phase difference is ~12 degree from the northern boundary. Therefore, the tidal currents are southwestward during ebb tides and northeastward during flood tides. The boundary conditions include two spring tides and one neap tide. The boundaries adjacent to the lagoons and the river boundary are closed. The river discharge is based on the daily mean observation of Hsin-Chung Station (1720H025) in which the observed surface elevation shows non-tidal variation. The significant wave heights and the peak period observed at the Anping station (16 m depth, see Figure 1 for the location) are applied to the boundary of the spectral wave model SWAN. The observed wave heights ranged from 0.5 to 3.1 m. The direction of observed waves is from the northwest during winter and from the southwest during summer. The Coriolis parameter is calculated using the latitude of 23.5 N.

4. Model Validation and Statistic Error Analysis

The model accuracy is assessed with the correlation coefficient and the root mean square error. The correlation coefficient (R) for continuous data ranges from −1 to 1 in which the higher linear association can be represented as the nearer of scatter points to the regression line. Considering that M_n and C_n are the measured data and the computed data, respectively, at N discrete points, and ${\mathsf{\sigma }}_{\mathrm{M}}$ and ${\mathsf{\sigma }}_{\mathrm{C}}$ are the standard deviation of the measured and the computed data, the formulas are given by

$$\mathrm{R}=\frac{\frac{1}{\mathrm{N}}{{\displaystyle \sum }}_{\mathrm{n}=1}^{\mathrm{N}}\left({\mathrm{M}}_{\mathrm{n}}-\overline{{\mathrm{M}}_{\mathrm{n}}}\right)\left({\mathrm{C}}_{\mathrm{n}}-\overline{{\mathrm{C}}_{\mathrm{n}}}\right)}{{\mathsf{\sigma }}_{\mathrm{C}}{\mathsf{\sigma }}_{\mathrm{M}}}$$

(15)

The root mean square error (RMSE) is used to measure the spreading of the observed data about the predicted values. The values of RMSE are between 0 and 1 in which the smaller value presents the good performance of the numerical model. RMSE is frequently used by researchers due to the easier interpretation of results by converting the error metric back into similar units. The formulas are given by

$$\mathrm{RMSE}={\left[\frac{1}{\mathrm{N}}{\displaystyle \sum }_{\mathrm{n}=1}^{\mathrm{N}}{\left(M_n-{\mathrm{C}}_n\right)}^2\right]}^{1/2}$$

(16)

The model results are compared with observations from the Sicao tidal gauge (year-round) and bottom-mounted wave gauges/current meters during 21–28 November 2018 and 23 August–4 September 2019 (see Figure 2 for the location). The measured current profiles are depth-averaged to compare with model results. The model accurately predicts the neap-spring variation of surface elevation at SiCao tidal station (Figure 5 and see Figure 1 for the location), at the main section of the river channel (station A, C and E, see Figure 6a as an example), and at the side of the main section of the river channel (station B and D).

Table 1. Model performance of surface elevation, significant wave heights, and current speed for all the sensors (see the location in Figure 2a) during 21–30 November 2018.

Station	R			RMSE
	$\mathit{\eta }$	Hs	Speed	$\mathit{\eta }$	Hs	Speed
A	0.96	0.489	NaN	0.096	0.18	NaN
B	0.96	0.441	0.39	0.098	0.13	0.157
C	0.96	0.498	0.236	0.098	0.139	0.11
D	0.96	0.459	0.55	0.093	0.157	0.161
E	0.96	0.501	0.275	0.098	0.084	0.195

and

Table 2. Model performance of surface elevation, significant wave heights, and current speed for all the sensors (see the location in Figure 2b) during 23 August–4 September 2019.

Station	R			RMSE
	$\mathit{\eta }$	Hs	Speed	$\mathit{\eta }$	Hs	Speed
A	−0.43	0.602	0.077	1.123	0.334	0.232
B	0.664	0.293	0.332	0.836	0.237	0.109
C	0.874	0.265	0.039	0.553	0.191	0.102
D	0.907	0.552	NaN	0.152	0.203	NaN
E	0.926	0.126	0.292	0.136	0.148	0.171
F	0.887	0.149	0.085	1.591	0.18	0.1
G	0.89	0.416	0.262	0.239	0.389	0.131
H	0.934	0.754	NaN	0.125	0.292	NaN

show model performance using the correlation coefficient and root mean square error during winter and summer, respectively. Overall, the correlation coefficient is around 0.96 and the root mean square error is less than 0.1 m in most of the stations in the winter survey. The correlation coefficient is lower in summer due to the storm surge on the 25th of August during the passage of typhoon was not captured by the model simulations based on astronomical tidal boundary conditions.

The model also predicts the observed period, direction, and the variation of significant wave heights reasonably well. The simulated wave period is around 3 to 4 s in winter and around 6 to 8 s in summer. During the winter survey waves are consistently propagated from the southwest and during the summer survey waves are mostly propagated from the southwest or the south direction to the river mouth. During winter, the observed wave direction is similar to the observation near the Anping Harbor but different from the data from the Qigu wave buoy in the north of the river mouth owing to the topographic refraction and the shedding effects of coastal geometry (see Figure 1 for the location). The RMSE of waves varies from 0.08 m to 0.18 m (see Table 1) and the Normalized Root Mean Square Error (NRMSE) are around 30%. Figure 6b shows the example of model data comparison of significant wave heights at station C (water depth 1.705 m) in the main section of the river channel during 21–28 November 2018. The correlation coefficient is around 0.5 and the root mean square error is less than 0.14 m at station C. In the summer survey, the model generally captures the stormy wave heights (more than 2 m) during the passage of the storm but not able to predict the trend during the condition with weaker wave energy (offshore wave heights less than 1 m) in summer. The correlation coefficient is high in the offshore locations, for example station A and H, but the correlation coefficient is low near the entrance of the river, for example at station E (water depth 0.961 m), where observed wave heights are random and smaller (less than 0.5 m) compared to the offshore stations due to wave breaking processes. The breaking wave heights on shallower shoals are also relatively small (for example station B, 2.837 m depth), so the correlation coefficient between observations and simulations are relatively low.

Figure 6c,d shows the model-data comparison of north–southward velocities and east–westward velocities. The model also predicts the observed north–southward velocity with a correlation coefficient around 0.66 and RMSE around 0.11 at station C. The correlation coefficient of current speed is up to 0.55 at the main section of the river channel (station D). Both observations and model data show these locations are dominant by tidal currents with more significant tidal variation in the north–southward velocity. The correlation coefficient of current speed reduces to 0.236 due to the lower correlation between the observed and modeled east–westward component of the current velocity. The data model comparison of the current speed during the summer survey shows a lower correlation coefficient with RMSE <0.2 m/s in most of stations due to the complex interactions of waves, tides, and riverine discharge. However, the model can predict the magnitude of observed flood discharge near GuoXing Bridge (see Figure 1 for the location), and the magnitude of current speed with the root mean square error less than 0.2 m/s in most of the stations near the river mouth.

5. Model Result and Discussion

5.1. Wave and Flow Field under Different Tidal Hours

The simulation results of surface elevation, waves and flow velocities in the estuary have reasonable agreements with observations. Hence, model results are used to provide insights into the patterns of flow residual for a range of tidal forcing and wave conditions (Figure 7) over the complex bathymetry. The wave and current field during the boundary conditions of (a) the maximum ebb with energetic waves (water level −0.7 m and significant wave heights 1.5 m), (b) flood tides with moderate waves (water level 0.1 m and significant wave heights 1.0 m), (c) the maximum flood tide with mild waves (water level 0.5 m and significant wave heights 0.55 m), and (d) ebb tide with mild waves (water level −0.1 m and significant wave heights 0.52 m) are discussed to understand the relevant importance of waves and tidal currents in the study site.

To understand the wave and tidal characteristics of the coastal regions and near the river mouth under different tidal conditions, the simulated nearshore wave, current fields are shown in Figure 8 and Figure 9. The location of wave breaking is modulated by the water level (the ratio of wave height and water depth is set as 0.78 in SWAN) and the background current velocity. The simulated wave field during different tidal hours shows that the breaking zone moves offshore during ebb tides and onshore during flood tides. The simulated wave field is more intense during the maximum ebb tide when the water level of boundary condition is −0.7 m with offshore significant wave heights 1.5 m (Figure 8a). During the slack tide with energetic waves (the water level is 0.1 m and offshore significant wave heights 1.0 m, Figure 8b), the simulated significant wave heights are up to 0.8 m near the 5 m contour line and decrease gradually over the slopping bathy of the shallow shoals in the river mouth.

Although the amount of riverine discharge is insignificant during 21–28 November 2018, the flow velocity in the river mouth is up to 0.5 m/s because of the effect of tides and waves (Figure 9a). The formation of ebb tidal jet is caused by the pressure gradient due to the water level difference between the river channel and the open ocean. Model results further demonstrate complex flow patterns due to wave-induced currents and the ebb tidal jet. The intensity of nearshore circulations is further enhanced by waves during more energetic wave conditions during ebb tides. Especially, the simulated current velocity on the shallower shoal outside of the river mouth is also more intense during ebb tides. Under similar wave conditions (significant wave heights ~ 0.5 m), the flow intensity near the river mouth is more intense during ebb tides (Figure 9d) but nearly zero during the maximum flood tide (Figure 9c). These flow patterns are absent when waves are excluded from the boundary of model simulations. Model simulations suggest that the flow field near the river mouth is strongly influenced by tides and waves.

5.2. The Seasonal Variations of Wave and Flow Field

The analysis of observed wave direction shows that waves are consistently propagated from the southwest direction to the river mouth in winter. In summer, observations at the station G and H (see Figure 2b for the locations) show that mostly waves are propagated from the southeast or the south direction to the river mouth but observations at station B and D show that waves are from the southwest due to the refraction of shallower shoals. During summer, the variation of observed waves at station A and H are similar to data provided by Qigu wave buoy in the north of the river mouth and the observation near the Anping Harbor. Model results during the stormy wave condition in summer are also discussed to provide insights into the flow patterns in the coastal regions and the effect on the sediment transport and morphological evolution. When the water level of boundary condition is around 0.2 m and offshore significant wave heights are 2.0 m, the simulated significant wave heights are up to 1.8 m near the 5 m contour line and decreases to 0.5 m on the shallow shoals in the river mouth. The flow velocity on the shallower shoals in the river mouth is around 0.5–0.8 m/s (Figure 10). The formation of nearshore circulations is caused by the interactions between wave-driven currents and riverine inflows during stormy wave conditions. However, the current velocities around the river mouth are low during milder wave conditions (<1 m) most of the time in summer. The flow velocities are small on the nearshore regions but more intense in the river channel due to the tidal effects.

Comparing with the observations and model simulations during the northeast monsoon, the wave heights with smaller magnitude and more significant tidal variations are shown at station B and D during most the time in the summer (exclude few days during the passage of typhoon). The observed and simulated significant wave heights are less than 0.5 m near the river mouth (at station C, F, and E, see Figure 2b for the location). In both seasons, waves are normal incident to the river mouth (at the station E) due to bathymetric refraction. The observed and simulated wave heights are relatively small at the station E so the observed current speed near the river mouth is mainly dominant by tidal currents. However, the observed and simulated currents at station B, and C are affected by both tidal currents and wave-induced longshore currents (see Figure 2a for the location).

5.3. The Residual Sediment Transport and Morphological Evolution

The simulated residual sediment flux and the resulting morphological evolution are further discussed based on the total load transport formula of Soulsby (1997) [26]. The mean grain size (D₅₀) is 0.110–0.172 mm during the winter experiment and 0.135–0.155 mm during the summer experiment according to the sample from the field experiment (see the sampling location in Figure 2). Therefore, the critical shear stress for erosion defined by Van Rijn (1984) [26] varies between 0.21 m/s and 0.30 m/s in the water depth of 1 m to 4 m based on the given D₅₀ = 0.1–0.2 mm. The bed roughness z_o in Equation (12) is set as 0.001 m which is approximated as D₅₀ = 0.2 mm over the entire domain. The magnitude of residual transport and the depth of erosion changes with the given parameters. The patterns of sediment fluxes and the locations of erosion spots are consistent through all test cases. The simulated nearshore circulations on the shallower shoals are greater than the estimated threshold when the significant wave heights near the river mouth are greater than 1.0 m. Therefore, in the study area, the transport of sediment on shallow shoals near the river mouth can be triggered by the relatively strong current velocities during the most time in the winter and the stormy wave condition in summer.

Our data analysis also shows that the observed sediment flux (turbidity multiply current speed, sediment concentration is not available in this field experiment) is highly correlated with the cubic of observed velocity due to the tidal currents and wave-driven currents. The variation of sediment fluxes in the main section of the river channel (station A, C, and E) has significant tidal variations. The maximum values of sediment fluxes are shown in ebb tides due to the interaction of ebb tidal jet and wave-driven longshore currents. The analysis shows that the sediment flux outside the breaker zone (in station A) significantly increases during more energetic wave conditions between the 23rd and 24th of November. However, the transport of sediment at shallower shoals is mainly modulated by tidal currents (Station C and E) or wave-driven currents (Station B and D). Model results show that the simulated current speed is in reasonable agreement with the observation. Therefore, the model is further utilized to study sediment transport. Model results indicate that the magnitude of sediment transport rate is more persistent under the northeast monsoon condition during the winter season. Therefore, the simulated residual sediment flux and the resulting morphological evolution during northeast monsoon are further discussed based on the total load transport formula of Soulsby. The simulated nearshore circulations in the river channel are greater than the estimated threshold during the maximum ebb and flood tide so the sediment can be transported landward to the river mouth and seaward to the adjacent coastal regions, respectively. To understand the net transport due to tidal effects, the residual flow, Ur, commonly known as the Eulerian-averaged residual flow, is useful for evaluating possible sediment transport and bottom stress patterns [26]. The tidally averaged velocity can be obtained directly by time-averaging the flow velocity vector over tidal cycles $U_{\mathrm{r}}=u$, where “< >” represents time-average over a tidal cycle.

Similar to the pattern of the simulated residual flow, the simulated residual sediment flux near the river mouth is influenced by tides and wave-induced currents. The simulated sediment fluxes also show that the transport of sediment converges on the shallower shoal near the river mouth (indicated by yellow vectors in Figure 11). The other hotspot for a high sediment transport rate is near the seawall with the newly constructed groins. The transport direction changes from alongshore to offshore due to the changing of coastal geometry. This pattern near the seawall is absent when waves are excluded in the model boundary. Model experiments indicate that residual tidal currents play an important role in net transport offshore in the river channel and the wave-driven circulations may further redistribute a significant amount of sediment from the riverine deposition near the mouth.

The simulated morphological evolution displays significant changes in the depth along the navigation channel during both seasons. The model results of morphological change after two weeks of simulation with a morphological factor of 12 are shown in Figure 12. The location with significant sea bed change is consistent with the spot with severe erosion in recent years. The offshore residual transport causes erosion at the northern part of the river mouth, and the onshore residual transport causes accretion in the ebb tidal shoals around the center of the river mouth. The simulated morphological changes are less significant (Figure 12b) when the newly constructed groins near the river mouth (Figure 1) are included in the model bathymetry. Model simulations suggest that the construction of groins can slow down the flow velocities and reduce the sediment transport rate at the northern corner of the river mouth.

6. Conclusions

Understanding hydrodynamics and sediment transport in estuaries and the adjacent coastal area is critical but challenging owing to the highly nonlinear interaction among tides, waves, and bathymetry. Locally intense circulations can be generated under the interaction of tidal currents wave-induced currents and complex bathymetry. A diagnostic study combing field and numerical methods is carried out in order to understand the dominant mechanisms causing the resulting complex flow pattern and sediment transport. An array of co-located wave gauges, ADCPs, and turbidity meters were deployed throughout the channels and ebb tidal shoals to obtain the time series and spatial distribution of hydrodynamic and sediment transport conditions at the river mouth in winter and summer. A quasi-3D nearshore community model is applied for the integrated observational and modeling study. NearCoM-TVD, couples SWAN and SHORECIRC, reproduces water levels, waves, currents observed at the river mouth reasonably well. The study indicates NearCoM-TVD can predict comprehensive physical processes over complex bathymetry and provide the information of sea sates for engineering applications in the coastal ocean. Model results suggest that flood flows over deltas funnel into the river mouth radially from all sides of the entrance, whereas ebb flows leave the entrance in a concentrated jet. The validated model is used to examine the complex hydrodynamic patterns of the channel–delta system. The circulation patterns at the mouth of Zengwen River may differ from prior model simulations of New River Inlet, NC [8,9,10], Ria de Ribadeo, Spain [35], and Willapa Bay Inlet, WA [7], owing to differences in the river mouth geometries and the complex channel–delta bathymetries. However, similar to the condition in those ebb dominant tidal inlets, tidally averaged flows on the ebb shoal were seaward directed owing to this asymmetry in flood versus ebb. The ebb-dominant tidal jet results in offshore transport in the main section of the river channel, whereas the circulation patterns resulting from waves and wave–current–bathymetry interactions may carry sediment landward on the ebb tidal delta. The model result shows that the breaking zone moves offshore during ebb tides and onshore during flood tides. Wave-break processes generate cross-shore/long-shore currents, which affect sediment transport, and morphological evolution. Similar to the assumption in the previous study based on an idealized inlet [7], our observational and numerical study in the Zengwen river demonstrates that the shallower shoal at the river mouth modulate the location of wave breaking and produce locally intense circulations due to the interactions of tidal currents and wave-induced currents. The simulated residual sediment transport patterns show that the transport process at the river mouth is dominated by the interaction of ebb tidal jet and wave-induced longshore current during low riverine flow conditions. The rate of residual flow and the residual sediment transport are higher in the location with significant topographic changes. The wave-induced current causes erosion at the northern part of the river mouth. The modeled morphological evolution suggests that the construction of groins can slow down the flow velocities and reduce the sediment transport rate at the northern corner of the river mouth.