Morphodynamic Evolution of a Nourished Beach with Artiﬁcial Sandbars: Field Observations and Numerical Modeling

: Beach nourishment, a common practice to replenish an eroded beach face with ﬁlling sand, has become increasingly popular as an environmentally friendly soft engineering measure to tackle coastal erosion. In this study, three 200 m long offshore submerged sandbars were placed about 200 m from the shore in August 2017 for both coastal protection and beach nourishment at Shanhai Pass, Bohai Sea, northeastern China. A series of 21 beach proﬁles were collected from August 2017 to July 2018 to monitor the morphological changes of the nourished beach. Field observations of wave and tide levels were conducted for one year and tidal current for 25 h, respectively. To investigate the spatial-temporal responses of hydrodynamics, sediment transport, and morphology to the presence of three artiﬁcial submerged sandbars, a two-dimensional depth-averaged (2DH) multi-fraction sediment transport and morphological model were coupled with wave and current model and implemented over a spatially varying nested grid. The model results compare well with the ﬁeld observations of hydrodynamics and morphological changes. The tidal range was around 1.0 m and the waves predominately came from the south-south-east (SSE) direction in the study area. The observed and predicted beach proﬁles indicate that the sandbars moved onshore and the morphology experienced drastic changes immediately after the introduction of sandbars and reached an equilibrium state in about one year. The morphological change was mainly driven by waves. Under the inﬂuences of the prevailing waves and the longshore drift toward the northeast, the coastline on the leeside of the sandbars advanced seaward by 35 m maximally while the rest adjacent coastline retreated severely by 44 m maximally within August 2017–July 2018. The model results demonstrate that the three sandbars have little effect on the tidal current but attenuate the incoming wave signiﬁcantly. As a result, the medium-coarse sand of sandbars is transported onshore and the background silt is mainly transported offshore and partly in the longshore direction toward the northeast. The 2- and 5-year model simulation results further indicate that shoreline salient may form behind the sandbars and protrude offshore enough to reach the sandbars, similar to the tombolo behind the breakwater.


Introduction
Beach nourishment, also known as beach replenishment, beach feeding, or beach recharge, is a beach restoration process to fill the eroded beach with the sand from somewhere else [1,2]. Compared with hard engineering coastal defense such as seawall and breakwater, a soft engineering alternative such as beach nourishment takes less time to construct, is more compatible with the environment, and provides ecosystem services. However, under the action of wave, wind, and tidal forcing, the filling sediment will be swept away so the beach nourishment needs to be repeated regularly [3]. The hybrid scheme with both hard and soft engineering coastal defenses takes the advantages of both approaches, and therefore has become popular recently [4].
Recent studies show that adding the filling sand in the form of artificial sandbar in the surf zone may serve the dual purpose of attenuating incident waves and restoring beach [5,6]. This new technique is also known as submerged berm, submerged mound, or shoreface nourishment as it is under the water surface most of the time. Unlike hard engineering defense, an artificial sandbar will be eroded and thereby feeds the beach with its constituent sand at the cost of self-sacrifice. The resulting migrations and morphological changes of sandbars would modify the nearby hydrodynamics which in turn would influence the sediment transport, coastline and beach morphology [7]. A number of recent studies have been dedicated to the beach morphological changes after the artificial sandbar construction [8][9][10][11][12].
Unlike artificial sandbars, the generation of natural sandbars has been attributed to forced response [13][14][15] or self-organization [16][17][18][19]. According to the former school of thought, the generation of natural sandbars is the passive response to forcing conditions, and therefore the morphological variations can be predicted solely by forces [14,15]. Wave breaking is one of the most important driving forces for sandbar formation in terms of sandbar location, size, and local depth [20]. The latter school of thought attributes the morphological changes to the feedback between hydrodynamics, sediment transport, and evolving morphology instead of external forcing [15][16][17][18][19]. The cross-shore migration of natural or artificial sandbar is of significance for understanding the mechanism of nearshore morphological evolution. Sandbars move shoreward under low-energy wave condition between storms [21], and seaward under high wave energy during storms [22]. Onshore sandbar migration under low-energy waves may be attributed to skewed fluid acceleration [23,24], onshore sheet flow in the boundary layer [21], Stokes drift [25], etc. Waves at the peak of a storm may produce a strong offshore undertow in the surf zone, which carries the sediment offshore and cause the sandbars to migrate seaward [26,27].
In the past decades, it has been a common practice to rely on empirical or semiempirical formulae to predict the beach nourishment profile, i.e., the eventual equilibrium profile after the nourishment project construction. The equilibrium beach profile is related to the local sand characteristics and nearshore hydrodynamics which is site specific, and therefore the empirical or semi-empirical profile equation is often an appropriate technique. Bruun [28] derived a two-thirds power law formula to express the field measurements of natural beach profiles at Monterey Bay, California, and the coast of Denmark. It was then improved by Dean [29] using further 502 observed beach profiles in the United States and was used in the demonstration of some important coastal processes [30][31][32]. However, it is also necessary to adopt both physical and numerical methods for complex profiles such as the presence of sandbars and troughs. Physical model experiments are usually performed to optimize the design scheme of a beach nourishment project. However, physical model experiments are expensive and have scaling effects. With the increasing computer power and improved modeling capability and accuracy, numerical modeling has become a popular tool to predict the morphological evolution after the construction of a coastal structure or beach nourishment. Morphodynamic models can be either process-oriented or behaviororiented [33][34][35]. In recent years, state-of-the-art process-based morphodynamic models such as Delft3D and Mike21 have gained popularity and widespread applications.
A lot of previous morphodynamic model studies of beach nourishment assumed a single-fraction sediment with uniform diameter [36][37][38][39][40]. During beach nourishments, however, it is a common practice to apply filling sand slightly larger than the original beach sand in order to slow down the erosion the introduced sand from the beach system [41]. Incorporating multi-fraction sediment transport model is expected to better capture the transport processes of sediment with mixed sand size, leading to more accurate predictions of geomorphological evolution of the nourished beach system [42][43][44]. Another advantage of multi-fraction sediment lies in the prediction of the sediment sorting process which is of significance for the design of beach nourishment sand and the assessment of the seabed appropriateness for the aquatic species. Broekema el al. [45,46] employed a two-dimensional vertically morphodynamical model with 8 fractions and considered the effects of short wave grouping and infragravity waves. Their results fit well with the observations of beach profile and bed composition development, reproducing the cross-shore sediment sorting properly. Huisman et al. [47] applied both 2DH and 3D based models with 5 sediment fractions to hindcast the spatiotemporal pattern of sediment size in a mega-nourishment area [48], which demonstrated that different kinetic characteristics of fine and coarse sand cause considerable D 50 spatial changes.
The study site of the Shanhai pass beach nourishment project is located in Qinhuangdao, a city in the middle of the west coast of Bohai Sea, in northeastern China. Nearly 45% of the beaches in Qinhuangdao, approximately 80 km long, are not up to the standard of tourism due to shortened beach width, steepened slope, and coarser sand arising from chronic erosion. To tackle the erosion, the local government implemented a series of beach nourishment projects [37,[49][50][51][52] together with submerged breakwaters, artificial headlands, submerged sandbars. However, too many coastal protection structures that calm the hydrodynamic environment too much may cause water quality degeneration. It is therefore important to optimize the combination of hard and soft engineering structures. At the Shanhai Pass site, shown in Figure 1, a shipyard and the coastline create a semi-closed region with low-energy waves and tidal currents, and therefore it is a suitable site to apply submerged sandbars to attenuate waves and feed the sheltered beach nearby. The objective of this study is to reveal both the short-term (a tidal cycle) and long-term (1-5 years) morphodynamic behavior of the nourished beach with submerged sandbars and examine the underlying physical mechanisms. Although the morphodynamic model with multi-fraction sediment has been applied successfully, it is still worthy of using it to study the sediment transport and morphodynamics together with their underlying driving forces after the beach nourishment in these unique semi-closed waters surrounded by the coastline and the Shanhai Shipyard. The model is driven by the tides and realtime waves in one year, during which no storms occurred, to reproduce generally annual morphodynamic changes. When analyzing the model results, the typical normal and strong waves during the real-time wave series are chosen, which are much weaker than a storm wave. The paper is organized in the following order: First, the field observation of tide, wave and cross-shore beach profiles from elevation −6 m to 2 m are described. Next, the wave-current coupling model and multi-fraction sediment transport and morphological model are set up to simulate the sea-bed profile evolution for one year after the installation of three submerged sandbars. Finally, both 2-and 5-year simulations are carried out to explore morphological characteristics in a longer term.

Study Area
The beach nourishment took place in the northeast of Qinhuangdao, China, which is close to the renowned tourist attraction, Shanhai Pass, the eastern end of the Great Wall. Qinhuangdao is a coastal city, 280 km east of Beijing, located at the coast of the northwestern Bohai Sea [37]. Qinhuangdao is one of the most famous sea resorts known for its plentiful high-quality sandy beaches including Shanhai pass. The local climate is characterized as semi-humid continental climate influenced by the monsoon circulation in the east coast of China and belongs to the eastern warm temperate monsoon zone. Under the great influence of the ocean, the local climate is also characterized by dry spring, mild summer, slightly cool but sunny autumn, and a long mild winter.
Based on the ratio between its diurnal and semidiurnal harmonics [53], tide is usually divided into four types: regular diurnal, irregular diurnal, regular semidiurnal, and irregular semidiurnal tide. The tide is regular diurnal with a tide ratio of 4.73 [54] in this coastal region. However, tidal current is characterized as semidiurnal current in Qinhuangdao. An amphidromic point of M 2 is located near Qinhuangdao coast, therefore, the tide range, with a maximum of 1.5 m and annual average of 0.74 m, is lower than that of most other coastal areas of China. According to the wave data of Qinhuangdao, shown as wave rose in Figure 1d derived from the measurements at the wave station located at 119 • 28'07" E and 39 • 47'08" N (W1 shown in Figure 2a) from March 2017 to December 2017, 79.33% of the wave height is mainly distributed between 0.1-0.5 m, and SSE, SE, ESE represent 54% of the prevailing wave directions of the time and approximately perpendicular to the coastline.
The coastline of Shanhai pass stretches roughly in the direction of WSW to ENE with a natural headland between the sandbars of S2 and S3. It consists of beaches with an average median grain size (D 50 ) of 0.055 mm. Shihe river (Figure 1a) is the major river connected to this coastal region, and supplied up to 9.03 × 10 4 tons of sediment per annum to the beach system before the construction of Shihe reservoir in 1975 [55]. Since the sediment carried by the river flow deposited at the estuary due to the channel expansion and the slowdown of flow, a river mouth bar formed here under the tide forcing. Driven by continuous sediment supply upstream, the mouth bar grew into an island called South Shihe Island and the river bifurcated into two river branches entering the sea. In 1975, however, Shihe reservoir was established to prevent flood and supply water resource for Qinhuangdao city. Thereafter, the sediment flux from Shihe dropped dramatically to 0.72 × 10 4 tons per annum and the average suspended sediment concentration near the estuary was reduced to less than 0.043 kg/m 3 [55]. In addition, other anthropogenic interferences such as construction of a rubber dam upstream and sediment dredging further decreased the sediment supply. As a result, the coastal erosion is exacerbated, leading to coarser sediment, shrinking beach width and deteriorating coastal environment. Since the wave and tidal energy is relatively low in the semi-closed bay study area and storm surge rarely attacks this region, it is worthwhile to investigate the coastal erosion over a long period of time such as 2-5 years in order to assess the feasibility and service life of the Shanhai Pass beach nourishment project under the relatively calm environmental forcing. In the past decades, the booming economy, growing population, and rapid urbanization in China has brought unprecedent influx of tourists and economic benefits to Qinhuangdao for its long stretched sandy beaches including Shanhai pass. To protect Shanhai pass from coastal erosion, a beach nourishment project was carried out at the end of August 2017. The main objective of this project is to replenish sand onto the beach berm and the three subaqueous artificial sandbars (S1-S3, shown in Figure 1a) along a 2 km long shoreline. Based on the geographical condition of the study area, the design elevation of beach berm was set to 1.6 m with respect to 85 datum (1985 National Elevation Datum of China). The intersecting profile [56] was applied as the designed profile with the slope less than 1:100 for the beach berm and the berm width was extended by 40-60 m along the coastline. The medium-coarse sand with a median grain size (D 50 ) of 0.42-0.61 mm was selected as the filling materials for the beach berm. Three artificial submerged sandbars were installed parallel to the shoreline approximately 200 m off the coast. Each sandbar was about 200 m long and 80 m wide at the bottom with the crest elevation/freeboard of −0.9 m below the still water level and the crest width of 50 m. The medium-coarse sand with D 50 of 0.5-2 mm was selected for the artificial sandbars. The total replenishment sediment of around 410,000 m 3 was used in the project.

Methodology
After the nourishment project was completed, field observations of the topographic and bathymetric profiles, tidal currents and waves were collected and used as the input for the wave-current model coupled with the sediment transport and morphological model based on Delft3D [57]. The tidal currents were measured by a direct-reading current meter continuously for 25 h at a sampling interval of 1 h and the wave data were measured by a wave buoy shown as W1 in Figure 2a. The beach profiles were surveyed differently on land and underwater as described in Section 3.1. Delft3D is an open-source software package, including FLOW, WAVE, SED, and MOR modules, which was developed by Deltares, the Netherlands for calculations of currents, waves, sediment transport and morphology in coastal, river and estuarine areas. Next, we will introduce the field measurement method of the beach profiles used for directly morphological analysis and the numerical model setup and validation for further analysis of the morphological modeling results.

Field Observations of Beach Profile
Beach erosion-accretion and source-sink of sediment may be extracted from regular beach profile surveys. A conventional beach profile survey comprises the offshore survey by a specialized boat and the inshore-land survey carried out using a level and transit survey equipment or a total station instrument. An accurate benchmark with definite coordinate and elevation is required to determine the position and elevation of the measuring points along the beach profile. A beach profile surveys from inshore to land is typically conducted in the cross-shore direction perpendicular to the shoreline. Sometimes, one must wade or even swim to measure the elevation of the upper beach along the profile into deep water as far as possible to overlap the offshore measurements, so the measurement is usually done at low tides. The offshore beach profile data are usually collected using a survey vessel equipped with a fathometer and GPS system so that the position of the boat can be correlated with the depth measurements. The beach profile surveys extend offshore into a water depth where there is little change in the profile, i.e., depth of closure, in order to provide a reliable estimate of sediment gains and losses between surveys.
In this study, the survey staff used the state-of-art measuring equipment with RTK (real-time kinematic) technique instead of level and transit survey equipment or a total station. After the beach nourishment, they carried out both the inshore-land and offshore surveys 5 times over 21 cross-shore sections (P1-P21, shown in Figure 1a) from August 2017 to July 2018. RTK applies a carrier-phase tracking method to achieve a centimeter accuracy. From inshore to on-land, the survey was operated by RTK mobile station to determine the height, while for offshore survey the combined technique of RTK and sounder with the centimeter accuracy was applied on a survey vessel. The measuring points are set at a spacing from 5 to 350 m. The beach profile changes after the beach nourishment are then analyzed and applied to validate the model.

Model Description
In order to investigate the morphological evolution mechanism (under the concurrent action of tides and waves) after the beach nourishment in Shanhai Pass, a morphodynamic model based on Delft3D [57,58] is employed. The hydrodynamic module Delft3D-FLOW simulates 2DH unsteady flow and transport phenomena under multiple forces. This module comprises three basic governing equations, i.e., the momentum equation, the continuity equation, and the transport equation. Delft3D-FLOW is the cornerstone of this program software suite and can be fully coupled with other modules such as Delft3D-WAVE, Delft3D-SED, and Delft3D-MOR to simulate morphdynamics. Delft3D is particularly suitable for modeling hydrodynamics in shallow water environments where the horizontally spatial scale is remarkably larger than the vertical scale, so that the vertical momentum equation is reduced to the hydrostatic pressure relation neglecting the vertical acceleration. Delft3D-WAVE can be one-way or two-way coupled with Delft3D-FLOW by the so-called communication file containing wave model results (i.e., RMS wave height, peak spectral period, wave direction, mass fluxes) from Delft3D-WAVE and hydrodynamic model results (i.e., water level and flow vector). Given the shallow depth of the study area, the 2DH model used in this study compromises the computational accuracy and efficiency to reproduce the hydrodynamic climate, sediment transport and morphological variations from 1 to 5 years. The detailed descriptions of the Delft3D-FLOW and Delft3D-WAVE modules can be referred to in other documents [57,58], and the Delft3D-SED and Delft3D-MOR for the sediment transport and morphological modeling are described in Appendix A.

Model Setup
Two grids, one nested in the other, were set up to resolve the regional scale of Bohai sea domain and local scale of Shanhai pass ( Figure 2) in order to balance the computational accuracy and efficiency for ocean circulation modeling. The regional grid domain (Figure 2a Figure 2b) ranges from 18 m to 239 m in the longshore direction and from 19 m to 247 m in the cross-shore direction with higher resolution in the study area and lower resolution in the outer area to capture the detailed processes near the coastline and coastal structures in the study area. The bathymetry of the Bohai coarse grid was derived through an interpolation of the latest bathymetric data from the admiralty chart in the past decade and the bathymetry of the nested fine grid at Shanhai pass was measured by the Eighth Geological Brigade of Hebei Geological Prospecting Bureau in 2017. The grid system for wave simulation employed two grids, i.e., the larger one covering the same area as the Qinhuangdao domain and the nested one covering the same area as the Shanhai Pass domain. The computational time steps of the Bohai regional model and the nested local model of Shanhai pass were set to 60 s and 30 s respectively to keep the Courant number under 10 to achieve both computational stability and accuracy.
The linear interpolation of the measured tidal levels at Dalian and Yantai tidal stations provides the offshore boundary condition of tidal level for the Bohai regional model. The model results of the Bohai are used as the boundary conditions to drive the Shanhai Pass local model. As for wave modeling, the wave boundaries of the larger grid are extrapolated from the measured wave data of W1 wave station and its computational results provide the wave boundary conditions for the nested grid. The wave model was two-way coupled with the flow model at the Shanhai pass grid with the real-time wave conditions. The bottom roughness is represented by the Manning coefficient which varies with the water depth and ranges from 0.012 to 0.018 m −1/3 ·s in the Bohai model. The Manning coefficient for Shanhai Pass model is specified as a constant of 0.0145 m −1/3 ·s following Kuang [59]. Since sediment transport properties such as incipient motion and settling velocity are dependent on grain sizes, it is not advisable to assume uniform sediment grain size in the model runs. According to the filling sand and the original beach sediment characteristics of the study area, the morphological model adopted 3 sediment fractions to represent the background sediment and the filling sand for the beach berm and the artificial sandbars. Both cohesive sediment and non-cohesive sediment are incorporated in the model, the former performs as the suspended load transport solely and the latter performs as simultaneously suspended load and bedload transport under hydrodynamics. The coarse sand with the median grain size D 50 of 0.515 mm and the very coarse sand D 50 of 1.25 mm are specified in the model to represent the non-cohesive filling sand at the beach berm and the sandbars, respectively. The coarse silt with D 50 of 0.055 mm is specified as the cohesive background sediment for the original condition of the beach. Since D 50 of the background sediment represents the sediment grain size of the entire study area at Shanhai pass which covers a large range of grain size, it is up to an order of magnitude smaller than D 50 of the filling sand to represent the larger proportion of fine sediment in the offshore region. The settling velocity is defined as 0.003 m −1 ·s for 0.055 mm sediment fraction, 0.073 m −1 ·s for 0.515 mm sediment fraction, and 0.134 m −1 ·s for 1.25 mm sediment fraction. The critical shear stress for erosion and deposition of 0.055 mm sediment fraction are calibrated as 0.04 N m −2 and 0.02 N m −2, respectively, with the erosion coefficient of 0.0001 kg·m −2 ·s −1 . The underlayer model uses the approach of layered bed stratigraphy with a transport layer on the top and three underlayers underneath it. The thickness of the transport layer is determined as 0.1 m referring to Huisman et al. [47] for field case modeling. In order to represent different sediment fractions at the beach berms, the sandbars and the background beach, the thickness of the top underlayer is defined as the nourished beach berm height at the beach berms with the beach berm sediment fraction and the thickness of other area is set to zero. Similarly, the thickness of the median underlayer is defined as the artificial sandbar height at the sandbars with the sandbar sediment fraction and the thickness of other area is set to zero. The thickness of the bottom underlayer is defined uniformly as 2 m throughout the whole area with the background sediment fraction. During the bed stratigraphy modeling, it is assumed that the erosion rate is proportional to the availability of the sediment fraction in the top-most layer of the bed stratigraphy. Sediment fractions of different sizes influence each other by means of burial and exposure: fine sediments buried among coarse sediments are partly shielded from the main flow while the coarser sediments are more exposed to the flow than they would be among other sediments of the same size. This effect is taken into consideration by increasing and decreasing the effective critical shear stress for fine and coarse sediments, respectively. For a longer morphological prediction of the study area in 2 and 5 years, the morphological time scale factor is implemented to accelerate the bed-level changes.

Model Validations
The morphodynamic model was validated by the measured tidal levels, currents, waves and the beach profiles. The Nash-Sutcliffe efficiency (NSE) method, based on the ratio between the model-data residual variance and the measured data variance [60], was employed to quantify the model performance: where Y obs i is the ith observed value, Y sim i is the ith simulated value, Y mean is the mean of the observed values, and n is the total number of observed values.
When NSE > 0.75, >0.36, and <0.36, the simulation results are deemed good, satisfactory, and unsatisfactory, respectively, according to Motovilov et al. [61]. Figure 2 shows the time history of the observed tidal level at two tidal stations at Qinhuangdao (T1) and Shanhai Pass (T2). The observed current velocity and direction at Shanhai Pass tidal station measured by the current meter (direct reading with a mechanical propeller) are shown in Hence, it is demonstrated that the numerical model is capable to capture the observed current, wave climates, sediment transport and morphological evolution and thus can be adopted to predict the hydrodynamics and morphological responses to the beach nourishment.

Field Observation and Model Results
The field measurement was carried out to obtain the beach profiles of the study area to capture the morphological evolution. Three typical profiles (i.e., P7, P12, and P17) traversing the correspondingly three sandbars (S1-S3) were selected to analyze the morphological evolution along the profiles with the emphasis on the movement and the deformation of the three artificial sandbars. They were measured 5 times between August 2017 and July 2018, during which no sizeable storms occurred. Then, the verified hydrodynamic and morphodynamic model based on Delft3D was run for one year to investigate the responses of tide, wave and seabed with and without the beach nourishment implementation at both beach berm and submerged shoreface. The combination of model and observed results are used to pinpoint the key mechanisms behind the morphological changes at the site due to the beach nourishment.

Observed Beach Profiles
In Figure 6a, beach profile P7 passing through S1 was selected to analyze the temporal and spatial evolution of the morphology of S1 in one year. At the beginning of the beach nourishment in August 2017, the profile of S1 was shaped roughly as a symmetrical trapezoid. There was a small sandbar about 50 m away from S1 towards the shore with smaller height and slope than S1. There was also a small sandbar with a lower crest height on the offshore side of S1. In October 2017, the crest height of sandbar S1 was reduced by approximately 0.8 m (at Distance = 1840 m), which suggests that the sediment on the sandbar crest was transported shoreward and deposited on the onshore side of the sandbar as illustrated later by the sediment transport in Figure 11, causing the sandbar crest to move about 10m towards the shore. The profile of S1 evolved from a symmetrical form to an asymmetrical form with its seaside face flattened and its shoreside face steepened. The small sandbar on the onshore side of S1 also experienced a small amount of erosion on the crest, while the small sandbar on the offshore side of S1 disappeared. In March 2018, S1 continued to be eroded, and the shape of sandbar was flattened further with the crest lowered by about 0.1m from the original level. The small sandbar on the onshore side of S1 also disappeared and the two troughs on its both sides accreted, leading to a smoother topography behind S1. Profile P7 after March 2018 hardly changed, which suggested that P7 became basically stable at this time and reached an equilibrium state.
Similarly, beach profile P12 passing through sandbar S2 was used to derive the temporal and spatial morphological evolution of S2 in one year. In the early stage of the project in August 2017, the morphology of S2 and S1 was similar and both roughly trapezoid. A small sandbar was located on the offshore side of S2, while there was only a trough on the onshore side of S2. In October 2017, the height of sandbar S2 was reduced by about 0.7 m (at Distance = 2195 m) and the crest moved nearly 20 m towards the shoreline. However, the small sandbar on the offshore side of S2 disappeared completely. In March 2018, the erosion at crest of S2 and accretion on the onshore side of the sandbar continued, so that the sandbar was widened from 50 m to nearly 150 m. Similar to P7 after March 2018, there was little change in morphology, which suggests that P12 seems to have reached an equilibrium state.
The evolutions of P7 and P12 are consistent with each other to some extent. In contrast, the morphological changes of P17 passing through S3 show a different behavior. Since P17 is within the semi-enclosed waters surrounded by the Shanhai Pass shipyard (shown in Figure 1) where the hydrodynamics are weak, S3 was designed to be smaller and further away from coastline than S1 and S2 considering the cost and construction period. As a result, the beach berm in the wave shadow zone of S3 is eroded more than P7 and P12 in the intertidal zone (between mean high water (MHW) and mean low water (MLW)) due to less wave energy dissipation by S3. In July 2018, the crest height S3 was decreased by 0.7 m and merged with its landside morphology.
The migration of the sandbars under mild waves is similar to the barred beach in Senigallia, Italy [62]. The Senigallia beach is also located near a river estuary and a rigid jetty. However, unlike the Shihe river, the river in Senigallia supplied a considerable amount of sediment, causing the nearby sandbars to increase. The measurements from 2016 to 2019 [63] show that the migration of the sandbars is sensitive to seasonal changes in wave climate.  Figure 6b shows the significant wave height predicted by the wave model, along P7, P12, and P17 under the typical strong wave condition in September 2017 shortly after the start of the beach nourishment with SE direction and significant wave height of 0.65-0.75 m at 400-500 m off the shore. The waves propagating over S1 and S2 broke first at the sandbar and again on the beach face until the swash zone. S1 and S2 remarkably attenuated the significant wave height by remarkable 36% and 52%, and the residual wave energy was dissipated on the beach. However, S3 with lower crest did not effectively attenuate the waves, which resulted in the severe erosion on the nearby beach berm. Figure 6c presents the sandbar crest elevations of S1, S2, and S3 (at Distance = 1840 m, 2195 m, and 2196 m respectively) in one year and the corresponding change rate. The trends of S1, S2, S3 crest change rates are similar, and they had the largest change rate (−0.42, −0.32, and −0.21 m/month) from August 2017 to October 2017 while remained almost unchanged (less than 0.1 m/month) from October 2017 to July 2018. The sandbar rear (i.e., the onshore side of the sandbar) elevations of S1, S2 and S3 (at Distance = 1795 m, 2160 m and 2165 m respectively) in one year and the corresponding change rate are shown in Figure 6d. The elevations of S1 and S2 underwent the process of first increasing and then decreasing illustrating the migration of the sandbar crest to the original sandbar rear position. Similar to the change rate of the sandbar crest elevation, the sandbar rear elevation change rate at S1 and S2 were the largest (0.98 and 0.59 m/month) from August 2017 to October 2017 and remained almost unchanged (less than 0.15 m/month) from October 2017 to July 2018. Unlike S1 and S2, S3 rear elevation changed little throughout this one year. It can be concluded that there is little change in morphology of artificial sandbars one year after the installation, reaching an equilibrium state. It can be seen that the profiles of S1 and S2 are similar in both shape and size. As illustrated by Eichentopf et al. [64], different initial beach profiles evolve towards the same final beach configuration under the same wave condition. It explains well the similarity of the quasi-equilibrium profiles in July 2018 for S1 and S2 due to the approximately identical wave forcing acting on them. However, the profile of S3 in July 2018 is different from that of S1 and S2 because the wave condition changes at S3 due to the refraction and reflection effect of the headland between S2 and S3.
The coastline changes relative to the baseline in August 2017 in Figure 7 indicate the influence of sandbars on the coastline changes. In fact, the residual current, shown later in Figure 12c, represents a northeast longshore current. The coastline on the upstream side of the longshore current from x = 0 m to x = 400 m of S1 was eroded severely up to~20 m landward before reaching an equilibrium due to the northeast longshore sediment transport driven by the dominant SE-SSE waves (Figure 1d). The coastline in the upstream side from x = 0 m to x = 400 m of S1 initially eroded severely up to~20 m landward before reaching an equilibrium due to the northeast longshore sediment transport driven by the dominant SE-SSE waves (Figure 1d). The coastline between S1 and S2 (from about x = 700 m to x = 800 m) where the hydrodynamic forcing is low-energetic under the sheltering effect of the two sandbars, migrated seaward continuously with time and widened the beach up to 25 m due to the sediment captured from the upstream longshore sediment transport and the sediment feeding from the two sandbars under the dominating SE-SSE waves. Similarly, the stretch of shoreline between x = 1000 m and 1200 m migrated seaward at first while retreated afterwards. The coastline at x > 1200 m retreated significantly in the same magnitude as the 0-400 m stretch, which is exposed to the SE-SSE dominant waves combined with the wave convergence around the natural headland between S2 and S3. Throughout the whole coastline stretch, the coastline mostly reached the equilibrium state except for the segment between x = 1000-1200 m. The unit-width volumes of profiles through three sandbars in five measurements (i.e., the area under the profile line shown in Figure 6a) in one year are presented in Table 1. The final profile volume changes of P7 and P17 after one year decreased while that of P12 increased, leading to the erosion on both sides and accretion in the middle of the triple-sandbar system. Standard deviations of profile volumes along P7, P12, and P17 are 14.85, 31.81, and 38.82 m 3 /m, respectively, which are relatively small, demonstrating the sediment volume conservation of the three profiles.

Effects of Sandbar on Tidal Currents and Waves
The bathymetric change due to the construction of the Shanhai Pass project would alter the spatial and temporal patterns of tidal current and waves, which would further influence the sediment transport and morphological evolution. In order to investigate the impact of the beach nourishments on the hydrodynamics, a flood-ebb cycle of spring tide (12 October 2017) was selected as the analysis period. The computational current field of the study area with or without the construction of the sandbars and their difference are shown in Figure 8. At the maximum flood of spring tide, the current fields are approximately the same in Figure 8a,c. The velocity is about 0.4 m/s about 3-4 km offshore, and the velocity decreases to less than 0.1 m/s about 1 km off the shoreline. Due to shallow water depth and the wave-induced current caused by the wave breaking, the peak velocity in the southern area near the South Shihe island (shown in Figure 1) increased above 0.4 m/s. The Shanhai Pass shipyard and the coastline to the north of Shihe estuary forms the semi-closed waters where the nourishment project was implemented. The current direction out of the semi-closed waters at the maximum flood of spring tide is west-southwest, roughly parallel to the coastline. At the inlet of the semi-closed waters the current deviates into another branch in the direction of north with the decreased velocity. The current field of the semi-closed waters is low-energetic on the whole with a low average velocity magnitude (0.05 m/s). The current fields resemble each other in Figure 8b,d at the maximum ebb of spring tide. The current direction out of the semi-closed waters is northeast opposite to that at the maximum flood. Unlike the maximum flood, near the north branch of Shihe river and S1, there exists a counterclockwise circulation in the maximum ebb with the current flowing south off S1. When the ebb current flows northeast against the long breakwater of the shipyard, it is blocked by the breakwater and partly flows north into the semi-closed waters. The increase of the water volume in the semi-closed waters drives the longshore current in the southwest direction to form the observed circulation.
The tidal velocity differences with and without the beach nourishment are illustrated in Figure 8e,f. At the maximum flood of spring tide, the velocity increases by 0.011 m/s, 0.01 m/s, and 0.006 m/s at sandbar S1, S2, and S3, respectively, due to the water depth reduction by the bathymetry uplift. The tidal current is weakened slightly behind the sandbars. The current velocities at all three sandbars at the maximum ebb of spring tide (Figure 8f), are comparable to that at the maximum flood of spring tide. Similarly, the current velocity is reduced by about 0.006 m/s averagely behind the sandbars. These model results suggest that the current under both tidal and wave forcing in this study area is small because the sandbars are located within the calm semi-closed waters surrounded by the coastline and the Shanhai Pass shipyard.
To investigate the wave climate changes due to the beach nourishment, the typical normal wave and strong wave during September to October in 2017 were selected to obtain the normal and strong wave height and their difference in the study area, as shown in Figure 9. Under the normal wave condition, the incident wave height of~0.7 m with the direction of SE-SSE offshore is reduced to~0.55 m at the sandbars, while under the strong wave condition, the height of the offshore incident waves of~1.2 m with the direction of SE drops to~0.65 m at the sandbars (in Figure 9a,b). It can be seen in Figure 9e that the wave height decreases by about 0.05 m, 0.15 m, and 0.1 m at S1, S2, and S3, respectively, with the average wave attenuation coefficient of 20% under the normal wave condition. Figure 9f shows that the wave height decreases by about 0.15 m, 0.25 m, and 0.2 m at S1, S2 and S3 respectively with the average wave attenuation coefficient of 30% under the strong wave condition. In both normal and strong wave conditions, the three sandbars effectively protected the beach by reducing the wave energy significantly.

Morphological Responses to Beach Nourishment
The beach nourishment of beach berm and three artificial sandbars along the coast have caused significant change in topography and bathymetry of the study area. The measured beach profiles indicate that the major erosion and accretion are concentrated in the sandbar area, therefore it is necessary to carry out a further investigation focusing on this area. Using the verified morphological model, the morphological erosion and accretion without and with the beach nourishment and the elevation difference between with and without beach nourishment in 1, 2, 3, 6, and 12 months after the beach nourishment are shown in Figure 10.    Figure 10a shows the general morphological processes in the nearshore area without the beach nourishment in one year. On the whole, less than 0.1 m erosion occurred in the first two months. Then, severer erosion took place near the headland with the maximum value of 0.78 m in one year. Without the beach nourishment, the beach was eroded significantly especially near the headland. However, one month after the completion of the beach nourishment (shown in Figure 10b), the beach berm deposited slightly (white patch in Figure 10b representing accretion from 0-0.4 m) with its lower side regions eroded, which means that the sediment was transported onshore. The sandbars started to erode with the landward sediment migration, which is consistent with the measured profiles. One month later the beach berm accreted more, and the sandbars was eroded further. Three months after the beach nourishment, the accretion of the beach berm and the erosion of its offshore side were more significant, and the artificial sandbars eroded more. The maximum erosion depth at the crest of S2 was up to 1.2 m, and a triangular sedimentation area was formed in the rear area just adjacent to S2. The cross-shore width of the sandbar was extended to more than 120 m and the sedimentation thickness was between 0.2-1 m. The erosion was severe at both ends of S1 while mild in the middle. The erosion depth was about 0.8 m at the crest of S1 and about 0.4-0.8 m of deposition occurred at the onshore side of the sandbar. Similar to S2, the cross-shore width of S1 was also expanded to its sheltered area by more than 100 m, which almost joined the beach. Half year after the beach nourishment, the erosion and accretion of the three sandbars began to slow down. Compared with three months after the beach nourishment, the bathymetry in the sandbar area half year after the beach nourishment did not change much, and the major changes were concentrated on the shoreward side of sandbars with continuing deposition. One year after the beach nourishment, the erosion and accretion of S1 and S2 changed very little. Figure 10c show the elevation difference between scenarios with and without the beach nourishment, in which the difference was reduced at the sandbars while increased near the coastline with time. It means that the sandbars were eroded, and its morphology was close to the natural condition while the sediment from them fed the beach behind. In summary, after one year, a relatively strong erosion occurred at the sandbars, and the maximum scour depth was more than 1.2 m. The maximum deposition occurred in the sheltered area behind the sandbars, with the maximum deposition depth of over 1.2 m. The resulting change in the water depth would alter the wave transformation across the sandbar which in turn affects the sediment transport and morphological change in this area. The bathymetry of the study area achieved equilibrium in about one year after the beach nourishment. The numerical modeling of morphodynamics with and without the beach nourishment was carried out in this section to figure out the erosion hotspots in the nourished area for the subsequent adjustment in the future project. In fact, physical models performed by Atkinson et al. [65] were also applied to optimize the beach nourishment project as well.

Driving Force for Transport of Multi-Fraction Sediment
In the sediment transport and morphological model, three sediment fractions were configurated to investigate their different transport mechanisms. The coarse sand with median diameter D 50 of 0.515 mm and the very coarse sand D 50 of 1.25 mm are specified in the model to represent the filling sand at the beach berm and the sandbars respectively, and the coarse silt with D 50 of 0.055 mm is specified as the background sediment for the origin condition of the beach. In order to investigate the transport of each of these three types of sediment, the predicted net sediment flux of each of them in a tidal cycle involving the strong waves (the same as in Figure 9b) are shown in Figure 11. The sediment transport at the beach berm close to shoreline is mainly shoreward. Smaller offshore and longshore transport in the northeast direction is evident in some areas as well. Due to the reduced wave height and current by the sandbars indicated in Figures 8 and 9, the net sediment flux of beach berm sediment is relatively small with the maximum value of 0.012 kg/(m·s). The shoreward sediment transport at the sandbars accounts for the major amount of the sediment transport and the offshore transport is negligible. Due to the larger wave and stronger wave-breaking induced current at the sandbars, the onshore net sediment flux of sediment in the sandbars is almost a magnitude larger than that at the berm, with the maximum value above 0.12 kg/(m·s). The net transport of the background silt in the study area is mainly in the northeast direction along the coast driven by the longshore current. In the northwest vicinity of S2 and S3, there is a large sediment transport along the coast, with the maximum net sediment flux up to 0.018 kg/(m·s). Overall, the sediment transport at sandbars is larger than that at beach berms. On the offshore side of the sandbars, the primary sediment transport of silt is in the northeast direction; at the sandbars and their sheltered areas, the main direction of sediment transport of medium-coarse sand is onshore, while the silt is transported in the northeast alongshore direction.
In order to investigate the driving force of multi-fraction sediment transport, two simulations carried out: one was run by real-time coupled tides and waves, i.e., tide-wave simulation, another was run solely by real-time waves i.e., wave-only simulation. A period of tidal cycle on 27 September 2017 is selected for analysis, during which the strong wave the same as that in Figure 9 occurred. The wave energy at the strong wave moment is much larger than that at the other moments in the selected period. The instantaneous current fields for the tide-wave simulation and the wave-only simulation at the strong wave moment during the selected tidal cycle period are shown in Figure 12a,b, respectively.
The residual current fields for the tide-wave simulation and the wave-only simulation during the selected tidal cycle period are shown in Figure 12c,d, respectively.In order to investigate the driving force of multi-fraction sediment transport, two simulations carried out: one was run by real-time coupled tides and waves, i.e., tide-wave simulation, another was run solely by real-time waves i.e., wave-only simulation. A period of tidal cycle on 27 September 2017 is selected for analysis, during which the strong wave the same as that in Figure 9 occurred. The wave energy at the strong wave moment is much larger than that at the other moments in the selected period. The instantaneous current fields for the tide-wave simulation and the wave-only simulation at the strong wave moment during the selected tidal cycle period are shown in Figure 12a,b, respectively. The residual current fields for the tide-wave simulation and the wave-only simulation during the selected tidal cycle period are shown in Figure 12c,d, respectively.   Figure 12b shows the breaking waves generates strong NW currents at S1 and S2, leading to the increase of water volume on the shoreward side of them which cause the offshore currents in the south vicinity of S1 and the trough between S1 and S2, and the northern alongshore current in the northward vicinity of S2. The current field of Figure 12a presents a similar pattern to Figure 12b in the nearshore area behind the sandbars, which means that the nearshore current under the coupled tides and waves is dominated by the wave-driven current. The magnitude of the nearshore current behind the three sandbars in Figure 12a is smaller than that in Figure 12b, possibly attributed to the high water lever driven by the tide at this moment. Figure 12c shows that the residual current in the direction of NE is consistent with the transport direction of background sediment (coarse silt), which is mainly transported along with the current as suspended load. However, the direction of the residual current is different from that of the beach berm and the sandbar sediment transport. The beach berm and the sandbar sediment transport directions are in line with the strong wave direction (shown in Figure 9d). The sediments of the beach berm and the sandbars are coarse sand and very coarse sand, which are more inclined to move as bedload than as suspended load. In fact, as shown in Equations (A9) and (A10) in the Appendix A, the direction of bedload transport is determined by the combination of the current and wave direction. Therefore, the similarity between direction of the beach berm and the sandbar sediment transport and the strong wave direction demonstrates that the strong wave plays a major role in the beach berm and sandbar sediment transport. At the three sandbars, the filling sand movement is driven mainly by the SE strong waves (Figure 9b), which is consistent with sediment transport direction in NW. At the beach berm, the majority of the beach berm sand migrates onshore and offshore in the direction perpendicular to the shoreline, which is consistent with the nearshore wave direction due to refraction. The onshore transport of the berm sand is mainly driven by waves as the wave-related bedload shown in Equations (A9) and (A10). The offshore transport is mainly driven by currents as the current-related bedload shown in Equations (A9) and (A10), and additional suspended load due to D 50 of berm sand is much finer than that of sandbars. Figure 12c,d show that the residual current field is generally in accord with the wave-induced current field, which implies that the residual current in the study area is dominated by the wave-induced current.

Regional Sediment Budget
To quantify the sediment transport and the spatial variation of erosion and accretion, six elements are defined which contain 3 sandbars and 3 rear sections directly behind the sandbars (shown in Figure 13a). The numbers and arrows in Figure 13 indicate the net sediment transport volume and direction between elements within a year after the beach nourishment and the value less than 0.1 m 3 are not marked in the figure. It can be seen from Figure 13b that the beach berm sediment is mainly transported alongshore whereas the quantity is relatively low compared with the sandbar sediment and background sediment. The sandbar sediment is mainly transported shoreward with the volume flux out of S1, S2 and S3 of 549.65 m 3 , 901.44 m 3 and 296.27 m 3 respectively, and the sediment transports from the seaside into Element2, Element4, and Element6 are 73.83 m 3 , 99.38 m 3 and 172.29 m 3 respectively (shown in Figure 13c). The background sediment is mainly transported along the coastline with a small portion transported offshore. In sum, the accretion behind the sandbars is caused by the shoreward transport of sediment from the sandbars. The erosion at the beach berm is mainly caused by the loss of background sediment, and the filling beach berm sediment is hardly changed. The background sediment in the nearshore zone is transported offshore, and the filling sand at the beach berm and the sandbars is concentrated in the replenishment site and its surrounding area.
The cumulative erosion/deposition in one year, erosion/deposition rate and average erosion/deposition thickness per unit area of each element are shown in Table 2. The total erosion volume of the sandbar area in one year after the beach nourishment is 5415 m 3 , the average erosion rate is about 15 m 3 /d, and the average erosion thickness is less than 0.02 m. The average deposition thickness of Element1 is about 5 mm/m 2 and the average erosion thickness of Element3 and Element5 is about −14 mm/m 2 and −11 mm/m 2 , which is less than the average value of total 6 elements (i.e., about −16 mm/m 2 ). Hence, it can be concluded that the beach in the study area has been effectively nourished and protected within one-year interval.

Long-Term Morphological Change and Model Improvement
The beach nourishment project consisting of sand supplement to the beach berm along the coast and three artificial submerged sandbars offshore has caused a significant change in the bathymetry of the study area. The morphological changes will inevitably disrupt the sediment transport balance, and therefore reconstruct the seabed until it reaches a new equilibrium. In addition, the study of the long-term evolution of the morphology is critical to determine the service life of the beach nourishment and the subsequent maintenance after the beach nourishment. The measured profiles indicate the hot spots of erosion and deposition are confined to near the sandbars and the shielded area behind it, which will be discussed later.
According to the verified sediment transport and morphological model based on the hydrodynamic condition of August 2017 to July 2018 without the introduction of storms, the erosion and accretion in scenario with the beach nourishment and the elevation difference between scenarios with and without the beach nourishment in 2 and 5 years after the completion time of the project were obtained and shown in Figure 14. Though the measured profiles of P7, P12, and P17 and the one-year morphological modeling of the study area suggest that the topography becomes stable one year after the beach nourishment, the erosion and accretion occur alternately in some locations in 2 and 5 years after the beach nourishment. From the predicted results of the morphological model in Figure 14a,b, the seabed of the offshore area is hardly changed, the three sandbars continue to move towards the shore, and the subaerial topography of the beach berm is hardly changed. It is consistent with the typical migration pattern of the sandbar under normal conditions (in calm seasons other than storms). As shown in Figure 14c,d, the elevation at the sandbar and beach berm area is still larger than the natural topography 5 years after the beach nourishment, indicating its long-standing effect. In the area shielded by the sandbars, salient came into existence, which were accreted by over 3m behind S2 5 years after the beach nourishment, protruding from the coastline into the sea due to the combination effect of the captured longshore transport and the feeding sand from the artificial sandbars. Furthermore, the area at sandbars showed obvious erosion in 5 years, and the maximum erosion depth is up to 2 m at S2.
As mentioned before, many nourishment projects have been implemented in Qinhuangdao to conserve the beach resources. Artificial sandbar, however, becomes more and more popular in beach nourishment for its eco-friendly nature. However, there is a lack of in-depth studies to interpret its deformation and migration mechanisms under various hydrodynamic conditions. On the basis of the previous studies, Kuang et al. [37] established a hydrodynamic-sediment-morphological model driven by the real-time series of tide and wave boundary condition to model the morphological processes of two artificial submerged shore-parallel sandbars in a nourishment project at Tiger-Rock Beach in Qinhuangdao. Although the model compares well with the measured tide and wave and reveals some important processes related to the migration of the artificial submerged shore-parallel sandbars, the model failed to capture the evolution of the backshore without consideration of the medium size nourished sand due to the limitation of MIKE MT module designed for simulating the transport of the cohesive sediment with grain size less than 63 µm and fine non-cohesive sediment with grain size of 63-125 µm. In the present 2DH model, three sediment fractions including coarse sand, very coarse sand, and coarse silt are incorporated in the morphodynamic model based on Delft3D to better simulate the sediment transport with different grain sizes. Despite that the present 2D model reproduce the onshore migration of sandbar under normal waves, it may not be adequate during storm conditions when the vertical variation of wave-induced current becomes important such as undertow due to wave breaking which contributes to either generation or offshore migration of the sandbar. In this case, to resolve the vertical variation of velocity and the associated sediment transport, a 3D or a qusi-3D modeling is required. The field observations by Elgar et al. [23] suggest that onshore sandbar migration, observed when breaker-driven mean flows are weak, may be related to the skewed fluid accelerations associated with the orbital velocities of nonlinear surface waves. Significant evolution of wave skewness and asymmetry has been observed on natural beaches by Elgar et al. [66] and over low-crested structures such as sandbars and breakwaters [67] and due to bottom slope and friction [68]. Using 1D wave-resolving eddy-diffusive vertical profile model of water and suspended sediment motion in the bottom boundary layer forced by measured pressure gradients during Duck 94 experiment, Henderson et al. [69] found that seaward (shoreward) bar migration was driven by a maximum in the current-generated (wave-generated) flux over the sandbar and the wave-generated momentum fluxes and that the Stokes drift were essential to the predictions of shoreward bar migration. A phase-resolving 2DV/3D hydrodynamic and morphological model such as that in Peng et al. [70] will be implemented to resolve these processes during a storm in the future study.

Conclusions
In this study, both field measurements and numerical modeling are employed to investigate the morphodynamic evolution of the nourished beach with submerged sandbars. Five field surveys were carried out along each of the 21 cross-shore sections of the beach from August 2017 to July 2018 from water depth −4 m to 2 m. The Delft3D modeling framework, including fully two-way coupled flow-wave model and sediment transport and morphological model, were verified by the field observations to simulate the flow and wave climate, sediment transport and morphological evolution at Shanhai Pass, Qinghuandao, located at the northwest coast of the Bohai sea.
The field observations of beach profiles indicate that the morphological change by the beach nourishments is negligible in the offshore area about 200 m away from the sandbars and that the most significant morphological change occurs around the sandbars and theirs onshore sheltered area. The submerged sandbars are gradually flattened and smoothened over time as sediment is entrained on the crest and the offshore side of the sandbars and deposit on the onshore side of the sandbars. One year after the beach nourishment, there are little changes in morphology which indicate the beach system has attained a new equilibrium state.
According to modeling results by Delft3D, the installation of the sandbars reduces the wave height by up to 20~30%. The predicted morphological change is consistent with the field observed severe erosion/deposition nearshore and stable bathymetry far offshore, and the bathymetry seems to reach a new equilibrium one year after the nourishment. It was found that the filling sand of sandbars is mostly transported onshore by wave-breaking induced current under strong waves, while the background sediment is transported by the residual flow in the alongshore direction, which is dominated by the wave-induced flow. The sandbar sediment is transported onshore with a larger amount than that of beach berm. The erosion and deposition pattern of sandbars in the sheltered zone in one-year after the beach nourishment is similar to that in two-and five-year after the beach nourishment but with a slightly larger magnitude. In the onshore areas sheltered by the sandbars, salients were formed protruding from the coastline into the sea due to the combination longshore transport and cross-shore transport of the feeding sand from the artificial sandbars, similar to the tombolo behind the shore parallel breakwaters.

Appendix A. Governing Equations of Sediment Transport and Morphological Model
In the sediment transport module, the components of cohesive sand and non-cohesive sand can be considered separately or together. They can be distinguished by particle size, and the critical particle size is approximately 0.063 mm. The sediment smaller than the critical particle size is taken as cohesive and transported in the form of suspended load; the sediment larger than that value is taken as non-cohesive and transported in the form of both bed load and suspended load. In the two-dimensional version of Delft3D model, the transport of suspended sediment is based on the advection-diffusion equation of average water depth: where c i is the concentration of the ith sediment component; ε s,x and ε s,y are the horizontal diffusion coefficient in the x and y directions; E i and D i are the erosion flux and deposition flux of the ith sediment component respectively, representing the sediment exchange between the water body and the bed surface. The formulations of E i and D i strongly differ for cohesive and non-cohesive sediment. For cohesive sediment fractions, the fluxes between the water phase and the bed are calculated with the well-known Partheniades-Krone formulations [71]: where M i is the user-defined erosion parameter; S(τ cw , τ cr,e,i ) is the erosion step function; ω s,i is the fall velocity; c b,i is the average sediment concentration in the near bottom computational layer; S(τ cw , τ cr,d,i ) is the deposition step function; τ cw is the maximum bed shear stress due to current and waves; τ cr,e,i is the critical erosion shear stress; τ cr,d,i is the critical deposition shear stress; the subscript i represents the ith sediment component. For non-cohesive sediment fractions, the transfer of sediment between the bed and the flow is modelled using sink and source terms acting on the near-bottom layer referred to as the reference layer. For non-cohesive sediment, E i is calculated based on upward diffusion and D i is calculated based on sediment settling, the equations of which are introduced in detail in Delft3D-Flow Manual [57].
The reference height a [72] is introduced for non-cohesive sediment. The sediment above a is treated as suspended load which is calculated following Equation (A7), while the sediment below a is treated as bedload. The reference height a is given by: where A k is a user-defined proportionality factor, k s is a user-defined current-related effective roughness height, ∆ r is the wave-induced ripple height, and h is the water depth. The magnitude of bedload sediment transport is calculated using the formula developed by Van Rijn [72] given as: where S b,i is the bedload transport, η is the relative availability of the sediment fraction in the mixing layer, ρ s,i is the dry density, ω s,i is the sediment settling velocity, D 50,i is the median grain size, M i is the sediment mobility number due to waves and currents, and M e,i is the excess sediment mobility number of the ith sediment component.
The direction of the bedload transport vector is defined by two parts, one driven by current (S bc,i ) in the direction of the near-bed current and one driven by waves (S bw,i ) in the direction of wave propagation. The magnitudes of these two parts are defined as follows: S bw,i = r S bc,i , (A7) where U on is the near-bed peak orbital velocity in the direction of wave propagation, v cr is the critical depth-averaged velocity for initiation of sediment motion, and v R is the magnitude of an equivalent depth-averaged velocity computed from the velocity in the bottom computational layer by assuming a logarithmic velocity profile. The bedload transport is the combination of the above two parts and an estimation of the suspended sediment transport due to wave asymmetry effects within about 0.5 m above the bed (S sw,i ). The bedload transport vector is defined as follows: where S bx,i and S by,i are the bedload transport components in x and y direction respectively, u bx and u by are the bottom-layer flow velocity components in x and y direction respectively, → u b is the bottom-layer flow velocity vector, f BED , f BEDW and f SUSW are the user-specified calibration factors, and φ is the angle between the direction of wave propagation and the computational grid.
In the morphological module, the morphological evolution is calculated based on the conservation equation of sediment mass, i.e., the balance between sediment transport and bed elevation change. The formula [73] is: where ρ s,i is the dry density, z b,i is the thickness of seabed, S x,i and S y,i are the sediment transport flux (involving both suspended load and bedload) of ith sediment component in the x and y directions. The total bed level change is the sum of the sediment thickness change of each component. The total bed level is updated in every time step of hydrodynamic computation cycle.