Next Article in Journal
Assessment of the Environmental Impact of Acid Mine Drainage on Surface Water, Stream Sediments, and Macrophytes Using a Battery of Chemical and Ecotoxicological Indicators
Previous Article in Journal
Research on the Coastal Marine Environment and Rural Sustainable Development Strategy of Island Countries—Taking the Penghu Islands as an Example
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
Ocean College, Zhejiang University, Zhoushan 316021, China
2
Ocean Research Center of Zhoushan, Zhejiang University, Zhoushan 316021, China
*
Author to whom correspondence should be addressed.
Water 2021, 13(10), 1435; https://doi.org/10.3390/w13101435
Submission received: 5 April 2021 / Revised: 13 May 2021 / Accepted: 14 May 2021 / Published: 20 May 2021
(This article belongs to the Section Oceans and Coastal Zones)

Abstract

:
In this paper, a computationally efficient shallow water model is developed for sediment transport in the Yangtze Estuary by considering mixed cohesive and non-cohesive sediment transport. It is firstly 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 coefficient). 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 coefficient 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 find wide applications in estuarine regions.

1. Introduction

The Yangtze Estuary has been intensively influenced by human activities. The suspended sediment discharge to the Yangtze Estuary has reduced from 4.2 × 109 t·yr−1 during the period 1953–2002 to less than 1.0 × 109 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 high-resolution 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] projected evolution of the Jiuduansha Shoal for 20 years. Using MF = 5, Kuang et al. [10] forecasted the 20-year evolution of Nanhui Tidal Flat after the closure of the Three Gorges Dam. Using MF ranging from 11 to 60, Luan et al. [11] conducted 20-year forecast modeling of the Yangtze Estuary with decreasing river inputs and relative sea-level rise. Using MF ranging from 100 to 400, Guo et al. [12,13] and Zhou et al. [17] investigated the role of river flow, tidal asymmetry and salinity in the long-term (600–4000 years) evolution trends of idealized estuaries. These studies greatly improved our understanding of the long-term evolution trends of the Yangtze Estuary. Nevertheless, it is still interesting to conduct intermediate-term (e.g., 10–30 years) simulations without resorting to the MF method. Such interests are motivated by both the great progress of efficient numerical schemes as well as the improvement of computing hardware. There has been great improvement in both the numerical accuracy and numerical efficiency. For example, for the regional oceanic modeling system (ROMS), Shchepetkin and McWilliams [18] proposed a new family of time-stepping algorithms that combine forward–backward feedback with the best known synchronous algorithms, allowing an increased time step due to the enhanced internal stability without sacrificing its accuracy. In addition, for finite-volume shallow water hydro-sediment-morphodynamic models, the numerical accuracy of such model builds on the use of Riemann solvers to capture shock waves and contact discontinuities [19,20,21]. The high computational cost of such models has been largely reduced by the recently proposed hybrid Local Time Step/Global Maximum Time Step (LTS/GMaTS) method [22,23]. Specifically, high computational efficiency is achieved by implementing the LTS to solve equations governing sediment-laden flows (i.e., the hydro-sediment part), and implementing the GMaTS to solve equations governing bed materials (i.e., the morphodynamic part). The simulation of a 2-year topographic evolution in the Taipingkou Channel of the Yangtze River shows that the calculation cost can be reduced by up to 90%. However, such models have not been applied to estuarine regions as they were developed for non-cohesive sediments. In this paper, the computationally efficient shallow water model [23] is improved by considering mixed cohesive and non-cohesive sediment transport, and the effects of salinity, sediment concentration and sediment diameter on the flocs settling velocity. The improved model was calibrated against field data (i.e., water level, tidal current velocity and SSC) in the Yangtze Estuary and Hangzhou Bay. Comparative numerical studies of key factors are conducted, including the bed resistance, erosion and deposition parameters and the initial bed sediment composition, as well as the Deepwater Navigation Channel (DNC). Finally, the computational efficiency of the improved model in the Yangtze Estuary is evaluated. This paper is structured as follows. In Section 2, the mathematical formulations are described. In Section 3, model setup and validations are introduced, Section 4 analyses the computational efficiency of the improved model, and Section 5 presents the conclusion.

2. Mathematical Formulations

2.1. Governing Equations and Empirical Relations

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.
U t + F x + G y = S
U = [ h h u h v h c k h c s a ] , F = [ h u h u 2 + g h 2 / 2 h u v h u c k h u c s a ] , G = [ h v h u v h v 2 + g h 2 / 2 h v c k h v c s a ] , S = [ ( E T D T ) / ( 1 p ) g h ( S b x S f x ) + f v h g h ( S b y S f y ) f u h E k D k 0 ]
z b t = E k D k 1 p
( δ f a , k ) t = E k D k 1 p f s , k η t
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 s p s , N s p s is the number of sediment size class; C s a = 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 b x =   z b / x and S b y = 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.
S f x = n 2 u u 2 + v 2 h 4 / 3 , S f y = n 2 v u 2 + v 2 h 4 / 3
D k = { max ( 0 , K s K s a K d k ω k c k ( 1 τ / τ d ) ) d k 0.03 mm α k c k ω k ( 1 α k c k ) m k d k > 0.03 mm
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
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 s a K d k ω k is the effective settling velocity of the flocs (formed by cohesive sediments with the diameter less than 0.03 mm); K s , K s a and K d k are the correction factors accounting for the effect of sediment concentration ( K s : Equation (9)), the effects of salinity ( K s a : Equation (10)) and the effect of sediment diameter ( K d k : 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.45 Re p 0.1 , Re p = ω k d k / υ is the Particle Reynolds number; υ = 1.2 × 10−6 m2·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 = g h S f x 2 + S f y 2 ; τ 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).
ω k = ( 13.95 υ d k ) 2 + 1.09 ρ s ρ w ρ w g d k 13.95 υ d k
K s = { 1 + 50 c T 1.3 × 2650 0 < c T c p ( 1 0.008 c T × 2650 1 0.008 c p ) 4 ( 1 + 50 c p 1.3 ) c T > c p
K s a = { ( c s a / c s a p ) 1 0.53 c s a , min < c s a c s a p c s a > c s a p
K d k = { ( d r / d k ) 1.8 1 d k d r d k > d r
c e , k = f a , k 1 20 ρ s [ ( u 2 + v 2 ) 3 / g h ω k ] 1.5 1 + [ ( u 2 + v 2 ) 3 / 45 g h ω k ] 1.15
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 s a p and c s a , min are threshold salinities: c s a p = 28 ppt and c s a , 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.

2.2. 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: Δ T = max i = 1 , N c ( Δ t L i ) or Δ T = 2 max i = 1 , N c ( m i * ) Δ t min , where N c is the total number of cells.
U i * * = U i * Δ t L i A i j = 1 3 E n i j * Δ L i j + Δ t L i S i
( z b ) i t 0 + Δ T = ( z b ) i t 0 + S c = 1 N p Δ t L i [ ( D T ) i S c ( E T ) i S c ] 1 p o
( δ f a k ) i t 0 + Δ T = ( δ f a k ) i t 0 + S c = 1 N p Δ t L i [ ( D k ) i S c ( E k ) i S c ] 1 p o f s k , i S c = 1 N p Δ t L i [ ( D T ) i S c ( E T ) i S c ] 1 p o
where A i is the area of the cell (i); E n i j = ( F n x + G n y ) i j is the numerical flux crossing the j-th edge of the cell (i) with j = 1 , 2 , 3 ; n i j = ( n x , n y ) i j represents the normal unit outward direction of the j-th edge of the cell (i); Δ L i j 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 , N c ( Δ 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 n j = F HLLC ( U j L , U j R ) , where U j L , U j R 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 = F n x + G n 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.
Δ t a m i = C r min j = 1 , 2 , 3 ( R i j u i j 2 + v i j 2 + g h i )   with   i = 1 , 2 , 3 , ~ , N c
where C r is the Courant number; Δ t a m i is the allowable maximum (subscript ‘am’) time step for cell (i); u i j and v i j 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 i j 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 a m i 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 a m i ) i = 1 , 2 , ~ N c
Then, the potential graded level m i is calculated for each cell:
m i = min ( int ( log ( Δ t a m i / Δ t min ) log ( 2 ) ) , m u s e r )
where m u s e r is a predefined value that represents the specified maximum graded level. Setting m u s e r = 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:
m f j = min ( m j L , m j R ) ,         j = 1 , 2 , 3 , ~ , N f
where m j L and m j R 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:
m i * = min ( m i , m i 1 , m i 2 , m i 3 )         i = 1 , 2 , 3 , ~ , N c
where m i 1 , m i 2 , m i 3 represents the potential graded level of the three neighboring cells of cell (i) (Figure 1).

3. Quantitative Accuracy

3.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 “”). 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 d50 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 m3·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:
R M S E = [ 1 L l = 1 L ( S l O l ) 2 ] 1 / 2
The skill score:
S S = 1 l = 1 L ( S l O l ) 2 l = 1 L ( O l O ¯ ) 2
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.
R E D ( S ) = 1 N c i = 1 N c [ ( S L T S / G M a T S ) i ( S G M i T S ) i ] 2
where S L T S / G M a T S is the hydrodynamic variables simulated by the present model using the hybrid LTS/GMaTS method, and S G M i T S is the hydrodynamic variables calculated by the traditional model using the GMiTS method.

3.2. 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.
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].

3.3. 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.
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.

4. 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 u s e r is a user-defined value, which limits the grade exponent, see Equation (18)). Setting m u s e r = 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 u s e r 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 u s e r . It is obvious that the application of the hybrid LTS/GMaTS method leads to significant reductions in the computational cost. Generally, larger m u s e r leads to a greater reduction. A reduction as high as 90% was obtained for m u s e r = 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 u s e r leads to larger relative difference, the maximum relative difference of water level, current velocity and direction for m u s e r = 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.

5. 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.

Author Contributions

Conceptualization, P.H.; Methodology, P.H., J.T., A.J. and Z.H.; Resources, P.H., W.L. and Z.H.; Investigation, W.L.; Software, J.T., A.J. and W.L.; Formal analysis, J.T.; Validation, J.T. and A.J.; Writing—original draft, J.T.; Writing—review & editing, P.H. and Z.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Key Research and Development Project (No. 2017YFC0405400), the National Natural Science Foundation of China (No. 11772300; 11872332), the Zhejiang Natural Science Foundation (No. LR19E090002), and the HPC Center OF ZJU (ZHOUSHAN CAMPUS).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Yang, S.; Milliman, J.D.; Li, P.; Xu, K. 50,000 dams later: Erosion of the Yangtze River and its delta. Glob. Planet. Chang. 2011, 75, 14–20. [Google Scholar] [CrossRef]
  2. Li, P.; Yang, S.; Milliman, J.D.; Xu, K.; Qin, W.; Wu, C.; Chen, Y.; Shi, B. Spatial, temporal, and human-induced variations in suspended sediment concentration in the surface waters of the Yangtze estuary and adjacent coastal areas. Estuaries Coasts. 2012, 35, 1316–1327. [Google Scholar] [CrossRef] [Green Version]
  3. Yang, S.; Milliman, J.D.; Xu, K.; Deng, B.; Zhang, X.; Luo, X. Downstream sedimentary and geomorphic impacts of the Three Gorges Dam on the Yangtze River. Earth Sci. Rev. 2014, 138, 469–486. [Google Scholar] [CrossRef]
  4. Luan, H.; Ding, P.; Wang, Z.; Ge, J.; Yang, S. Decadal morphological evolution of the Yangtze Estuary in response to river input changes and estuarine engineering projects. Geomorphology 2016, 265, 12–23. [Google Scholar] [CrossRef]
  5. Zhao, J.; Guo, L.; He, Q.; Wang, Z.; van Maren, D.S.; Wang, X. An analysis on half century morphological changes in the Changjiang Estuary: Spatial variability under natural processes and human intervention. J. Mar. Syst. 2018, 181, 25–36. [Google Scholar] [CrossRef] [Green Version]
  6. Zhu, C.; Guo, L.; van Maren, D.S.; Tian, B.; Wang, X.; He, Q.; Wang, Z. Decadal morphological evolution of the mouth zone of the Yangtze Estuary in response to human interventions. Earth Surf. Landf. Process. 2018, 44, 2319–2332. [Google Scholar] [CrossRef]
  7. Roelvink, J.A. Coastal morphodynamic evolution techniques. Coast. Eng. 2006, 53, 277–287. [Google Scholar] [CrossRef]
  8. Roelvink, J.A.; Reniers, A.J.H.M. A Guide to Coastal Morphology Modeling. Advances in Coastal and Ocean Engineering; World Scientific Publishing Company: Singapore, 2011; Volume 12. [Google Scholar]
  9. Hu, K.; Ding, P.; Wang, Z.; Yang, S. A 2D/3D hydrodynamic and sediment transport model for the Yangtze Estuary, China. J. Mar. Syst. 2009, 77, 114–136. [Google Scholar] [CrossRef]
  10. Kuang, C.; Liu, X.; Gu, J.; Guo, Y.; Huang, S.; Liu, S.; Yu, W.; Huang, J.; Sun, B. Numerical prediction of medium-term tidal flat evolution in the Yangtze Estuary: Impacts of the three gorges project. Cont. Shelf Res. 2013, 52, 12–26. [Google Scholar] [CrossRef]
  11. Luan, H.; Ding, P.; Wang, Z.; Ge, J. Process-based morphodynamic modeling of the Yangtze Estuary at a decadal timescale: Controls on estuarine evolution and future trends. Geomorphology 2017, 290, 347–364. [Google Scholar] [CrossRef]
  12. Guo, L.; van der Wegen, M.; Roelvink, D.J.A.; He, Q. The role of river discharge and tidal asymmetry on 1D estuarine morphodynamics. J. Geophys. Res. Earth Surf. 2014, 119, 2315–2334. [Google Scholar] [CrossRef] [Green Version]
  13. Guo, L.; van der Wegen, M.; Roelvink, D.J.A.; Wang, Z.; He, Q. Long-term, process-based morphodynamic modeling of a fluvio-deltaic system. Part I: The role of river discharge. Cont. Shelf Res. 2015, 109, 95–111. [Google Scholar] [CrossRef]
  14. Lesser, G.R.; Roelvink, J.A.; Van Kester, J.A.T.M.; Stelling, G.S. Development and validation of a three-dimensional morphological model. Coast. Eng. 2004, 51, 883–915. [Google Scholar] [CrossRef]
  15. Van der Wegen, M.; Roelvink, J.A. Long-term morphodynamic evolution of a tidal embayment using a two-dimensional, process-based model. J. Geophys. Res. Ocean. 2008, 113. [Google Scholar] [CrossRef] [Green Version]
  16. Van der Wegen, M.; Wang, Z.B.; Savenije, H.H.G.; Roelvink, J.A. Long-term morphodynamic evolution and energy dissipation in a coastal plain, tidal embayment. J. Geophys. Res. Earth Surf. 2008, 113. [Google Scholar] [CrossRef] [Green Version]
  17. Zhou, Z.; Chen, L.; Tao, J.; Gong, Z.; Guo, L.; van der Wegen, M.; Townend, I.; Zhang, C. The role of salinity in fluvio-deltaic morphodynamics: A long-term modelling study. Earth Surf. Process. Landf. 2020, 45, 590–604. [Google Scholar] [CrossRef]
  18. Shchepetkin, A.F.; McWilliams, J.C. The regional oceanic modeling system (ROMS): A split-explicit, free-surface, topography-following-coordinate oceanic model. Ocean. Model. 2005, 9, 347–404. [Google Scholar] [CrossRef]
  19. Toro, E.F. Shock-Capturing Methods for Free-Surface Shallow Flows; John Wiley & Sons Ltd.: Chichester, UK, 2001. [Google Scholar]
  20. Cao, Z.; Pender, G.; Wallis, S.; Carling, P. Computational dam-break hydraulics over erodible sediment bed. J. Hydraul. Eng. 2004, 130, 689–703. [Google Scholar] [CrossRef]
  21. Hu, P.; Cao, Z. Fully coupled mathematical modelling of turbidity currents. Adv. Water Resour. 2009, 32, 1–15. [Google Scholar] [CrossRef]
  22. Hu, P.; Lei, Y.; Han, J.; Cao, Z.; Liu, H.; Yue, Z.; He, Z. An improved local-time-step for 2D shallow water modeling based on unstructured grids. J. Hydraul. Eng. ASCE 2019, 145, 06019017. [Google Scholar] [CrossRef]
  23. Hu, P.; Lei, Y.; Han, J.; Cao, Z.; Liu, H.; He, Z. Computationally efficient hydro-morphodynamic modelling using a Hybrid local-time-step and the global maximum-time-step. Adv. Water Resour. 2019, 127, 26–38. [Google Scholar] [CrossRef]
  24. Hu, P.; Han, J.; Li, W.; Sun, Z.; He, Z. Numerical investigation of a sandbar formation and evolution in a tide-dominated estuary using a hydro-sediment-morphodynamic model. Coastal Eng. J. 2018, 60, 466–483. [Google Scholar] [CrossRef]
  25. Partheniades, E. Erosion and deposition of cohesive soils. J. Hydraulic. Div. ASCE 1965, 91, 105–139. [Google Scholar] [CrossRef]
  26. Hu, P.; Cao, Z.; Pender, G.; Liu, H. Numerical modelling of riverbed grain size stratigraphic evolution. Int. J. Sediment. Res. 2014, 29, 329–343. [Google Scholar] [CrossRef]
  27. Zhang, R.; Xie, J. Sedimentation Research in China: Systematic Selections; China and Power Press: Beijing, China, 1993. (In Chinese) [Google Scholar]
  28. Lin, Q.; Wu, W. A one-dimensional model of mixed cohesive and non-cohesive sediment transport in open channels. J. Hydraul. Res. 2013, 51, 506–517. [Google Scholar] [CrossRef]
  29. Wu, W.; Qi, H.; Wang, S. Depth-averaged 2-d calculation of tidal flow, salinity and cohesive sediment transport in estuaries. Int. Sediment. Res. 2004, 19, 1–10. [Google Scholar]
  30. Wu, W.; He, Z.; Lin, Q.; Sanchez, A.; Marsooli, R. Non-equilibrium sediment transport modeling-extensions and applications. Sediment. Transp. Monit. Modeling Manag. 2013, 21, 179–212. [Google Scholar]
  31. Zhang, R.; Xie, J.; Wang; River., M. Sediment Dynamics; China and Power Press: Beijing, China, 1989. (In Chinese) [Google Scholar]
  32. Hou, J.; Liang, Q.; Simons, F.; Hinkelmann, R. A 2d well-balanced shallow flow model for unstructured grids with novel slope source term treatment. Adv. Water Resour. 2013, 52, 107–131. [Google Scholar] [CrossRef]
  33. Liang, Q.; Marche, F. Numerical resolution of well-balanced shallow water equations with complex source terms. Adv. Water Resour. 2009, 32, 873–884. [Google Scholar] [CrossRef]
  34. Audusse, E.; Bouchut, F.; Bristeau, M.; Klein, R.; Perthame, B.T. A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows. SIAM J. Sci. Comput. 2004, 25, 2050–2065. [Google Scholar] [CrossRef]
  35. He, Z.; Wu, T.; Weng, H.; Hu, P.; Wu, G. Numerical investigation of the vegetation effects on dam-flows and bed morphological changes. Int. J. Sediment. Res. 2017, 32, 105–120. [Google Scholar] [CrossRef]
  36. Hu, P.; Li, Y. Numerical modeling of the propagation and morphological changes of turbidity currents using a cost-saving strategy of solution updating. Int. J. Sediment. Res. 2020, 35, 587–599. [Google Scholar] [CrossRef]
  37. Lu, C.; Jia, X.; Han, Y.; Bai, Y. Numerical simulation of sudden silting in the Yangtze Estuary deepwater channel by the wave of typhoon. Adv. Water Sci. 2018, 29, 696–705. (In Chinese) [Google Scholar]
  38. Song, D.; Wang, X. Suspended sediment transport in the Deepwater Navigation Channel, Yangtze River Estuary, China, in the dry season 2009: 2. Numerical simulations. J. Geophys. Res. Ocean. 2013, 118, 5568–5590. [Google Scholar] [CrossRef]
  39. Yang, H.; Yang, S.; Meng, Y.; Zhu, Q.; Wu, C.; Shi, B. Sediment distribution patterns in the Yangtze Estuary and comparison of particle size measurement methods. Shanghai Land Resour. 2017, 38, 75–79. (In Chinese) [Google Scholar]
  40. Willmott, C. On the validation of models. Phys. Geogr. 1981, 2, 184–194. [Google Scholar] [CrossRef]
  41. Warner, J.C.; Geyer, W.R.; Lerczak, J.A. Numerical modeling of an estuary: A comprehensive skill assessment. J. Geophys. Res. Ocean. 2005, 110. [Google Scholar] [CrossRef]
  42. Allen, J.I.; Somerfield, P.J.; Gilbert, F.J. Quantifying uncertainty in high-resolution coupled hydrodynamic-ecosystem models. J. Mar. Syst. 2007, 64, 3–14. [Google Scholar] [CrossRef]
  43. Dou, X.; Li, L.; Dou, G. Numerical model study on total sediment transport in the Yangtze Estuary. Hydro Sci. Eng. 1999, 2, 136–145. (In Chinese) [Google Scholar]
  44. Chen, G.; Niu, Y. Marine Atlas of Bohai Sea, Yellow Sea and East China Sea; China Ocean Press: Beijing, China, 1993. (In Chinese) [Google Scholar]
  45. Hu, K.; Ding, P.; Zhu, S.; Cao, Z. 2-D current field numerical simulation integrating Yangtze Estuary with Hangzhou Bay. China Ocean. Eng. 2000, 14, 89–102. [Google Scholar]
  46. Ge, J.; Ding, P.; Chen, C.; Hu, S.; Fu, G.; Wu, L. An integrated East China Sea-Changjiang Estuary model system with aim at resolving multi-scale regional-shelf-estuarine dynamics. Ocean. Dyn. 2013, 63, 881–900. [Google Scholar] [CrossRef]
  47. Yao, Z.; He, R.; Bao, X.; Wu, D.; Song, J. M_2 tidal dynamics in bohai and yellow seas: A hybrid data assimilative modeling study. Ocean. Dyn. 2012, 62, 753–769. [Google Scholar] [CrossRef]
  48. Ge, J.; Shen, F.; Guo, W.; Chen, C.; Ding, P. Estimation of critical shear stress for erosion in the Changjiang estuary: A synergy research of observation, GOCI sensing and modeling. J. Geophys. Res. Ocean. 2015, 120, 8439–8465. [Google Scholar] [CrossRef] [Green Version]
  49. Zhu, Q.; Prooijen, B.V.; Wang, Z.; Yang, S. Bed-level changes on intertidal wetland in response to waves and tides: A case study from the Yangtze River delta. Mar. Geol. 2017, 385, 160–172. [Google Scholar] [CrossRef]
  50. Cao, Z.; Wang, Y. Numerical Modeling of Hydrodynamics and Sediment; Tianjin University Press: Tianjin, China, 1994. (In Chinese) [Google Scholar]
Figure 1. Sketches of the unstructured triangular meshes: a cell surrounded by three cells.
Figure 1. Sketches of the unstructured triangular meshes: a cell surrounded by three cells.
Water 13 01435 g001
Figure 2. Sketch of the numerical structure (revised from Hu and Li [36]).
Figure 2. Sketch of the numerical structure (revised from Hu and Li [36]).
Water 13 01435 g002
Figure 3. (a) Model domain embedded with a set of meshes; (b) part of the bed topography. Also included in Figure 3b are stations for the available measured data; (c) local meshes considering the DNC engineering; (d) local meshes without considering the DNC engineering.
Figure 3. (a) Model domain embedded with a set of meshes; (b) part of the bed topography. Also included in Figure 3b are stations for the available measured data; (c) local meshes considering the DNC engineering; (d) local meshes without considering the DNC engineering.
Water 13 01435 g003
Figure 4. (a) Spatial distribution of bed sediment composition (using d50 as an example) in part of the Yangtze Estuary, which was obtained by interpolation using the available measured data in (b). (b) Indications for the available measured surficial bed sediment compositions in the Yangtze Estuary (data in Area I were measured in January 2016, whereas data in Area II were measured in July 2016).
Figure 4. (a) Spatial distribution of bed sediment composition (using d50 as an example) in part of the Yangtze Estuary, which was obtained by interpolation using the available measured data in (b). (b) Indications for the available measured surficial bed sediment compositions in the Yangtze Estuary (data in Area I were measured in January 2016, whereas data in Area II were measured in July 2016).
Water 13 01435 g004
Figure 5. Basic details of water discharge, sediment concentration, and composition. (a) Water discharge of the Yangtze Estuary during July 2016; (b) the sediment discharge at Datong station during July 2016; (c) the averaged riverine suspended sediment compositions at Datong in July from 2004 to 2010.
Figure 5. Basic details of water discharge, sediment concentration, and composition. (a) Water discharge of the Yangtze Estuary during July 2016; (b) the sediment discharge at Datong station during July 2016; (c) the averaged riverine suspended sediment compositions at Datong in July from 2004 to 2010.
Water 13 01435 g005aWater 13 01435 g005b
Figure 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).
Water 13 01435 g006
Figure 7. Verification of depth-averaged tidal velocity and direction (dots denote the observed data; line denotes the simulation data).
Figure 7. Verification of depth-averaged tidal velocity and direction (dots denote the observed data; line denotes the simulation data).
Water 13 01435 g007
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.
Water 13 01435 g008
Figure 9. Comparison of simulated and measured SSC under the different critical shear stresses and erosion coefficients at 8 stations shown in Figure 3b. The graphs on the left and right sides of break-line symbol show the variation of SSC during spring tide and neap tide, respectively.
Figure 9. Comparison of simulated and measured SSC under the different critical shear stresses and erosion coefficients at 8 stations shown in Figure 3b. The graphs on the left and right sides of break-line symbol show the variation of SSC during spring tide and neap tide, respectively.
Water 13 01435 g009aWater 13 01435 g009bWater 13 01435 g009cWater 13 01435 g009d
Figure 10. Comparison of simulated and measured SSC during spring tide under the different bed sediment composition, and subfigures (ah) represent the comparison results at 8 stations for SSC shown in Figure 3b, respectively.
Figure 10. Comparison of simulated and measured SSC during spring tide under the different bed sediment composition, and subfigures (ah) represent the comparison results at 8 stations for SSC shown in Figure 3b, respectively.
Water 13 01435 g010
Figure 11. Comparison of simulated and measured SSC during spring tide under the influence of Deepwater Navigation Channel (DNC), and subfigures (ah) represent the comparison results at 8 stations for SSC shown in Figure 3b, respectively.
Figure 11. Comparison of simulated and measured SSC during spring tide under the influence of Deepwater Navigation Channel (DNC), and subfigures (ah) represent the comparison results at 8 stations for SSC shown in Figure 3b, respectively.
Water 13 01435 g011
Figure 12. (a) Calculation efficiency varies with m u s e r ; (b) the relative difference of tide level varies with m u s e r ; (c) the relative difference of the tidal current velocity changes with m u s e r ; (d) the relative difference of the tidal current direction changes with m u s e r .
Figure 12. (a) Calculation efficiency varies with m u s e r ; (b) the relative difference of tide level varies with m u s e r ; (c) the relative difference of the tidal current velocity changes with m u s e r ; (d) the relative difference of the tidal current direction changes with m u s e r .
Water 13 01435 g012
Table 1. Model performance (RMSE and SS) in simulating hydrodynamics when using different roughness estimations.
Table 1. Model performance (RMSE and SS) in simulating hydrodynamics when using different roughness estimations.
StationsHydrodynamicsn = 0.01 + 0.01/hn = 0.013 + 0.012/hn = 0.016
RMSESSRMSESSRMSESS
SDKwater
level
0.363 m0.8120.214 m0.9370.297 m0.871
JGJ0.286 m0.9050.372 m0.7950.412 m0.782
NCD0.157 m0.9920.210 m0.9420.326 m0.865
BCZ0.218 m0.9520.179 m0.9760.268 m0.921
YSG0.229 m0.9660.279 m0.9120.358 m0.836
BL0.229 m0.9460.258 m0.9270.279 m0.891
DS0.259 m0.9220.288 m0.8830.331 m0.857
LH0.179 m0.9750.351 m0.8370.458 m0.753
SS0.145 m0.9810.196 m0.9550.287 m0.889
ZH0.218 m0.9440.188 m0.9680.275 m0.916
Mean value0.228 m0.9400.254 m0.9130.329 m 0.858
CS9STidal
current
velocity
0.287 m·s−10.8770.316 m·s−10.8470.375 m·s−10.815
NGN4S0.186 m·s−10.9650.198 m·s−10.9470.227 m·s−10.935
NCH60.305 m·s−10.8620.315 m·s−10.8510.371 m·s−10.821
Mean value0.259 m·s−10.9010.276 m·s−10.8820.324 m·s−10.857
CS9STidal
Current direction
13°37′0.94013°58′0.92514°26′0.907
NGN4S24°4′0.83724°15′0.82724°32′0.813
NCH615°20′0.88915°48′0.88216°5′0.876
Mean value17°20′0.88917°74′0.87818°36′0.875
Table 2. The calculation time and relative difference of different cases.
Table 2. The calculation time and relative difference of different cases.
MethodmuserMax Time Step (s)Computational Cost (h)Reduction in Computational CostRelative Difference
Water Level (m)Tidal Current Velocity
(m·s−1)
Tidal Current Direction (°)
GMiTS-0.42289.10----
LTS+ GMaTS00.42289.100---
10.85211.5826.81%0.00061.2 × 10−40°08′
21.70116.0359.87%0.00072.1 × 10−40°10′
33.4068.7576.22%0.00082.5 × 10−40°18′
46.8146.7283.84%0.00156.0 × 10−40°45′
513.6235.1087.86%0.00231.1 × 10−31°30′
627.2431.8888.97%0.00402.1 × 10−31°45′
754.4728.9889.98%0.00703.0 × 10−31°56′
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Hu, P.; Tao, J.; Ji, A.; Li, W.; He, Z. A Computationally Efficient Shallow Water Model for Mixed Cohesive and Non-Cohesive Sediment Transport in the Yangtze Estuary. Water 2021, 13, 1435. https://doi.org/10.3390/w13101435

AMA Style

Hu P, Tao J, Ji A, Li W, He Z. A Computationally Efficient Shallow Water Model for Mixed Cohesive and Non-Cohesive Sediment Transport in the Yangtze Estuary. Water. 2021; 13(10):1435. https://doi.org/10.3390/w13101435

Chicago/Turabian Style

Hu, Peng, Junyu Tao, Aofei Ji, Wei Li, and Zhiguo He. 2021. "A Computationally Efficient Shallow Water Model for Mixed Cohesive and Non-Cohesive Sediment Transport in the Yangtze Estuary" Water 13, no. 10: 1435. https://doi.org/10.3390/w13101435

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop