A Computationally Efﬁcient Shallow Water Model for Mixed Cohesive and Non-Cohesive Sediment Transport in the Yangtze Estuary

: In this paper, a computationally efﬁcient shallow water model is developed for sediment transport in the Yangtze Estuary by considering mixed cohesive and non-cohesive sediment transport. It is ﬁrstly shown that the model is capable of reproducing tidal-hydrodynamics in the estuarine region. Secondly, it is demonstrated that the observed temporal variation of suspended sediment concentration (SSC) for mixed cohesive and non-cohesive sediments can be well-captured by the model with calibrated parameters (i.e., critical shear stresses for erosion/deposition, erosion coefﬁcient). Numerical comparative studies indicate that: (1) consideration of multiple sediment fraction (both cohesive and non-cohesive sediments) is important for accurate modeling of SSC in the Yangtze Estuary; (2) the critical shear stress and the erosion coefﬁcient is shown to be site-dependent, for which intensive calibration may be required; and (3) the Deepwater Navigation Channel (DNC) project may lead to enhanced current velocity and thus reduced sediment deposition in the North Passage of the Yangtze Estuary. Finally, the implementation of the hybrid local time step/global maximum time step (LTS/GMaTS) (using LTS to update the hydro-sediment module but using GMaTS to update the morphodynamic module) can lead to a reduction of as high as 90% in the computational cost for the Yangtze Estuary. This advantage, along with its well-demonstrated quantitative accuracy, indicates that the present model should ﬁnd wide applications in estuarine regions. the on Comparisons of simulated and observed tidal currents and SSC the hydrodynamic and is the transport in the Yangtze Estuary. Sensitivity analysis of key factors (i.e., Manning roughness coefﬁcient, critical shear stresses for erosion and erosion coefﬁcient)


Introduction
The Yangtze Estuary has been intensively influenced by human activities. The suspended sediment discharge to the Yangtze Estuary has reduced from 4.2 × 10 9 t·yr −1 during the period 1953-2002 to less than 1.0 × 10 9 t·yr −1 in the past decade [1]. Not surprisingly, bed erosion and reduction in suspended sediment concentration (SSC) have been observed in various regions of the Yangtze Estuary [1][2][3][4]. However, it has been suggested that there might be a morphological lag response to riverine sediment supply changes of about 10-30 years in seaward regions [5,6]. In this regard, it is important to conduct highresolution 10-30 years prediction of the response of the hydro-sediment-morphodynamic system in the Yangtze Estuary. However, high-resolution numerical prediction of field scale numerical cases is computationally demanding. The idea of morphological accelerating factor (MF), which updates the bed level at each hydrodynamic time step by increasing sediment erosion and deposition fluxes (thus resultant bed level changes) using a constant MF, has found wide applications in long-term morphodynamic predictions [7][8][9][10][11][12][13][14][15][16][17]. One basic assumption of the MF approach is that the bed level changes in one tidal cycle are small so that "nothing irreversible happens within an ebb and flood phase, even when all changes are multiplied by the factor" [7,8]. Numerically, using MF = 10, Hu et al. [9] The governing equations include the mass and momentum conservation equations for the water-sediment mixture and the mass conservation equations for suspended sediment, the salinity and the bed materials [22][23][24]. They are as follows.   (4) where U = vector for the physical variables; F, G = the advective flux vectors in the two Cartesian directions; S is the source term vector; t = time; x, y = horizontal directions in the Cartesian system; h = water depth; u, v = depth-averaged velocities in x and y directions; g = 9.8 m·s −2 is the gravitational acceleration; c k is the depth-averaged volumetric sediment concentrations of the k-th sediment class with a mean diameter of d k , where the subscript k = 1, 2, 3, . . . N sps , N sps is the number of sediment size class; C sa = depth-averaged salinity concentration; z b = bed elevation; S f x and S f y are the friction slopes in the x and y directions, which are determined by Equation (5) following [23]; S bx = −∂z b /∂x and S by =−∂z b /∂y are the bed slopes in the x and y directions; f is the Coriolis force coefficient; E T = ∑ E k and D T =∑ D k are the total sediment erosion and deposition fluxes, respectively; E k and D k are the size-specific sediment erosion and deposition fluxes, which are determined by Equations (6) and (7). For cohesive sediment (d k ≤ 0.03 mm), the erosion and deposition fluxes are calculated applying the Partheniades-Krone formulations [25], whereas the erosion and deposition fluxes of non-cohesive sediment (d k > 0.03 mm) follow the approach of Hu et al. [23]; p = bed sediment porosity, which is set empirically as p = 0.42; f a,k and f s,k are the sediment fractions within the bed active layer and at the interface between the active layer and those below the active layer (see details in [26]), respectively; η = z b − δ is the bottom elevation (i.e., the interface) of the active layer and δ = thickness of the bed active layer.

∂U ∂t
E k = max(0, f a,k M(τ/τ e − 1)) d k ≤ 0.03 mm f a,k α k c e,k ω k (1 − α k c e,k ) m k d k > 0.03 mm (7) where n is the Manning roughness coefficient, which is determined in Section 3.2; ω k is the settling velocity of the k-th sediment class calculated by Zhang's formula [27] in Equation (8); K s K sa K dk ω k is the effective settling velocity of the flocs (formed by cohesive sediments with the diameter less than 0.03 mm); K s , K sa and K dk are the correction factors accounting for the effect of sediment concentration (K s : Equation (9)), the effects of salinity (K sa : Equation (10)) and the effect of sediment diameter (K dk : Equation (11)), respectively [28][29][30]; ω k (1 − α k c e,k ) m k and ω k (1 − α k c k ) m k are the effective settling velocities of non-cohesive sediments considering the influence of sediment concentration, m k = 4.45Re −0.1 p , Re p = ω k d k /υ is the Particle Reynolds number; υ = 1.2 × 10 −6 m 2 ·s −1 is the kinematic viscosity of water; α k is the ratio of near-bed to depth-averaged concentration, which is here set as α k = 1; c e,k is the sediment transport capacity of the k-th size class of sediment determined using the formula of Zhang et al. [31], see Equation (12); τ = ρu * 2 is the bed shear stress; u * is bed shear velocity with u 2 * = gh S 2 f x + S 2 f y ; τ d is the critical bed shear stress for deposition (N·m −2 ); τ e is the critical bed shear stress for erosion (N·m −2 ); M is the erosion coefficient (kg·m −2 ·s −1 ).
where ρ s = 2650 kg·m −3 is the density of sediment, ρ w = 1000 kg·m −3 is the density of fresh water; c T = ∑ c k is the total sediment concentration; c p is a threshold sediment concentration above which the settling velocity of the flocs decreases with sediment concentration, and c p = 1.5 kg·m −3 is used; c sap and c sa,min are threshold salinities: c sap = 28 ppt and c sa,min = 2.8 ppt [28,29]; the reference diameter d k has a range from 0.011 to 0.022 mm, d r = 0.0215 mm is used. The above governing equations and empirical relations differ from Hu et al. [23] in the following aspects. Firstly, the Coriolis force is considered because it plays a crucial role in the dynamics of fluid and the morphological evolution of large-scale water body. Secondly, salinity movement is considered because it may cause fine-grained suspended sediments to flocculate. Thirdly, both cohesive and non-cohesive sediment transport are considered. These improvements (e.g., consideration of the Coriolis force, salinity effects, cohesive sediment transport, etc.), as compared to the previous model by Hu et al. [23], enable the present model the potential to be applied in estuarine regions, where these effects are important. See Section 3 for numerical calibration.

Finite Volume Discretization Using the Hybrid LTS/GMaTS Method
Equations (13)- (15) are the discretized forms of the present governing equations using the finite volume method on unstructured triangular meshes. Figure 1 shows sketches for triangular meshes, including (a) cell (i) and its three surrounding cells, and (b) edge (j) and its two neighboring cells. Physical parameters shown in Figure 1 are introduced when they appear below. Bed elevations are firstly defined at nodes. The arithmetic average bed elevation of the three constituent nodes is used for bed elevation at a given cell, which is then updated during the simulation. Other physical variables are defined directly at the cell center. The inverse distance interpolation method is applied to obtain physical variables at nodes if necessary. As shown in Equations (13)- (15), the hybrid Local Time Step/Global Maximum Time Step (LTS/GMaTS) method is used for variable updating. Specifically, the hydro-sediment part is updated by the graded local time step: ∆t L−i = 2 m * i ∆t min , where the subscript L − i indicates local time step for the cell (i), m * i is the LTS level of the cell (i), ∆t min is the globally minimum time step. In contrast, the morphodynamic part is updated by the global maximum time step: is the total number of cells.
where A i is the area of the cell (i); E nij = (Fn x + Gn y ) ij is the numerical flux crossing the j-th edge of the cell (i) with j = 1, 2, 3; n ij = (n x , n y ) ij represents the normal unit outward direction of the j-th edge of the cell (i); ∆L ij is the length of the j-th edge of the cell (i); the superscripts * and ** represents two consecutive sub-time levels (the temporal interval is ∆t L−i ) between the two synchronized time levels t 0 and t 0 + ∆T; the interval between the two synchronized time levels t 0 and t 0 + ∆T is termed a full cycle; N p = ∆T/∆t min is the maximum number of sub-cycles in the full cycle; ∆t min = min i=1,Nc (∆t L−i ); S c = 1, 2, ∼, N p is used to indicate the sequence of sub-cycles. Figure 2 shows the numerical structure for the present model. Within a full cycle, the hydro-sediment-morphodynamic system will be updated from one synchronized time level to the next. To complete such an update, the morphodynamic part at all cells is updated only once (Equations (14) and (15)), whereas the times that the hydro-sediment part has to be updated at a specific cell (i) is equal to the ratio ∆T/∆t L−i (Equation (13)). If ∆t L−i = ∆t min , the times that the hydro-sediment part at cell (i) is updated are N p ; If ∆t L−i = 2 1 ∆t min , the times that the hydro-sediment part at cell (i) is updated are N p /2, indicating that hydro-sediment part in cell (i) will be updated every two sub-cycles. In a specific sub-cycle S c , the implementation of the hydro-sediment part is activated if this inequality m * i < l s (S c ) is satisfied, where l s is a function of the sequence S c (see Hu et al. [23] for the function). If the hydro-sediment part is to be updated in the S c sub-cycle, the source terms would also be estimated; otherwise, the source terms in the S c sub-cycle take zero-value. Specifically, the source terms for sediment erosion/deposition and the Coriolis force are evaluated directly using empirical relations with flow variables at the corresponding sub-cycle as input; the bed slope source term is evaluated using the slope flux method [32] with flow variables at the sub-cycle but bed elevations at the initial synchronized time level t 0 as input; the bed friction source term is estimated using the splitting point-implicit method [33] with flow variables at the corresponding sub-cycle as input. Estimation of numerical fluxes in a specific sub-cycle depends on whether MOD((S c − 1)/2 m f j ) = 0, where m f j is the graded LTS level for edge (j). Numerical fluxes at internal edges are estimated by the Harten-Lax-van Leer Contact (HLLC) approximate Riemann solver using physical variables of the cells at the left and right side of the edge (i.e., E nj = F HLLC (U jL , U jR ), where U jL , U jR are the two Riemann states of the physical variables at the left and right side of the edge (j)). Before calling the HLLC subroutine, the physical variables of the cells at the left and right side of the edge, as input, have to be modified to ensure non-negative water depth reconstruction (Audusse et al. [34]; He et al. [35]; Hu et al. [23]). Numerical fluxes at external edges (i.e., boundary edges) are estimated by the definition (i.e., E n = Fn x + Gn y ) using values of the physical variables (i.e., Equation (2)). For the subcritical inflowing boundary edge, the unit discharge and unit-width sediment transport rate (or the sediment concentration) must be given. The tangential velocity to the edge is set to zero. The water depth and the normal flow velocities are estimated by the method of characteristics. For subcritical outlet boundary edge, the water level must be given. The velocity tangential to the edge is also set to zero. All other physical variables are computed by the method of characteristics. For supercritical outlet boundary edge and a wall boundary edge, all physical variables at the edge are set equal to values at the neighboring cell. The estimation of the globally minimum time step ∆t min and the graded LTS levels m * i and m f j for cell (i) and edge (j) are as follows. Firstly, calculate the locally allowable maximum time step for each cell (Equation (16)). Secondly, the globally minimum time step of all cells is computed (Equation (17)). Thirdly, the potential graded LTS level of each cell is computed (Equation (18)). Fourthly, the actual graded level for cell edges and for cells are computed (Equations (19) and (20)). The locally allowable maximum time step of cell (i) is determined using the Courant-Friedrich-Lewy (CFL) condition.
where Cr is the Courant number; ∆t ami is the allowable maximum (subscript 'am') time step for cell (i); u ij and v ij are the depth-averaged flow velocities of the j-th edge of cell (i) (referred to as edge (ij) in the following section); h i is the water depth at cell (i); and R ij is the distance from the center of cell (i) to edge (ij). If the water depth of a specific cell is smaller than 10 −6 m, ∆t ami is set to the maximum value of those cells with water depth larger than 10 −6 m. Based on Equation (13), the globally minimum time step is computed as follows ∆t min = min(∆t ami ) i=1,2,∼Nc (17) then updated during the simulation. Other physical variables are defined directly at the cell center. The inverse distance interpolation method is applied to obtain physical variables at nodes if necessary. As shown in Equations (13)-(15), the hybrid Local Time Step/Global Maximum Time Step (LTS/GMaTS) method is used for variable updating. Specifically, the hydro-sediment part is updated by the graded local time step: x y ij n n = n represents the normal unit outward direction of the j-th edge of the cell (i); ij L Δ is the length of the j-th edge of the cell (i); the superscripts * and ** represents two consecutive sub-time levels (the temporal

Study Area and Numerical Setting
The computational domain is shown in Figure 3a, which covers a large part of the East China Sea, Yellow Sea, Bohai Sea and the Yangtze Estuary, the Hangzhou Bay as well as the Zhoushan area. Specifically, it ranges from 117.6524 E   [36]). Then, the potential graded level m i is calculated for each cell: where m user is a predefined value that represents the specified maximum graded level. Setting m user = 0 indicates that the graded level of all the cells is zero, which will make the present model equivalent to existing models that adopt the approach of the global minimum time step (GMiTS). The actual grade level of edge (j) is computed by: where m jL and m jR are the potential graded levels of two neighboring cells of edge (j). N f is the total number of edges. The potential graded time step level of each cell is finally computed as: where m i1 , m i2 , m i3 represents the potential graded level of the three neighboring cells of cell (i) (Figure 1).

Study Area and Numerical Setting
The computational domain is shown in Figure 3a, which covers a large part of the East China Sea, Yellow Sea, Bohai Sea and the Yangtze Estuary, the Hangzhou Bay as well as the Zhoushan area. Specifically, it ranges from 117.6524 • E to 127.3857 • E in longitude and from 25.4775 • N to 40.9230 • N in latitude. Embedded in Figure 3a is a set of meshes used in this paper, of which the total grids are 135,473 and the cell sizes vary from~19 m within the Deepwater Navigation Channel (DNC) to~38,934 m near the offshore boundary. The minimum cell size of about 19 m (see Figure 3c for more details) is in the DNC project region, which was completed in 2011 and created a 92 km long channel with a water depth of 12.5 m along the North Passage and South Channel of the Yangtze Estuary [37,38]. In this regard, this paper has also used another set of meshes without refining the meshes of the DNC region; see Figure 3d. The resultant cell sizes vary from about 205 m to 38,396 m, of which the total cells are 106,733. The initial topography is compiled from the following sources: for the whole South/North Branch, the South/North Channel and the South/North Passage of the Yangtze Estuary, the topography in February 2016 is adopted, which is measured using GPS RTK (Real-time kinematic) system; for the Hangzhou Bay, the adjacent coastal regions and the large part of the adjacent East China Sea, the bed topography is calculated from the electronic nautical chart in 2015; the topography for the rest of the area uses ETOP1 terrain data provided by NOAA (National Oceanic and Atmospheric Administration) (https://ngdc.noaa.gov/mgg/globalglobal.html (accessed on 5 July 2020)). Part of the initial topography is shown in Figure 3b. In Figure 3b, stations with available measured data are also indicated, including eight stations for SSC (e.g., NG3, CS9S, CS6S, CSWS, CS3S, NCH1, NCH4, NCH9; indicated as " "), and three stations for the tidal current (e.g., NGN4S, CS9S, NCH6; indicated as " Figure 2. Sketch of the numerical structure (revised from Hu and Li [36]).

Study Area and Numerical Setting
The computational domain is shown in Figure 3a Figure 3c for more details) is in the DNC project region, which was completed in 2011 and created a 92 km long channel with a water depth of 12.5 m along the North Passage and South Channel of the Yangtze Estuary [37,38]. In this regard, this paper has also used another set of meshes without refining the meshes of the DNC region; see Figure 3d. The resultant cell sizes vary from about 205 m to 38,396 m, of which the total cells are 106,733. The initial topography is compiled from the following sources: for the whole South/North Branch, the South/North Channel and the South/North Passage of the Yangtze Estuary, the topography in February 2016 is adopted, which is measured using GPS RTK (Real-time kinematic) system; for the Hangzhou Bay, the adjacent coastal regions and the large part of the adjacent East China Sea, the bed topography is calculated from the electronic nautical chart in 2015; the topography for the rest of the area uses ETOP1 terrain data provided by NOAA (National Oceanic and Atmospheric Administration) (https://ngdc.noaa.gov/mgg/globalglobal.html (accessed on 5 July 2020)). Part of the initial topography is shown in Figure 3b. In Figure 3b, stations with available measured data are also indicated, including eight stations for SSC (e.g., NG3, CS9S, CS6S, CSWS, CS3S, NCH1, NCH4, NCH9; indicated as " "), and three stations for the tidal current (e.g., NGN4S, CS9S, NCH6; indicated as " "). They were all measured by Acoustic Doppler Current Profiler (ADCP) during July 2016. The remaining 10 stations for the water level (e.g., SDK, BCZ, NCD, JGJ; YSG, DS, LH, SS, ZH, BL; indicated as " ") are measured by GPS RTK. The water level in SDK, BCZ, NCD and JGJ are measured from 10 July to 10 August 2016, whereas the water level in ZH, YSG, SS, LH, DS "). They were all measured by Acoustic Doppler Current Profiler (ADCP) during July 2016. The remaining 10 stations for the water level (e.g., SDK, BCZ, NCD, JGJ; YSG, DS, LH, SS, ZH, BL; indicated as "•") are measured by GPS RTK. The water level in SDK, BCZ, NCD and JGJ are measured from 10 July to 10 August 2016, whereas the water level in ZH, YSG, SS, LH, DS and BL are measured from 29 June to 9 July 2015. Multiple sediment fractions were considered, including two non-cohesive sediment fractions-i.e., the fine sand (0.1 mm) and the coarse silt (0.062 mm)-and two cohesive sediment fractions-i.e., the fine silt (0.03 mm) and and the clay (0.004 mm). Figure 4a presents the spatial distribution of surficial bed sediment composition (using d 50 as an example) in part of the Yangtze Estuary, which is obtained by interpolation using the available measured data (see Figure 4b). For other regions, the spatial distribution of surficial bed sediment composition refers to Yang et al. [39]. There are three open boundaries, including the boundary of the Yangtze River at the Sanjiangying, the boundary of Hangzhou bay at the Qiantang River Bridge and a seaward boundary. The open sea boundary was driven by nine tidal constituents (i.e., M2, S2, N2, K2, K1, O1, P1, Q1 and M4), which are obtained from TPXO 7.2; SSC at the seaward boundary was set to zero since the open boundaries are far away from the Yangtze Estuary and mostly deeper than 100 m; saline concentration at the seaward boundary is interpolated from the HYCOM Global Reanalysis salinity data (http://apdrc.soest.hawaii.edu/data/data.php (accessed on 6 July 2020)). For the landward boundary in the Yangtze River, the time-series of SSC and water discharge measured at Datong (e.g., Figure 5a,b) are used, whereas the sediment composition uses the averaged value over 2004-2010 (Figure 5c). At the Qiantang River bridge, constant water discharge of 1000 m 3 ·s −1 and sediment concentration of 0.07 kg·m −3 are prescribed. Saline concentrations at the two landward boundaries are zero. The spatial distributions of initial salinity in the whole region are interpolated from the HYCOM Global Reanalysis salinity data. The code of present model is written using Intel (R) FORTRAN. The model has realized parallel computing by using open multiprocessing (OPENMP), whereas the graphic processing unit (GPU)-acceleration and message passing interface (MPI) parallel computing is in the testing stage. Two parameters are introduced to quantify the performance of the model [40,41], including the RMSE (Equation (21)), and the SS (Equation (22)).    The root mean squared error: The skill score: where S l and O l are simulation results and observation results respectively; O is the mean values of the observation data; L is the number of observed data in situ. The value of SS represents the model performance, which is classified as follows: SS > 0.65 excellent; 0.5 < SS < 0.65 very good; 0.2 < SS < 0.5 good; SS < 0.2 poor [42]. The RED by Equation (23) is applied to quantify the relative difference between the hybrid LTS/GMaTS method and the GMiTS method.
where S LTS/GMaTS is the hydrodynamic variables simulated by the present model using the hybrid LTS/GMaTS method, and S GMiTS is the hydrodynamic variables calculated by the traditional model using the GMiTS method.

Model Performance for Tidal Hydrodynamics
The bed resistance, as represented by the Manning roughness coefficient, greatly differs for different modeling practices. Here, n = 0.01 + 0.01/h [43], n = 0.013 + 0.012/h [37] and n = 0.016 [10] were tested, where h is water depth. Tidal flows during July-August 2016 and June-July 2015 are simulated and compared against the available measured data. Table 1 summarizes the statistics of RMSE and SS for these simulations at different stations. From Table 1, when the relation n = 0.01 + 0.01/h is used, the maximum values for the RMSE of water level, tidal current velocity and direction are 0.363 m, 0.305 m·s −1 and 24°4′, respectively, which are consistently smaller than those for the relation n = 0.013 +

Model Performance for Tidal Hydrodynamics
The bed resistance, as represented by the Manning roughness coefficient, greatly differs for different modeling practices. Here, n = 0.01 + 0.01/h [43], n = 0.013 + 0.012/h [37] and n = 0.016 [10] were tested, where h is water depth. Tidal flows during July-August 2016 and June-July 2015 are simulated and compared against the available measured data. Table 1 summarizes the statistics of RMSE and SS for these simulations at different stations. From Table 1, when the relation n = 0.01 + 0.01/h is used, the maximum values for the RMSE of water level, tidal current velocity and direction are 0.363 m, 0.305 m·s −1 and 24 • 4 , respectively, which are consistently smaller than those for the relation n = 0.013 + 0.012/h and n = 0.016. Moreover, the values of SS for 0.01+0.01/h show that it is better than another two Manning roughness coefficients while simulating the hydrodynamics (i.e., SS of 0.812-0.992, 0.862-0.965 and 0.837-0.940 for water level, tidal current velocity and direction, respectively). Therefore, the Manning roughness coefficient n = 0.01 + 0.01/h is used in the following. The comparison between simulated and observed water level as shown in Figure 6, indicating that the simulated results are in good agreement with the measured data. Figure 7 presents the comparison between simulated and observed tidal velocity and direction, for which there are some deviations between the simulated tidal velocity and measured data, especially at NCH6. The reason for this discrepancy may be that the locations of these three stations are easily disturbed by the incident flow and reflected flow from the coasts and channel. The averaged values of SS for tidal velocity and direction are 0.901 and 0.889 (n = 0.01 + 0.01/h), respectively, which are greater than 0.65. Overall, the above validation generally shows the present model's good capability in the reproduction of tidal hydrodynamics.  6. Verification of water level (dots denote the observed data; line denotes the simulation data). Figure 6. Verification of water level (dots denote the observed data; line denotes the simulation data).
Using n = 0.01 + 0.01/h and the numerical setting as shown in Section 3.1 (using the meshes considering the DNC), the tidal characteristics over the whole East China Sea, Yellow Sea and the Bohai Sea are simulated. T_tide program is applied to tidal harmonic analysis. Specifically, the time interval is set to be 1 h and the starting time is 9:00 on 1 February 2016. Harmonic analysis is carried out for each grid node. Figure 8 presents the spatial distribution of co-amplitude and co-phase of tide constituents, including M2 (Figure 8a), S2 (Figure 8b), K1 (Figure 8c) and O1 (Figure 8d) over the East China Sea, Yellow Sea and the Bohai Sea. It can be seen that interactions between incident waves and reflected waves of the semi-diurnal tides M2 and S2 lead to four amphidromic points: two in the Bohai Sea and the other two in the Yellow Sea. The distribution of diurnal tides K1 and O1 is very similar in the Yellow Sea and Bohai Sea. They have one amphidromic point in the Bohai Strait and the middle of the Yellow Sea, respectively. The co-phase lines radiate from the amphidromic point to the surrounding areas, and the co-amplitude lines present concentric ring distribution centered on the amphidromic point. The above simulation results are consistent with observations for the major tide constituents of M2, S2, K1 and O1 [44]. Furthermore, the largest M2 tidal range shown in Figure 8a appears along the west coast of the Korean peninsula, especially in the area near Kyunggi Bay, where its maximum tidal amplitude exceeds 2.0 m. Along the coast of China, the tidal amplitude is relatively stronger in the Hangzhou Bay, and reaches 1.9 m at the top of Hangzhou Bay, which exhibits good agreement with the results described by Hu et al. [45], Ge et al. [46] and Yao et al. [47]. Using n = 0.01 + 0.01/h and the numerical setting as shown in Section 3.1 (using the meshes considering the DNC), the tidal characteristics over the whole East China Sea, Yellow Sea and the Bohai Sea are simulated. T_tide program is applied to tidal harmonic analysis. Specifically, the time interval is set to be 1 h and the starting time is 9:00 on 1 February 2016. Harmonic analysis is carried out for each grid node. Figure 8 presents the spatial distribution of co-amplitude and co-phase of tide constituents, including M2 ( It can be seen that interactions between incident waves and reflected waves of the semi-diurnal tides M2 and S2 lead to four amphidromic points: two in the Bohai Sea and the other two in the Yellow Sea. The distribution of diurnal tides K1 and O1 is very similar in the Yellow Sea and Bohai Sea. They have one amphidromic point in the Bohai Strait and the middle of the Yellow Sea, respectively. The co-phase lines radiate from the amphidromic point to the surrounding areas, and the co-amplitude lines present concentric ring distribution centered on the amphidromic point. The above simulation results are consistent with observations for the major tide constituents of M2, S2, K1 and O1 [44]. Furthermore, the largest M2 tidal range shown in Figure 8a appears along the west coast of the Korean peninsula, especially in the area near Kyunggi Bay, where its maximum tidal amplitude exceeds 2.0 m. Along the coast of China, the tidal amplitude is relatively stronger in the Hangzhou Bay, and reaches 1.9 m at the top of Hangzhou Bay, which exhibits good agreement with the results described by Hu et al. [45], Ge et al. [46] and Yao et al. [47].

Model Performance for SSC and Sensitivity Analysis
For SSC simulation that involves cohesive sediment, the following three parameters are important: the critical shear stresses for erosion (τ e ), the critical shear stress for deposition (τ d ) and the erosion coefficient (M). Previous numerical simulations of SSC in the Yangtze Estuary have used a very wide range of values for these parameters [9][10][11]37,38,45,48,49]: the critical shear stress for erosion can vary from 0.02 N·m −2 to 3.5 N·m −2 ; the critical shear stress for deposition is usually linked to the critical shear stress for erosion using empirical relations (e.g., τ d = 0.69τ e in Zhang and Xie [27]; τ d = 0.5τ e in Zhu et al., [49]; and τ d = 4 × τ e /9 in Cao and Wang [50]); the erosion coefficient can vary in the range of 10 −6 -10 −3 kg·m −2 ·s −1 . Here, three values were selected for τ e : 0.1 N·m −2 , 0.4 N·m −2 and 0.8 N·m −2 ; three values are selected for M: 10 × 10 −5 kg·m −2 ·s −1 , 3.0 × 10 −5 kg·m −2 ·s −1 and 5.0 × 10 −5 kg·m −2 ·s −1 ; and the critical shear stress for deposition is set as τ d = 4 × τ e /9 following Cao and Wang [50]. This results in a total of nine cases. Time variations of the SSC during July 2016 (including a spring tide and a neap tide) are simulated using these nine numerical cases. Water 2021, 13, x FOR PEER REVIEW 15 of 27 Figure 8. Distributions of the co-amplitude (blue dashed lines; unit: meter) and co-phase (red solid lines; unit: degree) of the M2, S2, K1 and O1 tide constituents around the East China Sea, Yellow Sea and the Bohai Sea. Figure 8. Distributions of the co-amplitude (blue dashed lines; unit: meter) and co-phase (red solid lines; unit: degree) of the M2, S2, K1 and O1 tide constituents around the East China Sea, Yellow Sea and the Bohai Sea. Figure 9 presents the comparison of the computed and measured SSC for these nine cases at eight stations. For convenience of comparison, the presentation of the simulation results is grouped depending on the M value. From Figure 9, the following are observed. Firstly, erosion and deposition parameters greatly affect the results. A smaller threshold shear stress for erosion and a larger erosion coefficient leads to a larger sediment erosion flux (see Equation (7)) and thus higher SSC. Secondly, it is noted that at different stations, there is no one case that can always get the best agreements between the simulation and the observation. At best, using τ e = 0.4 N·m −2 and M = 3.0 × 10 −5 kg·m −2 ·s −1 results in good agreements at five stations (NG3, CS9S, CS6S, NCH1, NCH9; see Figure 3b for the specific positions of these stations). However, the results of CS3S, CSWS and NCH4 stations during spring tide from this case deviates more from the measured data, as compared to another five stations. This indicates that the specification of τ e and M should be site-dependent. This is understandable because the values of τ e and M are both functions of the bed density, porosity, composition, consolidation and evolution of the sediment under the complex and mixed effects of the physical and biological interaction process. These characteristics vary significantly in space and in time. As a compromise, these parameters (i.e., τ e = 0.4 N·m −2 , τ d = 0.18 N·m −2 and erosion coefficient M = 3.0 × 10 −5 kg·m −2 ·s −1 ), which give the satisfactory agreements at most stations, are used in the following.
In order to investigate the effects of initial bed composition on SSC, we set two cases with different bed compositions, in which the first case considers multiple sediment fraction (including the content of clay: 0.004 mm, fine silt: 0.03 mm, coarse silt: 0.062 mm and sand: 0.1 mm) and another case only includes single cohesive sediment (i.e., 0.01 mm). Other numerical settings are presented in Section 3.1. We ran the model and compared the simulated results with the available measured data during July 2016 (presenting results during spring tide as an example). As shown in Figure 10, the SSC considering multiple sediment fraction exhibits good agreement with the observations. By comparison, the given single cohesive sediment in the model cannot capture the observed temporal variation of SSC, and it results in a smaller SSC than the measured one. This is because a large amount of cohesive sediment in the bed composition reduces erosion. Generally, the more accurately the model describes the sediment properties (e.g., particle size, composition and viscosity), the more realistic the sediment settling velocity and erosion process will be. Therefore, multiple sediment fractions are necessary to be considered in the morphological model of the Yangtze Estuary.
The effect of DNC on SSC. Because some measured stations for SSC as shown in Figure 3b are located in the DNC engineering. We use two sets of meshes with and without the DNC (for mesh details and numerical setting, see Section 3.1). We ran the model and compared the simulated results with the available measured data during July 2016 (presenting results during spring tide as example). The results as shown in Figure 11 indicates that the SSC at two sites far away from the North Passage (i.e., NCH1, NCH9; see Figure 3b for the specific positions of these stations) do not show significant improvement after considering the DNC. However, for those sites in the DNC (i.e., CS9S, CS6S, CS3S, CSWS), the SSC considering the DNC is in better agreement with the measured data, and larger than those without DNC, especially obvious during spring tide. The main reason is that the construction of the DNC increases current velocity and decreases sediment deposition in the North Passage. Moreover, for the two sites near the DNC (i.e., NG3, NCH1), the simulated SSC considering the DNC also has been greatly improved. Therefore, the influence of the DNC on sediment transport and topography evolution cannot be ignored in the morphological model of the Yangtze Estuary.

Computational Efficiency
In this paper, the hybrid LTS/GMaTS method was adopted for efficient variable updating, in which an important parameter was involved (i.e., m user is a user-defined value, which limits the grade exponent, see Equation (18)). Setting m user = 0 means the graded LTS levels of all cells are zero, which will make the model equivalent to a traditional model that uses GMiTS. We use different values of m user to simulate the tidal dynamics and sediment transport in July 2016. The total grids are 135,473 considering the DNC engineering and cell size vary from~19 m within the DNC to~38,934 m near the offshore boundary (see Figure 3). Table 2 summaries the relative computational cost and difference (see Equation (23)) between the present model and traditional model. Figure 12a shows the variation of calculation cost and relative difference with m user . It is obvious that the application of the hybrid LTS/GMaTS method leads to significant reductions in the computational cost. Generally, larger m user leads to a greater reduction. A reduction as high as 90% was obtained for m user = 7. Further increase in the grade level leads to only mild reduction in the computational cost. The relative difference between the models with hybrid LTS/GMaTS and GMiTS on hydrodynamics as shown in Figure 12b-d indicate that although the larger m user leads to larger relative difference, the maximum relative difference of water level, current velocity and direction for m user = 7 are 0.007 m, 3.0 × 10 −3 m·s −1 and 1 • 56 , respectively, which are significantly less than the RMSE in Table 1.

Conclusions
This paper improved the computationally efficient shallow water model (Hu et al., [23]) by considering (1) the mixed cohesive and non-cohesive sediment transport, and (2) the effects of salinity and sediment concentration, as well as sediment diameter, on the flocs settling velocity. Comparisons of simulated and observed tidal currents and SSC demonstrate that the model, given reasonable parameters, performs well in reproducing the distribution of hydrodynamic processes, and is capable of representing the sediment transport in the Yangtze Estuary. Sensitivity analysis of key factors (i.e., Manning roughness coefficient, critical shear stresses for erosion and deposition, erosion coefficient) shows the following: (1) consideration of multiple sediment fraction is important for accurate modeling of SSC in the Yangtze Estuary; (2) the Deepwater Navigation Channel (DNC) project may lead to enhanced current velocity and thus reduced sediment deposition in the North Passage; and (3) the critical shear stress and the erosion coefficient are shown to be site-dependent, for which intensive calibration may be required. This means that these parameters are spatially variable even within the studied area and their values would be different elsewhere. Therefore, the presented paper is, in the context of suspended sediment concentration, strictly a case study. Finally, the improved computational efficiency is demonstrated by comparing the computational cost of the present model

Conclusions
This paper improved the computationally efficient shallow water model (Hu et al. [23]) by considering (1) the mixed cohesive and non-cohesive sediment transport, and (2) the effects of salinity and sediment concentration, as well as sediment diameter, on the flocs settling velocity. Comparisons of simulated and observed tidal currents and SSC demonstrate that the model, given reasonable parameters, performs well in reproducing the distribution of hydrodynamic processes, and is capable of representing the sediment transport in the Yangtze Estuary. Sensitivity analysis of key factors (i.e., Manning roughness coefficient, critical shear stresses for erosion and deposition, erosion coefficient) shows the following: (1) consideration of multiple sediment fraction is important for accurate modeling of SSC in the Yangtze Estuary; (2) the Deepwater Navigation Channel (DNC) project may lead to enhanced current velocity and thus reduced sediment deposition in the North Passage; and (3) the critical shear stress and the erosion coefficient are shown to be site-dependent, for which intensive calibration may be required. This means that these parameters are spatially variable even within the studied area and their values would be different elsewhere. Therefore, the presented paper is, in the context of suspended sediment concentration, strictly a case study. Finally, the improved computational efficiency is demonstrated by comparing the computational cost of the present model against that of a traditional model using GMiTS. For the present simulated cases, the maximum reduction of the computational cost is approximately 90%, and the maximum relative differences are just 0.007 m (for water level), 3.0×10 −3 m·s −1 (for tidal current velocity) and 1 • 56 (for tidal current direction). This advantage, along with its well-demonstrated quantitative accuracy of its capability to deal with mixed cohesive and non-cohesive sediments, provide a basis for long-term morphodynamic modeling in estuarine regions.