An Unstructured-Grid Based Morphodynamic Model for Sandbar Simulation in the Modaomen Estuary, China

: The Modaomen Estuary is the most important passageway in discharging ﬂood and sediment of the Pearl River Delta, which is one of the most complex estuarine systems in China. Due to the coupling effect among tidal currents, waves, and sediments, an immense sandbar area evolved in the outer subaqueous delta, impeding the ﬂood releasing during wet season, as well as salinity intrusion during the dry season. In this work, an unstructured-grid based morphodynamic model was proposed to simulate the sandbar evolution process in the Modaomen Estuary. The proposed model was constructed by using the two-dimensional shallow water equations for tidal ﬂow, the advection-diffusion equations for salinity and suspended sediment transport, and the non-equilibrium formulation of the Exner equation for bed evolution. To simulate the wave-induced longshore currents in the Modaomen Estuary, an adaptive time-stepping approach was proposed to couple the unstructured-grid based Simulating WAves Nearshore (SWAN) model and the shallow ﬂow model. An integrated solver is proposed for computing ﬂow, salinity, and sediment transport ﬂuxes simultaneously, and then the shallow water equations and advection-diffusion equations are jointly solved by a high-resolution, unstructured-grid Godunov method. Application of the model to the Modaomen Estuary, using calibrated parameter values, gives results comparable to the measured data. The butterﬂy shape of the sandbar in the Modaomen Estuary is considerably well simulated by the proposed model, which matches well with the measured topography.


Introduction
Sandbars are omnipresent patterns along river estuaries that interact with waves. Their local morphology is expressed in terms of an underwater barrier [1]. An estuary sandbar, which provides natural protection for coastal beaches by reducing inshore wave energy through depth-induced wave breaking [2], is one of the most important morphological units for deltas and estuaries. It can considerably reduce the potential extreme wave run-up during major storms, which would induce severe flood disaster in tidal zones [3]. Moreover, estuary sandbar would play important role in flood releasing during the wet season, as well as salinity intrusion during the dry season.
Sandbar behaves in significantly different ways on different timescales. On shorter timescales of a storm period, rapid offshore sandbar migration might occur due to intense wave breaking over the barrier crest [4]. On longer timescales of several years, sandbars might exhibit a cyclic progressive net offshore migration [5]. To simulate sandbar behavior on different timescales, various kinds of numerical

Governing Equations for the Circulation Model
The improved two-dimensional shallow water equations [18] in conservation formulation, further integrating the wave radiation stresses, are adopted as governing equations for the circulation model:

∂U ∂t
in which, U denotes the vector of conservation variables; E adv , G adv denote the vector of convective flux in the xand y-coordinates, respectively; E diff , G diff denote the vector of diffusion flux in the xand y-coordinates, respectively; and S represents the source terms: in which, h is the water depth; u, v are velocities in the x-, and y-coordinates, respectively; b is the bed elevation; ν t is the turbulent viscosity coefficient, and ν t = ακu * h, α = 0.2, κ = 0.41 is the von Karman's constant, where u * is the shear velocity; g is the gravitational acceleration; f is the Coriolis force coefficient; τ s = (τ s x , τ s y ) denotes the wind stress, and τ s x = C d U w U 2 w + V 2 w ρ a /ρ, τ s y = C d V w U 2 w + V 2 w ρ a /ρ, C d = min[0.0035, 0.001 × (0.8 + 0.065 U 2 w + V 2 w )] is the wind drag coefficient, where U w and V w are the wind speeds measured at 10.0 m above the sea level in the xand y-directions, respectively, ρ is the density of water, and ρ a is the density of air; S xx , S xy , S yx , S yy are wave radiation stresses; S f x and S f y are friction slopes in the xand y-directions, respectively; and S 0x and S 0y are bed slopes in the xand y-directions, respectively. The friction slopes are estimated by using the Manning formulae: (2) in which, n is the empirical Manning coefficient.

Governing Equations for Suspended Load
The following advection-diffusion equation is used as the governing equation for suspended load, in which, the superscript k denotes the k-th groups of sediment diameters; S denotes the sediment concentration; ν sx , ν sy denote the turbulent diffusion coefficients, which are calibrated constant parameters for simplification; and F s denotes the scour and deposit rate. The sediment-carrying capacity is used to represent the function of scour and deposit, as well as the equation of bed deformation, which is given by [19], in which, α 1 , α 2 denote the saturation recovery coefficients of sediment concentration and sediment-carrying capacity, respectively; α 3 denotes the settling probability; ω denotes the sediment settling velocity; S * denotes the sediment-carrying capacity; ρ s denotes the sediment density; and α 1 , α 2 and α 3 are calibrated constant parameters.

Governing Equations for Salinity
The following advection-diffusion equation is used as the governing equation for salinity, ) (5) in which, c is the salinity, which is measured and expressed by the chlorinity; and D x , D y are diffusion coefficients of the salinity, which are calibrated constant parameters for simplification. It should be noted that, in the proposed coupled model, salinity plays a role only in the evaluation of the sediment settling velocity. Since the settling velocity is an important factor for sediment transport, salinity must be considered by numerical models.

Governing Equations for Waves
The spectral action balance equation is used as the governing equation for waves, which is given by [20], σ (6) in which, N is the action density; σ is the relative radian frequency; θ is the wave propagation direction; C x , C y are the propagation velocities of wave energy in spatial x-, y-space, respectively; C σ , C θ are the propagation velocities in spectral σ-, θ-space, respectively; and S(σ, θ) is the source/sink term that represents all physical processes which generate, dissipate, or redistribute wave energy. The source/sink term represents the wave growth by the wind, nonlinear transfer of wave energy through three-wave and four-wave interactions and wave decay due to whitecapping, bottom friction and depth-induced wave breaking. Using the linear wave theory and the conversion of wave crests, the wave propagation velocities in spatial space within the Cartesian framework and spectral space can be obtained from the kinematics of a wave train [20]: in which, s is the space co-ordinate in the wave propagation direction of θ, and m is a co-ordinate perpendicular to s; → χ = (χ x , χ y ) is the wave speed; and U = (U x , U y ) is the tidal current velocity.

The Calculation of the Settling Velocity of Sediment
The Zhang Ruijin formula [21] is adopted for computing the settling velocity of sediment, which is given as, in which, ν denotes the kinematic viscosity coefficient (1.01 × 10 −6 m 2 /s); d denotes the sediment grain size; γ denotes the specific gravity of water; and γ s denotes the specific gravity of sediment.
With the increase of cohesive sediment concentration, the group settling velocity increases due to the flocculation. Therefore the Richardson & Zaki formula [22] is used to correct the settling velocity, in which, the superscript k denotes the k-th groups of sediment diameters; N s is the number of groups; ω s denotes the corrected settling velocity; ω 0 denotes the settling velocity computed by Equation (11), using the corrected kinematic viscosity coefficient ν m for muddy water; C µ is the correction coefficient of sediment concentration; S v is the concentration of muddy water in the form of the volume ratio, and S v = S/ρ s ; S vm is the limit concentration of muddy water in the form of the volume ratio; and d k , P k are the representative diameter and weight percentage of the k-th grain-size fractions, respectively. Furthermore, the salinity and sediment grain size have a direct impact on sediment flocculation and settling. Considering the flocculation effect of salinity, the settling velocity of sediment is corrected by the following equations [23], which are obtained from measured data in the Pearl River Estuary.
Study shows that there is the following relationship between the flocculation settling velocity, the sediment practical size and the deposition distance [23]: (16) in which, ω f denotes the settling velocity with flocculation; and R denotes the deposition distance. Based on Equation (16), if the chlorinity is lower than or equal to 20 psu, then [23] Otherwise, i.e., the chlorinity is higher than 20 psu, then [23] 0.89 (18) in which, ω f c is the corrected settling velocity by considering salinity, which is used in Equation (4) of the function of scour and deposit.

The Initiation Condition of Sediment Movement
The Zhang Ruijin formula [21] is adopted for computing the incipient velocity, 0.72 1 /2 (19) in which, the superscript k denotes the k-th groups of sediment diameters; and U c denotes the incipient velocity. The formula of wave height for the initiation of sediment movement developed by Dou Guoren [24] is adopted in this work, (20) in which, T denotes the wave period; L denotes the wave length; ∆ denotes the bed roughness height; δ = 2.3 × 10 −5 ; ε 0 denotes the bonding strength; β denotes the infilled aggregate coefficient in sediment; β w denotes the infilled aggregate coefficient in waves; and a, b denote the constant coefficients, which are determined by experimental data of sediment incipient motion under the wave.
Based on Equations (19) and (20), the initiation condition of sediment movement is given by, in which, U is the velocity of tidal current and H is the significant wave height.
So Equation (4) is equivalent to the following equations,

The Sediment-Carrying Capacity
Based on the analysis of the measured data in the Pearl River Delta, the sediment-carrying capacity of the tidal flow is proportional to the third power of the tidal flow velocity, and is inversely proportional to the settling velocity and water depth. In this work, the sediment-carrying capacity of the water flow proposed by Zhang Ruijin [21] is adopted, which is given as, in which, S * c denotes the sediment-carrying capacity of the water flow; and k, m denote the empirical parameters. The sediment-carrying capacity of the waves proposed by Dou Guoren [24] is adopted in this work, in which, S * w denotes the sediment-carrying capacity of the waves; α 0 , β 0 denote the empirical parameters; and H s denotes the significant wave height. The carrying capacity for a wave-current coexistent system is computed by, in which, S * is the total sediment-carrying capacity.

Adjustment of Bed-Load Gradation Composition
The concept of a mixing layer is introduced in this work to construct the relationship between the bed-load gradation composition and the erosion strength, which is given by [25], in which, H m denotes the thickness of the mixing layer; P b,k denotes the thickness ratio of the k-th grain-size fractions; B m denotes the width of the mixing layer; G k denotes the flux of sediment transport, which is the sum of the convection flux and diffusion flux of Equation (3); ε denotes the porosity of the bed material; and P b0,k denotes the thickness ratio of the bed material below the lower boundary of the mixing layer.
To adjust the bed-load gradation composition, the sea bed is divided into four layers. The top layer is the layer for sediment exchange between the bed material and the suspended sediment. The two middle layers are designed for the transition from the top layer to the bottom layer. Both the settlement of suspended sediment and the re-suspension of bed material will first induce the adjustment of bed-load gradation composition at the top layer. Following this, the thicknesses of the top layer and two middle layers are assumed to remain constant, so the bed-load gradation composition of the two middle layers and the bottom layer, as well as the thickness of the bottom layer, are adjusted by using the conservation law.
Assume that the initial thickness ratio of the bed material of the n-th layer is P erosion-sedimentation of the bed and the k-th groups of sediment are ∆Z b and ∆Z kb , then the thickness ratio of the bed material of the partial top layer, which is above the bottom surface of the former top layer, is updated by, (26) in which, P k1 is the updated thickness ratio of the bed material; H 1 is the thickness of the top layer. So for the case of sedimentation (∆Z b > 0), the bed-load gradation composition is adjusted by, For the case of erosion (∆Z b < 0), the bed-load gradation composition is adjusted by, in which, P kn is the updated thickness ratio of the bed material of the n-th layer; H 2 and H 3 are the thicknesses of the middle layers, and it should be noted that H 1 , H 2 , and H 3 are assumed to remain constant; and H 0 4 and H 4 are the initial and updated thickness of the bottom layer, respectively.

Numerical Method for Tidal Currents Model
The robust well-balanced finite volume method proposed by Song et al. [18,26,27] is adopted in this work for numerical resolving the governing equations of tidal currents. The method uses a high-resolution finite volume Godunov-type scheme with the HLLC (Harten, Lax, vanLeer, Contact) approximate Riemann solver [18]. The MUSCL (Monotone Upstream Scheme for Conservation Law) reconstruction technique is used to achieve second-order accuracy in space, and the second-order accuracy in time is obtained by using Hancock's predictor-corrector scheme [27]. Details about the numerical method can be found in the references [18,26,27]. Briefly, the process of time stepping is as follows: Step 1: Calculate the limited gradient by central values for each cell.
Step 2: Calculate the predictor results for each cell. The predictor results are given by, Water 2018, 10, 611 9 of 23 in which, η is the water level; ∆t is the computational time step controlled by the CFL (Courant-Friedrichs-Lewy) condition; ∂ x and ∂ y denote the limited gradient in the xand y-directions respectively; and the superscripts k and k + 1/2 denote the current time step and the predictor step, respectively.
Step 3: Calculate the reconstructed values for each edge.
Step 4: Calculate the numerical fluxes for each edge. The HLLC approximate Riemann solver is adopted to compute face fluxes, which is given by where U L and U R denote the reconstructed values; and s 1 , s 2 , s 3 are wave speeds. More details about the HLLC solver can be found in Song et al. [18].
Step 5: Calculate the source terms approximation for each cell.
Step 6: Calculate the solutions at the next time step for each cell. The solutions at the next time step k + 1 are given as, where m and L are the edge number and length respectively; F is the numerical flux; Ω is the area; and S is the source approximation. The superscripts k + 1/2 denote that the reconstructed variables are computed on the predictor results.

Numerical Method for the Coupled Tidal Current and Wave Models
To simulate the offshore wave propagation and coastal currents in the Pearl River Delta region, the two-dimensional tidal current model and the third generation wave model SWAN are adopted to develop a coupled hydrodynamic model. The dynamic coupling between the tidal current model and wave model is constructed by considering the wave radiation stress in the tidal current model, as well as the influence of tidal flow on waves in the wave model. The tidal current model and wave model are executed alternately. As the wave model uses a bigger fixed time step, while the time step of tidal current model is dynamically limited by the CFL condition, an asynchronous approach is adopted for the coupling of the tidal current model and wave model. So an adaptive time-stepping method is used for the tidal current model to trace the time line of the wave model, which is given as,

Numerical Method for the Advection-Diffusion Equation
Assume that the water flow fluxes are denoted as F WF = (F 1 , F 2 , F 3 ), in which F 1 denotes the flux of water volume, and F 2 and F 3 denote the momentum fluxes. The integrated solver proposed by Song et al. [28,29] is adopted for computing flow and mass transport fluxes simultaneously, which is given by, in which, F adv denotes the convection flux; and c L , c R denote the reconstructed concentration values from the left-cell and the right-cell, respectively. In Equation (40), F 1 > 0 means that water flows from the left-cell to the right-cell, and F 1 < 0 means that water flows from the right-cell to the left-cell.
By using the divergence theorem, the diffusion flux of mass transport is given by, in which, F dif is the diffusion flux; h Rec L , h Rec R are the reconstructed water depth from the left-cell and the right-cell, respectively; and n x , n y are the components of the unit outward normal vector.
Similar to the tidal current model, the MUSCL reconstruction technique is used to achieve second-order accuracy in space, and the second-order accuracy in time is obtained by using Hancock's predictor-corrector scheme. Details about the numerical method for resolving the advection-diffusion equation can be found in the references [28,29]. Briefly, the process of time stepping is as follows: Step 1: Calculate the limited gradient by central values for each cell.
Step 2: Calculate the predictor results for each cell. The predictor results are given by, in which, ∆t is the computational time step of tidal current model; ∂ x and ∂ y denote the limited gradient in the xand y-directions respectively; and the superscripts k and k + 1/2 denote the current time step and the predictor step, respectively.
Step 3: Calculate the reconstructed values for each edge.
Step 4: Calculate the numerical fluxes for each edge.
Step 5: Calculate the source terms approximation for each cell.
Step 6: Calculate the solutions at the next time step for each cell. The solutions at the next time step k + 1 are given as, where m and L are the edge number and length respectively; F is the numerical flux; Ω is the area; and S is the source approximation.
It should be noted that, the numerical method described above can be applied in the same way to solve Equation (3) of suspended load, as well as to solve Equation (5) of salinity.

The Coupling Approach
The conceptual framework for coupling the tidal currents, waves, and sediment models is shown in Figure 1. As presented in Figure 1, the tidal current model, wave model, salinity model, and the sediment model are tightly coupled with each other. and S is the source approximation.
It should be noted that, the numerical method described above can be applied in the same way to solve Equation (3) of suspended load, as well as to solve Equation (5) of salinity.

The Coupling Approach
The conceptual framework for coupling the tidal currents, waves, and sediment models is shown in Figure 1. As presented in Figure 1, the tidal current model, wave model, salinity model, and the sediment model are tightly coupled with each other.

The Study Area
To simulate the sandbar in Modaomen Estuary, the Lingding Sea is selected for this case study (see in Figure 2). The width of the study area is about 112 km, the length is 125 km, and the total waters area is about 7000 km 2 .
Water 2018, 10, x FOR PEER REVIEW 12 of 24

The Study Area
To simulate the sandbar in Modaomen Estuary, the Lingding Sea is selected for this case study (see in Figure 2). The width of the study area is about 112 km, the length is 125 km, and the total waters area is about 7000 km 2 .  Figure 3 shows the computational grids, and the bed elevation. The upstream boundaries include Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station (see in Figure 2). The western boundary reaches to Zhuhai airport, and the eastern boundary reaches to Hong Kong waters (see in Figure 2). P1-P2-P3-P4 are the tidal level boundaries (see in  Figure 3 shows the computational grids, and the bed elevation. The upstream boundaries include Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station (see in Figure 2). The western boundary reaches to Zhuhai airport, and the eastern boundary reaches to Hong Kong waters (see in Figure 2). P1-P2-P3-P4 are the tidal level boundaries (see in Figure 3). The number of computational grids is 51,800, of which the minimum edge length is 5 m, and the minimum grid area is 20 m 2 .   The correspondingly initial accumulative grading compositions, for example, are 32.0%, 56.5%, 76.5%, 85.8%, 92.7%, 97.1%, 99.6%, 100.0%, 100.0%, and 100.0%. It should be noted that there are different accumulative grading compositions, and the initial spatial distributions of the accumulative grading compositions for both the suspended load and the sea bed were set through practical experience in the Pearl River Delta. The sea bed is divided into four layers, and from the top layer to the bottom layer, the initial thicknesses of each layer are 0.05, 0.10, 0.5, and 2.0 m, respectively. Those thicknesses were also set through practical experience in the Pearl River Delta. The computational time step for the sediment model is equivalent to that of the tidal currents model. The diffusion coefficients of the suspended load vary from 5 to 20 m 2 /s. The empirical parameters k, m, which are used to compute the sediment-carrying capacity of the water flow, are set to 0.013, 0.24 when the water depth is smaller than 1.0 m; otherwise, if the water depth is smaller than 3.0 m the k, m are set to 0.08, 0.54, respectively; otherwise, the k, m are set to 0.10, 0.55, respectively. The empirical parameters α 0 , β 0 , which are used to compute the sediment-carrying capacity of the waves, are set to 0.023, 0.0004, respectively. To limit the sediment-carrying capacity, a maximum value of 4.5 kg/m 3 is used for the water flow, while a maximum value of 0.5 kg/m 3 is used for the waves. The settling probability, α 3 , is set to 1.0, and the saturation recovery coefficients of sediment concentration and sediment-carrying capacity, α 1 and α 2 , are set to 0.05 and 1.0 when S ≥ S * ; and are set to 0.5 and 1.0 when S < S * , respectively.

Model Parameters
The measured data, including tidal level, velocity, and sediment concentration, from 26 June 1998 20:00 to 27 June 1998 21:00 were used for parameter validation. Figure 4 presents the positions of the measuring points V1-V9. The measured data, including tidal level, velocity, and sediment concentration, from 26 June 1998 20:00 to 27 June 1998 21:00 were used for parameter validation. Figure 4 presents the positions of the measuring points V1-V9.   Figure 5b presents the comparison of the flow velocities and directions between the numerical results and the measured data; Figure 5c presents the comparison of the sediment concentration between the numerical results and the measured data. From Figure 5 it can be seen that the computed tidal levels, flow velocities and directions match well with those from the measured data, and the simulated sediment concentrations are consistent with that of the measured data, validating the reliability of the model parameters.  Figure 5a presents the comparison of the tidal levels between the numerical results and the measured data; Figure 5b presents the comparison of the flow velocities and directions between the numerical results and the measured data; Figure 5c presents the comparison of the sediment concentration between the numerical results and the measured data. From Figure 5 it can be seen that the computed tidal levels, flow velocities and directions match well with those from the measured data, and the simulated sediment concentrations are consistent with that of the measured data, validating the reliability of the model parameters. The measured data, including tidal level, velocity, and sediment concentration, from 26 June 1998 20:00 to 27 June 1998 21:00 were used for parameter validation. Figure 4 presents the positions of the measuring points V1-V9.  Figure 5a presents the comparison of the tidal levels between the numerical results and the measured data; Figure 5b presents the comparison of the flow velocities and directions between the numerical results and the measured data; Figure 5c presents the comparison of the sediment concentration between the numerical results and the measured data. From Figure 5 it can be seen that the computed tidal levels, flow velocities and directions match well with those from the measured data, and the simulated sediment concentrations are consistent with that of the measured data, validating the reliability of the model parameters.

Boundary Conditions
The combined conditions of the "99.7" flood season are used in this work for simulation. The "99.7" flood season denotes a flood that happened in July 1999. The time period was 15 July1999 23:00~23 July 1999 17:00, and the flood during this time was about 200 h. The "99.7" peak discharge of Sanshui station, which is an upstream hydrology control station in the Pearl River Delta, was 9200 m 3 /s and was close to the mean annual peak discharge of 9640 m 3 /s. Besides this, the "99.7" peak discharge of Makou station, which is another upstream hydrology control station in the Pearl River Delta, was 26,800 m 3 /s and was close to the mean annual peak discharge of 27,600 m 3 /s. Thus, it can be concluded that the combined conditions of the "99.7" flood season could be chosen as a typical dynamic condition for the flood season. Furthermore, the combined conditions of the "2001.2" dry season should be chosen as a typical dynamic condition for the dry season. In this work, we only present the results of the "99.7" flood season for simplification. Figure 6 presents the boundary conditions of the tidal currents model, including discharges at Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station, and tidal levels at P1-P4. From Figure 6, it can be seen that the computational period was about 200 h, and during this period the amount of water delivery of Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station are 368,937, 405,458, 206,921, 289,548, and 558,068 × 10 4 m 3 , respectively. It should be noted that the large difference in discharge was due to the tidal current; positive discharge indicates an ebb tide, while negative discharge indicates a flood tide. Figure 7 presents the boundary conditions of the sediment model, from which it can be seen that the maximum sediment concentration is 1.2 kg/m 3 , and the averages are 0.07, 0.34, 0.30, 0.29, and 0.37 kg/m 3 for Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station, respectively. A steady wind condition is used in this case, the wind speed is 3.6 m/s, and the wind direction is SE. This corresponds to the statistical value of the Dawanshan Station for the flood season. Since the wind wave is the dominant wave pattern in the Lingding Sea, the steady wind condition is appropriate for the wave and sediment simulation.

Boundary Conditions
The combined conditions of the "99.7" flood season are used in this work for simulation. The "99.7" flood season denotes a flood that happened in July 1999. The time period was 15 July1999 23:00~23 July 1999 17:00, and the flood during this time was about 200 h. The "99.7" peak discharge of Sanshui station, which is an upstream hydrology control station in the Pearl River Delta, was 9200 m 3 /s and was close to the mean annual peak discharge of 9640 m 3 /s. Besides this, the "99.7" peak discharge of Makou station, which is another upstream hydrology control station in the Pearl River Delta, was 26,800 m 3 /s and was close to the mean annual peak discharge of 27,600 m 3 /s. Thus, it can be concluded that the combined conditions of the "99.7" flood season could be chosen as a typical dynamic condition for the flood season. Furthermore, the combined conditions of the "2001.2" dry season should be chosen as a typical dynamic condition for the dry season. In this work, we only present the results of the "99.7" flood season for simplification. Figure 6 presents the boundary conditions of the tidal currents model, including discharges at Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station, and tidal levels at P1-P4. From Figure 6, it can be seen that the computational period was about 200 h, and during this period the amount of water delivery of Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station are 368,937, 405,458, 206,921, 289,548, and 558,068 × 10 4 m 3 , respectively. It should be noted that the large difference in discharge was due to the tidal current; positive discharge indicates an ebb tide, while negative discharge indicates a flood tide. Figure 7 presents the boundary conditions of the sediment model, from which it can be seen that the maximum sediment concentration is 1.2 kg/m 3 , and the averages are 0.07, 0.34, 0.30, 0.29, and 0.37 kg/m 3 for Dahu station, Nansha station, Fengmamiao station, Hengmen station, and Denglongshan station, respectively. A steady wind condition is used in this case, the wind speed is 3.6 m/s, and the wind direction is SE. This corresponds to the statistical value of the Dawanshan Station for the flood season. Since the wind wave is the dominant wave pattern in the Lingding Sea, the steady wind condition is appropriate for the wave and sediment simulation.   Figure 8 shows the simulated distribution of significant wave height and the wave period. From Figure 8 it can be seen that the significant wave height was about 0.5 m~1.0 m in Modaomen Estuary, while the wave period was about 4.5 s~5.5 s.   Figure 8 shows the simulated distribution of significant wave height and the wave period. From Figure 8 it can be seen that the significant wave height was about 0.5 m~1.0 m in Modaomen Estuary, while the wave period was about 4.5 s~5.5 s.  Figure 8 shows the simulated distribution of significant wave height and the wave period. From Figure 8 it can be seen that the significant wave height was about 0.5 m~1.0 m in Modaomen Estuary, while the wave period was about 4.5 s~5.5 s. Figure 9 shows the simulated distribution of sediment concentration at flood fast tide and at ebb fast tide. It should be noted that the moments of flood fast tide and ebb fast tide are confirmed with regard to the Modaomen Estuary, and there is a phase difference of tidal current between the Modaomen Estuary and other northern estuaries. From Figure 9 it can be seen that the sediment concentrations near river estuaries are much higher than those of other areas. Especially, near the Modaomen Estuary, the sediment concentration at ebb fast tide is much higher than that at flood fast tide. So it can be concluded that the tide is the dominant factor influencing the variation of suspended sediment concentration. Similar results were found in Deep Bay, Pearl River Estuary, China [30]. Besides, the diffusion phenomenon of suspended load is well reproduced in Figure 9.

Model Validation and Discussion
The high sediment concentration zone along the shoreline indicates the effect of "wave-lifting-sand, tidal-current-transporting-sand" in the Lingding Sea.  Figure 9 shows the simulated distribution of sediment concentration at flood fast tide and at ebb fast tide. It should be noted that the moments of flood fast tide and ebb fast tide are confirmed with regard to the Modaomen Estuary, and there is a phase difference of tidal current between the Modaomen Estuary and other northern estuaries. From Figure 9 it can be seen that the sediment concentrations near river estuaries are much higher than those of other areas. Especially, near the Modaomen Estuary, the sediment concentration at ebb fast tide is much higher than that at flood fast tide. So it can be concluded that the tide is the dominant factor influencing the variation of suspended sediment concentration. Similar results were found in Deep Bay, Pearl River Estuary, China [30]. Besides, the diffusion phenomenon of suspended load is well reproduced in Figure 9. The high sediment concentration zone along the shoreline indicates the effect of "wave-lifting-sand, tidalcurrent-transporting-sand" in the Lingding Sea.   Figure 9 shows the simulated distribution of sediment concentration at flood fast tide and at ebb fast tide. It should be noted that the moments of flood fast tide and ebb fast tide are confirmed with regard to the Modaomen Estuary, and there is a phase difference of tidal current between the Modaomen Estuary and other northern estuaries. From Figure 9 it can be seen that the sediment concentrations near river estuaries are much higher than those of other areas. Especially, near the Modaomen Estuary, the sediment concentration at ebb fast tide is much higher than that at flood fast tide. So it can be concluded that the tide is the dominant factor influencing the variation of suspended sediment concentration. Similar results were found in Deep Bay, Pearl River Estuary, China [30]. Besides, the diffusion phenomenon of suspended load is well reproduced in Figure 9. The high sediment concentration zone along the shoreline indicates the effect of "wave-lifting-sand, tidalcurrent-transporting-sand" in the Lingding Sea. From Figure 9 it can be concluded that the northwestern estuaries of the Lingding Sea, including Humen, Hengmen, Jiaomen, Hongqimen, and Modaomen, are the main inputs of the suspended sediment, since freshwater and sediments of the Pearl River basin are discharged by these estuaries into the Lingding Sea, and the Modaomen Estuary is the biggest gate with regard to the water delivery and sediment discharge. During the flood season, the suspended sediment concentration of the Pearl River is much higher than that of "wave-lifting-sand", so it can be concluded that these estuaries of the Pearl River provide most of the suspended sediment. The sediment concentration near those estuaries is obviously higher than that of other waters. Especially, the boundary (see in Figure 9) from Lingding shallows (LS) to Tonggu shallows (TS), by which the Lingdingyang Sea can be divided into two parts regarding the relatively higher/lower sediment concentration, is considerably well simulated. As a whole, the simulated result of sediment concentration distribution is in accordance with that of the actual situation [31], i.e., in the Lingdingyang Bay, under the influence of the strong tidal flow on the east-side, the sediment is transported southward along the west-side, and spread to the east-side. The sediment concentration distribution is higher on the westside and lower on the east-side, and higher on the north-side and lower on the south-side [31]. Moreover, the suspended sediment dynamics computed by the proposed model match well with other research results [32,33]. Figure 10a shows the simulated thickness of deposition and scour near the Modaomen Estuary; Figure 10b shows the annual average change of the measured terrain from 2005 to 2011. Figure 10c presents the 3D view of the measured terrain of the Modaomen Estuary in 2011. To compare the simulated thickness with the annual bed change, an amplification coefficient was used to compute the simulated thickness, i.e., the simulated thickness from the combined conditions of "99.7" was amplified by the coefficient. Since the net sediment transport during the "99.7" flood period was about 10% of the annual sediment transport through the Pearl River estuaries, the amplification coefficient was set to 10. To present quantitative comparison, we set seven blocks A-G in the (b) From Figure 9 it can be concluded that the northwestern estuaries of the Lingding Sea, including Humen, Hengmen, Jiaomen, Hongqimen, and Modaomen, are the main inputs of the suspended sediment, since freshwater and sediments of the Pearl River basin are discharged by these estuaries into the Lingding Sea, and the Modaomen Estuary is the biggest gate with regard to the water delivery and sediment discharge. During the flood season, the suspended sediment concentration of the Pearl River is much higher than that of "wave-lifting-sand", so it can be concluded that these estuaries of the Pearl River provide most of the suspended sediment. The sediment concentration near those estuaries is obviously higher than that of other waters. Especially, the boundary (see in Figure 9) from Lingding shallows (LS) to Tonggu shallows (TS), by which the Lingdingyang Sea can be divided into two parts regarding the relatively higher/lower sediment concentration, is considerably well simulated. As a whole, the simulated result of sediment concentration distribution is in accordance with that of the actual situation [31], i.e., in the Lingdingyang Bay, under the influence of the strong tidal flow on the east-side, the sediment is transported southward along the west-side, and spread to the east-side. The sediment concentration distribution is higher on the west-side and lower on the east-side, and higher on the north-side and lower on the south-side [31]. Moreover, the suspended sediment dynamics computed by the proposed model match well with other research results [32,33]. Figure 10a shows the simulated thickness of deposition and scour near the Modaomen Estuary; Figure 10b shows the annual average change of the measured terrain from 2005 to 2011. Figure 10c presents the 3D view of the measured terrain of the Modaomen Estuary in 2011. To compare the simulated thickness with the annual bed change, an amplification coefficient was used to compute the simulated thickness, i.e., the simulated thickness from the combined conditions of "99.7" was amplified by the coefficient. Since the net sediment transport during the "99.7" flood period was about 10% of the annual sediment transport through the Pearl River estuaries, the amplification coefficient was set to 10. To present quantitative comparison, we set seven blocks A-G in the Modaomen Estuary (see in Figure 10a,b). Table 1 presents the quantitative comparison of blocks A-G, from which it can be concluded that the simulated thickness of the deposition and scour in blocks A, B, C, D, G consist of the average change of the measured terrain from 2005 to 2011. However, due to human activity, namely sand extraction, block E was scoured, while the simulated results exhibited a deposition trend. Block F, which is the western branch (see in Figure 10c), was partially simulated by the proposed model. Besides, the eastern branch (see in Figure 10c) was not well simulated by the proposed model. The main reason for this is probably that the combined conditions of the "99.7" flood season are not adequately typical for the formation and evolution of the eastern branch. From Figure 10a it can be seen that there exists sediment deposition in the Modaomen Estuary, constructing a butterfly-shaped sandbar. Nevertheless, from Figure 10a,c it can be concluded that the butterfly shape of the sandbar was well simulated by the proposed model. Modaomen Estuary (see in Figure 10a,b). Table 1 presents the quantitative comparison of blocks A-G, from which it can be concluded that the simulated thickness of the deposition and scour in blocks A, B, C, D, G consist of the average change of the measured terrain from 2005 to 2011. However, due to human activity, namely sand extraction, block E was scoured, while the simulated results exhibited a deposition trend. Block F, which is the western branch (see in Figure 10c), was partially simulated by the proposed model. Besides, the eastern branch (see in Figure 10c) was not well simulated by the proposed model. The main reason for this is probably that the combined conditions of the "99.7" flood season are not adequately typical for the formation and evolution of the eastern branch. From Figure 10a it can be seen that there exists sediment deposition in the Modaomen Estuary, constructing a butterfly-shaped sandbar. Nevertheless, from Figure 10a,c it can be concluded that the butterfly shape of the sandbar was well simulated by the proposed model.  From Figure 10, it can be concluded that the formation and evolution process of the sandbar in the Modaomen Estuary is correctly simulated. The butterfly shape of sandbar in the Modaomen Estuary is considerably well simulated by the proposed model, which matches well with the measured terrain. The main reason for this would be that the combined conditions of the "99.7" flood season are typical for the formation and evolution of the sandbar.
As a whole, results show that the proposed model can accurately simulate the process of "wavelifting-sand, tidal-current-transporting-sand" in the Lingding Sea, and the sandbar evolution process in Modaomen Estuary is well simulated.

Conclusions
In this work, a two-dimensional coupled tidal currents-waves-salinity-sediment model is proposed to simulate the sandbar evolution process in the Modaomen Estuary on an unstructured grid. The tidal currents, waves, and sediment are governed by the two-dimensional shallow water equations, the action density balance equation, and the advection-diffusion equation, respectively. To simulate the coupled effect between tidal currents and waves, the radiation stress of waves is considered in the source terms of the shallow water equation, and the effect of current on waves is computed by dynamically taking the water levels and velocities from the tidal current model to the wave model. An integrated solver is proposed for computing flow and sediment transport fluxes From Figure 10, it can be concluded that the formation and evolution process of the sandbar in the Modaomen Estuary is correctly simulated. The butterfly shape of sandbar in the Modaomen Estuary is considerably well simulated by the proposed model, which matches well with the measured terrain. The main reason for this would be that the combined conditions of the "99.7" flood season are typical for the formation and evolution of the sandbar.
As a whole, results show that the proposed model can accurately simulate the process of "wave-lifting-sand, tidal-current-transporting-sand" in the Lingding Sea, and the sandbar evolution process in Modaomen Estuary is well simulated.

Conclusions
In this work, a two-dimensional coupled tidal currents-waves-salinity-sediment model is proposed to simulate the sandbar evolution process in the Modaomen Estuary on an unstructured grid. The tidal currents, waves, and sediment are governed by the two-dimensional shallow water equations, the action density balance equation, and the advection-diffusion equation, respectively. To simulate the coupled effect between tidal currents and waves, the radiation stress of waves is considered in the source terms of the shallow water equation, and the effect of current on waves is computed by dynamically taking the water levels and velocities from the tidal current model to the wave model. An integrated solver is proposed for computing flow and sediment transport fluxes simultaneously, and then a two-dimensional coupled flow-sediment transport model based on unstructured Godunov-type scheme is developed. The MUSCL-Hancock scheme combined with the slope limiter is adopted to achieve the high-accuracy and high-resolution property. The effect of waves on sediment is considered by the sediment-carrying capacity and the condition of sediment incipient motion.
For the case study, the simulated result of sediment concentration distribution is in accordance with that of the actual situation in the Lingding Sea. The boundary of sediment concentration from Lingding shallows to Tonggu shallows is considerably well simulated. Besides, the butterfly shape of sandbar in the Modaomen Estuary is considerably well simulated by the proposed model, which matches well with the measured terrain. The performance of the proposed model is validated through an application to the case study of sandbar simulation in the Modaomen Estuary, China. The comparison between the simulated result and the measured terrain demonstrates that the complex coupling effects among tidal currents, waves, salinity and sediments are well simulated by the proposed model. So, it can be concluded that the proposed coupled numerical model can be used in the estuary regulation project, and thus has a bright application prospect. Since the finite-volume scheme used in this work is explicit, the computational time steps of the tidal currents model, salinity model, and sediment model are limited by the CFL condition, and the computational efficiency would be lower than the traditional implicit schemes if we use serial computation. To improve the computational efficiency of the proposed model, GPU-based high-performance computing method [34] will be used in the near future. Furthermore, the proposed morphodynamic model, which combines the SWAN model and Godunov-type tidal current model which uses the finite volume scheme as the solver, can be used to simulate the morphodynamic processes under sharp hydrodynamic gradients. It is worth assessing the sensibility of the proposed morphodynamic model under extreme hydrodynamic forcings [35,36].