*2.1. Governing Equations*

Depth-averaged 2D shallow water equations (SWEs), with Coriolis terms, are used as the governing equations for the HDM, which are given by

$$\frac{\partial \eta}{\partial t} + \frac{\partial (l\mathbf{u})}{\partial \mathbf{x}} + \frac{\partial (l\mathbf{w})}{\partial y} = 0 \tag{1}$$

$$\frac{\partial \mu}{\partial t} + \mu \frac{\partial \mu}{\partial \mathbf{x}} + v \frac{\partial \mu}{\partial y} = f \upsilon - g \frac{\partial \eta}{\partial \mathbf{x}} + \frac{\tau\_{\mathbf{x} \mathbf{x}}}{\rho h} - g \frac{n\_m^2 \mu}{h^{4/3}} \frac{\sqrt{\mathbf{u}^2 + \upsilon^2}}{h^{4/3}} + \nu\_l \left( \frac{\partial^2 \mu}{\partial \mathbf{x}^2} + \frac{\partial^2 \mu}{\partial \mathbf{y}^2} \right) \tag{2a}$$

$$\frac{\partial v}{\partial t} + u \frac{\partial v}{\partial x} + v \frac{\partial v}{\partial y} = -fu - g \frac{\partial \eta}{\partial y} + \frac{\tau\_{sy}}{\rho h} - g \frac{n\_m^2 v \sqrt{u^2 + v^2}}{h^{4/3}} + \upsilon\_l \left(\frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right) \tag{2b}$$

where *h*(*<sup>x</sup>*, *y*, *t*) is the water depth, (m); *<sup>u</sup>*(*<sup>x</sup>*, *y*, *t*) and *<sup>v</sup>*(*<sup>x</sup>*, *y*, *t*) are the components of depth-averaged velocity in the horizontally in the *x*- and *y*-directions, respectively, (m/s); *t* is the time, (s); *g* is the gravitational acceleration, (m/s2); η(*<sup>x</sup>*, *y*, *t*) is the water level measured from an undisturbed reference water surface, (m); υ*t* is the coefficient of the horizontal eddy viscosity, (m<sup>2</sup>/s); *f* is Coriolis factor; *nm* is Manning's roughness coefficient, (m−1/<sup>3</sup> s); ρ is the water density, (kg/m3); and τ*sx* and <sup>τ</sup>*sy* are the wind stress in the *x*- and *y*-directions, respectively, (N/m2).

The above equations construct a set of equations for *u, v* and η. Their forms are invariable in the rotating frame of unstructured grids. The wind stress is imposed as per [12,13]. For a given location (*<sup>x</sup>*, *y*) of the Yangtze Estuary, the Coriolis factor *f* is given by

$$f = 2\Omega \sin\left(\frac{\pi}{180}\phi + \frac{y - y\_c}{6357.0 \times 1000}\right) \tag{3}$$

where Ω (7.29 × 10−<sup>5</sup> rad/s) is the angular velocity of rotation of the Earth; φ (31.38724◦) is the latitude of the reference location (*xc*, *yc*) which is shown in Figure 2.

The annual bed-load quantity transported through the outlets of the Yangtze Estuary is about 500–1000 × 10<sup>4</sup> tons, accounting for 1–2% of the total sediment load [15]. The bed-load transport therefore contributes little to the horizontal circulations of global water–sediment fluxes in the Yangtze Estuary, and is not solved by the present model. The suspended sediment is regarded to be nonuniform and is described by a fraction method. The vertically averaged 2D advection–diffusion equation, with a source term describing sediment exchange between flow and riverbed, is used to describe the transport of nonuniform suspended load:

$$\frac{\partial(\text{hC}\_{k})}{\partial t} + \frac{\partial(\text{uhC}\_{k})}{\partial \mathbf{x}} + \frac{\partial(\text{vhC}\_{k})}{\partial y} = \frac{\nu\_{l}}{\sigma\_{c}} \left[ \frac{\partial^{2}(\text{hC}\_{k})}{\partial \mathbf{x}^{2}} + \frac{\partial^{2}(\text{hC}\_{k})}{\partial y^{2}} \right] + \text{aw}\_{\text{sk}}(\text{S}\_{\star k} - \text{C}\_{k}) \tag{4}$$

where *k* is the index of the sediment fraction, *k* = 1, 2, ... , *Ns* (*Ns* is the number of fractions); *Ck* and *S\*k* = sediment concentration and the sediment-carrying capacity of flows for the *k*th fraction of the nonuniform suspended load, respectively, kg/m3; *wsk* = settling velocity of sediment particles for the *k*th fraction of the suspended load, m/s; α = sediment recovery coefficient, which is set to 1.0 and 0.25, respectively, in case of erosion and deposition [16].

According to particle size and physical/chemical property, the nonuniform sediment is divided into four fractions. The size ranges of fractions 1–4 are, sequentially, 0–0.031, 0.031–0.125, 0.125–0.5, and >0.5 mm. In real applications, researchers often determined the settling velocity (*w*s) of the fine particles according to field data, experiments or their experience [11,12,17–20]. In the present model, the *w*s of fraction 1 is set according to field data in the Yangtze Estuary, while the primitive settling velocity is directly used for other fractions.

Zhang's formula [21], which is widely used in evaluating the sediment-carrying capacity of flows in real applications, is used in our model and given by

$$S\_{\ast k} = \mathcal{K} \Big[ \mathcal{U}^3 / (\mathcal{g}hw\_{\ast k}) \Big]^m \tag{5}$$

where *U* is a vertically averaged velocity (*U* = √*u*<sup>2</sup> + *v*2); *m* is an exponent and set to 0.92 in our model; *K* is sediment-carrying coefficient and determined by calibrations with field data. In the model, Zhang's formula [21], with the help of the method in [22], is used to determine the fractional sediment-carrying capacities of flows for the nonuniform sediment.

Corresponding to Equation (4), riverbed deformation induced by the transport of the *k*th fraction of the nonuniform suspended load is described by

$$
\rho' \frac{\partial z\_{bk}}{\partial t} = \alpha w\_{sk} (\mathbb{C}\_k - S\_{\ast k}) \tag{6}
$$

where *zbk* = riverbed deformation caused by the *k*th fraction sediment, m; ρ = dry density of bed materials, kg/m3. The gradation state of the bed materials is also updated using the method of [22].

The coefficient ofManning's roughness, *nm*, in the HDM and the coefficient of the sediment-carrying capacity, *K*, in the STM are determined by calibration tests with field data. Because the Yangtze Estuary is large and includes various regions with different characteristics of flows and sediment transports (e.g., river reach, tidal reach, coast sea area, and sea region), non-constant model parameters are used in different regions.
