Simplified Spectral Model of 3D Meander Flow

Most 2D (two-dimensional) models either take vertical velocity profiles as uniform, or consider secondary flow in momentum equations with presupposed velocity profiles, which weakly reflect the spatio-temporal characteristics of meander flow. To tackle meander flow in a more accurate 3D (three-dimensional) way while avoiding low computational efficiency, a new 3D model based on spectral methods is established and verified in this paper. In the present model, the vertical water flow field is expanded into polynomials. Governing equations are transformed by the Galerkin method and then advection terms are tackled with a semi-Lagrangian method. The simulated flow structures of an open channel bend are then compared with experimental results. Although a zero-equation turbulence model is used in this new 3D model, it shows reasonable flow structures, and calculation efficiency is comparable to a depth-averaged 2D model.


Introduction
A meandering channel pattern is the most common fluvial channel type, and numerical simulation concerning meander flow is an efficient technology to analyze hydrodynamic processes with sediment transport and morphological evolution. The secondary flow (Prandtl's first kind of secondary flow, or helical secondary flow) at the bend cross section generates transversal momentum transport from inner bank to outer bank [1], as a result of the local imbalance of the centrifugal force induced by flow curvatures and the transversal pressure gradient. Accounting for the effects of secondary flow has been a subject of continuous research. For an idealized mild bend of uniform curvature, Rozovskii [2], Engelund [3], and Kikkawa et al. [4] have given analytical solutions to the secondary flow at the centerline. Spatial variation of the flow field governs sediment flux, and accordingly, bed deformations and flow separation [5]. In return, effects of topographical steering are also important for the flow structure, such as pool-riffle sequences [6] and bars on the inside of the bends [7]. Smith and McLean [8] have shown that both the advective acceleration associated with the curved flow path and bed topography play a decisive role. In many experimental channels, flow field cannot adapt instantaneously to changes in curvature, and secondary flow lags substantially [1,9]. Rozovskii [2], Odgaard [10], Ikeda and Nishimura [11], and Johannesson and Parker [1] established semi-heuristic relaxations for the phase lag between the vertical flow structure and curvature changes at the centerline of bends, which Ikeda et al. [12] later extended throughout the flow width.
On the other hand, secondary flow shows different characteristics according to channel curvatures. Many studies have focused on the secondary flow at mild and moderate bends where secondary flow monotonically increases with channel curvature [1,3,11,13]. However, these so-called linear models neglect the feedback of the secondary flow to the primary flow, which leads to an over-prediction of secondary flow for sharp bend cases. For laboratory experiments, Blanckaert [14] reported that there are capacity constraints for secondary flow development, and the saturation of secondary flow further contributes to energy loss and turbulence. There is a smaller counter-rotating circulation near the The present paper focuses on bend flow simulation in a simple 3D manner. To avoid vertical discretization and calculations in the vertical direction, this paper proposes a 3D model based on a horizontal 2D mesh combined with a vertical spectral method. The vertical spectral method expands velocity profiles into polynomials. The velocity profiles can be calculated dynamically in response to local boundary conditions, without numerical discretization. Although a zero-equation turbulence closure model is used in this new 3D model, it shows reliable flow structures compared with experimental data. Balanced on computational efficiency and precision, this 3D model is practically available.

Governing Equations
A set of 3D shallow water equations based on curvilinear coordinates is used for the water velocity calculation. The curvilinear coordinates include two horizontal axes and one vertical axis, which make the structure of the equations simple and easy to calculate. For convenience, equations in Cartesian coordinates are illustrated in the model, as in Equations (1) ∂u ∂t ∂v ∂t where t is time; H is water surface elevation; h is water depth; g is acceleration due to gravity; x, y, and z are horizontal and vertical Cartesian coordinates; u, v, and w are velocity components in the curvilinear coordinate system; ρ is fluid density; ν t = αhu × ζ(1−ζ) is the eddy viscosity coefficient with constant α = 0.2; u * is shear velocity ν t is depth-averaged eddy viscosity (αhu * ); and ζ = (z−z b )/h is relative elevation from 0 to 1.
The assumption of a hydrostatic pressure distribution for meander flow simulation is reasonable here. The hydrostatic assumption has been employed by many researchers to produce acceptable results [12,19]. The vertical acceleration and vertical velocity are significant and cannot be ignored near the inner and outer banks [32]. However, the hydrostatic pressure assumption causes large errors near lateral walls for simulation, which requires further discussion.

Spectral Method
The spectral method can be viewed as a special case of the weighted residual methods or the Galerkin method, where the trial and basis functions are the same orthogonal polynomials. We assume that the vertical profiles of velocity components are well approximated by a finite sum of basis functions. The basis functions can be polynomial or trigonometric functions of relative elevation ζ. The Legendre polynomials of degree from zero to any positive integer N for horizontal velocity components, and of degree N + 2 for vertical velocity, are chosen as: where c i , d i , e i and these notations with prime are polynomial and Legendre polynomials coefficients of degree i, and p i is a Legendre polynomial of degree i. Here Legendre polynomials are defined in (0, 1) rather than the original (−1, 1).
To determine the coefficients, two horizontal momentum Equations (2) and (3) are multiplied by the ith Legendre polynomial and integrated with respect to ζ on the interval (0, 1) with weighting function equaling 1. Since the integral variable ζ is independent of x and y, the integral equations are performed as Equations (8) and (9).
Substituting the Equations (5)-(7) into integral Equations (8) and (9), all integral terms except the advection terms can be explicitly reduced in relation to polynomial coefficients. The final equation sets of coefficients to be solved are written as (10) and (11).
where B N+1, N+1 is the transform matrix (12) from general polynomials of degree N to Legendre polynomials of degree N.
For the u expression in Equation set (5), B N+1, N+1 meets the condition:

Boundary Conditions
The inner shear stress close to the bed is equal to the frictional stress. A finite number of polynomials cannot accurately represent the large velocity gradient in the near-bed region. The velocity gradients and shear stress calculated from the polynomial profile deviate greatly from reality. Hereafter, by the integration of vertical diffusion terms as in Equations (14) and (15), the momentum equations have surface and bottom conditions, which can be explicitly shown as expression (18) where τ bx and τ by are contravariant components of bed shear stress are expressed using the velocity magnitude near bed u b .
where C f = gn 2 h −1/3 is the coefficient of bed shear force, and n is a resistance coefficient to be calibrated. Note that the Manning resistance coefficient refers to the depth-averaged velocity. The velocity components perpendicular to all fixed boundaries are set to zero. As the vertical velocity is obtained from the vertical integration of the continuity equation from the bottom, only the vertical velocity component at the bottom is needed, as: The vertical velocity component at the water surface is numerically satisfied as At the lateral boundaries, the slip boundary condition is assumed, and a quadratic bottom stress formulation can be applied. Here, the tangential shear stress at the lateral boundaries is neglected.
For the downstream boundary conditions, all velocity components are supposed to be uniform. Downstream water surface elevation is determined by the upstream discharge.
To make the upstream velocity components adjustable with the discharge, local surface slopes, and bottom slopes, the velocity magnitude and shape per upstream boundary grid has to be adjusted to keep consistent with the adjacent grid downstream while maintaining the total discharge after every time step. Thus, the upstream velocity can be kept uniform after a few time steps, and therefore, the initial velocity profiles can be arbitrarily set under the prescribed discharge.

Advection Terms
Under the Gauss-Legendre quadrature rule, integrals of polynomial degree no larger than 2N + 1 can be equivalently substituted by a weighted sum of function values at the N + 1 roots of the Nth Legendre polynomial as shown in Equations (23) and (24), where ζ j is the jth root of the Nth Legendre polynomial. To solve advection terms of a hyperbolic nature, Eulerian-Lagrangian methods (ELMs) or the arbitrary Lagrangian-Eulerian method [33] have been extensively investigated [34][35][36]. ELMs are attractive because the advection term is equivalent to the temporal change of one particle moving along the characteristic line (backtracking) within one time step. However, a Courant number that does not exceed 1, as a constraint of Eulerian computational power for the appearance of advection expression, would only be taken as a reference for time step scale. Advections in Equations (23) and (24) at ζ = ζ j are changed as in Equations (25) and (26).
∂uu ∂x ∂uv ∂x where u n+1 and v n+1 are velocity components solved from non-advection terms for the n+1 th time step, and u * +1 and v *n + 1 are velocity components at the root of the characteristic line after backtracking in the velocity field of u n + 1 and v n + 1 . The backtracking process can be solved either by simple one-step tracking, such as the CIP scheme [35] for a small time step, or by multi-step tracking [37] for a large time step. In this paper, we set up an ELM scheme with multi-step tracking. Divide time step ∆t into m sub-steps. Suppose a particle located at point P *k at sub-step k, originally from a certain vertex (i, j), moves back to point P *k + 1 at sub-step k + 1. The displacement is equal to the velocity magnitude at P *k multiplied by sub-step ∆t/m. This backtracking process stops at the location of P *m . Velocity components and other physical qualities at P *k are obtained by bilinear interpolation from grid vertexes. The velocity field at sub-step k + 1 is also a temporal linear interpolated value between n and n + 1.

Water Surface Elevation
The integration of mass conservation Equation (1) from the riverbed (z = z b ) to the water's surface (z = H) results in water depth Equation (27).
where u and v with bar are the depth-averaged u and v, Once the horizontal velocity components are solved, water surface elevation is calculated by Equation (27) with (28) and (29).

Numerical Algorithm
Water surface elevation and polynomial coefficients for velocity components are variables to be solved. This numerical model is set on a horizontal staggered grid in a curvilinear coordinate system. The numerical algorithm has a flow diagram ( Figure 1) and follows the sequence: 1.
In the beginning of every time step, calculate the first three terms on the r.h.s. for Equations (10) and (11) using an explicit scheme; 2.
The vertical diffusion terms (the 6th, 7th, and 8th terms on the r.h.s. of Equations (10) and (11)) is implicitly discretized and moved to the l.h.s. of the equations for the sake of stability. The water surface elevation in Equation (27) and other terms in Equations (10), (11), (23) and (24) are also calculated explicitly, but using the iteration method. The criterion for the end of iteration in each time step is that the total water surface elevation change is smaller than a given value; 3.
Solve the advection terms using the ELM scheme; 4.
Solve the vertical velocity component using Equation (1).

Laboratory Flume Experiment Verification
The sharp bend experiment by Blanckaert [13,38] was used for model evaluation in this paper. The flat bed flume consists of a straight inflow of 9 m, a 193° curved bend, and Compute diffusion and pressure gradients using explicit scheme Compute water level and pressure gradients using semi-implicit scheme

Laboratory Flume Experiment Verification
The sharp bend experiment by Blanckaert [13,38]   The simulation was run on a staggered grid with 130 × 20 cells with inlet and outlet reaches both shortened to 2.6 m (Figure 2b). The calculation time step was 0.02 s. The water surface level was calculated at the outlet, using uniform flow conditions which employ the overall slope along the bend centerline. Fixed discharge (0.089 m 3 /s) was set at the upstream boundary. The resistance coefficient was set as 0.02 and no friction was considered in the simulation on the lateral walls for the sake of simplicity. The simulation was run on a computer with intel i7 3770 CPU and 4 GB RAM. Simulation for a duration of 400 s cost a CPU time of about 252 s, 1021 s, 2085 s, and 2877 s for N = 0, 2, 4, and 6 respectively. As N = 0 is the normal depth-averaged 2D model case, the model was efficient on the 2D model level. Here only the results of N = 2 are shown. Grid independency was tested on another grid with 130 × 40 cells. For the same experiment, van Balen et al. [39] previously ran the simulation using an LES-based approach on a grid with 1260 × 192 × 24 cells, and Zeng et al. [40] employed RANS with 127 × 101 × 35 cells. The computational efficiency of this model is obviously high, and is especially obvious considering the cell number.
We chose the S90 and S180 cross sections for the velocity profile comparison. Van Balen et al. [39] also included these two cross-section data sets for comparison with the LES simulation. As shown in Figures 3 and 4, the model shows acceptable agreement for both stream-wise and transversal velocity profiles at the S90 and S180 cross sections. The averaged errors of stream-wise and transversal velocity were 0.053 m/s and 0.034 m/s for the S90 cross section, with a mean velocity of 0.419 m/s. The averaged errors of stream-wise and transversal velocity were 0.046 m/s and 0.039 m/s for the S180 cross-section. The dip phenomena observed in the measurement were successfully simulated. The transversal velocity, which is a good reference for secondary flow strength, again showed a good agreement both in magnitude and shape near the centerline. The simulation was run on a staggered grid with 130 × 20 cells with inlet and outlet reaches both shortened to 2.6 m (Figure 2b). The calculation time step was 0.02 s. The water surface level was calculated at the outlet, using uniform flow conditions which employ the overall slope along the bend centerline. Fixed discharge (0.089 m 3 /s) was set at the upstream boundary. The resistance coefficient was set as 0.02 and no friction was considered in the simulation on the lateral walls for the sake of simplicity. The simulation was run on a computer with intel i7 3770 CPU and 4 GB RAM. Simulation for a duration of 400 s cost a CPU time of about 252 s, 1021 s, 2085 s, and 2877 s for N = 0, 2, 4, and 6 respectively. As N = 0 is the normal depth-averaged 2D model case, the model was efficient on the 2D model level. Here only the results of N = 2 are shown. Grid independency was tested on another grid with 130 × 40 cells. For the same experiment, van Balen et al. [39] previously ran the simulation using an LES-based approach on a grid with 1260 × 192 × 24 cells, and Zeng et al. [40] employed RANS with 127 × 101 × 35 cells. The computational efficiency of this model is obviously high, and is especially obvious considering the cell number.
We chose the S90 and S180 cross sections for the velocity profile comparison. Van Balen et al. [39] also included these two cross-section data sets for comparison with the LES simulation. As shown in Figures 3 and 4, the model shows acceptable agreement for both stream-wise and transversal velocity profiles at the S90 and S180 cross sections. The averaged errors of stream-wise and transversal velocity were 0.053 m/s and 0.034 m/s for the S90 cross section, with a mean velocity of 0.419 m/s. The averaged errors of stream-wise and transversal velocity were 0.046 m/s and 0.039 m/s for the S180 cross-section. The dip phenomena observed in the measurement were successfully simulated. The transversal velocity, which is a good reference for secondary flow strength, again showed a good agreement both in magnitude and shape near the centerline.  In the S90 cross section, there was an obvious difference at the region close to the surface of profile 0.27B from the inner (left) bank. In the S180 cross section, there was an obvious difference in the upper region of the profiles. The measured velocity always achieved a maximum near the water surface, while the simulated maximum velocity appeared at the middle for each profile. For the transversal velocity, a weaker agreement appeared only at two profiles close to the outer bank at the S90 cross section. At the S180 cross-section, only the upper region of the profiles close to the inner bank had an obvious difference. This model apparently overestimated the secondary flow strength at the outer bank of the S90 cross section. Figure 5 gives a general comparison for the other cross sections. From cross sections S15 to P25, each cross section had an uneven distribution of velocity, and the higher-velocity core moved from the inner bank to the outer bank. The simulated distribution showed the same trend as the measured distribution. One obvious difference is that the measured higher-velocity core was located closer to the water surface.  In the S90 cross section, there was an obvious difference at the region close to the surface of profile 0.27B from the inner (left) bank. In the S180 cross section, there was an obvious difference in the upper region of the profiles. The measured velocity always achieved a maximum near the water surface, while the simulated maximum velocity appeared at the middle for each profile. For the transversal velocity, a weaker agreement appeared only at two profiles close to the outer bank at the S90 cross section. At the S180 cross-section, only the upper region of the profiles close to the inner bank had an obvious difference. This model apparently overestimated the secondary flow strength at the outer bank of the S90 cross section. Figure 5 gives a general comparison for the other cross sections. From cross sections S15 to P25, each cross section had an uneven distribution of velocity, and the higher-velocity core moved from the inner bank to the outer bank. The simulated distribution showed the same trend as the measured distribution. One obvious difference is that the measured higher-velocity core was located closer to the water surface. In the S90 cross section, there was an obvious difference at the region close to the surface of profile 0.27B from the inner (left) bank. In the S180 cross section, there was an obvious difference in the upper region of the profiles. The measured velocity always achieved a maximum near the water surface, while the simulated maximum velocity appeared at the middle for each profile. For the transversal velocity, a weaker agreement appeared only at two profiles close to the outer bank at the S90 cross section. At the S180 cross-section, only the upper region of the profiles close to the inner bank had an obvious difference. This model apparently overestimated the secondary flow strength at the outer bank of the S90 cross section. Figure 5 gives a general comparison for the other cross sections. From cross sections S15 to P25, each cross section had an uneven distribution of velocity, and the higher-velocity core moved from the inner bank to the outer bank. The simulated distribution showed the same trend as the measured distribution. One obvious difference is that the measured higher-velocity core was located closer to the water surface. . From top to bottom, the cross sections are M05, S15, S30, S60, S90, S120, S150, S180, P05, P15, and P25, as shown in Figure 1a.
The measured and predicted horizontal distributions of depth-averaged velocit illustrated in Figure 6. The simulated velocity distribution was in good agreement the measured result. Just downstream of the bend inlet, the higher-velocity core peared at the inner bank region, while the low-velocity region appeared closer to outer bank. The higher-velocity core slowly moved to the outer bank downstream cause of secondary flow in the transversal direction. At the bend outlet, the core wa cated close to the outer bank. There is only a marginal difference between Figure proving that the model had relatively strong grid independence. The core shifting . From top to bottom, the cross sections are M05, S15, S30, S60, S90, S120, S150, S180, P05, P15, and P25, as shown in Figure 1a.
The measured and predicted horizontal distributions of depth-averaged velocity are illustrated in Figure 6. The simulated velocity distribution was in good agreement with the measured result. Just downstream of the bend inlet, the higher-velocity core appeared at the inner bank region, while the low-velocity region appeared closer to the outer bank. The higher-velocity core slowly moved to the outer bank downstream because of secondary flow in the transversal direction. At the bend outlet, the core was located close to the outer bank. There is only a marginal difference between Figure 6b,c, proving that the model had relatively strong grid independence. The core shifting process across the bend is illustrated in Figure 7. The transverse distance of the core from the inner bank at each cross section along the stream-wise direction is depicted. The measured core location gradually shifted from the inner band at the bend inlet to the outlet. The present simulation results show an acceptable agreement with the measured data, though the simulated shifting became slower at upstream regions of the outlet than in the measured data.
Water 2021, 13, x FOR PEER REVIEW 12 of 15 cess across the bend is illustrated in Figure 7. The transverse distance of the core from the inner bank at each cross section along the stream-wise direction is depicted. The measured core location gradually shifted from the inner band at the bend inlet to the outlet. The present simulation results show an acceptable agreement with the measured data, though the simulated shifting became slower at upstream regions of the outlet than in the measured data.

Discussion
As the vertical velocity is calculated from a continuity equation, no vertical velocity effect on horizontal equations is considered except vertical advection terms. This kind of hydrostatic assumption is not applicable to sharp bends, where the vertical velocity is not negligible. On the other hand, the polynomials' degree is only two. When the profile is more complicated, low-degree polynomials are not enough, and the systematic error is very large. Turbulent viscosity is another mechanism for the velocity profile. As a larger velocity gradient brings higher viscosity, this would suppress the velocity gradient in return. Therefore, a parabolic distribution of eddy viscosity is not appropriate for bend open channel flow. All of these factors would introduce numerical errors to the simulation.
For a given discharge, water surface elevation depends on boundary resistance and water turbulence. As a zero-equation turbulence model is employed and turbulence decouples with flow structure in this model, adjustment of the bed resistance would not affect the flow strength effectively. Here we ignore the water surface elevation comparison. cess across the bend is illustrated in Figure 7. The transverse distance of the core from the inner bank at each cross section along the stream-wise direction is depicted. The measured core location gradually shifted from the inner band at the bend inlet to the outlet. The present simulation results show an acceptable agreement with the measured data, though the simulated shifting became slower at upstream regions of the outlet than in the measured data.

Discussion
As the vertical velocity is calculated from a continuity equation, no vertical velocity effect on horizontal equations is considered except vertical advection terms. This kind of hydrostatic assumption is not applicable to sharp bends, where the vertical velocity is not negligible. On the other hand, the polynomials' degree is only two. When the profile is more complicated, low-degree polynomials are not enough, and the systematic error is very large. Turbulent viscosity is another mechanism for the velocity profile. As a larger velocity gradient brings higher viscosity, this would suppress the velocity gradient in return. Therefore, a parabolic distribution of eddy viscosity is not appropriate for bend open channel flow. All of these factors would introduce numerical errors to the simulation.
For a given discharge, water surface elevation depends on boundary resistance and water turbulence. As a zero-equation turbulence model is employed and turbulence decouples with flow structure in this model, adjustment of the bed resistance would not affect the flow strength effectively. Here we ignore the water surface elevation comparison.

Discussion
As the vertical velocity is calculated from a continuity equation, no vertical velocity effect on horizontal equations is considered except vertical advection terms. This kind of hydrostatic assumption is not applicable to sharp bends, where the vertical velocity is not negligible. On the other hand, the polynomials' degree is only two. When the profile is more complicated, low-degree polynomials are not enough, and the systematic error is very large. Turbulent viscosity is another mechanism for the velocity profile. As a larger velocity gradient brings higher viscosity, this would suppress the velocity gradient in return. Therefore, a parabolic distribution of eddy viscosity is not appropriate for bend open channel flow. All of these factors would introduce numerical errors to the simulation.
For a given discharge, water surface elevation depends on boundary resistance and water turbulence. As a zero-equation turbulence model is employed and turbulence decouples with flow structure in this model, adjustment of the bed resistance would not affect the flow strength effectively. Here we ignore the water surface elevation comparison.
For 3D shallow water equations, Legendre polynomials were selected with the domain of relative elevation variable from 0 to 1. The depth gradient was small enough in this test case, and relative elevation was almost uniform horizontally in the local region. Otherwise, it could cause an obvious system error. That is, inhomogeneous depth would bring in extra items for momentum equations under non-orthogonal vertical and horizontal coordinates.
By using the semi-Lagrangian method, advections of momentum were solved. There was no grid in the vertical direction, and so an upwind scheme was impossible for the velocity profile calculation. We transferred the profile calculation into a weighted sum of calculation at the roots of Legendre polynomials. This transformation was accurate for this model.
In the test case, no wall friction was considered for the simulation. This led to an overestimation of velocity close to the lateral wall. A rough zero-equation turbulence model was employed for flow structure. Vertical momentum was neglected. The Legendre polynomials' degree was only 2. These simplifications led to a deviation of the simulation from the measurement data. In contrast to the measured data, there was no second circulation cell near the outer bank in the simulation. The zero-equation turbulence model is not strictly valid in the vicinity of the banks. Moreover, it will not reproduce the second circulation cell near the outer bank, as the cell is a result of the complex interaction between turbulence stresses and centrifugal effects [41]. To reproduce the second circulation cell, it is necessary to solve the vertical momentum equation without reliance on the assumption of hydrostatic pressure distribution and the use of an isotropic turbulence model.
There is no physical meaning for variable domains of ζ < 0 and ζ > 1. The momentum flux across the water surface and bed surface should be zero, which is not automatically satisfied in the spectral method. In this model, resistance from bed friction is expressed by setting non-zero flow velocity numerically on the bed boundary. The hydrostatic pressure assumption is applicable if the vertical velocity or its gradient is assumed to be small enough. However, for the flow at a sharp bend, a vertical velocity is very large and a non-hydrostatic pressure model is appropriate.

Conclusions
The authors propose a new 3D numerical model using a 2D plane model combined with a vertical spectral method. The present model was applied to simulate open channel flow around a sharp bend, and the simulated flow structures are compared with the experimental results, while a zero-equation turbulence model is employed. The obvious deviation of the simulated velocity from the measured velocity may have been caused by the presumed viscosity distribution. Another factor is that the Legendre polynomials employed were of degree from 0 to 2. Both the vertical profiles at the selected cross sections and the horizontal distribution of the depth-averaged velocity showed rather good agreement between the experimental and simulated data. We conclude that this 3D model has the ability to simulate the complicated flow structure of a sharp bend channel. The present computational results suggest that the proposed model is a powerful tool to predict and analyze open channel flows around sharp bends.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data are available on request from the authors.