Computational investigation on hydrodynamic and sediment transport responses influenced by reclamation projects in the Meizhou Bay, China

Abstract

Reclamation projects are the main method of coastal exploitation, and the hydrodynamic environmental effect, together with the sediment transport response of the reclamation project, is important to the project’s site selection and environmental protection. Herein, a 3D numerical model based on the finite volume community ocean model (FVCOM) is applied to simulate the changes in the Meizhou Bay’s hydrodynamic environment and sediment transport after a reclamation project. The reclamation project greatly alters the shape of the shoreline and narrows the bay, leading to a significant change in its hydrodynamic environment and sediment transport. After the project, the clockwise coastal residual current in the corner above the Meizhou Island gradually disappears. An obvious counter-clockwise coastal residual current emerges around the rectangular corner. The tidal prism decreases by 0.65 × 10⁹ and 0.44 × 10⁹ m³ in the spring and neap tides, respectively. The residence time presents a major increase. These changes lead to the weakening of the water exchange capacity and the reduction of the self-purification capacity of the bay. Currents in the tidal channel weaken, whilst currents in the horizontal channel strengthen. The strength and scope of particle trajectories around the horizontal channel and the Meizhou Island enhance. The suspended sediment concentration (SSC) increases in the majority of the Meizhou Bay but decreases in the lateral bay. The eastern corner of Z2 shows a tendency to erode. The western region of the Meizhou Island, the upper portion of the rectangular corner and the western corner of Z4 show a tendency to deposit. The reclamation project increases the maximum storm surges by 0.06 m and decreases the maximum significant wave heights by 0.09 m.

Keywords

Meizhou Bay / FVCOM+SWAN / reclamation project / hydrodynamic environment / SSC / typhoon

Introduction

Meizhou Bay is located in China’s south-eastern coast and covers an area of approximately 458 km². It has a shallow water system including a 200 km² area of mudflat with an average depth of 1.5 m; along its tidal channel, the maximum depth can reach up to 50 m, causing strong tidal currents (Fig. 1). The tide in the bay is classified as an anomalistic semidiurnal tide, and it is largely controlled by the M₂ harmonic constituent. The perennial average tidal ranges vary from 4.65 m at the entrance of the bay to 6.52 m at the top of the bay. Connecting to the Taiwan Strait, Meizhou Bay is of considerable significance to the city’s political affairs and economy. Similar to other coastal cities, the innate shortage of land has constrained the Meizhou’s economic development and population growth. To acquire more land resources and stimulate its local development, the Meizhou government has carried out a large number of reclamation projects. However, reclamation projects considerably narrow the estuary and cause obvious shoreline relocation and seabed deformation, presumably impacting the economic sustainable development and the natural ecosystem equilibrium of the bay. Meanwhile, vulnerability to extreme conditions, such as typhoons, is another inherent challenge to coastal cities. Every summer and autumn, Meizhou Bay suffers enormous losses from storm surge, strong wind, and inland flooding (^{Xu et al., 2015}). When a cyclone generated above the Pacific Ocean moves toward the south-eastern region of China, its strength gradually intensifies and it eventually becomes upgraded into a typhoon before landfall. As the typhoon approaches the estuary mouth, the low barometric pressure at the typhoon center and the strong surface wind push a large volume of seawater into the estuary, causing remarkable increase in water level and great threats to the bay (^{Chen et al., 2012}).

The impacts of reclamation projects and typhoons have been the focus of recent research. ^{Guo et al. (2009}, ²⁰¹²) established a 3D model based on the FVCOM to reproduce the storm surges generated by Typhoon Agnes and studied the effects of reclamation projects and cyclonic parameters on the hydrodynamics of the Hangzhou Bay. ^{Liu and Huang (2009)} modelled sediment resuspension and transport induced by storm winds in Apalachicola Bay; the simulation results demonstrated that storms could cause intense resuspension and transport activities. Research on numerical simulations in the Meizhou Bay are not abundant and are largely limited to tidal process and barotropic mode. ^{Liu et al. (2009)} applied a 2D numerical model to simulate and analyze the characteristics of tidal currents in the Meizhou Bay; results indicated that the 50% and 80% water exchange periods in the bay are 5 and 15 days, respectively. ^{Zhu et al. (2014)} applied a finite volume model with unstructured grids to simulate the distributions and variations of the temperature field in the Meizhou Bay; in this case, however, the salinity process was ignored. Simulations of typhoon and typhoon-induced storm surges and waves in the Meizhou Bay are rare. Most typhoon simulations are conducted over a large domain. ^{Xu et al. (2015)} simulated the characteristics of storms on the coastal region of Jiangsu Province with a coupled current and wave model; the results of this work showed that the maximum storm surges would occur in the south-eastern coast. Sediment simulation in the Meizhou Bay is another research hotspot. However, studies related to this issue largely deal with the analysis of sediment measurement, and the numerical simulations of sediment transport are limited. ^{Chen (1988)} pointed out that the sediment resources in the Meizhou Bay are the transport process from the ocean and the suspension process above the mudflat.

A 3D numerical model that considers current, baroclinic process, wave, and sediment within the Meizhou Bay is necessary to provide a comprehensive simulation of the hydrodynamic environment and sediment transport. This study aims to: (i) analyze the characteristics of the hydrodynamic environment and sediment transport, including particle trajectory, tidal prism, residence time, residual current, typhoon-induced storm surges and waves, the advection-diffusion process, and seabed deformation in the Meizhou Bay; and (ii) analyze how reclamation projects affect these characteristics. The conclusions of this paper can provide a basis for decisions on engineering design, environmental protection, pollutant emission control, and plans for dredging projects.

Numerical model description

The FVCOM model developed by ^{Chen et al. (2003}, ²⁰⁰⁷), which has been widely applied to hydrodynamic simulations, is used in this study. It is based on 3D primitive equations and uses the finite volume approach to ensure a better conservation of mass and momentum. It also uses a sigma coordinate in the vertical direction and a triangular mesh in the horizontal plane, thereby enabling it to capture rough seabed shapes and meandering shorelines. The Mellor-Yamada 2.5 turbulence closure (^{Mellor and Yamada., 1982}) is applied to the FVCOM to calculate the coefficient of vertical mixing. The mode separation method largely increases the computational efficiency, and the point drying/wetting treatment greatly enhances the computational stability of the model.

The third-generation wave model SWAN developed at the Delft University of Technology is used to compute random, short-crested, wind-generated waves in coastal regions; it is based on the discrete spectral action balance equation (^{Hasselmann et al., 1973}). Transformation from the primitive structured grid version to an unstructured grid version has been made, and the SWAN model has been fully coupled to the FVCOM. The fully coupled current-wave-sediment module in the FVCOM plays an important role in simulating hydrodynamic environment and sediment transport under various conditions, such as tide, current, and typhoon process.

Larger model setup and verification

To enable a better simulation of the hydrodynamic environment of the Meizhou Bay, this study implements a two-nested grid method in which the larger grid model provides the information of interest, such as water level, salinity, and temperature to the smaller model. Therefore, variable motivations along the entrance (the open boundary) of the Meizhou Bay, as well as the initial field for salinity and temperature, can be obtained. In addition, simulation of the typhoon process also requires a larger domain covering the majority of the typhoon-affected area so that the boundary wave input can be neglected. Consequently, a larger model is adopted within the Taiwan Strait to meet these considerations.

Model domain and mesh

The computational domain covers the vast majority of the Taiwan Strait, which has an extensive span of 115–120 °E in longitude and is confined at the north by 26 °N, the south by 22 °N. The triangular grids in the horizontal plane have 2203 nodes and 4406 elements, with an extensive grid size of 5000 m in the open boundary and a narrow grid size of approximately 100 m in the low-lying areas (Fig. 2). Three levels are applied in the vertical direction, with a maximum depth of 500 m in the southern boundary and a minimum depth of 0.1 m around the mudflat. The topography in the majority of the computational domain is obtained based on the global topography database ETOPO2, whilst the topography within the Meizhou Bay is acquired from high-resolution field measurements.

Initial and boundary condition

The model starts with the mean sea level and zero flow velocity. The coastlines are regarded as closing boundaries with no-flow condition. Zhangzhou, Jinjiang, and Hanjiang are set as river inputs, with their monthly average discharges applied based on the multi-year database analysis. Along the open boundary, tidal level forcing is adopted by combining four main tidal constituents with four minor constituents. Water temperature and salinity along the open boundary are also included in the model, which are obtained based on the marine atlas and the HYbrid Coordinate Ocean Model (HYCOM) database; the atlas provides a monthly average value, and the HYCOM database provides a 3-hourly variation tendency. Atmospheric forcing, including wind speed, surface net flux, evaporation, and pressure perturbation, is imposed on the model as sea-surface conditions. Data for these parameters are acquired from the European Centre for Medium-Range Weather Forecasts database with a resolution of 0.75° × 0.75°.

Parameter settings

Simulation begins at 00:00 (Beijing time) on 1 January 2005 and ends at 00:00 on 31 December 2006, with a 5 s time step for the internal mode and a 60 s time step for the external mode. Bottom roughness coefficient applied to this model is considered as a function of depth and topography gradient. The basic coefficients for the minimum and maximum depths are 0.03 and 0.005, respectively, and the high-order coefficients for the minimum and maximum topography gradients are 0.0 and 0.01, respectively. Therefore, the total coefficient is equal to the basic coefficient plus the high-order coefficient; these two coefficients are obtained from the linear interpolation among the ranges mentioned above according to the local depth and its local gradient. The point drying/wetting treatment scheme is applied with a minimum depth of 0.05 m; thus, when the depth of an element is lower than 0.05 m, it is locked and is not considered into calculation.

Model verification

To evaluate the accuracy of the results, comparisons are made between the simulated results and the observed data. Two stations, Meizhou and Shantou, are chosen for comparisons of tidal level according to the availability of data from 00:00 on October 6, 2005 to 00:00 on November 5, 2005. As shown in Fig. 3, the simulated tidal levels are in good agreement with the real-time data in terms of tidal phase and range. The root mean square errors (RMSEs) between the simulated results and the observations are 0.11 and 0.13 m at the Meizhou and Shantou stations, respectively, and the correlation coefficients for these stations are 0.95 and 0.93, respectively. Table 1 presents the errors of the amplitudes and phases of the four main tidal constituents at various stations. The maximum absolute mean errors of the amplitudes of M₂, S₂, K_1, and O₁ are 8.3, 5.4, 5.2, and 2.5 cm, respectively, which indicates that the model can capture the main characteristics of tidal process in the Taiwan Strait.

The Shantou, Meizhou, and Gaoxiong stations are chosen to compare the observed salinity and temperature with the simulated results. As Fig. 4 shows, the modelling results are in good agreement with the monthly average data. The RMSEs for salinity at the Shantou, Meizhou, and Gaoxiong stations are 0.5, 0.3, and 0.6 PSU, respectively, and their correlation coefficients are 0.90, 0.92, and 0.89, respectively. The Shantou and Gaoxiong stations present a similar salinity development pattern, i.e. lower salinity in summer due to frequent rainfall and higher salinity in winter due to a small amount of freshwater supplementation. However, the Meizhou station shows an opposite tendency due to the strong summer monsoon, which brings high-salt sea water into the bay. The RMSEs for temperature at these stations are 0.8°C, 1.2°C, and 0.3°C, respectively, and their correlation coefficients are 0.91, 0.85, and 0.93, respectively. All stations reveal a similar temperature development mode, i.e. lower in winter due to the long distance between the earth and sun and higher in summer because of their relative proximity. Figure 5 shows the spatial distribution of the salinity and temperature around the Meizhou Bay and the Xiamen Bay in summer and winter. In summer, the temperature is higher around the estuary but lower in the outer sea caused by the high-temperature water supplementation from the Zhangzhou River. By contrast, in winter, the temperature in the estuary is lower, whilst that in the ocean is higher. As for salinity, given the continuous supplementation of fresh water from river, both seasons show similar distributions.

In conclusion, the characteristics of the tidal process and the salinity and temperature development modes in the Taiwan Strait, especially around the Meizhou Bay, are in good agreement with the observed data and empirical theory. Thus, the larger model can provide reliable information to the smaller model and lay the foundation for accurate typhoon simulation.

Smaller model setup and verification

Model setup

The computational domain covers the entire portion of the Meizhou Bay, with a small span of 118.8–19.4°E in longitude and is confined at the north by 25.4°N, the south by 24.9°N. Triangular grids in the horizontal plane have 4096 nodes and 8871 elements, with an extensive grid size of 800 m along the open boundary and a narrow grid size of approximately 30 m at the top of the bay (Fig. 6). Ten levels are used in the vertical direction, with a maximum depth of 54 m in the tidal channel and a minimum depth of 0.1 m above the mudflat.

Similar to the larger model, the smaller model adopts the cold-start method, together with the initial temperature and salinity fields. The shoreline and island are regarded as solid boundaries with no-flow condition. Along the open boundary, tidal level and baroclinic motivations are adopted. Atmospheric forcing is imposed on the model by the sea-surface conditions. All of the information required, such as tidal level, salinity, and temperature, are obtained from the interpolation of the simulated results in the larger model.

Configurations for the sediment model are described as follows: Under normal circumstances, SSC slightly fluctuates around 0.05 and 0.03 g/L during the spring and neap tides within the Meizhou Bay, respectively. Compared with the suspended-load movement, the bed-load movement is relatively small; thus, only suspended sediment transport is considered here (^{Guo et al., 2014}). Meanwhile, the cohesive sediment type is chosen in the sediment model (^{Hou et al., 2017}). The main supplies of sediment are the suspension process above the mudflat and the replenishment from the ocean (^{Chen, 1988}). The suspension process is realized with a critical shear stress about 0.13 N/m² and a settling velocity about 0.007 m/s (^{Luo et al., 2007}). The open boundary motivation of SSC is set as a time-variable series changing from 0.03 to 0.08 g/L based on the observations.

Simulation begins at 00:00 on 1 April 2006 and ends at 00:00 on 8 May 2006, with a 5 s time step for the internal mode and a 60 s time step for the external mode. The bottom roughness coefficients are considered as spatially variable values changing from 0.035 to 0.039 based on the different level of local depth. For depths above 30, 30–10 and below 10 m, the bottom roughness coefficients are 0.035, 0.037 and 0.039, respectively. The minimum depth is 0.05 m as the larger model does.

Model verification

Tidal level, current, salinity, temperature, and SSC are used for verification over the period from 00:00 on April 28, 2006 to 23:00 on April 29, 2006. Figure 1 presents the locations of the observation stations. T1, T2 and T3 are the tidal level stations, whilst C1–C6 are the current, salinity, temperature, and SSC stations.

Verification of tidal level and current

As shown in Fig. 7, the simulated tidal levels are in good agreement with the observed data in terms of tidal phase and range. The RMSEs for tidal level at T1, T2, and T3 are 0.07, 0.09 and 0.11 m, respectively, and their correlation coefficients are 0.94, 0.91 and 0.89, respectively. The tidal range increases from the entrance to the top of the bay owing to the shoreline convergence effect. Meanwhile, the simulated results of current speed and direction agree well with the observed data. During the flood tide stage, the current runs into the bay with an average direction of about 200 °C along the tidal channel. During the ebb stage, however, the current runs out of the bay with an average direction of about 300 °C along the tidal channel.

Verification of temperature and salinity

Given the limited volume of runoff and other fresh water supplements, the temporal and spatial variations in temperature and salinity demonstrate a strong tidal characteristic. As shown in Fig. 8, alteration processes of temperature and salinity are similar to that presented in the tidal level process except for their amplitudes. During a tidal period, temperature and salinity change with the flood and ebb currents. For temperature, the RMSEs for it at the surface, middle, and bottom of C3 are 0.13°C, 0.12°C, and 0.09°C, respectively, and their correlation coefficients are 0.90, 0.91, and 0.94, respectively. For salinity, the RMSEs for it at the surface, middle, and bottom of C3 are 0.12, 0.13, and 0.21 PSU, respectively, and their correlation coefficients are 0.92, 0.93, and 0.85, respectively. Salinity clearly fluctuates more intensely than temperature, which may result from the moderate temperature fluctuation along the open boundary obtained from the larger model.

Verification of SSC

As shown in Fig. 9, the simulated SSC generally agree with the observed data, and it is positively correlated with the tidal current velocity. Figure 10 presents the average distribution of SSC during the spring and neap tides. During the spring tide, the stronger tidal current and the higher bottom shear stress strengthen the suspension and transport processes, thus leading to a higher SSC. However, during the neap tide, the average SSC is lower due to the weak suspension and transport processes. Tidal channel presents various SSC distributions in these periods. In the spring tide, most of the channel has sufficient sediment supplementation, and the whole channel presents a higher SSC distribution; by contrast, in the neap tide, this supplementation is weakened, and the whole channel presents a lower SSC distribution. Mudflat in either of these periods presents a higher SSC distribution due to the strong erosion process at shallow water situation.

Typhoon model setup and verification

Typhoon selection

On 10 September 2010, Typhoon Fanapi (1011) formed as a tropical cyclone above the Pacific Ocean. It then moved northwest with increasing intensity and transformed into a tropical storm on September 14. On September 15, it was upgraded to a typhoon with a maximum wind speed of approximately 50 m/s and a central pressure of approximately 946 hPa. Fanapi eventually made landfall in the coastal region of the Zhangpu Town, Fujian Province, on 20 September. Information of Fanapi before and after its land fall is present in Tab. 2. Detailed information on this typhoon can be obtained from the Chinese typhoon weather website.

Typhoon model description

Applying an accurate atmospheric pressure and wind field model is crucial in typhoon simulation. To obtain an optimal model, a number of parametric equations with reasonable settings have been adopted. The radius of the maximum wind speed is calculated by the following empirical formula (^{Willoughby and Rahn, 2004}):

(1)R_{\max }=\mathrm{51.6}EXP(-\mathrm{0.0223}{V}_{\max }+\mathrm{0.0281}\varphi ）,

where R_{\max } represents the radius of maximum wind speed, \varphi represents the latitude of the typhoon center (deg), and V_{\max } represents the maximum wind speed.

A number of formulas, such as the Fujita formula (^{Fujita, 1952}), the Myers formula (^{Myers, 1957}), and the Jelesnianski formula (^{Jelesnianski, 1965}), have been developed to calculate the pressure and wind field during the typhoon period. Based on the recent research related to the typhoon simulation around the Taiwan Strait, the improved Holland-Bformula is chosen to calculate the pressure and wind field in this simulation, expressed as follows (^{Holland, 1980}, ²⁰⁰⁸):

(2)P\left(r\right)=P_c+(P_n-P_c){(-{\displaystyle \frac{R_{\max }}{r}})}^B;

V\left(r\right)=

(3)\sqrt{(P_n-P_c){\displaystyle \frac{B}{{\rho }_a}}{\left({\displaystyle \frac{R_{\max }}{r}}\right)}^{B}EXP{(-{\displaystyle \frac{R_{\max }}{r}})}^B+{\left({\displaystyle \frac{rf}{\mathrm{2}}}\right)}^{\mathrm{2}}}-{\displaystyle \frac{rf}{\mathrm{2}}},

where r is the distance from the calculation point to the typhoon center, P\left(r\right) is the pressure at a distance r from typhoon center, P_c is the cyclone center pressure, P_n is the ambient pressure with a constant value about 1013 hPa, V\left(r\right) is the velocity at a distance r from the typhoon center, f is the Coriolis force parameter, and B is an index that affects the intensity and kurtosis of a typhoon , given by the following formula (^{Vickery et al., 2000}):

(4)B=\mathrm{1.881}-\mathrm{0.00557}{R}_{\max }-\mathrm{0.01295}\varphi .

Typhoon model verification

The computational mesh used in the typhoon model is identical to that used in the large model. Simulation of Typhoon Fanapi begins at September 10, 2010 with a time step of 1 s in the coupled FVCOM-SWAVE module. Water elevations, consisting of the astronomical tide and the typhoon-induced water surface variation based on the ‘inverted barometer’ estimation, are used in the open boundary forcing. Other parameters, such as bottom roughness coefficient, are set similar to those in the larger model. As Fig. 11 shows, the simulated tidal level is slightly lower than the observation, which may be due to the differences between the symmetrical typhoon model and the actual typhoon field. The RMSEs for the tidal level at the Pingtan, Chongwu, and Xiamen stations are 0.05, 0.13, and 0.11 m, respectively. Storm surge deriving from the differences between the tidal level in the typhoon model and that in the modelling with only astronomical tide is used to show the effect of typhoon. The storm surge development mode among these stations is similar (Fig. 12). The maximum values of the surges in these stations are 41, 62, and 110 cm, respectively. Among these stations, Xiamen has the maximum surge because it is located closest to the typhoon center trace. Pingtan shows the maximum surge at the earliest time compared with other stations due to its northern location, which is where the typhoon first passed within the Taiwan Strait.

Effects of the reclamation project on the hydrodynamic environment and sediment transport

A sketch map of the reclamation project in the Meizhou Bay is shown in Fig. 13. The area of the reclamation project is about 137.32 km² and covers the majority of the bay’s mudflat. Modifications have been made to the well-verified smaller model to account for the new shape of the shoreline after the reclamation project. The smaller and modified smaller models are applied to simulate and analyze the effects of the reclamation project on the hydrodynamic environment and sediment transport of the bay. The typical spring and neap tides in summer are acquired based on the analysis on the multi-year tidal database. By applying these tides to the simulation, issues such as particle trajectory, tidal prism, residence time, and residual current can be further analyzed. A real-time simulation from 00:00 on 1 May 2006 to 00:00 on 8 May 2006 is conducted to analyze the effects of the reclamation project on the advection-diffusion process. For sediment transport analysis, long-term simulation is conducted until the stable distributions of seabed height and tidal period averaged SSC are achieved. To investigate the effect of typhoons, a 7-day simulation of the Typhoon Fanapi beginning at 00:00 on 18 September 2010 is imposed on the smaller model with the configurations and boundary conditions obtained from the typhoon model. Differences between the results of the typhoon model and those obtained in the modelling with barely astronomical tide are used to indicate the effects of a typhoon.

Effects on particle trajectory

Particle trajectory can provide vivid representations of the way the hydrodynamic field behaves. Based on the Lagrangian particle-tracing method, 48 particles are released in the Meizhou Bay driven by the typical spring and neap tides. As Fig. 14 shows, the tidal ranges of the typical spring and neap tides are 6.48 and 3.26 m, respectively. The colour alteration from blue to green indicates the time spans from 0 to 24.83 h. Given this colour alteration, the particle trajectory in Fig. 15 can also shows the time course, where the blue point is the start location, whilst the green point is the end location.

Figures 15(a) and 15(b) present the particle trajectories before the reclamation project. Particle in the tidal channel has a linear and narrow trajectory shape, which is in accordance with the rectilinear stream observed by the Meizhou Marine Office. Particle connecting to the ocean presents a curvilinear trajectory shape and can be classified as the rotary current, as proposed by the related study (^{Zhu et al., 2014}). Different tidal currents can result in various hydrodynamic fields, thus leading to the diversity of the characteristics of the particle trajectory. Under stronger tide, the motion strength and coverage of particles are enhanced. Meanwhile, particle trajectories around the Meizhou Island present a major difference. These trajectories form a circular shape around the island in the spring tide, but the circular shape disappears in the neap tide, which means the flow around the island is more intense in the spring tide than that in the neap tide.

To investigate how the reclamation project affects the particle trajectory, Figs. 15(c) and 15(d) present the particle trajectories after the reclamation project. In general, a reclamation project does not remarkably influence the particle trajectory. The interior and exterior portions of the bay remain the rectilinear and rotary streams, respectively, but the motion strength and coverage of these particles decrease. However, the strength and scope of the particle trajectories around the horizontal channel and the Meizhou Island are slightly enhanced, thus indicating that regions around these places present a stronger flow after the reclamation project. Moreover, the more compact shoreline shape narrows the trajectory along the tidal channel and changes the trajectory in the horizontal channel from a curvilinear shape into a linear one.

Effects on tidal prism

The tidal prism refers to the volume of tidal water that can be accepted by a bay, and the size of it reflects the self-purification capacity of the bay. The tidal prism can be calculated via the general formula (^{Jia et al., 2018}):

(5)W={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(S_{\mathrm{1}}+S_{\mathrm{2}})H,

where S_{\mathrm{1}} represents the water area at the time of high tide, S_{\mathrm{2}} represents the water area at the time of low tide, and H represents the tidal range. For computational simulation, this formula can be expressed as:

(6)W={\displaystyle \sum _{i=\mathrm{1}}^n{S}_i(h_{\mathrm{1}i}-h_{\mathrm{2}i})},

where S_i represents the area of elements, h_{\mathrm{1}}{}_i represents the high-tide level of the elements, h_{\mathrm{2}}{}_i represents the low-tide level of the elements, and n represents the total number of elements.

Tidal prisms in the typical spring and neap tides are calculated based on Eq. (6). The calculated results and comparisons are presented in Table 3. In general, tidal prism in the spring tide is roughly twice as much as that in the neap tide, which is mainly caused by the differences in the tidal strength (^{Liu et al., 2009}). Reclamation project decreases the tidal prism in the typical spring and neap tides by 0.65 × 10⁹ and 0.44 × 10⁹ m³, respectively. This decline means the water exchange volume between the bay and the ocean decreases, thus leading to the exacerbation of the water pollution. Compared with the spring tide, the neap tide presents a larger relative loss in the tidal prism. This implies that the water exchange capacity largely weakens during the neap tide than the spring tide.

Effects on residence time

Residence time is the time required by a virtual water parcel to escape from a given area. The particle-tracking method using Lagrangian floats is applied to calculate this timescale. To ensure computational efficiency, the main study domain is divided into five zones of interest, Z1–Z5 (Fig. 13), and particles within the zones of interest are selected for calculation. Residence time of a certain zone is calculated from the average residence times of all the particles within this designated zone. Solid boundaries, such as coastlines and islands, are treated as non-penetrable boundaries. The open boundary is regarded as a threshold signaling when the simulation should cease. When a particle reaches the open boundary, the computation stops, and the time taken by this particle is recorded as its residence time. Considering the particle may stay in the bay for a relatively long time, 30-day is selected as the upper boundary of the simulation time. The typical spring and neap tides are imposed on the open boundary to provide periodic motivation. Table 4 presents the statistics and comparisons of residence times in five zones within the Meizhou Bay.

Generally speaking, residence time is highly associated with the boundary restriction and the hydrodynamic environment. As Z1 and Z4 are narrower than other zones, the former reveal higher residence times than the latter. Z2 and Z3 are located at the tidal channel and present a small residence time owing to their powerful current velocity. Z5 features a good advection-diffusion process and is nearest to the open boundary; thus, the residence time of it is the lowest. As expected, residence time in the spring tide is generally smaller than that in the neap tide, which is probably caused by the stronger tidal current in the former. However, Z2 has a smaller residence time in the neap tide period than that in the spring tide period. After the reclamation project, residence time presents a substantial growth owing to the reduced tidal current and weakened tidal prism. Z1 presents the maximum growth, while Z5 presents the minimum growth.

Effects on residual current

Tidal movement is the main dynamic process in the Meizhou Bay, and the residual current has a significant influence on the material transport, ecology, and pollutant diffusion. Based on the time-averaged velocity components of the typical spring and neap tides, the Euler residual current distributions in the Meizhou Bay before and after the reclamation project are obtained, as shown in Fig. 16 (^{Jia et al., 2018}).

Because of the complicated shoreline shape and significant shoreline contraction of the bay, its residual current presents a relative intricate pattern. The counter-clockwise coastal residual current in the junction between Z1 and Z2 is in agreement with the finding in the related research (^{Luo et al., 2007}). Two powerful clockwise coastal residual currents form around the western part of Z2. Two obvious coastal residual currents flow in the opposite direction at the bottom of Z3, creating a low-velocity zone., A powerful current flows from east to west in the horizontal channel, causing the formation of a clockwise coastal residual current in the corner above the Meizhou Island. Some feeble clockwise coastal residual currents form around the western corner of Z4, the southern corner of Z2, and the middle portion of Z1. Some counter-clockwise coastal residual currents form in the east of Z4 and the south of Z1.

A reclamation project greatly alters the shape of coastline and narrows the Meizhou Bay, thus leading to a significant change in the residual current. After the project, an obvious counter-clockwise coastal residual current emerges around the rectangular corner. The clockwise coastal residual current in the corner above the Meizhou Island disappears due to the smoother shoreline after reclamation. Some clockwise coastal residual currents in the western corner of Z4 and the southern corner of Z2 are weakened owing to the major mudflat loss. The low-velocity zone shows a major reduction due to the intense flow caused by the narrower coastline.

Effects on the advection-diffusion process

Characteristics of the advection-diffusion process are presented by the concentration tracer method. The main study domain is divided into five zones of interest, and each scenario is studied with an initial tracer concentration of 1 g/L in the selected zone and 0 g/L in the rest of the study domain. Concentrations of tracer at the highest and lowest water levels are used to study the influence of the reclamation project on the advection-diffusion process within the bay.

Figures 17(a)–17(d) present the advection-diffusion process in Z1 before and after the reclamation project. The reclamation project does not obviously change the entrance width of Z1, which largely determines its local hydrodynamic process; thus, the reclamation project exerts little impact on this zone. However, the advection-diffusion process in this zone is relatively weakened with the reduction of water-bearing areas, leading to the decrease in the range that concentration diffuses. Meanwhile, the relative area of the high-concentration region of this zone increases, mainly because of the narrower coastline after the project. Alterations in Z2 are limited, as shown in Figs. 17(e)–17(h). The relative area of the high-concentration region decreases, which means the strength of the flow around the inner portion of Z2 slightly increases. Z3 is located in the tidal channel, where a strong tidal current exists. As shown in Figs. 17(i)–17(l), the reclamation project obviously affect the advection-diffusion process in Z3. After the project, a rectangular corner is formed, which locally inhibits water escape. This feature weakens the advection-diffusion and increases the concentration around this corner. Meanwhile, the distance that the concentration travels along the east-west direction increases, which is highly associated with the enhancement of the horizontal flow. The reclamation project largely decreases the area of Z4, and the weakened advection-diffusion process caused by the decreased water-bearing capacity is clearly presented in Figs 17(m)–17(p). The distance the concentration travels declines, and the relative area of the high-concentration region increases. Compared with other zones, Z5 has a superior advection-diffusion condition owing to its stronger flow and wider exchange boundary. After the reclamation project, the concentration travels farther in the east-west direction on account of the reduced excursion of the horizontal current and its accompanying enhanced current strength. In the north-south direction, however, the travel distance is slightly weakened by the reduced water-bearing capacity (Figs. 17(q)–17(t)).

Effects on sediment transport

By subtracting the average SSC before the reclamation project from that obtained after the reclamation project, distribution of the variation in average SSC during a tidal period is presented in Fig. 18(a). Values above zero reflect the increase in SSC after the project, whilst values below zero reflect the corresponding decrease. Given the weakened north-south current, the suspension process of sediments is reduced, thus leading to the decrease in SSC in most portion of the Meizhou Bay. Meanwhile, the mudflat is a major supplier of sediment, and the decrease in the coverage area of the former may contribute to the decrease in the latter. The SSC decreases largely in the western portion of Z4 by up to –0.2 g/L, which can be most attributed to the mudflat loss. However, SSC in some areas presents a tendency to increase. Z1 shows an obvious increase in SSC of up to 0.2 g/L. The hydrodynamic environment here is relatively weak, which is helpful for sediment accumulation. The western portion of the Meizhou Island reveals a high SSC distribution after the reclamation project resulting from the strengthened horizontal current, which brings more suspended sediments to this area. The rectangular corner demonstrates a complex change pattern in which the SSC in its upper and lower portions increases but that in its middle part decreases. This pattern is highly associated with the alterations in the corner’s hydrodynamic environment. In the middle part, where the current is enhanced locally, more sediment is transported to the upper and lower parts with flow. By contrast, in the upper and lower areas, flow is restricted by the narrow shoreline. Thus, the current is unable to flow freely, and sediment is locked in these areas. The increase in SSC in the upper part of the lateral bay may also result from the intensified horizontal current, which brings more sediment after the reclamation project.

The erosion-deposition process is complex and highly associated with the surrounding hydrodynamic environment, topography, SSC, and other factors. To illustrate the erosion-deposition alteration caused by the reclamation project, Fig. 18(b) presents the seabed deformation within the Meizhou Bay. Here, values above zero indicate further deposition, whereas values below zero indicate further erosion. The corner of Z2 presents erosion that is mainly caused by the locally enhanced flow, which scours the seabed. Moreover, the north-eastern corner of the horizontal channel has a tendency to erode. The enhanced horizontal flow around this corner initiates and enhances the motion of sediment. Together with the good advection-diffusion condition, more sediment is moved out of this corner. The peripheral region of the Meizhou Island shows slight erosion mainly due to the intensified current around there. However, the western portion of this island shows a tendency to deposit, which is highly associated with its outstretched crescent shape; such a shape promotes the collection of sediment and leads to the deposition. The upper portion of the rectangular corner also presents deposition, which is largely caused by the outstretched shoreline shape too. Deposition is also found in the western corner of Z4 and is likely due to the locally reduced tidal current and tidal prism. The scouring process is weakened under these conditions, thus leading to an increase in the seabed height after the reclamation project.

Effects on storm surge and waves

To illustrate how the reclamation project affects storm surges and waves, the maximum storm surges and maximum significant wave heights before and after the project are presented in Table 5. In general, the reclamation project does not exert a significant impact on them. It increases the maximum storm surge by 0.06 m likely because of the reduced tidal prism and the mudflat loss; the amphidromic point may be shifted to a dissipation-enhanced location and consequently alter the entire region, which also contributes to this increase (^{Pelling and Mattias Green., 2013}). It decreases the maximum significant wave height by 0.09 m; this negative effect is consistent with findings in a previous study (^{Shen et al., 2018}) and may be caused by the reduced tidal prism owing to shoreline contraction.

Figure 19 presents the spatial distribution of the maximum storm surges and the maximum significant wave heights before and after the reclamation project in the Meizhou Bay. In general, the reclamation project exerts an inconspicuous effect on the spatial distribution of these features. The maximum storm surges of the inner portions show higher values than the exterior area, which connects to the ocean. This distribution is similar to that found in a related research (^{Xu et al., 2015}) and may be caused by the narrow shoreline restriction and shallow water condition, both of which help accumulate water. Deeper regions tend to have a higher wave height than shallower ones, and areas close to the typhoon track demonstrate a higher wave height than those located farther from typhoon threat. However, areas near the reclamation project region present major changes. After the project, the newly formed rectangular corner shows higher storm surges than before. Storm surges around the inner portions of Z2 and Z3 decrease, and significant wave heights around the north of the Meizhou Island slightly increase.

Conclusions

In the present paper, a 3D numerical model coupled with current, baroclinic process, wave, and sediment is implemented to simulate the hydrodynamic environment and sediment transport within the Meizhou Bay. After analysing the simulated results before and after the reclamation project, the main conclusions are as follows:

1) The reclamation project decreases the tidal prism in the spring and neap tides by 0.65 × 10⁹ m³ and 0.44 × 10⁹ m³, respectively, thus leading to the weakening of the water exchange capacity and the self-purification capacity of the bay. After the reclamation project, the motion strength and coverage of the particle trajectory decrease. Trajectory along the horizontal channel is translated from a curvilinear shape into a linear shape. Current along the tidal channel weakens, whilst current along the horizontal channel strengths. Flow in the eastern corner of Z2 slightly intensifies, whilst flow within Z4 slightly weakens. The residence time presents a major increase. An obvious counter-clockwise coastal residual current emerges around the rectangular corner. The clockwise coastal residual current in the corner above the Meizhou Island disappears.

2) After the reclamation project, the SSC declines in the main area of the Meizhou Bay. The SSC in Z4 shows a major decrease of up to –2 g/L. The western region of the Meizhou Island reveals a high SSC distribution. The rectangular corner has a complex change pattern, where the SSC in its upper and lower portions increases but that in its middle part decreases. The corner of Z2 shows a tendency to erode. However, the western portion of the Meizhou Island, the upper portion of the rectangular corner, and the western corner of Z4 show a deposition tendency.

3) The reclamation project increases the maximum storm surge by 0.06 m and decreases the maximum significant wave height by 0.09 m. After the reclamation project, the rectangular corner shows an increased storm surge, but storm surges around the inner portion of Z2 and Z3 decrease. Significant wave heights around the north of the Meizhou Island slightly increase.

References

[1] Chen J C (1988). A primary study on sedimental resource in Meizhou Bay. J Fujian Normal U (Natural Science), 4(3): 104–111

[2] Chen C, Liu H, Beardsley R C (2003). An unstructured, finite-volume, three-dimensional, primitive equation ocean model: application to coastal ocean and estuaries. J Atmos Ocean Technol, 20(1): 159–186

[3] Chen C, Huang H, Beardsley R C, Liu H, Xu Q, Cowles G (2007). A finite volume numerical approach for coastal ocean circulation studies: comparisons with finite difference models. J Geophys Res, 112(C3): C03018

[4] Chen W B, Liu W C, Hsu M H (2012). Computational investigation of typhoon-induced storm surges along the coast of Taiwan. Nat Hazards, 64(2): 1161–1185

[5] Fujita T (1952). Pressure distribution within typhoon. Geophys Mag, 23(4): 437–451

[6] Guo Y K, Zhang J S, Zhang L X, Shen Y M (2009). Computational investigation of typhoon-induced storm surge in Hangzhou Bay, China. Estuar Coast Shelf Sci, 85(4): 530–536

[7] Guo Y K, Wu X G, Pan C H, Zhang J S (2012). Numerical Simulation of the Tidal Flow and Suspended Sediment Transport in the Qiantang Estuary. J Waterw Port Coast Ocean Eng, 138(3): 192–202

[8] Guo X F, Wang C, Chen C H, Tang J J (2014). Numerical simulation of the transport diffusion of suspended matter during the construction of Fengwei reclamation project at Meizhou Bay. J Appl Oceanography, 33(1): 125–132 (in Chinese)

[9] Hasselmann K, Barnett T P, Bouws E, Carlson H, Cartwright D E, Enke K, Ewing J A, Gienapp H, Hasselmann D E, Kruseman P, Meerburg A, Müller P, Olbers D J, Richter K, Sell W, Walden H (1973). Measurements of Wind-Wave Growth and Swell Decay during the Joint North Sea Wave Project (JONSWAP). Deutches Hydrographisches Institut

[10] Holland G J (1980). An analytic model of the wind and pressure profiles in hurricanes. Mon Weather Rev, 108(8): 1212–1218

[11] Holland G J (2008). A revised hurricane and pressure-wind model. Mon Weather Rev, 136(9): 3432–3445

[12] Hou Q Z, Zuo L Q, Lu Y J, Wang Z L, Mo S P, She X J (2017). Hydrodynamic environment response to human interventions in a Macro Tidal Bay: the example of the Quanzhou Bay. J Basic Sci Eng S, 25(6): 1124–1138 (in Chinese)

[13] Jelesnianski C P (1965). A Numerical calculation of storm surges induced by a tropical storm impinging on a Continental shelf. Mon Weather Rev, 93(6): 343–358

[14] Jia H, Shen Y M, Su M R, Yu C X (2018). Numerical simulation of hydrodynamic and water quality effects of shoreline changes in Bohai Bay. Front Earth Sci, 12(3): 625–639

[15] Liu X H, Huang W R (2009). Modeling sediment resuspension and transport induced by storm wind in Apalachicola Bay, USA. Environ Model Softw, 24(11): 1302–1313

[16] Liu H, Pan W R, Luo Z B, Shi X F (2009). Numerical study on the tidal regime in Meizhou Bay. J Oceanogr Taiwan, 28: 118–122 (in Chinese)

[17] Luo Z B, Pan W R, Zhang G R (2007). The sediment properties and transportation rules of the Xiuyu Sea-route in Meizhou Bay. J Xiamen U Natural S, 46(3): 382–385

[18] Myers V A (1957). Maximum hurricane winds. Bull Am Meteorol Soc, 38(4): 227–228

[19] Mellor G L, Yamada T (1982). Development of a turbulence closure model for geophysical fluid problems. Rev Geophys Space Phys, 20(4): 851–875

[20] Pelling H E, Mattias Green J A (2013). Sea level rise and tidal power plants in the Gulf of Maine. J Geophys Res Oceans, 118(6): 2863–2873

[21] Shen Y M, Jia H, Li C P, Tang J (2018). Numerical simulation of saltwater intrusion and storm surge effects of reclamation in Pearl River Estuary, China. Appl Ocean Res, 79: 101–112

[22] Vickery P J, Skerlj P F, Steckley A C, Twisdale L A (2000). Hurricane wind field model for use in hurricane simulations. J Struct Eng, 126(10): 1203–1221

[23] Willoughby H E, Rahn M E (2004). Parametric representation of the primary hurricane vortex. Part I: Observations and evaluation of the Holland (1980) model. Mon Weather Rev, 132(12): 3033–3048

[24] Xu Y, Chen D W, Li L, Li P L (2015). Numerical modeling of strom surge in Taiwan Strait caused by Typhoon FANAPI. Periodical Ocean U C, 45(2): 008–017 (in Chinese)

[25] Zhu Z C, Hu Z J, Zhang Y Q, Chi W Q, Liu J Q, Bian S H, Xiong C B (2014). Numerical simulation and observations of vertical non-uniform of water temperature in the Meizhou Bay influenced by warm water discharge. Adv Marine Sci, 32(2): 230–238 (in Chinese)

RIGHTS & PERMISSIONS

Higher Education Press and Springer-Verlag GmbH Germany, part of Springer Nature