1. Introduction
The Vietnam Mekong Delta (VMD) comprises 13 provinces and cities in the south of Vietnam (
Figure 1), with a total area of 3.95 million ha (12% of areas of Vietnam) and a population of 19 million [
1]. It provides 50%, 90%, and 60% of Vietnamese food, rice production (second largest exporter of rice worldwide), and seafood, respectively. The delta is also the most concentrated fish biodiversity per unit area and is ranked second only after the Amazon in overall biodiversity [
2].
In recent decades, the Vietnam Mekong Delta Coastal Zone (VMDCZ) has undergone severe coastal erosion due to sediment depletion, inappropriate land use, accelerated land subsidence, sand mining, climate change, and sea-level rise. Dinh C.S. et al. [
3], using remote sensing, show that, within 25 years (from 1990 to 2015), the VMDCZ expanded unevenly. For the first 15 years of this period, the erosion rate was about 85 ha.year
−1, while, in the last decade, it reached approximately 325 ha.year
−1 on average. In certain coastal zones, such as Ca-Mau and Go-Cong, severe erosion has caused coastlines to retreat about a hundred meters in recent years, threatening destruction of sea-dyke systems. Every year, coastal erosion causes approximately 300 ha of agricultural land loss. Thus, studies on sustainable measures for protecting the VMDCZ from erosion are needed.
In the Asian Southeast region, Project Demak in Centre Java, Indonesia would be considered as an excellent example of building with nature in protecting against coastal erosion. The location in Demak concerns a tropical muddy mangrove coast. The key lessons learned from this project are that a combination of a thorough understanding of the biophysical, socio-economic, and governmental system and early stakeholder involvement results in higher vital benefits, reduces costs, and provides the setting for sustainable design solutions. Based on this system understanding, temporary permeable structures made from local material have been chosen to create wave-sheltered areas that stimulate the settlement of sediment and create a habitat favourable to mangrove recolonisation. Once the mangrove forest is fully grown, it will provide protection against waves. It will also provide other ecosystem services, such as food provisioning, tourism, a nursery habitat for fishery production, and CO
2 storage [
4].
In the past, along the coastlines of the Go-Cong Coastal Zone (GCCZ), there has existed a mangrove belt, the width of which varies from several hundred up to 1000 m. The mangrove belt has protected the sea dyke system, 580,000 people, and 63,000 ha of natural lands, including 43,000 ha of agricultural lands, from sea waves and coastal floods. However, according to the observations in situ, in recent decades, the mangrove belt has been destroyed by waves, storm surges, and erosions. It is observed that, every year, the mangrove belt loses about 8–10 m of width. There are locations where the width of the mangrove belt rests several meters only. This positions the sea dyke just in front of the sea and threatens it to be fully destroyed by waves and storm surges.
In the framework of the MARD (Ministry of Agriculture and Rural Development) program, Le Manh Hung et al. [
5] studied the flow regime and distribution of coastal sediments from the Soai-Rap estuary to the Tieu mouth and proposed solutions to prevent erosion of the Go-Cong sea dyke. Physical models have been performed to test the efficiency of mangrove belts in attenuating wave energies. The results showed that at least a mangrove belt of width 200–300 m is needed to efficiently protect the Go-Cong coastlines from erosion. In order to create and maintain such a belt, however, some hard or soft structures must be built even temporarily.
In this context, the European Union (EU), Agence Française de Developpement (AFD), and Southern Institute of Water Resources Research (SIWRR) signed the “partnership research agreement” in 2016 named “Study on the erosion process and the measure for protecting the Vietnam Mekong Delta Coastal Zones from erosion”. The main objective of this agreement is to define the protection measures for the Go-Cong Coastal Zone (
Figure 1). As we all know, mangrove forest belts can be considered as a nature-based solution to protect coastal zones from erosion sustainably. The local authorities of the Go-Cong Coastal Zone would like to restore such mangrove belts. However, other soft and/or hard measures should be implemented first.
In the literature, there is much research on the impacts of emerged or submerged breakwaters on attenuation of wave energy. Koley et al. [
6] have provided a fairly complete summary of this experimental and numerical research. Yet, there are very few references relative to the impacts of protective measures on the evolution of seabed morphology once construction has been completed. Farhad Sakhaeea and Fatemeh Khalili [
7] have used MIKE21 ST to study seabed changes in the pool of Nowshahr port (the southern part of the Caspian Sea) due to construction of breakwaters. They have used the nesting technique to obtain boundary conditions for coupling three subdomains: the smallest one includes an area of 47 hectares, which includes only the pool of Nowshahr port. Areas two and three include larger areas. Luca Carpi et al. [
8] have posed the question “Detached breakwaters, yes or not?” since water stagnation, rip current development in the gaps between multiple structures, and obstruction of the natural panorama are well-recognized issues. They propose a modelling approach to obtain a preliminary evaluation of a theoretical removal intervention and discuss its consequences. The numerical results show that every intervention involving detached breakwater removal is wrong. Indeed, even though a natural panorama is always desirable, a coastal urban settlement without protection from storms is not acceptable. Moreover, the results show that beach nourishment must be considered as a real winning tool for coastal protection. Ali Temel and Mustafa Dogan [
9] have studied wave-induced local scour at the toe of rubble mound breakwaters placed on an erodible seabed both experimentally and numerically. Time-dependent scour depth and wave characteristics in a laboratory flume were observed with ultrasonic methods. To observe and compare scour mechanisms, an equivalent numerical wave channel was used. Numerical simulations using Flow3D provide insight into the hydrodynamics and scour/deposition pattern in front of the breakwater trunk under non-breaking regular wave impacts.
This paper presents the development of a set of numerical models for assessing the impacts of different soft and hard measures to protect a tropical coastal zone from erosions, which undergo complex (semi-diurnal and mixed) tidal regimes with strong monsoon winds. The development comprises adapting two numerical codes, MIKE21-FM and TELEMAC-2D, to the natural and geophysical conditions of the Go-Cong Coastal Zone (GCCZ). A nesting technique has been used to provide the boundary conditions for the different models: regional, local, and study ones (
Figure 2). T-shape breakwaters as hard measures and sandbars as soft ones have been considered. The present paper aims to: (i) assess the performance/effectiveness of T-shape breakwater and sandbars against coastal erosion in the GCCZ, and (ii) compare the capacity of MIKE21-FM and TELEMAC-2D in modelling the erosion/accretion process with and without protection measures. To our knowledge, in the literature, there is no study that compares TELEMAC-2D and MIKE21 yet.
2. Methodology
In the VMD, there are two hydrological seasons: wet and dry, associated with SW (southwest) and NE (northeast) monsoon climates, respectively. In the wet season, the Mekong River discharges into the sea about 80% of its water and 90% of its sediments, and river flows prevail in the estuarine and coastal zones. In the dry season, the tidal influence is dominant. In the Sai Gon–Dong Nai (SG-DN) River, the main-stream flows are significantly regulated by the upstream hydropower reservoirs. Thus, they do not affect the estuaries, and tidal forcing is also dominant.
The Go-Cong Coastal Zone (GCCZ) is situated between the Soai-Rap Estuary (Vam-Co and Dong-Nai Rivers) and the Cua-Tieu one (Mekong Delta). This zone is oriented in the south –north direction and extends over 15 km along the coastline (see
Figure 1).
GCCZ is subject to monsoon winds: from November to March, they blow toward the northeast, and from May to October, they blow toward the southwest. Under the norteast wind, the outflow from the Soai-Rap estuary lightly sediment laden is directed towards the south, combined with wave action seriously erodes the Go-Cong coastline. From May to October, under the southwest wind and the region of freshwater influence (ROFI) of the VMD [
10], the very turbid coastal flow coming from the Cua-Tieu and Cua-Dai estuaries creates accretions offshore rather than near the coast in the GCCZ.
Some well-known protection measures have been investigated, analysed, and modified to apply to the Go-Cong Coastal Zone. The proposed solution is a combination of solid revetment on the shoreline slope, local material fences in the middle, and porous breakwaters offshore. The porous breakwaters attenuate waves and currents and promote accretion. They also protect local material fences, which are installed inside the sheltered area to survive longer and trap fine sediments for mangrove restoration. The dimension of these structures depends on the site-specific conditions of waves, currents, sediments, etc. [
11]. We assess the impacts of two measures: hard structure—the “porous breakwaters” and soft structure—the sandbars. Because MIKE21 cannot simulate the porosity of breakwaters, whereas TELEMAC-2D can do it, henceforth, we consider that the breakwaters are not porous in both two numerical models. The sandbar measures were imported from (Building with Nature Indonesia, Project Demak, Demak City, Central Java, Indonesia), as this project has similar natural conditions to the VMDCZ.
In another study, as an alternative, physical models were used to check wave attenuation and exchanges of suspended sediments through porous breakwaters. As well, the deformation of sandbars (cross-shore processes) caused by waves is studied experimentally in the hydrodynamic laboratory of the Southern Institute of Water Resources Research (SIWRR) [
12]. The wave transmission and reflection coefficients versus the relative crest freeboard values, as well as sediment volumes trapped in the sheltered zones, have been measured and analysed. Limited by the flume width, these experiences could not show how the horizontal arrangement, i.e., the gap dimension between T-shape breakwaters or sandbanks, will play in reducing erosion in the studied coastal zone of Go-Cong. This horizontal arrangement will be one subject of this paper.
2.1. Numerical Models Description
Two numerical models, MIKE21-FM (Danish Hydraulics Institute-DHI, Hørsholm, Denmark, version 2012) and TELEMAC-2D (version 7.2.2), are used simultaneously. TELEMAC-2D is open-source software, while MIKE21-FM is commercial software that prevents user interference. Both models are the solvers of the Saint-Venant equations, which are derived from the Navier-Stokes equations using the hydrostatic hypothesis. In TELEMAC-2D, the Saint-Venant equations are decomposed by a projection method [
13]. This method allows the reconstruction of a Poison’s equation, which is implicitly solved to get water levels. Thus, velocity components can be determined once the water levels are known. In MIKE21 FM, the Saint-Venant equation has been solved by explicit numerical schemes, based on a cell-cantered finite volume technique. TELEMAC-2D is based on a finite element technique. Both use unstructured meshes. Implicit numerical schemes used in TELEMAC-2D have the 1st accuracy order in time and space and allow the use of a high Courant number, i.e., a large time step. Using MIKE21 FM, the users have the option of 1st or higher accuracy order in space. Explicit numerical schemes used in MIKE21-FM restrict the Courant number to less than 1, i.e., the time step should be relatively small. MIKE21-FM is coupled with MIKE21 SW to compute wave propagation, and with MIKE21 MT for mud transport. TELEMAC-2D is coupled with the module TOMAWAC for wave computations and with SISYPHE for sediment transport simulations. It is needed to keep in mind that both TOMAWAC and MIKE21 SW cannot take into consideration the reflection of waves. However, TOMAWAC can take into consideration the breaking waves that would form just upstream of the breakwaters. These breaking waves will promote the resuspension of bottom sediments and erode the seabed at the breakwater toe. SYSIPHE and MIKE21-MT are the solvers of the sediment transport equation.
The parameters used in SYSIPHE and MIKE21-MT are as follows: Mean diameter d50 of suspended sediment is 60 µm; Settling velocity is 0.2 mm.s−1. Sediments are cohesive (muddy). Fine cohesive sediments are transported in suspension (no bed-load) and transport processes strongly depend on the state of flocculation of the suspension and consolidation of the bed. The erosion rate mainly depends on the consolidation of the sediment bed, while the settling velocity depends on the state of flocculation and aggregate properties. Cohesive sediment beds are generally non-uniform. As a result of self-weight consolidation, the bed becomes stratified, with density increasing with distance from the surface. The top layer is generally made of freshly deposited soft mud, while the bottom consolidated layers present higher resistance to erosion. The vertical stratification of the cohesive bed is therefore a key issue, which controls the amount of material to be put in suspension. In SYSIPHE and MIKE21-MT, a cohesive sediment bed can be represented by a fixed number of layers. Each layer is characterized by its concentration and resistance to erosion. The potential depth to be eroded is estimated using the iterative procedure. For each layer j and at each time step, the erosion rate Ej is calculated as a function of the difference between the applied bed shear stress and the critical bed shear strength, using the classical Partheniades formula. If Ej > 0, layer j is eroded. At each time step, the top layer is first eroded. Once it is empty, the next layer is eroded, etc. The erosion/accretion process has been simulated using 3 layers under the bottom.
As the finality is to obtain the numerical simulations for assessing the impacts of the selected hard and soft measures to protect the Go-Cong Coastal Zone from erosions, a nesting technique has been used. Therefore, three models have been developed: Regional, Local, and Study, the one nests the other successively (
Figure 2). The regional model that covers the South China Sea provides the boundary conditions to the local model. The local model occupies the whole Lower Mekong Delta Coastal Zone. This model gives the boundary conditions to model the Go-Cong Coastal Zone.
Table 1 shows the number of elements, maximal and minimal element sizes, and CPU (Central Processing Unit) time (in hours) needed for a one-month simulation using 64 processors, for regional, local, and study domains. First order schemes in space and in time have been used in MIKE21-FM and TELEMAC-2D.
CPU time consumed by a model depends in a large part on the computational mesh. The larger the number of elements, the smaller the size of the computational elements, and the greater the required CPU time.
Table 1 shows that, compared to MIKE21 FM, TELEMAC-2D has a higher speed-up, which allows the use of finer computational meshes. Indeed, for the study model (
Figure 2), for example, MIKE21 FM using a mesh of 48,077 elements, the smallest size of which is 10 m, consumes 60 h (64 Processors) for a one-month simulation, whereas under the same conditions, TELEMAC-2D only needs 18 h using a mesh of 101,126 elements (more two times of elements number) with the smallest size of 8 m.
2.2. Models’ Calibration and Validation
Two field surveys, lasting 15 days each, were performed in 2016′s southwest (SW) from 16 October to 1 November, and 2017’s northeast (NE) monsoons from 24 February to 12 March (see
Figure 3 for the location of 183 mobile samplings and two fixed in situ stations in Go-Cong and U-Minh Coastal Zones). Data obtained from these field surveys are used to calibrate and validate the models.
The surveys include the following river and coastal stations (see
Figure 3):
- ✓
2 fixed river stations in Mekong (at My Thuan) and Bassac (at Can Tho) rivers (coincided with National Hydrology Stations) for measuring discharges (Q) and suspended sediment concentration (SSC) by ADCP (Acoustic Doppler Current Profiler).
- ✓
2 fixed coastal stations at Go-Cong and U-Minh to measure waves (height, period, and direction), and SSC;
- ✓
Marine sampling stations from 2 mobile ships cruising along the east and west seas are:
- +
Water level (hourly at the 2 fixed stations);
- +
Vertical distribution of velocity;
- +
Vertical distribution of salinity (5 points for each vertical line);
- +
Vertical distribution of sediments (5 points for each vertical line);
- +
Waves (height, period, and direction).
In the Regional Model, tidal elevations are imposed at open boundaries, such as Taiwan, Luzon (Philippines), and Malacca Straits (Malaysia–Indonesia) (
Figure 2). The 71 evenly distributed point data of the FES-2014 global tidal model and height, direction, and period of Wave Watch III data (6 points) are used for calibrating the Regional Model. Atmospheric pressure, wind speeds, and direction data of NOAA (National Oceanic and Atmospheric Administration) Global Forecast System are used as surface forcing. These forcing are also imposed on the Local and Study models.
The Local Model is set up using a computational domain with My Thuan on the Tien River, and Can Tho on the Hau Rivers, as upstream open boundaries at about 100 km from the sea. It covers over 774 km of the VMD coastline, from Soai-Rap Estuary (on the east coast) to Ha-Tien (on the west coast). The sea open boundaries are about 130 km far from the coastlines (
Figure 2). We extract the tide and wave boundaries from the Regional Model. The calibration has been performed for flows, tide levels, waves, and SSC (suspended sediment concentration) using data from two fixed in situ stations, Go-Cong and U-Minh, and 11 tide gauging stations (set up by the Vietnamese General Directorate of Hydro-Meteorology—VGDHM), 183 in situ (mobile) samplings performed in 2017’s NE and 2016’s SW monsoons. At Can-Tho and My-Thuan as the upstream boundaries, discharge and SSC values are imposed. At Vam Co, Soai-Rap, Dong Tranh—Rivers connecting to the Mekong Delta, discharges and SSC measured from VGDHM gauging stations and simulated from a 1D model are used as the other upstream boundaries of the Local Model.
The Study Model covers the Go-Cong Coastal Zone connecting to the Soai-Rap Estuary, into which the Saigon-Dong Nai and Vam-Co Rivers discharge (
Figure 2). The conditions on discharges and SSC at the upstream boundaries, as well on tides, waves, and SSC at the downstream ones, are extracted from the Local Model.
Figure 4,
Figure 5 and
Figure 6 show the validation of both MIKE21 and TELEMAC-2D on tide elevation, significant wave height, and suspended sediment concentration (SSC). We note that the numerical tidal elevations provided by MIKE21 and TELEMAC-2D cannot be distinguished. They coincided with the measured values (
Figure 4) at different stations. Concerning significant wave heights,
Figure 5 shows that the significant wave heights obtained from both models are in agreement with the measured values. TELEMAC2D-TOMAWAC gave values that more approach to the measured values than MIKE21-SW, in particular for the first 8 days.
Figure 6 shows the depth-averaged values of SSC obtained from TELEMAC-2D and MIKE21 compared with SSC values measured on the surface and bottom. We remark that TELEMAC-2D’s results are in good agreement with the measured ones, correctly varying between SSC values measured on the surface and bottom. MIKE21 gives the values, which are still acceptable, although they are only slightly greater than the ones measured on the bottom.
In calibration analysis, the values of Nash–Sutcliffe efficiency (NSE) are computed and shown in
Table 2.
The
NSE values are provided by [
14] as follows:
NSE indicates how well the plot of observed versus simulated values fits.
NSE ranges between −∞ and 1.0 (1 inclusive), with
NSE = 1 being the optimal value. Values between 0.0 and 1.0 are generally viewed as acceptable levels of performance, whereas values < 0.0 indicate that the mean observed value is a better predictor than the simulated value, which indicates unacceptable performance [
14].
Overall, all models are calibrated for discharges and tide levels (
NSE = 0.65 ÷ 0.98) and significant wave heights (0.61–0.73). As the calibration for SSC is often very difficult, there is no accurate correlation between
NSE and SSC.
Figure 6 shows that the numerical model correctly simulates the variability of SSC even if the values measured in the field and those simulated by the numerical model are different and their correlation changes over time. Settling and flocculation of cohesive sediments are very complex processes that change the size, shape, and colour of particles during a long time of sedimentation and deposition. Many factors and sediment properties strongly affected them.
No comparison between computed and observed morphological change is possible because there is no such observation in situ. We suggest monitoring should be implemented once protection measures have been completed.
3. Numerical Model Results and Discussion
Two protection structures have been chosen: T-shape breakwaters and sandbars with different configurations. No protective measure is considered as a baseline.
The 2014 data have been used to simulate the impacts of the chosen structures because 2014 has been considered as a mean year on waves and river discharges.
The results of typical two-month simulations in the NE and SW monsoons provided by both models, MIKE21 FM and TELEMAC, are discussed herein.
3.1. No Protection Measures (Baseline)
Both TELEMAC-2D and MIKE21 FM have been used to simulate the baseline scenario in the study domain, which is affected by currents coming from the Soai-Rap Estuary (Dong-Nai and Vam-Co Rivers) in the northeast and the Cua-Tieu Estuary (Mekong Delta) in the southwest, and particularly by waves and tides coming from the East Sea.
Numerical results of NE (January 2014) and SW monsoon seasons (September 2014) at four sampling points P1 to P4 (see
Figure 7 for their location) are extracted and analysed.
According to the numerical results, currents in the Soai-Rap and Cua-Tieu Estuaries and their adjacent coastal zone are strongly dominated by tides. In the Go-Cong Coastal Zone, they can reach up to 1.35 m.s−1 oriented toward the north and northeast in the rising tide and 1.65 m.s−1 toward the south and southeast in the falling tide. The magnitude of currents in the SW monsoon is higher than that in the NE.
Figure 8 and
Figure 9 show wave rises at P1, P2, P3, and P4 in the NE and SW monsoons of 2014. Note that the significant wave heights in the NE monsoon are much higher than in the SW monsoon. In the SW monsoon, the probability of significant wave heights < 0.1 m is over 90%, while, in the NE monsoon, it is only about 10% to 25%. The probability of significant wave heights of 0.2 m is less than 2% in the SW monsoon, while it can reach up to 40% in the NE monsoon. The differences between NE and SW conditions cannot be attributed to river flow because two main forcings in this study are tides and monsoon winds. River flow influence is weak compared with these two forcings.
3.2. Impact of Breakwaters
As above mentioned, The GCCZ has suffered serious erosions in NE monsoons. Longshore currents driven by waves and tides cause erosions. The beach profile of the GCCZ can be divided into three segments with different slopes (
Figure 10c): 3% over the first 100 m from the coastline, nearly 0.00% on the 2nd segment of the next 100 m, and 1.10% for the last segment extending from 200 m to 500 m.
Breakwaters have been proposed as a protection measure. In this sub-section, we will assess the capacity of breakwaters to attenuate wave energy and trap sediments to restore mangrove belts. TELEMAC and MIKE21 systems will model wave propagation, tide- and wave-driven longshore currents, sediment transport, and morphological evolution, with and without protection measures.
Three scenarios of T-shape breakwaters are investigated to predict their impacts on currents, waves, and morphology in the study area (
Figure 10), which extends until 300 m (Y = 300 m) offshore. T-shape breakwaters attenuate wave energy, facilitate trapping fine sediments and thus promote mangrove rehabilitation. The dimensions of T-shape breakwaters (
Ls,
Lg; see
Figure 10) are defined by the national guideline named “Hydraulic structures—Requirements for sea dike design—TCVN 9901: 2014) and based on the physical models [
12]: the length
Ls = 600 m and the gap between the two successive breakwaters
Lg = 30 m, 50 m, and 70 m (denoted as KB1, KB2, and KB3 scenarios and KB0 as no-breakwater or baseline, respectively). Note that the results from physical models have suggested a crest level of +2.2 m above MSL. Therefore, the “groin” parts (cross-shore parts) have a crest level of +0.50 m to prevent just a part of longshore currents. These groin parts can be installed with local material (such as bamboo or melaleuca fences) to reduce the costs and can be adjusted if needed. Additional scenarios (KBi-75%) are proposed in assuming a sediment-flux reduction by 75% due to sediment trapping by upstream dams. It is necessary to keep in mind that, because of grid resolution, in which the smallest size varies within 7–15 m, the numerical results for KB1 (
Lg = 30 m) should be taken with caution.
3.2.1. Tidal Exchanges in Pumping Mode between the Sheltered Zones Protected by Breakwaters and the Sea
In the area studied, currents were generated by two forcings: tides and waves. Fisher et al. [
15] pointed out that the tidal exchange between a bay and the ocean is often described by an ebb–flood asymmetry, in which the ebb flow exits from the bay as a jet, while the return flow during the flood is aspirated approximately as a radial sink, a process often referred to as “tidal pumping”.
The tidal pumping should be effective when the offshore length scale of the jet,
Ljet, significantly exceeds the radial scale,
Lsink, of the withdrawal zone. According to Chadwick and Largier [
16], we can calculate the ratio of two length scales by:
Further, Vp = 2ηA is the tidal prism in the bay; η is the tidal amplitude, A is the surface area of the bay, Am is the area of the mouth cross-section, and b0 is the bay mouth width, w is the fraction of a complete circular sink occupied by the withdrawal zone.
The sheltered zone protected by two neighbouring T-shape breakwaters can be considered as a small bay. For KB1, b0 = 30 m, water depth of 3 m, Am = 90 m2, A= 171,000 m2, η = 1.50 m, Vp = 513,000 m3, w is about 10,000 m2, then Lx = 5700 m. The ratio Ljet/Lsink = 1727 is large. The “tidal pumping” could be effective in our study case. However, it should be kept in mind that wave-driven longshore currents can largely modify this mechanism.
Tidal exchanges by pumping modes combined with wave actions at the gap have been numerically modelled for September 2014 (SW monsoon).
Figure 11 shows the mean unit discharges at ebb and flood periods of KB1 (gap size 30 m) and KB3 (gap size 70 m), computed by TELEMAC-2D. The unit discharges have been averaged over half of the tidal cycle corresponding to ebb and flood periods, respectively. We note that even “tidal pumping” is effective, but, under longshore currents, the mechanism of “tidal pumping” has been modified. Although the flood flow takes the form of a jet inside the sheltered zone, the withdrawal zone is not radial, and the ebb flow is not ejected as a jet. Because of a small gap, the velocity values in KB1 are strong at the gap, with a maximum velocity of 0.53 m·s
−1, whereas, in KB3, the velocity at the gap is smaller and can reach a maximum value of 0.29 m·s
−1. The water depths at the gap vary from 1.20 to 4.30 m. The flood discharges for KB1 and KB3 are 15.74 m
3·s
−1 and 16.15 m
3·s
−1, respectively. Thus, KB3 allows a larger flood discharge, i.e., bringing more suspended sediment than KB1. The flood currents produce two large eddies that are asymmetrical inside the sheltered zones. The eddies in KB1 are stronger than those in KB3. The strong eddies in KB1 bring sediments following a circular trajectory, which is anti-cyclonic in the north part and cyclonic in the south part, and return them to the east side. The weak eddies of KB3 transport a great deal of suspended sediments to the west side of the sheltered zone, and few sediments arrive on the east side. This supports the remarks made later in the next section: the larger the gaps, the more accretion for zone 1 (west side). The results are the opposite in zone 2 (east side); i.e., the smaller gap, the higher the accretion there.
3.2.2. Assessment of Different Scenarios of the Breakwaters
During our study, we note that the breakwater impact on currents and waves provided by MIKE21 FM and TELEMAC-2D is almost the same. Therefore, we only present MIKE21 FM’s one.
In the sheltered area (such as at P5; see
Figure 7 and
Figure 10a for its position), one notes a wave height reduction both in intensity and duration.
Figure 12 shows how the gap size could affect wave energy reduction. We can remark that, the smaller the gap (
Lg), the greater the reduction in the wave energy. Clearly, scenario KB1 reduces the most energy compared to KB2 and KB3.
Figure 13 and
Figure 14 illustrate the wave reduction at P5 in the NE and SW monsoons for different scenarios. The probabilities of calm waves (Hs < 0.1 m) are increased from 17.42% (KB0) to 78.17, 70.36, and 61.42% for KB1, KB2, and KB3, respectively, in the NE monsoons. The results are similar in the SW monsoons, with the calm probabilities increasing from about 91% to 99% for all scenarios.
Considering KB3, in the shelter area, we note that Hs > 0.6 m is about 50 m from the breakwater toward the coastlines. To reduce Hs < 0.4 m (for mangrove rehabilitation purposes), Lg should be lower than 70 m if no additional measure (such as bamboo fences) would combine with breakwaters.
The morphological changes in KB0 ÷ KB3 scenarios provided by MIKE21 after 1 month in NE monsoon (January 2014) are presented in
Figure 15.
In order to assess the capacity of breakwaters in wave reduction and in sediment trapping, we divide the sheltered zone (behind the breakwater line) into two zones, ones 1 and 2. Each zone is 13 km long. Zone 1 begins from the shore and has a width of 150 m. Zone 2 expands 150 m from zone 1′s border seaward and ends at the breakwater line, 300 m from the shoreline. Zone 3 begins from the breakwater line and ends at 1000 m offshore. Therefore, the three zones have an area of 2.087, 2.236, and 9.395 m
2, respectively (see
Figure 11b). Looking at the evolution of the thickness of erosion/accretion layers in different zones and the reduction in wave rises, we can deduce how breakwaters and sandbars react to incident wave energy.
Table 3 summarises the numerical results. In January during the NE monsoon, all the zones would be eroded if no protective measure is used (KB0, baseline) or in zone 3, beyond the breakwater line, outside the protected area. The numerical results approve the protective role of breakwaters for the studied zones. No erosion has been observed; inversely, accretions have occurred in zone 1 and zone 2. In the shelter area, the smaller
Lg, the bigger the accretion volumes; in zone 1, the net accretion volumes of KB1, KB2, and KB3 are 0.046, 0.036, and 0.031 Mm
3 (equivalent to an average thickness of 2.2, 1.7, and 1.5 cm), respectively; in zone 2, the net accretion values of KB1, KB2, and KB3 are 0.068, 0.060, and 0.060 Mm
3 (equivalent to the average thickness of 3.0, 2.7, and 2.7 cm), respectively. However, in zone 3, the smaller the
Lg values, the larger the eroded volumes. Indeed, KB1 creates the largest erosion volume.
In September during the SW monsoon, in all the zones, as mentioned in
Section 2, there are only accretions. The net accretion volumes increase sea-ward, ranging from zone 1 to zone 3 for all the scenarios. In the sheltered area protected by breakwaters, the gap between the two successive breakwaters (
Lg) can affect accretion: in zone 1, the net accretion volumes of KB1, KB2, and KB3 are 0.038, 0.031, and 0.025 Mm
3 (equivalent to an average thickness of 1.6, 1.5, and 1.2 cm), respectively; in zone 2, the net accretion of KB1, KB2, and KB3 is 0.091, 0.087, and 0.087 Mm
3 (equivalent to the average thickness of 4.1, 3.9, and 3.9 cm), respectively.
Lg does not have any effect on the accretion volume in zone 3 beyond the breakwater line: the net accretion volume is 0.123, 0.091, 0.093, and 0.096 Mm
3, respectively, for KB0, KB1, KB2, and KB3 (equivalent to the average accretion thickness of 1.30, 1.00, 1.00, and 1.00 cm).
3.2.3. Impacts the Breakwater on Morphology with 75% Sediment-Flux Reduction
Construction of 22 upstream reservoirs on the Mekong River will trap about 75% of current sediment fluxes [
17]. The question that will be naturally posed is how the morphology will be changed in the Go-Cong Coastal Zone. To answer this question, the simulations at the sediment-flux reduction by 75% are performed. The morphology response for the three scenarios was presented in
Table 3. KB1 (
Lg = 30 m) reacts better to 75% sediment reduction with the highest accretion in zones 1 (2.2 cm) and zone 2 (2.5 cm) compared to the other scenarios.
Figure 16 presents the morphology changes provided by TELEMAC-2D in the study area for one month of the NE monsoon (January 2014). The accretion/erosion volumes and average thickness were compared for different scenarios in
Table 4. Without protective measures (KB0), all the zones would be eroded. With breakwaters, we can see the accretion in the sheltered zones (behind the breakwaters). Zone 1 receives higher accretion thickness than that in zones 2 and 3, beyond the breakwater line, which is eroded in all the scenarios. This means that breakwaters can trap sediment in the sheltered zone but do not prevent erosion offshore. The larger the gaps, the more the accretion for zone 1 and the less the erosion for zone 3. However, the results are the opposite in zone 2 (the smaller the gap, the higher the accretion). KB3 is the most appropriate scenario because of its highest accretion and smallest erosion.
3.2.4. The Selected Breakwater Scenario
To select the best scenario, we suggest two criteria. They are (1) the highest accretion volume, and (2) the lowest erosion thickness.
According to the MIKE21’s results (
Table 3), KB1 (gap of 30 m) is better than other scenarios (gaps of 50 m and 70 m) in a net accretion in the sheltered area, but it creates the largest eroded volume in zone 3. TELEMAC’s results (
Table 4) provide a slight advantage to KB3, with the highest accretion and smallest erosion for all zones.
Based on the criteria of highest accretion volume, lowest erosion thickness, and relatively better reaction to sediment reduction of 75% due to upstream reservoirs, KB3 (gap of 70 m) has been chosen.
3.3. Impacts of Sandbars
Wave reduction by sandbars is mostly due to bed frictions and depth-induced breaking along the width of sandbars. The configuration of the sandbar should be properly determined. Based on the experiences from a similar existing project (Demak, Central Java, Indonesia) and on the physical models [
12], two sandbars of the width of 50 m (SB50) and 70 m (SB70), with a crest elevation of −2.4 m (MSL), a length (
Ls) of 1000 m, and a gap (
Lg) of 200 m, installed at 500 m offshore (Y = 500 m), are investigated (
Figure 17).
Under wave energy, while playing the role of wave attenuation, the sandbars will be deformed from the initial shape. This deformation should be taken into consideration in simulation of wave propagation and sediment transport over sandbars. In the EU-ADF Project LMDCZ (Lower Mekong Delta Coastal Zone, [
18]), the deformation of sandbars has been investigated by using a physical model [
12]. Thus, the simulation of sandbars will be performed in two steps. In the first step, we validate the numerical models in simulating the deformation of the sandbar from the initial shape, ignoring morphological changes in the study area. In the second step, we use the already validated models to simulate the wave propagation and sediment transport over the sandbars and their impact on morphological changes in the study area.
As the Go-Cong Coastal Zone extends over 15 km in a north to south direction, we study the deformation of the sandbar of 15 km by decomposing it into three 5 km segments, which are South, Central, and North, respectively.
Table 5 shows that, after January and September 2014, the sandbar lost 11.64%. As the hydrodynamics of the Soai-Rap Estuary are more active than those of the Cua-Tieu Estuary, the lost sand volumes and the erosion rate decrease southward. The lost volume in January (9.85%) is much higher than that in September (1.8%). Stronger winds/waves of the NE monsoon can explain it. A one-year simulation is needed to estimate the lost volume correctly. We roughly evaluate it at about 20% per year. Thus, this sediment lost volume can contribute to the accretion behind the sandbars, but it cannot be considered as a permanent solution.
Because the morphological impacts of two scenarios of crest elevation of −2.4 m are insignificant, an additional scenario with a crest elevation of 0.0 m (MSL) and width of 70 m is considered (denoted as SB70(0.0)).
Figure 18 illustrates the morphological changes by using breakwaters KB3 in comparison with sandbars SB70(0.0), predicted by MIKE21 FM. Moreover,
Table 6 shows that, in January during NE winds, the sandbars SB70(0.0) have fewer impacts on erosion/accretion than breakwaters KB3 in all zones. In September during SW winds, we have the same remark for Zone 1. However, in zones 2 and 3, sandbars make two and three times more accretion than breakwaters. Sandbars react to sediment reduction (−75%) better than breakwaters. The net accretion volumes are 0.075, and 0.049 Mm
3 in zone 2 and 3, for breakwaters KB3 (
Lg = 70)-75%, while they are 0.126, and 0.083 Mm
3 for sandbars SB70 (0.0)-75%.
We illustrated the impact of sandbars provided by TELEMAC in January 2014 (NE winds) in
Figure 19.
Table 7 compares the net accretion volume/thickness engendered by sandbars SB70 (0.0) and breakwater KB3 (
Lg = 70). Without protective measures (KB0), all the zones would be eroded. This result agreed with MIKE21’s one. With the presence of KB3, there are important accretions in zones 1 and 2 (accretion thickness is 15.3 and 10.0 cm, respectively), but it provokes serious erosion in zone 3 (eroded thickness is 5.3 cm). Contrarily, sandbars prevent wave reflection and trap higher sediment volumes in zone 3 (2.9 cm of accretion thickness). However, its role in reducing erosion in the sheltered zone is weak. Zone 2 is still eroded by −3.9 cm with SB70 (0.0).
3.4. Comparison between MIKE2-FM and TELEMAC in Simulating the Impacts
To compare the morphological impacts predicted by MIKE21-FM and TELEMAC, the accretion/erosion volumes in Zone 1+2 (sheltered area) and zone 3 are extracted for January and September 2014.
Table 8 shows that the erosion/accretion trends of the two models are very similar. Without protecting measures (KB0—baseline), both models confirm that the study area should be strongly eroded. According to the numerical results from these two models, the proposed protection measures, i.e., breakwaters and sandbars, cause the accretion in zone 1+2, with the breakwaters being more efficient in generating more accretion volumes than the sandbars. In zone 3, breakwaters KB3 (
Lg = 70) promote erosion by wave breaking. The eroded volumes computed by TELEMAC-2D are even twenty times bigger than that provided by MIKE21 (−0.495 vs. 0.029 Mm
3). In the same zone, sandbars prevent this phenomenon and trap sediments to generate accretion.
For breakwaters, the TELEMAC-2D yields higher accretion in zone 1+2 and higher erosion in zone 3 than MIKE21 FM. For sandbars, MIKE21-FM produces higher accretion than TELEMAC-2D in zone 2, but the result in zone 3 is the opposite.
4. Conclusions and Recommendations
Very good calibration is obtained for tidal levels and currents. It is good enough for wave heights and sediments. This concerned both MIKE21 FM and TELEMAC-2D models. The water and SSC exchanges through the gap between two successive T-shape breakwaters or across the sandbars were taken into consideration by the models. Thus, the effect of tidal exchanges by pumping modes at the gap has been numerically examined for September 2014 as an illustration.
The numerical results of the protection measure impact (breakwaters and sandbars) in the Go-Cong Coastal Zone have approved wave attenuation and sediment accretion trends inside the sheltered area for two typical months of NE and SW monsoons.
For protecting the Go-Cong Coastal Zone, two solutions should be considered: (i) breakwaters installed at 300 m offshore, with a unit length of 600 m, the gap between two consecutive units is 70 m, and crest elevation is 2.20 m above MSL; (ii) sandbars installed at 500 m offshore, with a unit length of 1000 m, the gap between the two consecutive units is 200 m, and crest elevation is 0.00 m above MSL.
Sandbars considerably reduce erosion in the sheltered area in comparison with the baseline. The positive impacts of sandbars have been observed in zone 3 seaward with a higher accretion thickness. T-shape breakwaters firmly protect the sheltered areas, i.e., Zones 1 and 2. As supported by the pumping mode concept and numerically shown, the larger the gaps, the more accretion for zone 1. However, the results are the opposite in zone 2; i.e., the smaller the gap, the higher the accretion. Erosion can occur in zone 3 because of the wave breaking and reflection effect. Therefore, if this solution is chosen, a rock and gravel backfilter should be set up at the breakwater toes to avoid scouring there.
Accompanying these protection measures, T-shaped fences made from local materials should be installed behind these protective structures. These fences enable better trapping of fine sediments for mangrove rehabilitation.
Down drift erosion due to breakwaters was not found in the results of MIKE21-FM and TELEMAC-2D models. It could be explained by the fact that the end part of the breakwater system is directly connected to the Cua-Tieu River.
Calibration and comparison of morphological changes require many in situ measurements, which are very expensive. Only one project cannot cover such calibration and comparison. In this study, we would like to provide an initial idea regarding the capacity of breakwaters and sandbars. We will need future studies to be more conclusive. The appropriate solutions should be defined by long-term simulations, which should last at least the full NE and SW monsoon seasons consecutively, especially for estimating the lost volume of sandbars.
TELEMAC-2D and MIKE21-FM provide the same trends for the erosion/accretion process, with a slight advantage in the calibration accuracy of TELEMAC-2D compared with MIKE21-FM. However, in this study, TELEMAC-2D needs a CPU time of three times less than that of MIKE21-FM.
The commercial DHI MIKE toolkit is popular because of interface convenience. However, it is a black box that does not allow the users to intervene when necessary. As its computational speed-up is low, it consumes exorbitant CPU time. Conversely, the worldwide well-known open-source EDF TELEMAC toolkit has no interface and is not convenient to use. However, it has very high computational speed-up with perfect computing parallelisation, consuming little CPU time. It allows using very fine meshes to obtain the results exhibiting physical phenomena correctly. Moreover, users can change the code sources based on their needs. As the data format is open, with TELEMAC, users can use a variety of public software programs or create their tools (using FORTRAN, C++, Python, Perl, etc.) to process data and visualise results. That is why TELEMAC would be preferable for an adaptation for waves, currents, and sediment transport modelling, which require long-term simulations and some development.