A Three-Dimensional Flow and Sediment Transport Model for Free-Surface Open Channel Flows on Unstructured Flexible Meshes

Three-dimensional (3D) hydrostatic-pressure-assumption numerical models are widely used for environmental flows with free surfaces and phase interfaces. In this study, a new flow and sediment transport model is developed, aiming to be general and more flexible than existing models. A general set of governing equations are used for the flow and suspended sediment transport, an improved solution algorithm is proposed, and a new mesh type is developed based on the unstructured polygonal mesh in the horizontal plane and a terrain-following sigma mesh in the vertical direction. The new flow model is verified first with the experimental cases, to ensure the validity of flow and free surface predictions. The model is then validated with cases having the suspended sediment transport. In particular, turbidity current flows are simulated to examine how the model predicts the interface between the fluid and sediments. The predicted results agree well with the available experimental data for all test cases. The model is generally applicable to all open-channel flows, such as rivers and reservoirs, with both flow and suspended sediment transport issues.


Introduction
Numerical models are gaining popularity for solving a wide range of environmental fluid mechanics problems.When dealing with free-surface flows and sediment transport processes in open channel cases (e.g., rivers and reservoirs), one-dimensional (1D) and two-dimensional (2D) numerical models have been widely used.1D models are cross-section averaged, while 2D models are mostly depth averaged.Some of the widely used 1D flow and sediment models include HEC-RAS [1], MIKE11 [2], CCHE1D [3], and SRH-1D [4].2D flow and sediment transport models include the research models of Wu [5], Hung et al. [6], and Huang et al. [7], and widely used public-domain models, such as CCHE2D [8], TELEMAC-MASCARET [9], UnTRIM [10,11], Delft3D [12], and SRH-2D [13,14].1D and 2D models will remain useful, particularly for applications with a relatively large spatial extent or over a long period of time.
Three-dimensional (3D) numerical models are less used in river engineering but start to gain attention owing to the rapid advancement of computer hardware.There are many instances where 1D and 2D models are insufficient for environmental modeling.Flow examples include rivers with sharp bends, impacts of in-stream structures, and lake stratification; sediment transport examples are the local scour in rivers, sediment sluicing at reservoirs, and turbidity currents in lakes and estuary, ROMS-CSTMS is a 3D hydrostatic model developed and used in Europe, and it has the circulation and sediment transport modeling capability [27].The model used also a curvilinear orthogonal mesh in the horizontal and a stretched, terrain-following mesh in the vertical.The model has been applied to many cases; it was applied to the Waipaoa River continental shelf offshore of the North Island in New Zealand [28].Satisfactory results were reported in comparison with a 13-month field data set.
Most of the 3D hydrostatic models adopted a restrictive mesh type in the horizontal plane and the sigma or Z mesh in the vertical direction.The structured, curvilinear mesh is the most widely used horizontal mesh type in the existing models.These meshes are inflexible in representing the complex geometry of environmental flows.A structured mesh, for example, is difficult to generate as there are requirements on the number of mesh points and the restricted shape and quality of the mesh cells.The vertical sigma or Z mesh has its own issues.With the sigma mesh, the governing equations are transformed from the physical space to the computational space.The benefit is that the vertical coordinate ranges from −1 to 0, which does not change with the free surface and bed movement.A drawback, however, is that the vertical mesh points are inflexible, which may lead to inadequate resolution around a density interface [29].Significant errors were reported in the areas with steep bottom topography [30].For flows with density interfaces, a significant number of vertical points may be needed to reduce the numerical error.The Z mesh was developed to overcome the problem facing the sigma mesh; it did not make the coordinate transformation and kept the physical coordinates so that vertical points may be easily moved adaptively.Another benefit of the Z mesh is that different numbers of vertical points may be used in different zones.The horizontal mesh planes, however, needed to be orthogonal to the vertical coordinate; this restriction caused the mesh layer to not align on the density interface, leading to large numerical errors [31,32].Further, the Z mesh is not terrain-following, so the bed is represented only by the staircase approximation.Such a zig-zag representation was found to cause inaccuracy in the bed shear stress computation and near-bed horizontal advection [31,32].
In this study, a new mesh type is proposed-an unstructured, physical-coordinate sigma mesh.In the horizontal plane, the unstructured polygonal mesh is generated; in the vertical direction, the physical-coordinate based sigma mesh is used.The polygonal mesh is the most general mesh type; other types are merely special cases.For example, the hybrid mesh of mixed triangular and quadrilateral cells may be used by the proposed new mesh.The hybrid mesh was shown to be a flexible and accurate type with the 2D depth-averaged hydraulic models [14].With the vertical physical coordinate mesh, the mesh point distribution may be arbitrarily assigned.Points can freely follow interfaces and remove the restrictions of both the traditional sigma mesh and Z mesh.Key features of the proposed new mesh type, different from the traditional meshes, are listed below:

•
The horizontal mesh uses the polygonal cells;

•
The vertical mesh points may be distributed arbitrarily conforming to free surface, bed, or interface changes; and

•
The physical-Cartesian-coordinate based governing equations are solved directly without the coordinate transformation.

Flow Equations
The present model makes the following assumptions: (a) the vertical pressure distribution is hydrostatic; (b) the density variation impacts only the buoyancy force (Boussinesque assumption); and (c) the flow-sediment mixture is incompressible.We concern primarily with rivers and reservoirs so that special costal processes, such as ocean wave-generated processes and Coriolis force, are not included.The above assumptions lead to the following flow equations:

+
∂V ∂y ∂U ∂t ∂V ∂t In the above, t is time; (x, y, z) are the Cartesian coordinates (x and y are in the horizontal directions; z is the vertical direction opposite to the gravity); ρ is the mixture density; U, V and W are the mean velocity components along the Cartesian coordinates x, y, z, respectively; g is the acceleration due to gravity; ζ is the free surface elevation; T xx , T xy , T yx , and T yy are the horizontal turbulence stresses; ν V is the vertical turbulent viscosity.Note that the horizontal and vertical turbulences are treated separately due to the quite different turbulent characteristics caused by the different spatial scales of typical open channel flows.It is a practice widely adopted by almost all hydro-static assumption models [12].The horizontal turbulence stresses in Equations ( 2) and ( 3) are computed by: where ν is the kinematic fluid viscosity, ν H is the horizontal turbulent viscosity, and k is turbulent kinetic energy.
Various turbulence models may be used to compute the horizontal and vertical viscosities.In this study, the large eddy simulation (LES) model is used for the horizontal turbulence and k-ε model is used for the vertical turbulence.The horizontal LES is based on the Smagorinsky subgrid scale model, following the formulation of [33].The horizontal eddy viscosity is computed by: In the above, A 2D is the 2D horizontal cell area and α is a model constant (normally taken to be 0.11 but can be a user calibration parameter).
The k-ε model is used to compute the vertical viscosity; the formulation of [26] is adopted as: Fluids 2019, 4, 18 5 of 19 In the above, k is the turbulence kinetic energy and ε is the turbulence dissipation rate.The 3D transport equations for k and ε are expressed as: ∂ε ∂t where G is the turbulence generation due to vertical velocity gradient and B is turbulence generation due to buoyancy force.The generation term G is computed by: The buoyance production represents energy conversion between turbulent kinetic energy and potential energy in a stratified flow; it is computed by: In general, σ ρ = 1 is assumed for sediment induced stratification cases.Many turbulence model constants have been proposed and the following standard values [34] are used: C µ = 0.09 and C ε3 = 0 for stable stratification; C ε3 = 1 for unstable stratification.

Sediment Equations
A general set of sediment governing equations are adopted.The suspended sediments are divided into a number of size classes so that the non-uniform processes are taken into account.A non-equilibrium and fully unsteady sediment transport equation is utilized for each size class.The sediment concentration is tightly coupled to the flow in that the flow dictates the sediment concentration distribution, while the sediment alters the flow simultaneously through the baroclinic term in the momentum equations.For size k, the following sediment transport equation is solved: In the above, the subscript k designates the variable associated with the sediment size class k; C k is the volumetric concentration; ω k is the fall velocity; D Hk and D Vk are the horizontal and vertical diffusivities computed by: In the above, σ Ck is the Schmidt parameter that is assumed to be a constant from 0.5 to 1.0 [35].The fall velocity is computed by the method of van Rijn [36].
The solution of Equation ( 13) needs a special boundary condition on the stream bed.The net sediment flux on the bed reflects the net sediment exchange rate between sediments in water and those on the bed; it is not zero unless the flow reaches equilibrium.The net sediment rate out of the water column is computed by [35]: Fluids 2019, 4, 18 6 of 19 where D k = ω k C k and E k are the sediment deposition and entrainment rate, respectively.The entrainment rate is computed by: In the above, entrainment rate is proportional to the local equilibrium concentration near the bed.The equilibrium concentration is determined by an empirical equation assuming unlimited sediment supply from the bed.In the numerical implementation, (D k − E k ) is the net sediment rate out of the water column and is added to bed as the depositional sediment.
The equilibrium concentration equation is often computed at a reference height δ above the bed, but the interpretation of the bed location (z = 0) varies among researchers [37].Theoretically, z = 0 should be the edge of the bed-load layer, but practically, such an edge is hard to measure as there is no clear boundary between the bedload and the suspended load.In this study, the z = 0 bed is chosen to be located at about ( 0.15 ∼ 0.35) d 50 below the top of the roughness element, as recommended by [38], where d 50 is the medium size of the roughness elements on the bed.
Seven equilibrium equations were reviewed in [39]; later a more comprehensive list were provided by [40].In this study, the equilibrium concentration is determined by the formula of [41] as follows: C * = 0.331(θ − 0.045) 1.75   1 + 0.331 0.46 (θ − 0.045) 1.75   (17a) In the above, u τ is the bed shear velocity, d is the representative sediment diameter such as d 50 , and g is the acceleration due to gravity.

Numerical Method
The above governing equations are first discretized on a suiTable 3D mesh; the discretized equation set is then solved with a robust solution algorithm.The 3D mesh generation consists of two steps in this study: a 2D horizontal mesh generation and an automatic vertical mesh point creation.The 2D mesh may assume polygonal shapes and the 2D nodes are assigned with the bed elevation.The vertical nodal distribution is carried out automatically by the model given: (a) the water surface elevation predicted by the model, (b) the number of vertical nodes, and (c) the distribution criterion.Vertical nodal distribution is performed every time step so that the free surface and bed changes are followed.Each 3D mesh cell consists of a closed set of faces-a set of vertical faces coinciding with the 2D mesh and the top and bottom faces.
The discretization of the governing equations is performed using the finite-volume method.The partial differential equations (PDE) are integrated over 3D mesh cells.The Gaussian integral is used to transform the volume integral to flow fluxes on mesh faces.They are described next.

Flow Discretization
The continuity Equation ( 1) is discretized for a 3D cell k as: Fluids 2019, 4, 18 7 of 19 In the above, subscript k refers to a variable associated with the cell k; ∀ is cell volume; The 3D mesh nodes are allowed to move in the vertical direction to conform to the free surface and bed elevation changes.The moving mesh fluxes, F gTk and F gBk , are computed with the method of [42].With this approach, a new mesh is first formed and the new cell volume is then computed.The moving mesh fluxes, F gTk and F gBk , are finally computed to satisfy the volume conservation constraint, i.e., The depth-averaged continuity Equation ( 4) is used to compute the free surface elevation in the flow solver.Instead of discretizing (4) directly, like most models did, we derive the discretized depth-averaged equation by summing up Equation (18a) vertically.Doing so over all vertical 3D cells at a 2D cell k leads to the following: In the above, → F f is the volume flow flux through the entire vertical face from bed to free surface; is the 2D velocity vector comprised of only the two horizontal velocity components.The momentum Equations ( 2) and ( 3) are rewritten in the following form: In the above is the velocity vector in the horizontal plane; ∇ is the 3D Kroneker Delta operator, ρ re f is a constant reference density, e.g., the clear water or average mixture density.Discretization of Equation ( 20a) is also carried out with the finite volume method.In the following, the subscript k is dropped for clarity.
The discretization of the unsteady term is performed using the first-order Euler scheme.The convection term is discretized as a second-order scheme: Fluids 2019, 4, 18 8 of 19 In the above, subscripts f, T, and B refer to the values at the centroids of the vertical, top, and bottom cell faces; −→ Φ f is obtained with a second-order averaging equation involving known values at the neighboring cells as [43]: In the above, stands for summation over all edges of the cell face; → n f ; and subscripts 1 and 2 denote values of the left and right cells.From now on, the above averaging of a variable from cell center to cell face is denoted as → Φ with "< >" as the averaging operator.Discretization of the diffusion term is more involved and follows the method of [43].A detailed derivation is omitted and the final second-order accuracy discretized equation for an arbitrary polygon is expressed as: In the above, ∑ F is summation over all faces of the cell (vertical, top, and bottom faces).
The water elevation gradient term is the same for all vertical 3D cells that share the same horizontal mesh (i.e., the 2D mesh).It is discretized as a second-order accuracy scheme by: In the above, ∑

F2D
is summation over all edges of the 2D mesh cell; A H is horizontal area of the cell, ζ F is free surface elevation at the center of the 2D cell edge, d F is edge distance of the 2D cell, and n F is unit normal vector on an edge of the 2D cell.The final discretized horizontal momentum equations, say at a cell P, may be assembled as: where ∑ nd is summation over all neighboring cells which share the same faces with cell P.

Flow Solution Procedure
A special procedure is used to obtain flow flux through vertical faces of a 3D mesh cell [43].The horizontal velocity at a cell face is computed by: where Equation ( 26) is vertically summed to obtain the volumetric flow flux vector as follows: insertion of Equation ( 27) into Equation (19a) leads to the following water elevation equation: and Equation ( 28) is used to solve the water surface elevation.
The iterative solution procedure proceeds as follows.First, the two momentum equations in Equation ( 25) are solved.That is, given all variables at time level 0, the 3D horizontal velocities are obtained (starred superscript denotes provisional values at the new time n): the flow flux at vertical faces is updated using Equation ( 27) as: The corrector step is next performed so that the new values of f are obtained that satisfy the continuity Equation (28).Expressed in the incremental form, the following equations are obtained: In the above, f .Equation (31c) is the elevation correction equation that is solved with the ILU preconditioned CGS solver [14].With the new water elevation obtained at time level n, the new velocity and flow fluxes are updated and the vertical velocity at time level n is computed using Equation (18a).This completes one cycle of the iterative solution of the governing equations.A number of iterations are usually performed to achieve the convergence within a time step.Note that the sediment transport equation is also included in the above iterative solution process if sediment transport is to be simulated.

Sediment Discretization
Discretization of the suspended sediment Equation ( 13) is carried out similarly to the momentum equations.The concentration Equation ( 13) is rewritten in the following form: The discretization of the unsteady term adopts the first-order Euler scheme and the second-order convection term is discretized as: where C f , C T , and C B are the concentration at the centroids of vertical, top, and bottom cell faces.Discretization of the diffusion term is second-order and is expressed as: The definition of D n and D edge C are discussed before for the momentum equations; ∑ F is a sum over all faces of the cell (vertical, top, and bottom faces).
The final discretized concentration equation, say at a cell P, is rearranged in the following form: where ∑ nb is summation over all neighboring cells which share the same faces with cell P.

Model Verifications
Flow and suspended sediment transport cases are selected to test and verify the new numerical model, and they are reported next.In the discussion, the proposed new 3D model is named SRH-3D, consistent with the 2D model SRH-2D developed in [13,14].

Diversion Channel Flow
Channel bifurcation occurs frequently in natural rivers.Flow features are complex at the juncture.A laboratory case studied in [44] is selected to verify the new flow solver.The case was also simulated by other 2D and 3D models in the past [13].
The solution domain consists of a main channel, 6.0 m in length and 0.3 m in width, and a side diversion channel normal to the main channel, 3.0 m in length and 0.3 m in width.Two horizontal meshes are used: one is the same as the 2D mesh of [13], which used pure quadrilaterals, and the other adopts the triangle-quadrilateral hybrid mesh.A total of 16 vertical cells are selected leading to 76,800 cells for the first mesh and 107,280 cells for the mixed mesh.A close-up view of the mixed mesh is shown in Figure 1a.
The flow conditions are as follows.The main channel has a discharge of 0.00567 m 3 /s (applied at X = −3.0m) and water surface elevation of 0.0555 m at the exit (X = 3.0 m); the water surface elevation is 0.0465 m at the exit of the diversion channel (Y = 3.3 m).The simulation with the two meshes produces essentially the same results, so only results from the mixed mesh are reported herein.The predicted flow pattern is shown in Figure 1b.Key flow features are predicted well, such as flow acceleration near the junction and into the diversion channel and the flow recirculation in the diversion channel.Simulated results are compared with the measured data of [44]; the water surface elevation is in Figure 2 and the depth-averaged velocity is in Figures 3 and 4 in both the main and side channels.Good agreements are obtained between the model results and the data.The SRH-3D and SRH-2D results are similar, which points to the fact that 2D models are applicable for such flows, at least for the case simulated.
predicted flow pattern is shown in Figure 1b.Key flow features are predicted well, such as flow acceleration near the junction and into the diversion channel and the flow recirculation in the diversion channel.Simulated results are compared with the measured data of [44]; the water surface elevation is in Figure 2 and the depth-averaged velocity is in Figures 3 and 4 in both the main and side channels.Good agreements are obtained between the model results and the data.The SRH-3D and SRH-2D results are similar, which points to the fact that 2D models are applicable for such flows, at least for the case simulated.acceleration near the junction and into the diversion channel and the flow recirculation in the diversion channel.Simulated results are compared with the measured data of [44]; the water surface elevation is in Figure 2 and the depth-averaged velocity is in Figures 3 and 4 in both the main and side channels.Good agreements are obtained between the model results and the data.The SRH-3D and SRH-2D results are similar, which points to the fact that 2D models are applicable for such flows, at least for the case simulated.

Flow in a Sharply Curved Bend
An open channel flow with a sharp 180° rectangular bend was experimentally studied in [45].A bend with a mean radius-to-width ratio of 3.0 and less is considered strongly curved and exhibits highly 3D flow characteristics [46].The Rozovskii bend has a mean radius-to-width ratio of 1.0, so it belongs to the very sharp bend category.The same case was subject to a number of numerical studies, such as a 2D depth-averaged modeling in [46], a 2D modeling explicitly taking into account the secondary flow in [47], and a 3D NS solver modeling in [48].
The model domain and bed geometry are shown in Figure 5a.The approach and exit channels are straight, 1.6 m in length, and 0.8 m in width.The bend has a radius of 0.4 m for the inner wall.The numerical modeling parameters are as follows.The entire channel has a flat smooth bottom.The discharge at the entrance is 0.0123 m 3 /s, resulting in an average velocity of 0.265 m/s, Reynolds number of 15,600, and Froude number of 0.11.The water elevation at the model exit is extrapolated from the experimental data.The horizontal mesh adopts the mixed triangles (two zones) and quadrilaterals (Figure 5a).A total of 20 vertical mesh cells are chosen, leading to a total of 165,280 3D mesh cells.
Predicted water depth along the inner and outer sidewall is compared with the experimental data in Figure 5b.It is seen that the water surface elevation is higher at the outer bank than at the inner bank-the expected super-elevation effect.The model results match well along the outer wall and the elevation is under-predicted along the inner wall.The SRH-3D model indeed has a better

Flow in a Sharply Curved Bend
An open channel flow with a sharp 180 • rectangular bend was experimentally studied in [45].A bend with a mean radius-to-width ratio of 3.0 and less is considered strongly curved and exhibits highly 3D flow characteristics [46].The Rozovskii bend has a mean radius-to-width ratio of 1.0, so it belongs to the very sharp bend category.The same case was subject to a number of numerical studies, such as a 2D depth-averaged modeling in [46], a 2D modeling explicitly taking into account the secondary flow in [47], and a 3D NS solver modeling in [48].
The model domain and bed geometry are shown in Figure 5a.The approach and exit channels are straight, 1.6 m in length, and 0.8 m in width.The bend has a radius of 0.4 m for the inner wall.The numerical modeling parameters are as follows.The entire channel has a flat smooth bottom.The discharge at the entrance is 0.0123 m 3 /s, resulting in an average velocity of 0.265 m/s, Reynolds number of 15,600, and Froude number of 0.11.The water elevation at the model exit is extrapolated from the experimental data.The horizontal mesh adopts the mixed triangles (two zones) and quadrilaterals (Figure 5a).A total of 20 vertical mesh cells are chosen, leading to a total of 165,280 3D mesh cells.
Predicted water depth along the inner and outer sidewall is compared with the experimental data in Figure 5b.It is seen that the water surface elevation is higher at the outer bank than at the inner bank-the expected super-elevation effect.The model results match well along the outer wall and the elevation is under-predicted along the inner wall.The SRH-3D model indeed has a better prediction than the 2D model, indicating that vertical variation of the velocity and the secondary flow produced by the sharp bend have some effects on the flow.Predicted depth-averaged velocity is compared with the experimental data in Figure 6.The 3D model results agree with the measured data well.It is noted that the 2D model results agree also well with the data, contrary to the finding of [47].The improvements of the 3D model over the 2D are mainly in the shift of the maximum velocity from the inner to outer sidewalls when flow comes out of the bend.This rapid shift can be explained by the falling transverse water elevation slope and release of the additional momentum by the secondary flow when the radius of curvature at the bend Predicted depth-averaged velocity is compared with the experimental data in Figure 6.The 3D model results agree with the measured data well.It is noted that the 2D model results agree also well with the data, contrary to the finding of [47].The improvements of the 3D model over the 2D are mainly in the shift of the maximum velocity from the inner to outer sidewalls when flow comes out of the bend.This rapid shift can be explained by the falling transverse water elevation slope and release of the additional momentum by the secondary flow when the radius of curvature at the bend exit abruptly changes to infinity.Finally, the ability of SRH-3D model to predict secondary flows is examined, since primary reasons to use a 3D model are to predict vertical velocity distribution and secondary flows.There is no secondary flow data with the present case, so a separate simulation using the NS solver U 2 RANS [43] is carried out.Comparison of the predicted secondary flows at two transects are compared in Figure 7 between SRH-3D and U 2 RANS.Very similar secondary flows are predicted by the two models.This shows the hydrostatic models, such as SRH-3D, are adequate in predicting the secondary flows, even in sharp bends.Predicted depth-averaged velocity is compared with the experimental data in Figure 6.The 3D model results agree with the measured data well.It is noted that the 2D model results agree also well with the data, contrary to the finding of [47].The improvements of the 3D model over the 2D are mainly in the shift of the maximum velocity from the inner to outer sidewalls when flow comes out of the bend.This rapid shift can be explained by the falling transverse water elevation slope and release of the additional momentum by the secondary flow when the radius of curvature at the bend exit abruptly changes to infinity.Finally, the ability of SRH-3D model to predict secondary flows is examined, since primary reasons to use a 3D model are to predict vertical velocity distribution and secondary flows.There is no secondary flow data with the present case, so a separate simulation using the NS solver U 2 RANS [43] is carried out.Comparison of the predicted secondary flows at two transects are compared in Figure 7 between SRH-3D and U 2 RANS.Very similar secondary flows are predicted by the two models.This shows the hydrostatic models, such as SRH-3D, are adequate in predicting the secondary flows, even in sharp bends.

Suspended Sediment in a Channel
3D modeling of suspended sediment transport is not trivial; correct implementation of even the bed boundary condition is not well understood [37].The flume results in [49] are used to test the sediment module.The experiment was carried out in a flume with rough bottom created by placing a single layer of pebble-sized stones (median diameter: 7 mm).The bed slope was 0.07252%.The water elevation from the roughened bottom was 62 mm after the flow was fully developed.A log-fit

Suspended Sediment in a Channel
3D modeling of suspended sediment transport is not trivial; correct implementation of even the bed boundary condition is not well understood [37].The flume results in [49] are used to test the sediment module.The experiment was carried out in a flume with rough bottom created by placing a single layer of pebble-sized stones (median diameter: 7 mm).The bed slope was 0.07252%.The water elevation from the roughened bottom was 62 mm after the flow was fully developed.A log-fit of the measured velocity profile led to the shear velocity of u τ = 0.021 m.The mean velocity was 0.22 m/s.The mesh has 62 cells along the 12.4 m flow direction and 31 vertical cells distributed non-uniformly.The cell center of the first cell near the bed is δ 1 = 2.4415 mm, so that the point is within the log law (y + = δ 1 u τ ν = 53.4).Upstream discharge is 0.01364 m 3 /s, downstream depth is 62 mm, and bed roughness height is k S = 13 mm.
The computed velocity and vertical eddy viscosity distribution are compared with the experimental data in Figure 8. Good agreement is obtained.Vertical eddy viscosity is important in predicting the suspended sediment transport.For the present case, two mechanisms are important for the sediment motion: settling velocity and turbulent resuspension.The well-cited experimental data sets of [50,51] are used for comparison of the turbulent viscosity.
Two model runs are made to compute the suspended sediment corresponding to two sediment sizes.The same runs were reported by [37].Run 1 has a sediment size of 0.0236 mm, sediment fall velocity of 0.0516 mm/s, and Rouse number of 0.06.Run 2 has the sediment diameter of 0.079 mm, fall velocity of 5.79 mm/s, and the Rouse number of 0.672.The sediment sizes are chosen, such that the Rouse parameters are significantly different.The simulation starts from clear water and stops once the equilibrium sediment concentration profiles are obtained.Predicted concentration is compared with the theoretical Rouse profile in Figure 9.It is seen that the Rouse profile is well predicted by the model.Note that the Rouse profile is only approximate and was derived by assuming the parabolic distribution of the vertical turbulence viscosity.

Lock Exchange: Intrusive Turbidity Current into a Two-layer Fluid
Intrusion of a turbidity current into a two-layer fluid is simulated next in an attempt to verify the ability of the model to simulate the sediment transport and capture the fluid-sediment interface.The experimental cases of [52] are chosen and the setup is illustrated in Figure 10.The case was also studied by [53] with an NS solver.The flume had a glass tank measuring 197.1 cm long, 19.9 cm wide, and 48.5 cm tall.The lock-length (l) behind the gate was fixed at 18.6 cm and the total water depth (H) was 20 cm.

Lock Exchange: Intrusive Turbidity Current into a Two-layer Fluid
Intrusion of a turbidity current into a two-layer fluid is simulated next in an attempt to verify the ability of the model to simulate the sediment transport and capture the fluid-sediment interface.The experimental cases of [52] are chosen and the setup is illustrated in Figure 10.The case was also studied by [53] with an NS solver.The flume had a glass tank measuring 197.1 cm long, 19.9 cm wide, and 48.5 cm tall.The lock-length (l) behind the gate was fixed at 18.6 cm and the total water depth (H) was 20 cm.

Lock Exchange: Intrusive Turbidity Current into a Two-layer Fluid
Intrusion of a turbidity current into a two-layer fluid is simulated next in an attempt to verify the ability of the model to simulate the sediment transport and capture the fluid-sediment interface.The experimental cases of [52] are chosen and the setup is illustrated in Figure 10.The case was also studied by [53] with an NS solver.The flume had a glass tank measuring 197.1 cm long, 19.9 cm wide, and 48.5 cm tall.The lock-length (l) behind the gate was fixed at 18.6 cm and the total water depth (H) was 20 cm.Two cases are simulated.The first is the symmetric flow turbidity current (SFTC), in which the depth of the two layers in the ambient fluids is equal (h 0 = h 1 = 10 cm).The densities of the two layers are: ρ 0 = 1000 kg/m 3 and ρ 1 = 1020 kg/m 3 .SFTC has the symmetrical flow type according to [52], as the density of the lock fluid is equal to the depth-weighted average of the upper and lower layers.The second case is an asymmetric flow turbidity current (AFTC), in which the depth of the two layers in the ambient fluids are set as h 0 = 17.5 cm and h 1 = 2.5 cm, respectively.The densities of the two layers are the same as in the first case.
In the numerical modeling, no variations are expected in the lateral direction (y), so only 3 lateral cells are used with the symmetry boundary condition specified in the lateral back and front boundaries.In the two other directions, 533 cells are in the longitudinal direction (x) and 50 cells are in the vertical direction (z).Initially, fluid is stationary.For the SFTC case, the initial concentration is C = 0.006061 in the lock, C = 0.0 and 0.012122 in the top and bottom of the ambient, respectively.For the AFTC case, C = 0.009091 in the lock, and C = 0.0 and 0.012122 in the top and bottom of the ambient, respectively.
The SFTC case results are shown in Figure 11 to visualize the temporal evolution of the intrusive gravity current.They are compared with the images taken from the laboratory experiments of [52].After the lock gate was removed, the fluid contained behind the lock gate collapsed symmetrically and propagated along the interface.The head already started to form and was visible at 2 s.The initial collapse began with rapid acceleration.As it propagated to the right end of the wall, the gravity current brought strong mixing, resulting in mass entrainment and dilution.These processes are predicted by the numerical model well.The results, however, confirm the finding of [54] that the hydrostatic models, such as the current one, failed to predict the formation of the Kelvin-Helmholtz billows.The present model otherwise is capable of predicting the overall turbidity current features, as well as the propagation speed of the current. respectively.
The SFTC case results are shown in Figure 11 to visualize the temporal evolution of the intrusive gravity current.They are compared with the images taken from the laboratory experiments of [52].After the lock gate was removed, the fluid contained behind the lock gate collapsed symmetrically and propagated along the interface.The head already started to form and was visible at 2 s.The initial collapse began with rapid acceleration.As it propagated to the right end of the wall, the gravity current brought strong mixing, resulting in mass entrainment and dilution.These processes are predicted by the numerical model well.The results, however, confirm the finding of [54] that the hydrostatic models, such as the current one, failed to predict the formation of the Kelvin-Helmholtz billows.The present model otherwise is capable of predicting the overall turbidity current features, as well as the propagation speed of the current.
Simulated results of the AFTC case are shown in Figure 12; the temporal evolution of the intrusive gravity current is compared between the model and the experiment of [52].Similar to the findings of the SFTC case, the hydrostatic model predicts the initial turbidity current formation and travel well.The model, however, is incapable of predicting the Kelvin-Helmholtz waves as well as the wave reflection process after the front reaches the end wall.Simulated results of the AFTC case are shown in Figure 12; the temporal evolution of the intrusive gravity current is compared between the model and the experiment of [52].Similar to the findings of the SFTC case, the hydrostatic model predicts the initial turbidity current formation and travel well.The model, however, is incapable of predicting the Kelvin-Helmholtz waves as well as the wave reflection process after the front reaches the end wall.

Conclusions
A new 3D hydrostatic numerical model is developed that may be used to simulate flow and suspended sediment transport in rivers and reservoirs with free surface and fluid-sediment interface.The new model is formulated in a general way so that it is applicable to a wide range of environmental fluid flows.A new solution algorithm is developed that derives the water elevation

Conclusions
A new 3D hydrostatic numerical model is developed that may be used to simulate flow and suspended sediment transport in rivers and reservoirs with free surface and fluid-sediment interface.The new model is formulated in a general way so that it is applicable to a wide range of environmental fluid flows.A new solution algorithm is developed that derives the water elevation equation from the vertical summation of the discretized continuity equation.The resultant elevation correction equation is consistent with the rest of the discretized equations, leading to a robust and stable iterative solver.A major contribution of the present study is the adoption of a new mesh-the horizontal polygonal mesh coupled with the arbitrarily distributed vertical mesh.The 3D mesh is terrain-following and movable vertically to changes in water surface elevation and bed elevation.
The flow solver is verified with two 3D cases having experimental data.The model performs well in stability and accuracy.The sharply curved bend flow modeling shows that the 3D hydrostatic model works well for flows in meander bends, as the predicted secondary flows match that from the non-hydrostatic NS solvers.In comparison with the 2D model, the 3D model improves upon the flow through the sharp bend.A primary advantage of the 3D model over 2D ones is that the vertical distribution and secondary flows are predicted.
The suspended sediment transport module is fully coupled to the flow solver, in that the flow dictates the sediment movement and the changes of sediment concentration alter the flow velocity.The coupling is achieved implicitly within the same time step.The sediment module is first tested for its ability to predict the vertical distribution of the sediment concentration, and then used to predict the turbidity current formation and travel.Despite the inability of the model to predict the waves associated with the Kelvin-Helmholtz instability, the model is adequate in predicting the sediment collapse, current formation and movement speed.
Future works include: (a) demonstration and application of SRH-3D to practical rivers; (b) development of the adaptive mesh movement capability to track density interface; and (c) extension to bedload sediment transport.

and→Vn
Bk are the flow velocity vector at the centroids of the vertical, top, and bottom faces; Bk are the unit normal vectors of the vertical, top, and bottom faces; A f k , A Tk , and A Bk are the areas of the vertical, top, and bottom faces; summation ∑ f aceV is over vertical faces only; F gTk and F gBk are volume fluxes at the top and bottom faces due to vertical mesh movement.In our solution procedure, Equation (18a) is used to compute the vertical flow velocity once the horizontal velocities are computed.

→Φ
edge is horizontal velocity at the center of the edge; δ → r edge is the edge distance vector; → r 1 and → r 2 are left-to-right distance vectors between the neighboring cell center to the face centroid ( → n f is used to decide left or right);

Figure 1 .Figure 2 .
Figure 1.A portion of the mixed mesh at the juncture (a) and the predicted surface flow velocity for the flow diversion case (b).

Figure 1 .
Figure 1.A portion of the mixed mesh at the juncture (a) and the predicted surface flow velocity for the flow diversion case (b).

Figure 1 .
Figure 1.A portion of the mixed mesh at the juncture (a) and the predicted surface flow velocity for the flow diversion case (b).

Figure 2 .
Figure 2. Comparison of water elevation along both walls of the main channel (a) and the side channel (b).Figure 2. Comparison of water elevation along both walls of the main channel (a) and the side channel (b).

Figure 2 .
Figure 2. Comparison of water elevation along both walls of the main channel (a) and the side channel (b).Figure 2. Comparison of water elevation along both walls of the main channel (a) and the side channel (b).

Figure 3 .
Figure 3.Comparison of x-velocity profiles at selected x locations in the main channel (Solid Blue Line: SRH-3D; Red Dashed Line: SRH-2D; Symbol: Experiment).

Figure 3 .
Figure 3.Comparison of x-velocity profiles at selected x locations in the main channel (Solid Blue Line: SRH-3D; Red Dashed Line: SRH-2D; Symbol: Experiment).

Figure 3 .
Figure 3.Comparison of x-velocity profiles at selected x locations in the main channel (Solid Blue Line: SRH-3D; Red Dashed Line: SRH-2D; Symbol: Experiment).

Figure 4 .
Figure 4. Comparison of V velocity profiles at selected y locations in the side channel (Solid Blue Line: SRH-3D; Red Dashed Line: SRH-2D; Symbol: Experiment).

Figure 4 .
Figure 4. Comparison of V velocity profiles at selected y locations in the side channel (Solid Blue Line: SRH-3D; Red Dashed Line: SRH-2D; Symbol: Experiment).

Fluids 2018, 3 ,Figure 5 .
Figure 5. Model domain and the horizontal mixed mesh (a) and a comparison of the water depth along the inner and outer sidewall for the bend case (b).

Figure 5 .
Figure 5. Model domain and the horizontal mixed mesh (a) and a comparison of the water depth along the inner and outer sidewall for the bend case (b).

Figure 5 .
Figure 5. Model domain and the horizontal mixed mesh (a) and a comparison of the water depth along the inner and outer sidewall for the bend case (b).

Fluids 2018, 3 ,Figure 8 .
Figure 8.Comparison of predicted results with the measured data for the Fuhrman case.(a) velocity; (b) turbulent viscosity.

Figure 8 .Figure 8 .Figure 9 .
Figure 8.Comparison of predicted results with the measured data for the Fuhrman case.(a) velocity; (b) turbulent viscosity.

Figure 10 .
Figure 10.Sketch and parameters for the intrusive gravity current cases.

Figure 9 .
Figure 9.Comparison of suspended sediment concentration vertical distribution for the Fuhrman case (Black line: the model; Red line: Rouse equation).(a) Run 1; (b) Run 2.

Figure 8 .Figure 9 .
Figure 8.Comparison of predicted results with the measured data for the Fuhrman case.(a) velocity; (b) turbulent viscosity.

Figure 10 .
Figure 10.Sketch and parameters for the intrusive gravity current cases.Figure 10.Sketch and parameters for the intrusive gravity current cases.

Figure 10 .
Figure 10.Sketch and parameters for the intrusive gravity current cases.Figure 10.Sketch and parameters for the intrusive gravity current cases.