1 Introduction

The bed morphology of a river bend is the result of the effect of its complex dynamic factors, and the relationship between the main factors and the bed response has been investigated by researchers all over the world. In recent decades, theoretical models and laboratory tests have been carried out by hydrodynamicists and geomorphologists to quantitatively depict the processes of a variety of bedforms. The results of the models and tests show that most of the morphological features are related to the mechanism of hydrodynamic instability, which is caused by the unstable interface between the flow and the erodible bed. The first investigator who mentioned this mechanism was Yalin^[1], who once stated that local disturbance could cause flow instability, and accordingly resulted in sand ripples. However, he did not give a theoretical model to explain it. After that, Johannesson^[2] and Johannesson and Parker^[3] established mathematical models for the erodible beds in river bends. They added small perturbation into the equations to investigate the mechanism that the flow shaped the bed morphology. Bai and Luo^[4] stated that the shear stress in the bedform caused by disturbance waves, coherent structures near the bed boundary, and the flow itself was greater than the Shields shear stress, and accordingly sand ripples appeared. They also studied the loss of laminar flow stability in open channels and the mechanism of sand ripple formation. Seminara^[5] and Pittaluga et al.^[6] established a linear model and a weakly nonlinear model, and studied the different responses of the bed in curved river. Martin and Jerolmack^[7] performed the flume experiment under the effect of unsteady flow rates, obtained the response of the sandy bed, and quantitatively described the relationship between different flow rates and bedforms. Nelson and Dube^[8], based on extreme floods, calculated the dynamic geomorphology of the channel in view of three conditions including the uniform channel width, non-uniform channel width, and actual river.

Xu and Bai^[9] built a theoretical model for the resonant triad interaction between the flow perturbation and the bedform, and analyzed the dynamic process and evolution of sandy ripples in straight river. Xu and Bai^[10] established a theoretical model for the flow in a narrow, deep channel with variable curvature, studied the unsteady flow characteristics, and discussed the relationship between the parameters of the curvature, wavenumber, and wave frequency and the dynamic instability of the flow with the Muller method and quick response method. Xu et al.^[11] investigated the erodible bed with fixed sides, and showed that the point bars and pools appearing in the channel bends had an obvious tendency of downstream migration and the sediment transport rate increased when the inlet flow rate increased while decreased when the channel sinuousness increased. Based on the narrow, deep channel bend with constant curvature, Bai et al.^[12] studied the characteristics of the dynamic stability of the flow, and obtained the theoretical formula for the flow velocity distribution and the critical Reynolds number. The results showed that, compared with straight channel (ψ = 0), stable neutral curve moved forward along the coordinate, the value of the unstable Reynolds number increased, the response range of the disturbance wavenumber decreased, and the flow pattern was much easier to maintain. Bai et al.^[13] theoretically analyzed the stability and adaptive feature of the turbulent coherent structure with constant curvature in the channel bends, assigned the turbulent coherent structure as a perturbation, calculated the perturbation growth rate of the coherent vortex and the response range of the wavenumber, gained the response of different scaled turbulent structures affected by the river bends, and explained the dynamic mechanism that river bends could maintain.

In this paper, we establish a theoretical model for meandering river by coupling the incompressible viscous fluid, sediment transport, and bed deformation equations in a curvilinear coordinate system, where the depth-to-width ratio of the channel is greater than 5. The flow characteristics and bed morphology of the river are analyzed with the perturbation method (see Section 2). The results are discussed in Section 3 in order to forecast the evolving trend of the bedforms under different flow conditions. Finally, conclusions are described in Section 4.

2 Theoretical model 2.1 Coordinate transformation

The flow field in a meandering river is referred to as an orthogonal reference system. The relationship between the Cartesian coordinate system and the orthogonal curvilinear coordinate system is shown in Fig. 1, where s is the longitudinal coordinate of the channel axis lying on a plane, n is the transverse coordinate defined along a horizontal axis orthogonal to s, and z is the coordinate of the axis orthogonal to s and n and pointing upwards. In the Cartesian coordinate system, x is the coordinate of the horizontal axis, y is the normal coordinate perpendicular to the horizontal axis, and b is the half width of the channel.

In Fig. 1, (x₀, y₀) is the centerline coordinate of the channel, and (x, y) is the arbitrary point coordinate.

According to Fig. 1, the conversion between the orthogonal curvilinear coordinate system and the Cartesian coordinate system is , and .

The curvature radius is

The Lame coefficients are h_s = 1−N, and h_n = 1, where h_s and h_n are the Lame coefficients in the s-and n-directions, respectively. When the shape of the channel bend is close to the sine or cosine function shape, these parameters can be expressed as follows:

where α is the wavenumber of the river, and ω is the propagating frequency of the river. For a fixed channel bend with constant curvature, ω = 0, and α = 0.

2.2 Dimensionless equations and parameterization

Half of the river width b, the water depth H, and the minimum radius of the curvature r_m are spatial independent variables. The average velocity U is an independent variable of velocity. The independent variable of time is b/U. Therefore, the following terms can be determined:

(ⅰ) The dimensionless terms of space are

(ⅱ) The dimensionless terms of velocity are

(ⅲ) The dimensionless terms of time are

In the above equations, the terms with ~ refer to their dimensionless forms. Define , and , where ψ is the sinuousness, which is the ratio of the half-width to the curvature radius in the channel centerline, and γ is the ratio of the half-width to the water depth. For the fixed channel with constant curvature, ψ and γ can be considered as small parameters.

The dimensionless continuity equation and Navier-Stokes (N-S) equations can be obtained as follows:

(1)

(2)

(3)

(4)

where

Re is the Reynolds number defined by Re = Ub/ν, and ρ is the density of water. u_s, u_n, and u_z are the velocity components in the s-, n-, and z-directions, respectively. f_s, f_n, and f_z are the mass forces in the s-, n-, and z-directions, respectively, f_s = gsinθ = gJ, f_n = 0, and f_z = −g(1 − cosθ). J = sinθ is the bottom slope of the river (streamwise). P is the pressure. ν is the kinematic viscosity, and ν = 0.077U_*H, where U_* is the shear velocity near the bed.

The dimensionless bed response equation is

(5)

where some dimensionless terms are given as follows:

In the above equations, z_b is the bed surface elevation, q_s and q_n are the sediment transport rates of the unit width in the s-and n-directions, respectively. D_s is the median diameter of the sediment, and D_s = 0.5 mm is adopted in this research from the natural river sand. The sediment porosity p' = 0.4^[14]. R_q is the submerged specific gravity of the sediment, and R_q = 2 650 kg/m³.

(ⅳ) The dimensionless sediment transport equations are

(6)

(7)

where the velocity index formula of the sediment transport rate is adopted from EngelundHansen. M is the empirical coefficient, and M = 5 (plane bed)^[2]. is the transverse velocity at the initial time. T(ξ) is the dimensionless function of the velocity distribution, where ξ = z/h, and h is the real water depth. γ is the ratio of the half-width to water depth. β is the empirical coefficient.

The empirical coefficient β^[18] can be expressed as follows:

where α_* is the ratio of the lift coefficient to the drag coefficient for a spherical sand particle placed on a rough bed, μ is the dynamic angle of the Coulomb friction, f_* is an order-one coefficient to be determined from the data, τ_C^* is the critical Shields share stress, and τ_G^* is the grain Shields share stress. The suggested values of Kikkawa^[19] are α_* = 0.85 and μ = 0.43, and the value of f_* is determined from Zimmermann and Kennedy^[20] as f_* = 1.19. It is assumed that all shear stresses are active in the bed-load process rather than just τ_G^*. Therefore, τ_G^* = τ^* (plane bed), and τ^* = /(ρRgD_s).

2.3 Perturbation analysis

In a meandering channel, the flow structure is affected by both the hydrodynamic instability itself and the curved boundary of the channel. According to the theoretical method for the calculation of the turbulent structure in fluid mechanics^[15-17], the flow structure in the channel with a slight curve can be solved with the perturbation method.

The solution of the flow structure is given as follows:

(8)

where , , and are the different velocity components in the s-, n-, and z-directions, respectively, , , and . is the pressure, and . , and are small disturbances acting on the flow of the channel. ε_T, representing the scale of the small disturbances, is the internal disturbance parameter, and is defined by the ratio of the amplitude of the disturbance velocity to the average velocity.

Since the channel width is much smaller than the radius of the curvature, ψ is a small parameter. The terms of the flow, i.e., , and (i = 1, 2, …), are the correction terms to the flow in the straight channel, which are caused by the curved configuration of the channel.

2.3.1 Analysis of the basic term and bending correction term

(Ⅰ) Solution of the ψ⁰ term

If ψ = 0, the terms (i≥1) in Eq. (8) will vanish. Then, Eqs. (8)–(11) will turn into the flow dynamic equations, which describe straight channels. Therefore, the corresponding terms , and actually represent the flow parameters of the straight channel.

For the steady uniform flow,

The boundary conditions are given as follows:

(ⅰ) There are no-slip bank conditions at = ±1, and (±1, z) = 0.

(ⅱ) There is a zero-shearing stress boundary condition at the water surface , and , where is the depth of the free water surface, and is the atmospheric pressure at the water surface.

(ⅲ) There are no-slip conditions on the bed at , and .

In view of a large depth-to-width ratio, i.e., γ = (b/H) 1, Eqs. (1)–(4) can be simplified into two-dimensional forms, and the following expressions can be obtained:

(9)

(10)

(11)

where Fr is the Froude number, and Re is the Reynolds number.

(Ⅱ) Solution of the ψ¹ term, i.e., bending modification

For the boundary conditions of banks at = ±1, = 0, and = 0. With the same procedure and the expressions shown in Eqs. (9)–(11), we can express the correcting solutions as follows:

(12)

(13)

(14)

2.3.2 Stability analysis of coherent disturbances

The curved geometry and fluctuating characteristics of river bends can further induce the instability of the flow. In view of the slightly bending channels with large depth-to-width ratios, the second and higher order terms of ψ can be neglected. Therefore, Eq. (8) can be simplified as follows:

(15)

The solution of the disturbance can be written as follows:

(16)

Substituting Eq. (15) into Eqs. (1)–(4) and subtracting the corresponding parts of the basic flow, we can obtain the amplitude equations with two-dimensional disturbance as follows:

(17)

(18)

(19)

Equations (17)–(19) can be solved with the boundary conditions = ±1 and , and the solution satisfies

(20)

where

In the above equations, .

If it is the temporal mode question, the wavenumber of the disturbance α is a real number, whereas the frequency of the disturbance ω = α_*c is a complex number, where c is the propagating speed of disturbance. The positive or negative sign of the imaginary part of ω represents the growth or decay of the disturbance.

Since Eq. (20) has nonzero solutions, the determinant of the coefficient matrix should be equal to zero, i.e.,

(21)

which is the Othmer-Stevens (O-S) equation depicting the stability characteristics of the flow in the channel bends. From Eq. (21), it can be observed that, different from the classical O-S equation, some bending parameters, e.g., , ψ, and γ, are added into the O-S equation, which become significant parameters for the hydrodynamic instability characteristics of the bending river.

After providing the parameters Re, α, ψ, , and γJ, the parameter ω can be gained. When the imaginary part of the disturbance ω_i is positive, the coherent disturbances develop over time, and the flow change becomes more and more unstable.

2.4 Solution to the bed response equation

Substituting Eqs. (15) into Eqs. (5)–(7) yields

(22)

where

Equation (22), which has been simplified by the author (the detailed derivation process is shown in Appendix A), is the second-order linear ordinary differential equation with variable coefficients. Then, we can easily solve the equation with boundary conditions.

The no-slip bank condition is at = ±1, = 0, = 0, = 0, and = 0. Thus,

(23)

When the discrete governing equation and boundary conditions are given, we can get the solution with the chasing method in the domain −1 ≤ ≤ 1.

After is determined, we can gain the equation of bed deformation as follows:

(24)

3 Calculation results and analysis 3.1 Flow field calculation

In the mathematic model used in this paper, the channel with the large depth-to-width ratio (H/B = 5) is considered. Therefore, the change in the flow velocity along the vertical direction can be ignored.

The velocity distribution of the basic flow along the transverse direction is shown in Fig. 2. The velocity of the basic flow presents a parabolic distribution, and it is a function of the parameter n. when n = 0, the flow velocity reaches the maximum. When n = ±1, the flow velocity is equal to 0, and it satisfies the no-slip boundary condition.

Considering the above velocity distribution, we now introduce the small disturbance of the streamwise and crosswise directions, and discuss the response of the bedform of the river with large depth-to-width ratio.

Figure 3 is the shape function of the flow velocity considering the turbulent coherent perturbation. From Figs. 3(a) and 3(b), we can see that, if the bending is close to zero (ψ = 0.00), the shape function of the streamwise velocity will be antisymmetric with respect to the centerline of the channel, while the shape function of the crosswise velocity is symmetric. With the increase in the sinuousness (ψ = 0.00, ψ = 0.05, and ψ = 0.10), the shape function of the streamwise disturbance velocity is skewed, and the velocity value is significantly enlarged in the concave bank. In Fig. 3(b), the value of decreases, and the distribution is close to symmetric.

3.2 Response of the bedform

Figure 4 is the shape function of the bed response in view of various sinuousnesses and Reynolds numbers (see Eq. (22)), where is the real part of the shape function . As shown in Fig. 4, the elevation of the bed response is approximately the sinusoidal distribution along the crosswise direction. Within a certain range, its amplitude decreases when the Reynolds number increases. The amplitude of the bed response climbs when the sinuousness increases (see Fig. 5). When the hydraulic gradient in the channel bends increases, the interaction between the flow and the banks becomes stronger, the concave bank experiences erosion, and the convex bank experiences deposition.

It is shown that the growth rate is the maximum at some value of the meander wavenumber^[5]. Moreover, the meandering river exists as typical nonlinear waves, and the selected wavenumber depends on the amplitude of the initial perturbation. Bai et al.^[13] found that multiple disturbance waves with different growth rates might occur in the turbulent flow of the channel under the same hydraulic conditions, and presented a series of characteristic models, in which the characteristic model with the maximum growth rate was defined as the least stable spectra. In this paper, we will discuss the least stable spectra and the demarcation point between the stable band and the unstable band of the characteristic model around α = 1.00 (Model Ⅰ) and α = 1.05 (Model Ⅱ). The unsteady waveband of Model Ⅰ, in which the imaginary part of the disturbance frequency is greater than 0.00, is studied in detail (see Figs. 6–8). It is shown that, in view of various Reynolds numbers, the disturbance wave group near the characteristic wavenumber α = 1.00 is discussed, and the characteristic spectra of the wavenumber-growth rate can be seen in Figs. 6–8 when ψ = 0.00, 0.05, and 0.10. With the same routine, we also discuss the unsteady waveband of Model Ⅱ, in which the imaginary part of the disturbance frequency is greater than 0.00. The corresponding characteristic spectra of the wavenumber-growth rate when the sinuousness ψ is 0.00 is presented in Fig. 9.

According to the foregoing characteristic spectra of the unable waveband in Figs. 6(b), 7(b), 8(b), and 9(b), which are parts of Figs. 6(a), 7(a), 8(a), and 9(a), respectively, the bed response will be affected by the disturbance wave of the turbulent. The disturbance frequencies, which corresponds to different Reynolds numbers and curvatures in Figs. 6–9, are shown in Table 1. From Table 1, we can see that the value of the disturbance frequency declines when the sinuousness increases, and the imaginary part of the disturbance frequency reaches its maximum when the growth rate reaches a peak value.

From Figs. 6–9, we can also see that, with increased bending, the growth rates of the coherent disturbance wave in the channel with ψ = 0.05 and ψ = 0.01 are significantly smaller than that of the channel with ψ = 0.00, and the growth rate increases gradually when the Reynolds number increases. Similarly, the turbulent flow is more drastic and unstable with a larger Reynolds number, and the response range of the disturbance wavenumber moves towards smaller wave numbers with a larger growth rate.

By searching the model and analyzing the characteristic spectrum of the wavenumber-growth rate, the unstable range for different Reynolds numbers can be obtained (see Figs. 6(b)–9(b)).In order to further explore the bed response, Figs. 6(c), 7(c), 8(c), and 9(c) are given, which are the characteristic spectra when the change in the Reynolds number is very small, where α = 1.00 and α = 1.05 are the characteristic wavenumbers, and ω_i = 0 is the basic line representing the growth rate of 0.

According to the results of Figs. 6–9, the bed response can be further calculated in view of the characteristic wavenumbers α = 1.00 and α = 1.05 in various sinuousnesses including ψ = 0.00, ψ = 0.05, and ψ = 0.10.

Using α = 1 as a typical example, the shape function characteristics of the bed response can be analyzed. In the micro-bend channel, the contour map of the shape function of the bed response in two wavelengths is shown in Figs. 10 and 11, which correspond to the sinuousnesses ψ = 0.05 and ψ = 0.10, respectively. As shown in Figs. 10 and 11, the amplitude of the bed response decreases when the Reynolds number increases. Since the coherent disturbance gets more unstable with the increase in the Reynolds number, it is difficult to form the vortex turbulence structure. Thus, the amplitude of the bed response tends to decrease. At the same time, the shape function of the amplitude demonstrates a wave morphology presenting antisymmetrical distribution with regard to the centerline channel. When the bending increases, the amplitude of the bed response shows an escalating trend.

The bed responses in two wavelengths are presented in Figs. 12 and 13, which corresponds to Fig. 8(c), where T is the period of sand ripples. The charts in Figs. 12 and 13 show that, when Re = 7 923, and Re = 25 152, the bed response decreases with time, while when Re = 7 924 and Re = 25 151, the bed response increases with time. These results prove that the method in this study can accurately discriminate the response of the bed morphology.

The characteristics of the bed response are analyzed, where ψ = 0.10. Figures 12(a) and 13(b) correspond to the bed response attenuation, and Figs. 12(b) and 13(a) correspond to the increasing bedform. Combined with the corresponding legend, it shows that the bed response varies with time. For the attenuation condition, the bed response eventually disappears. Figures 12(a) and 12(b) present the decay and growth of the bed response at Re = 7 923 and Re = 7 924, respectively. The results strongly establish that the method can accurately distinguish the increasing or decreasing bed response.

4 Conclusions

This research discusses the bed response in rivers with large depth-to-width ratios and constant curvature. The unstable range of the Reynolds number and the bed response are obtained by searching the model and analyzing the characteristic spectra of the disturbance wavenumber-growth rate.

(ⅰ) When curved geometry appears, the channel bend will affect the development of the disturbance wave. Then, with the increase in the sinuousness (ψ = 0.05 and ψ = 0.10), the growth rate of the disturbance wave will be obviously smaller than that of the straight channel (ψ = 0.00). This is consistent with the phenomenon that natural rivers always develop a specific curvature and tend to be relatively stable.

(ⅱ) For various types of bending and Reynolds numbers, the shape function of the bed response in the crosswise direction presents a sinusoidal distribution, and the value of bed response near the convex (concave) bank is positive (negative). A higher Reynolds number leads to a much more unstable disturbance, which implies that it is difficult to form a large turbulent vortex structure. Thus, the amplitude of the bed response decreases with increasing the Reynolds number, and the bed response develops toward the stable form.

(ⅲ) The shape function of the bed response amplitude is expressed as the antisymmetrical distribution of the channel centerline. The amplitude of the bed response expands with the increase in the sinuousness, while it reduces with the increase in the Reynolds number. When the sinuousness is uniform, the imaginary part of the disturbance frequency gets larger with the increase in the Reynolds number, which leads to rushing, tempestuous turbulence, and also unstable flow.

(ⅳ) This study obtains methods for identifying the trend of bed response development with the curvature, Reynolds number, disturbance wavenumber, and bed morphology gradient as important parameters.

Appendix A

For Eq. (23), a more detailed derivation process is shown as follows:

(A1)

Substituting Eq. (15) into Eqs. (5) and (6) and only considering the first-order term of the perturbation velocity, we have

(A2)

(A3)

Substitute Eqs. (A2) and (A3) into Eq. (A1), carry out integral on every term of Eq. (A1) in time, and consider the simplified conditions as follows:

Then, we have

(A4)

Equation (A4) can be simplified as follows:

(A5)

where