Numerical Simulation of Hydrodynamic Processes for East Shantou Reclamation Engineering

1. Introduction

Coastal reclamation is an important approach to effectively alleviate the shortage of land resources in coastal areas. While providing more development space and bringing enormous economic benefits to humans, coastal reclamation alters the topography and geomorphology of natural coastlines and adjacent marine areas, exerting certain impacts on hydrodynamic conditions such as tidal currents and waves in nearby waters [1]. Meanwhile, it also triggers a series of ecological problems and exerts negative effects on the natural environment[2,3,4,5].

Numerous scholars at home and abroad have adopted various methods to investigate the impacts of coastal reclamation on hydrodynamic environments and sediment transport. He and Xin[6] simulated the tidal current movement in Shenzhen Bay before and after engineering projects based on a two-dimensional tidal numerical model, revealing that the environmental impacts of the projects in Shenzhen Bay were mainly concentrated in the vicinity of the construction area. Zeng et al.[7] explored the effects of coastal reclamation projects on the hydrodynamic environment of Xiangshan Port during the periods of 1963-2003 and 1963-2010 by using the Princeton Ocean Model (POM). Based on the Finite Volume Coastal Ocean Model (FVCOM), Park et al.[8] and Gao et al.[9] demonstrated that coastline changes induced by coastal reclamation altered the reflection of tidal waves, which would reduce tidal amplitude after reclamation, thus revealing the evolution law of tidal dynamics caused by reclamation. Kuang et al.[10] established a three-dimensional hydrodynamic model with Delft3D to simulate the impacts of coastal reclamation projects on tidal circulation, flow velocity and water exchange. Maren et al.[11] analyzed the long-term effects of channel deepening, port construction and dredging works on suspended sediment concentration (SSC) via a three-dimensional Delft3D numerical model. Yao et al.[12] and Zou et al.[13] employed the MIKE21 numerical model to simulate the variations in tides, currents and sediments before and after the implementation of engineering projects, conducting a quantitative analysis of the impacts of reclamation on hydrodynamics, sediment diffusion and topographic erosion-siltation in surrounding marine areas.

For hydrodynamic simulations under wave-current coupling, the research findings of Fearghal and Scott[14] indicated that current-wave coupled models can effectively improve the simulation accuracy of tidal current fields and wave fields. Wu[15] adopted the MOCOE hydrodynamic model to simulate the tidal current fields and wave fields in the coastal waters of the Yamen River. Based on the Delft3D model, Vu et al.[16] demonstrated that under the combined action of waves and runoff, suspended sediment transport in the Mekong Delta is mainly governed by wave effects. Liu and Cai[17] constructed a three-dimensional wave-current-sediment coupled numerical model to evaluate the scenarios of suspended sediment distribution in the Pearl River Estuary under the influence of different driving forces. Elias et al.[18] employed the Delft3D model to simulate the hydrodynamic distribution in the Egmond coastal waters, and the results showed that the model can effectively capture the core characteristics of water levels, waves and current fields in the study area.

In the research on wave-current coupled models, most studies currently focus on hydrodynamic simulations under existing conditions. In terms of the impacts of reclamation projects on the hydrodynamic environment, the majority of relevant analyses only investigate the variations in hydrodynamic conditions before and after project construction from the perspective of tidal current fields. Moreover, there is a dearth of research literature on reclamation projects in areas with multiple distributary estuaries (e.g., the Hanjiang River Delta) and complex topographic conditions. Therefore, it is of great significance to adopt a wave-current coupled hydrodynamic mathematical model to study the impacts of the reclamation project in the Eastern Urban Economic Belt of Shantou on the local hydrodynamic environment.

The Coastal Estuary Regulation and Comprehensive Development Project for the Eastern Urban Economic Belt of Shantou is located in the shoal area of the Hanjiang River Delta in eastern Shantou City, Guangdong Province. Initiating from the sand barrier of Shantou Port, the project runs roughly parallel to the existing coastline and extends northeastward, passing through the estuaries of the Xinjin River and Waisha River before terminating at the tombolo of Laiwu Peninsula, with a reclamation area of approximately 20 km². After the project implementation, the coastline advances seaward by 1.5-2.4 km, showing a significant positional shift. Meanwhile, a shortage of about 6100×10⁴ m³ of reclamation materials is required for land creation, which needs to be addressed by sand mining from sand extraction pits near the project area. Given the enormous sand mining volume, the resulting sand extraction pits cover a relatively extensive area (see Figure 1). Owing to the substantial changes in the coastline and seabed topography within the project area after the reclamation project, coupled with the project involving multiple distributary estuaries of the Hanjiang River, the topographic and flow boundary conditions are extremely complex. Consequently, the traditional hydrodynamic simulation that only considers tide-level-driven changes can no longer meet the requirements. Therefore, to clarify the influence mechanism of the reclamation project on the hydrodynamic environment of the surrounding waters, this paper adopts the MIKE21 wave-current coupled hydrodynamic mathematical model to conduct numerical simulations of the hydrodynamic fields in the adjacent sea areas before and after the project. Further analysis is performed on the impacts and variation characteristics of the reclamation project on the hydrodynamic environment of the eastern sea area of Shantou. The research results of this paper will provide reference value for the feasibility evaluation of the construction project and the marine environmental risk management and control of the Eastern Urban Economic Belt of Shantou.

2. Materials and Methods

2.1. Research Methodology and Framework

To systematically clarify the implementation process of this study, the technical roadmap is presented as follows (Figure 2). This roadmap outlines the core steps from the definition of research objectives to the derivation of conclusions. With the core objective of analyzing the impacts of coastal reclamation projects on hydrodynamic environments, this study sequentially carries out data preprocessing, parameter setting, and model coupling and verification. Subsequently, comparative simulations of the flow fields before and after the project are conducted to extract the key hydrodynamic characteristics. Finally, the comprehensive impacts of the project on the hydrodynamic environment are summarized.

2.2. Model Introduction

2.2.1. Hydrodynamic Model

The two-dimensional unsteady shallow water equations in the Cartesian coordinate system are as follows[19]:

The continuous equation is

$${\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial h\overline{u}}{\partial x}}+{\displaystyle \frac{\partial h\overline{v}}{\partial y}}=hS$$

(1)

The momentum equation is

$${\displaystyle \frac{\partial h\overline{u}}{\partial t}}+{\displaystyle \frac{\partial h{\overline{u}}^2}{\partial x}}+{\displaystyle \frac{\partial h\overline{v}\overline{u}}{\partial y}}=f\overline{v}h-gh{\displaystyle \frac{\partial \eta }{\partial x}}-{\displaystyle \frac{1}{{\rho }_0}}h{\displaystyle \frac{\partial P_a}{\partial x}}+A_x+h{u}_{s}S$$

(2)

$${\displaystyle \frac{\partial h\overline{v}}{\partial t}}+{\displaystyle \frac{\partial h{\overline{v}}^2}{\partial y}}+{\displaystyle \frac{\partial h\overline{u}\overline{v}}{\partial x}}=-f\overline{u}h-gh{\displaystyle \frac{\partial \eta }{\partial y}}-{\displaystyle \frac{1}{{\rho }_0}}h{\displaystyle \frac{\partial P_a}{\partial y}}+A_y+h{v}_{s}S$$

(3)

In Equations (1) to (3), $t$represents time; $x$and$y$represent the horizontal and vertical coordinates respectively, $d$represents static water depth, $h=\eta +d$represents total water depth; $\eta$ is water level; $u$,$v$represent the components of the flow velocity in the$x$and$y$directions, respectively; ${\rho }_0$is the reference water density; $P_a$ represents local atmospheric pressure; $f=2\Omega \sin \varphi$ is parameter of Corioli ($\Omega$is the Earth’s rotation rate, $\varphi$is geographic latitude); $f\overline{v}$and$f\overline{u}$represent the acceleration caused by Earth’s rotation; $S$is source and sink terms, $u_s$, $v_s$are the flow velocity of the source sink term. The horizontal line represents the average value of depth.

The stress terms are provided by the stress model, including horizontal viscous stress, surface wind stress, bottom shear stress, and wave radiation stress. The corresponding equation is

$$A_x=-{\displaystyle \frac{1}{{\rho }_0}}\left({\tau }_{bx}-{\tau }_{sx}+{\displaystyle \frac{\partial S_{xx}}{\partial x}}+{\displaystyle \frac{\partial S_{xy}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial x}}\left(h{T}_{xx}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(h{T}_{xy}\right)$$

(4)

$$A_y=-{\displaystyle \frac{1}{{\rho }_0}}\left({\tau }_{by}-{\tau }_{sy}+{\displaystyle \frac{\partial S_{yx}}{\partial x}}+{\displaystyle \frac{\partial S_{yy}}{\partial y}}\right)+{\displaystyle \frac{\partial }{\partial x}}\left(h{T}_{xy}\right)+{\displaystyle \frac{\partial }{\partial y}}\left(h{T}_{yy}\right)$$

(5)

• Horizontal viscous stress

$T_{xx}$, $T_{xy}$ and $T_{yy}$ represent horizontal viscous stress (N/m²), mainly caused by viscous resistance, turbulent resistance, etc.

2.

Surface wind stress vector

${\overline{\tau }}_s$ represents surface wind stress vector, ${\overline{\tau }}_s=\left({\tau }_{sx},{\tau }_{sy}\right)={\rho }_a{c}_d\left|{\overline{u}}_w\right|{\overline{u}}_w$, ${\rho }_a$ is atmospheric density, ${\overline{u}}_w=\left(u_w,v_w\right)$ shows the wind speed vector at 10 m above the sea surface; $c_d$ is drag force coefficient of wind.

3.

Bottom shear stress

${\overline{\tau }}_b$ represents the bottom shear stress caused by water flow, which can be defined according to Newton’s law of friction as:

$${\overline{\tau }}_b=\left({\tau }_{bx},{\tau }_{by}\right)={\rho }_0{c}_f\left|{\overline{u}}_b\right|{\overline{u}}_b$$

(6)

In the formula, ${\rho }_0$ represents the density of water, ${\overline{u}}_b=\left(\overline{u},\overline{v}\right)$ represents the average velocity of depth, $c_f={\displaystyle \frac{g}{C^2}}$ represents the drag coefficient of water flow, and $C$ stands for Chézy coefficient.

4.

Radiation stress

$S_{xx}$, $S_{xy}$, $S_{yx}$, and $S_{yy}$ represent wave radiation stress component provided by wave models. The expression of radiation stress component is:

$$S_{xx}=E\left[\left(2n-1/2\right)-n{\sin }^2\theta \right]$$

(7)

$$S_{xy}=S_{yx}=En\sin \theta \cos \theta$$

(8)

$$S_{yy}=E\left[\left(n-1/2\right)+n{\sin }^2\theta \right]$$

(9)

where the wave energy $E=\rho g{H}^2/8$, $H$ represents wave height; $n=c_g/c$, where $c_g$ is wave group velocity and $c$ is wave velocity; $\theta$ is wave direction.

2.2.2. Wave Model

In the hydrodynamic model, the effect of waves is realized by coupling the radiation stress field, which is provided by the large-scale wave field calculated by the MIKE 21 SW wave model. The SW model is a wind-wave spectral model based on the principle of wave energy conservation, which can account for factors such as wave diffraction, refraction, bottom friction dissipation, whitecapping dissipation, wind energy input, and wave breaking. Its governing equations are as follows [20]:

$${\displaystyle \frac{\partial N}{\partial t}}+\nabla \cdot \left(\overrightarrow{v}N\right)={\displaystyle \frac{S}{\sigma }}$$

(10)

$$\left(c_x,c_y\right)={\displaystyle \frac{d\overrightarrow{x}}{dt}}=\overrightarrow{c_g}+\overrightarrow{U}$$

(11)

$$c_{\sigma }={\displaystyle \frac{d\sigma }{dt}}={\displaystyle \frac{\partial \sigma }{\partial d}}\left[{\displaystyle \frac{\partial d}{\partial t}}+\overrightarrow{U}\cdot {\nabla }_{\overrightarrow{x}}d\right]-c_g\overrightarrow{k}\cdot {\displaystyle \frac{\partial \overrightarrow{U}}{\partial s}}$$

(12)

$$c_{\theta }={\displaystyle \frac{d\theta }{dt}}={\displaystyle \frac{1}{k}}\left[{\displaystyle \frac{\partial \sigma }{\partial d}}{\displaystyle \frac{\partial d}{\partial m}}+\overrightarrow{k}\cdot {\displaystyle \frac{\partial \overrightarrow{U}}{\partial m}}\right]$$

(13)

where: N is the action density, t is time, $\overrightarrow{x}$ denotes the Cartesian coordinate system, $\overrightarrow{v}\left(c_x,c_y,c_{\sigma },c_{\theta }\right)$ is the wave group velocity, $\nabla$ is the differential operator, s is the wave propagation direction, $\theta$ and m are the directions perpendicular to s, ${\nabla }_{\overrightarrow{x}}$ is the two-dimensional differential operator in $\overrightarrow{x}$ space, and $S$ is the source term in the energy balance equation, which includes wind energy input, nonlinear wave-wave interactions, whitecapping, bottom friction, and breaking dissipation terms.

2.3. Model Construction

2.3.1. Selection of Calculation Parameters

• Calculation parameter selection

The average elevation difference between the sea level in the engineering sea area and the theoretical lowest tidal level is approximately 1.3 meters. To ensure the unification of the model calculation datum, all nautical chart water depths used shall be revised to the mean sea level.

2.

Computing domain and grid design

The calculation adopts a nested grid method, that is, for the entire northern South China Sea area, a larger grid is used for wind and wave field simulation. In the engineering area and nearby waters, a high-resolution grid is implemented, and the large-area simulation wind and wave results are used as open boundary conditions to calculate the wave current interaction and sediment movement process in the engineering area. The calculation mode adopts an unstructured triangular grid, with a total of 78,230 calculation nodes and 147,518 elements deployed during the validation stage. The minimum spatial step is about 3.15 m (for the breakwater section), and the calculation time step is adaptively adjusted. The minimum time step is 0.03 s, and the average time step is 0.06 s. The grid division for the large-scale and engineering areas is shown in Figure 3 and Figure 4.

3.

Boundary

The shore boundary adopts sliding and non-entry conditions, and the water boundary near the land is determined based on the measured tide level value. At the river mouth, this is determined based on the flow rate and input flux of runoff and sediment. The outer seawater boundary is calculated using the Global Tidal Model (TPXO10) [21], which establishes astronomical tide levels through 13 tidal components, including tidal components M2, S2, N2, K2, K1, O1, P1, Q1, Mf, Mm, M4, MS4, MN4, and is generally able to construct real astronomical tidal processes in deep waters of the open sea:

$${\xi }_0\left(x\right)={\zeta }_p\left(x\right)+{\displaystyle \sum _{i=1}^{13}{A}_i\left(x\right)}\cdot \sin \left({\omega }_{i}t+{\alpha }_i\left(x\right)\right)$$

(14)

where ${\xi }_0$ represents the tide level at the boundary, ${\xi }_p$ represents the static pressure water level at the boundary, i is in the range of 1 to 13, corresponding to the above tidal components, $A_i$ and ${\alpha }_i$ represent the amplitude and delay angle of the tidal component at the three boundaries, respectively, and ${\omega }_i$ is the angular frequency of tidal components. During the model calculation and debugging process, real-time adjustments are made based on the measured tide level values of some hydrological observation stations for the best possible fitting of the tide level hydrograph.

The wind stress calculation in the model comprehensively considers the wind field model calculation results as well as the measured wind fields at Shantou and Yun’ao Ocean Stations. The wave action is achieved by coupling the radiation stress field, and the wave boundary conditions are provided by the large-scale wave field calculated by the SWAN model.

2.3.2. Model Validation

The hydrodynamic model verification includes the comparison of tide levels, vertical average flow velocity, and flow direction of four voyages from September 2012 to January 2015, including spring-tide, mid-tide, and neap-tide stations. Here, taking Voyage 4 as an example, the results are shown in Figure 5 and Figure 6. It can be seen that the measured and calculated values of the tide level are in good agreement, and the peak and trough values of tidal level and the trend of flood and ebb tides are consistent. The values of tidal flow velocity and direction are similar to the measured data, and the variation process is also basically consistent, indicating an optimal verification effect. The established two-dimensional hydrodynamic mathematical model can better reflect the spatiotemporal distribution of the tidal current field in the engineering area and can be used for analyzing the changes in flow field in response to the project.

The wave elements (including effective wave height H_1/3 and effective wave period T_1/3) observed by ocean station and buoy were used to verify the calculation results of the wave model. The positions of verification points are shown in Figure 7, and the verification results are shown in Figure 8 and Figure 9. According to the comparison results, it can be seen that the results of simulation calculation are basically consistent with the measured data, which can meet the needs of subsequent model calculation.

3. Results and Discussion

This section takes Voyage 4 as an example to analyze the changes in tide levels and flow fields in the construction area before and after the reclamation project, as shown in Figure 10 and Figure 11. At this point, the shoreline of the reclamation area has stabilized and the seabed terrain is in a natural adjustment state. The coastline, terrain and calculation conditions used in the simulation process are as follows:

Pre-project simulation: The pre-project coastline, topographic data and the model calculation conditions of Survey Cruise 4 were adopted.

Post-project simulation: The coastline, topographic data and the model calculation conditions of Survey Cruise 4 for the project were adopted.

3.1. Tide Level Changes Due to the Reclamation Project

After the reclamation project in the enclosed and reclaimed area, the outward extension of the reclamation area to the open sea has occupied part of the original tidal water area of the sea area. Meanwhile, the runoff discharged from the estuaries of the Xinjin River and Waisha River, whose outflow directions have been altered, is subject to the backwater effect of the flood tide. As a result, the high tidal levels in the front of the Xinjin and Xinxi Districts have generally risen by 1 to 2 cm, and those on the east side of the sand-retaining dyke have generally risen by 2 to 3 cm. The topographic changes in the front of the Tagangwei District before and after the project are drastic, with the water depth increased by approximately 6 m compared with that before the project. Despite the supplement of runoff discharged from the nearby rivers during the flood tide, a water drawdown phenomenon still occurs, with the high tidal level reduced by 2~4 cm, as shown in Figure 12. During the ebb tide, due to the straight coastline of the reclamation area and the direct connection of the outflow river channels to the open sea, the ebb tide current velocity is faster than that before the project. Under the condition of unchanged ebb tide discharge, the lowest tidal level in the open sea outside the reclamation area is 1~2 cm lower than that before the project, as shown in Figure 13.

After the reclamation project, the high tide level of the Shantou Port waterway generally increased by 2~4 cm, while the low tide level generally decreased by 2~4 cm. Due to the greater changes in high and low tide levels compared to the water area on the right side of the sand barrier, and the fact that the Shantou Port waterway underwent channel excavation activities during the construction of the reclamation area, the channel dredging activities are the main reason for the significant changes in high and low tide levels of the Shantou Port waterway after the reclamation.

3.2. Flow Field Changes After Engineering

3.2.1. Changes of Flood Peak Current

Before the project, the estuaries of the Xinjin River and Waisha River were bell-shaped with a south or southwestward trend. A large number of sandy deposits such as sand bars and sand spits were developed in the vicinity of the two estuaries, with the most prominent development observed at the Waisha River Estuary. During the flood tide process, obstructed by the coastal sand bars and the sand bar at the entrance of the Waisha River Estuary, a portion of the flood current flowed counterclockwise along the coastline toward Jinshazui to the southwest. This current then converged with the flood current moving northwestward along the sand bar at the entrance of the Xinjin River Estuary before entering the Xinjin River together. For the Waisha River Estuary, the flood current was blocked by a massive triangular sand bar at its entrance and split into two branches, which bypassed the sand bar and flowed into the Waisha River channel. Meanwhile, the presence of the sand bar in this area formed a low-velocity zone with a flow velocity of 0.1~0.2 m/s and a northwestward flow direction between the sand bar on the left bank of the Waisha River Estuary and the Laiwu Peninsula, as shown in Figure 14.

After the project completion, the sand bars at the estuaries of the Xinjin River and Waisha River were excavated, and the two estuaries were diverted southeastward. The reclamation area protrudes seaward, forming a regular coastline with a northeast-southwest strike. During flood tide, the flood current now flows northwestward along the straight coastline into the Xinjin and Waisha Rivers, with the flow velocity at the estuary entrances reaching 0.2~0.3 m/s. Owing to the frictional resistance of the seabed, the flood current velocity in the nearshore zone has been attenuated, forming three low-velocity zones with a flow velocity of less than 0.1 m/s. These zones are respectively the newly formed sea area from the right bank of the Xinjin River Estuary to the sand barrier of Shantou Port, and the small water areas within the arcuate coastlines of the Xinxi and Tagangwei Districts. After the project, the flow velocity in the sea area from the Waisha River Estuary to the Laiwu Peninsula has increased from 0.1~0.2 m/s to 0.2–0.3 m/s, and the flow direction has also shifted from northwestward to northeastward, as shown in Figure 15.

3.2.2. Changes of Ebb Peak Current

Before the project, the ebb currents from the Xinjin and Waisha rivers spread outward at the open estuary entrances with a flow velocity of 0.4~0.5 m/s, which decreased to 0.1~0.2 m/s after the project. The ebb current of the Xinjin River flowed southeastward along the sand barrier into the open sea, whereas the ebb current of the Waisha River was blocked by the sand bar at its estuary entrance and split into two branches that flowed southward and southeastward into the sea, respectively. The ebb discharge from the two estuaries experienced a gradual velocity increase as it flowed seaward, reaching 0.3~0.4 m/s at the head of the sand barrier. Water bodies with a flow velocity of less than 0.2 m/s were primarily distributed in the nearshore waters between the Xinjin and Waisha rivers, as well as in the small embayments between the Waisha River and the Laiwu Peninsula. In contrast, water bodies with a flow velocity exceeding 0.5 m/s were concentrated in the tidal channels between the Xinjin River Estuary and the sand bars at the Waisha River Estuary, as shown in Figure 16.

After the project, the sand bars at the estuaries of the Xinjin and Waisha rivers were excavated, and the newly formed northeast-southwest coastline was nearly parallel to the ebb current direction. The ebb current velocity at the estuary entrances ranged from 0.2~0.3 m/s, an increase of approximately 0.1 m/s compared with that before the project. The ebb current now flows along the straight coastline of the Tagang, Xinxi and Xinjin Districts toward the sand barrier of Shantou Port, with a nearshore flow velocity of 0.2~0.3 m/s. The velocity diminishes to 0.1~0.2 m/s in the triangular area between the Xinjin District and the root of the Shantou Port sand barrier, after which the current flows southeastward along the sand barrier toward the open sea and gradually accelerates to 0.3~0.4 m/s. The ebb current velocity in the open sea reaches 0.4~0.6 m/s, with the flow direction gradually shifting from southwestward to southward. The ebb current velocity in the Shantou Port channel is approximately 0.6 m/s, an increase of about 0.1 m/s relative to the pre-project value. In the reclamation area, the nearshore ebb current velocity has decreased by approximately 0.03~0.05 m/s compared with the pre-project condition, equivalent to a reduction of about 25%. Low-velocity zones with a flow velocity of less than 0.1 m/s emerge in the triangular area between the coastline and the sand barrier in the Xinjin District, as well as in the small water body at the northeast corner of the Tagang reclamation area. Owing to the project, the ebb current velocity in the southeastern waters off the Laiwu Peninsula has increased by approximately 0.1 m/s with the flow direction remaining southwestward, as shown in Figure 17.

3.3. Wave Field Changes After Engineering

The wave field variations in the project area before and after reclamation are closely associated with the changes in coastline morphology, seabed topography, and the presence of sand barriers/shoals, with key characteristics analyzed as follows (Figure 18 and Figure 19):

3.3.1. Pre-Project Wave Field Characteristics

Before the project, the coastline in the study area was a naturally formed horn-shaped coastline of the Hanjiang Delta. A large number of sandbars and barrier shoals were developed in the front of the Xinjin River and Waisha River estuaries, resulting in an overall shallow water depth (corresponding to the isobath labels in Figure 18). When offshore waves propagated into the project area, their directions deflected toward the northeast-southwest extension of the natural coastline, and the wave direction lines tended to be perpendicular to the isobaths or coastline.

In the offshore waters far from the shore, waves transmitted energy normally without obstruction from shoals, leading to a gradual decrease in wave height as the distance to the shore shortened. In the nearshore areas, however, waves were affected by the flow constriction effect of the horn-shaped coastline and the blocking of barrier shoals. This caused the waves to undergo obvious refraction and breaking, resulting in rapid energy attenuation and a significant reduction in wave height. Specifically, the significant wave height in the nearshore estuary areas was only 0.0–0.2 m, while that in the offshore areas reached 0.4–0.8 m, as shown in Figure 18.

3.3.2. Post-Project Wave Field Characteristics

After the reclamation project, the reclamation areas (Xinjin, Xinxi, and Tagangwei areas) advanced 1.5–2.4 km seaward, and the coastline was transformed from a natural horn shape to a regular straight artificial coastline with a northeast-southwest strike. The water depth in the front of the reclamation areas increased compared with the pre-project period, and local sand extraction pits were formed in the areas near the Xinjin River and Waisha River estuaries due to sand mining from nearby sand sources (corresponding to the deep-water areas in Figure 19).

When offshore waves propagated into the project area after reclamation, the overall variation trend of wave height was similar to that before the project: waves in offshore areas transmitted energy normally, with wave height decreasing gradually as the distance to the shore shortened. In the nearshore areas, the increase in water depth of the sand extraction pits reduced wave energy loss compared with the pre-project period, making the reduction trend of wave height more gentle.

Notably, the front area of the reclamation coastline lacked the blocking effect of barrier shoals on waves. Combined with the drastic local topographic changes caused by large-scale sand extraction pits, waves refracted significantly toward both sides of the “deep pits” of the sand extraction areas. This led to obvious wave energy concentration in the area from the slope top of the sand extraction pits to the front of the reclamation coastline: the significant wave height increased by approximately 0.1–0.3 m compared with the pre-project period, with the most prominent increase observed in the estuary areas of the Xinjin River and Waisha River (Figure 19). In these estuary areas, the significant wave height increased from 0.0–0.2 m (pre-project) to 0.2–0.4 m (post-project), which was mainly attributed to the combined effects of reduced wave energy dissipation by sand extraction pits and enhanced wave refraction by artificial straight coastlines.

4. Conclusions

A wave-current coupled mathematical model was adopted in this paper to simulate the hydrodynamic variations in the coastal waters surrounding the reclamation project of the Eastern Urban Economic Belt in Shantou, both before and after the project execution. The results demonstrate that after the implementation of the reclamation project, the tidal level changes in the project sea area were negligible, with a variation range of merely 2~4 cm. During the flood tide, the flood current direction was northeastward across most areas, except for the estuaries of the Xinjin River and Waisha River, where the flow direction was northwestward. The current velocity in the waters ranging from the outer area of the Waisha River Estuary to Laiwu Peninsula increased to 0.2~0.3 m/s, and the flow direction shifted from northwestward to northeastward. This indicates that the reclamation project not only enhanced the current velocity in this region, but also altered the flood-tide sediment transport direction prevailing prior to the project.

Due to the dredging and deepening of the Shantou Port waterway, the current velocity in a small area west of the sand barrier head slightly decreased compared with the pre-project level. Nevertheless, the average current velocity in the channel center was maintained at 0.4~0.5 m/s, with no significant change in flow direction. During the ebb tide, the ebb current velocity at the estuaries of the Xinjin River and Waisha River increased relative to the pre-project conditions, while the along-shore ebb current velocity along the coastline of the reclamation area decreased. In the southeastern waters of Laiwu Island, the current velocity increased with the flow direction remaining unchanged, which suggests that the reclamation project exerted a velocity-enhancing effect on this region without altering the pre-project ebb-tide sediment transport direction. The current velocity and direction in the shoal area east of the Shantou Port waterway showed minimal differences from the pre-project state, indicating that the reclamation project had little impact on the current velocity of the Shantou Port waterway. Instead, the variations in current velocity within the waterway were mainly induced by dredging activities.

Regarding the wave field, before the project, the naturally formed horn-shaped coastline of the Hanjiang Delta and the extensive sandbars and barrier shoals at the estuaries of the Xinjin River and Waisha River led to obvious refraction and breaking of nearshore waves, resulting in rapid energy attenuation. The significant wave height in nearshore areas was only 0.0~0.2 m, while that in offshore areas reached 0.4~0.8 m. After the reclamation, the coastline was transformed into a regular straight artificial coastline with a northeast-southwest strike, and local sand extraction pits were formed near the estuaries due to sand mining. The increase in water depth of the sand extraction pits reduced wave energy loss, making the wave height reduction trend gentler. However, the absence of barrier shoals and drastic topographic changes caused by sand extraction pits led to significant wave refraction toward both sides of the pits, resulting in obvious wave energy concentration in the front of the reclamation coastline. The significant wave height increased by approximately 0.1~0.3 m compared with the pre-project period, with the most prominent increase observed in the estuaries of the Xinjin River and Waisha River (from 0.0~0.2 m to 0.2~0.4 m).

Overall, based on a wave-current coupled model, the research results of this paper comprehensively analyze the impacts of the coastal reclamation project on the hydrodynamic environment (tidal levels, flow fields, and wave fields) of the waters surrounding the Eastern Urban Economic Belt of Shantou, providing technical support for the feasibility evaluation of the project and the marine environmental risk management and control of the Eastern Urban Economic Belt of Shantou. However, these results are based on a two-dimensional wave-current coupled mathematical model, which cannot adequately reflect the three-dimensional vertical distribution of flow velocity and wave energy. Therefore, conducting further research by adopting a three-dimensional wave-current-sediment coupled mathematical model to explore deep-seated impacts such as suspended sediment transport and seabed erosion-siltation will be the focus of future work.