Deformable River Beds: A Numerical Study

1. Introduction

In most cases, river channels' position and depth change over time through channel deformations. Therefore, a forecast of channel deformations is required for high-quality design of relevant structures and work performance. In addition, most rivers are regulated due to the construction of reservoirs for energy and agricultural needs. There are also a large number of water intake structures that significantly affect the flow dynamics. All these rivers' special flow aspects significantly change the natural course of the riverbed process, and a forecast of rivers' bed change is required. Therefore, studying and developing the theory of channel processes and the dynamics of channel flows has always attracted the attention of researchers [1,2].

However, despite the abundant work devoted to forecasting channel deformations, its solution is still far from practical completion due to the complex and multifactorial nature of riverbed processes in space and time. Particularly great difficulties arise when designing various river structures in rivers, the bed of which, due to large bottom slopes, high flow velocities, and easy erosion of bottom sediments (represented by fine sandy soft soils), is subject to extremely complex and intense deformations. An example of such a river in Uzbekistan is the Amu Darya.

The behavior of channel processes in easily eroded soils is very complex. Due to such complexity and multifactorial causes that determine channel processes, as well as due to the lack of a rigorous theoretical solution to the problems of river hydraulics and the dynamics of channel flows, researchers resort to the use of methods of physical and numerical modeling of channel processes [3,4]. Both of these methods complement each other, and as a result, it is possible to obtain more reliable solutions to the problem posed by the deformation of the channel in a specific section.

Numerical or physical modeling is used to solve the problem of channel processes, which can give a specific forecast of channel deformations in the area of the structure. Despite such long practical experience, achieving an acceptable balance between the requirements of water transportation systems and the structures themselves, especially in the conditions of meandering channels composed of easily eroded soils, is a complex and pressing engineering task [5,6].

A typical research object of the channel processes is the Amu Darya River, which is easily eroded in the middle reaches. Numerical and physical are the most important directions in the channel process research. However, big problems with similar and distinctive features arise in both research directions.

The general problem is finding the basic physical laws and creating a mathematical model (closed system of equations) to describe the process accurately. In this finding, the mathematical model can be both stochastic and deterministic. After drawing up such a model, the problems of both modeling methods become completely different.

Physical modeling is based on finding the conditions for connecting the model and nature based on analyzing the system of equations under study [7,8,9]. This analysis is related to:

−

the use of a similarity transformation that would show that a process on a smaller scale (on a model) is equivalent to a process on a larger scale (in reality);

−

finding modeling criteria based on it;

−

establishing areas of self-similarity according to various criteria if they actually exist.

Ultimately, the rules for converting from model to actual are obtained. Very often, it turns out that a very large model scale or expensive materials is required, which is associated with large material costs.

Currently, hydraulics is experiencing a period of rapid development of numerical models. The numerical model represents a way of calculating quantities necessary for practice, described by the above system of equations. A numerical model can be either an exact or an approximate solution to a system of equations. In more or less complex cases, exact solutions are impossible, and therefore, it is necessary to use approximate numerical models. The most developed numerical models are based on discretizing temporal and spatial variables.

Moreover, the main disadvantage of such models is the unknown degree of approximation of the solution to the original equations. It is almost impossible to strictly prove a discrete problem's convergence to the original solution. It is necessary to demonstrate such convergence empirically by comparison with test problems. However, compared to physical modeling, numerical modeling requires significantly less time and labor, allowing for more multivariate studies and considering a greater number of factors influencing the process. Therefore, numerical modeling of physical processes and numerical modeling of channel processes play an increasingly important role.

For modeling channel flows, approaches based on the numerical solution of two-dimensional Saint-Venant equations have shown high efficiency and sufficient accuracy. The derivation of these equations, numerical integration algorithms, and examples of calculations are given, for example, in [5,6,7,8]. However, the Saint-Venant model is open-ended for deformable channels and requires addition, for example, with an equation or system of equations to find variables in the time and space of bottom marks.

Previous studies [10,14] have revealed that the process of deformation, including the flattening of an underwater slope with a rectilinear generatrix directed along the flow velocity, can be effectively described by a one-dimensional diffusion equation governing the bottom elevation within a watercourse:

$$(1-\mathrm{p})\frac{\partial {\mathrm{Z}}_{\mathrm{b}}}{\partial \mathrm{t}}=\frac{\partial }{\partial y}D\frac{\partial Z_b}{\partial y},$$

(1)

where:

- p is the soil porosity coefficient,

- Z_b (t,y) is the bottom elevation,

- t is the signifies time,

- y is the coordinate across the slope;(

- D is the diffusion coefficient characterizing the specific transverse sediment flow rate.

Furthermore, the studies [15] have established the relationship for the diffusion coefficient D. Analytical solutions to the initial boundary value problem associated with equation (1) were also obtained, specifically for the hypothetical case of deformation of an underwater slope with an initial profile in the form of a vertical step in a channel of unlimited width.

The authors in the work[16] facilitated the derivation of analytical solutions for a slot of finite width. Additionally, the model's performance was rigorously verified against a substantial dataset of field observations, demonstrating its effectiveness in solving practical problems:

$$(1-\mathrm{p})\frac{\partial {\mathrm{Z}}_{\mathrm{b}}}{\partial \mathrm{t}}=-\frac{\partial U{S}_{*}h}{\partial x}-\frac{\partial V{S}_{*}h}{\partial y},$$

(2)

where:

−

U,V is the components of the depth-average water velocity vector along the x and y axes, respectively,

−

h (x,y,t) is the flow depth

−

S (x,y,t)* is the saturation turbidity, denoting the vertical average volume concentration of sediment in an equivalent uniform flow.

Equation (2) has been instrumental in predicting channel transformations by assuming that saturation turbidity, or equilibrium sediment concentration, primarily depends on local flow characteristics and is related to them like uniform flow conditions. Using equation (2) in conjunction with the two-dimensional Saint-Venant equations has yielded satisfactory results in predicting channel transformations for various scenarios [17,18].

Equation (2) presents limitations for describing planned deformations in channels that feature bottom sections with relatively steep slopes. Such slopes are dredging slots, quarries, or steep banks at bends. In a straight channel with an underwater slope, where the generatrix aligns with the velocity vector (and both the velocity vector and flow depth remain constant along the longitudinal coordinate x), equation (2) implies that Z_b=const. Consequently, equation (2) does not account for a well-recognized phenomenon: the alteration of an underwater slope with a generatrix parallel to the flow velocity vector. The non-accounting arises from the assumption in equation (2) of the collinearity of water velocity and solid phase flow vectors, which is not true in all cases.

Several researchers have attempted to enhance equation (2) by introducing terms that enable it to account for the Several researchers have attempted to enhance equation (2) by introducing terms that enable it to account for the effect of noncollinearity of water velocity and solid phase flow vectors. In works [19,20], the theoretical justification for incorporating diffusion terms into equation (2) was established. The diffusion terms' incorporating led to the derivation of an expression for the diffusion coefficient of the bottom mark, represented as a spherical tensor, with proportionality to the longitudinal specific sediment flow rate in the direction of the velocity vector:

$$(1-\mathrm{p})\frac{\partial {\mathrm{Z}}_{\mathrm{b}}}{\partial \mathrm{t}}+\frac{\partial U{S}_{*}h}{\partial x}+\frac{\partial V{S}_{*}h}{\partial y}=\frac{\partial }{\partial x}D\frac{\partial Z_b}{\partial x}+\frac{\partial }{\partial y}D\frac{\partial Z_b}{\partial y}$$

(3)

where $D=ahS{}_{*}\left|\overrightarrow{U}\right|$, $\left|\overrightarrow{U}\right|=\sqrt{U^2+V^2}$, and $\alpha$ is proportionality factor.

An alternative model was introduced in [21,22], which diverges notably from the prior one. The key distinction lies in formulating the diffusion term exclusively in the direction orthogonal to the flow velocity vector. Additionally, an extra term is incorporated to consider the streamline curvature's impact on bottom elevations' alterations. The diffusion coefficient is proportionate not to the flow velocity but to the value of the non-shearing flow velocity. Applying this alternative model has yielded remarkable results in predicting channel reformations. The previously mentioned deformation models [21,22] all operate assuming that sediment particle concentration in the flow closely approximates equilibrium conditions. It has been established in [10] that the specific sediment flow from the bottom into the flow thickness is proportionate to the value ($S-S_{*}$).

The current study aims to amalgamate the strengths of various previously mentioned models and propose a mathematical model for sediment transport in turbulent and dynamic river flows. It is recommended to conduct numerical simulations of river flows characterized by deformable bottoms, as proposed in [8], to employ a mathematical model governed by a system of two-dimensional Saint-Venant equations [17] incorporating partial derivatives and closing relations.

The authors in works [25,26,27] determined the dependence and influence of the channel shape on the ratio of the friction coefficient and the Reynolds number for laminar flow in an open channel. In the articles [28,29,30], the researchers conducted an experimental analysis of the rotation of carbon fibers flowing in the discharge channel at an angle of 90°. The velocity and velocity gradients inside the channel are calculated using computer hydrodynamics simulation. It was found that the fibers are affected by a relatively high local shear rate and that they propagate together with the flow, which indicates a parallel alignment in the flow direction while simultaneously demonstrating a 90° rotation in the bend angle. In studies [31,32,33], the authors performed numerical simulations of turbulent flow in a channel with superhydrophobic surfaces. The contribution of turbulence to the volumetric average velocity was found to remain almost unchanged in the absence of sliding due to the negative coherent component of the Reynolds shear stress. In articles [34,35,36,37], the thermal characteristics of a compressible laminar flow of natural convection induced at a high-temperature difference in a vertical channel with an open end were studied by optimizing the distance between the channel plates using numerical modeling. As a result, the authors presented a correlation of the optimal aspect ratio with the Rayleigh number, which maximizes heat transfer inside the channel. In [38,39,40,41,42,43], the authors investigate stable, incompressible, laminar flow in a channel bounded by rough and/or permeable walls. The study found a closed solution of the Navier–Stokes equations for flow in a channel with conditions at each boundary linking sliding velocities with shear stress and pressure gradient along the flow.

Simulations of the hydrodynamic and dynamic changes in the relief of the riverbed were carried out, and the models for accurate calculation of changes in the layer level and the area of deposition and erosion were proposed [44,45,46,47]. In the work [48,49,50], a hydraulic experiment was conducted to study the hydraulic phenomena of the dam, comparing hydraulic surges and flow characteristics. It was found that although the sluice gates generated hydraulic surges similar to those in stationary dams, their supercritical flow increased downstream, which ultimately lengthened the overall hydraulic surge.

Based on the above analysis, the numerical research method is identified as the main research method for this work. According to the purpose of the study, the main objectives of the study are outlined in the following interpretation:

-

selection of the basic flow equation;

-

conducting numerical studies of the reformation of inclined walls of a channel wall with a moving bottom;

-

conducting experimental studies

-

comparison of the obtained results with the results of experimental studies.

2. Methods

2.1. Analytical studies

The analytical approach employed in this work involves:

−

Choice of Flow Equation: The selection of the equation governing flow movement in the river bed is represented by a system of two-phase hydrodynamic equations of Saint-Venant. These equations are complemented by including equations related to sediment balance, transport, and the Bagnold equation.

−

Numerical and Experimental Studies: Numerical and experimental studies are conducted under conditions of a moving bottom by utilizing the chosen equations.

−

Comparison of Results: The obtained results from both numerical and experimental studies are thoroughly compared and analyzed.

This analytical approach aimed to comprehensively understand and model the dynamics of riverbed processes and deformations.

2.1. Experimental studies

The laboratory experiments were conducted to validate the model. The experiments focused on a slope's erosion by a current, where the current's velocity vector aligns with the slope's generatrix. These experiments were carried out in the Hydraulics Laboratory of Moscow State University of Civil Engineering, Moscow, Russia. The experimental setup and procedures were consistent with those previously documented in references [10,11,12,13].

The laboratory experiments were conducted within a channel featuring variable slope dimensions, with the following primary specifications: a length of 18 meters, a width of 2 meters, and a wall height of 0.8 meters. The channel's capacity allowed for a flow rate of up to 4 cubic meters per second, and the slope ranged from 0 to 0.1. Water flow measurements were facilitated through a triangular weir positioned at the channel's outlet. A specialized pocket was integrated at the end of the flume's working section to collect sediment to accommodate the experiments. During the experiments, a swinging shield was also employed to regulate the channel's flow level. The longitudinal slope maintained during the experiments was 0.0027.

Figure 1 illustrates the eroded model. Table 1 shows the sand's granulometric composition. Specifically, the sand had a mean diameter of 0.24 millimeters.

The cross-sectional shape of the model closely resembled half the cross-section of a trapezoidal channel with a slope coefficient represented by m=2. Before commencing water flow, a thorough horizontal survey of the bottom was performed. Water was incrementally introduced into the flume, simultaneously from the upper and lower tails, to ensure minimal bottom disturbances before the start of each experiment. Throughout the experiments, measurements were taken at defined intervals to record the water's bottom surface and free surface in fixed sections. These measurements were then averaged to provide comprehensive data for analysis.

3. Results

3.1. Mathematical model

As recommended in [8], a mathematical model based on a system of two-dimensional Saint-Venant equations [17] incorporating partial derivatives and closing relations to perform numerical calculations of river flows characterized by deformable bottoms is employed. This mathematical model serves as a foundational framework for understanding and predicting the complex dynamics of river flows:

$$\frac{\partial hS}{\partial t}+\frac{\partial USh}{\partial x}+\frac{\partial VSh}{\partial y}=-K(S-S_{*})$$

(4)

$$(1-\mathrm{p})\frac{\partial {\mathrm{Z}}_{\mathrm{b}}}{\partial \mathrm{t}}=K(S-S_{*})+\frac{\partial }{\partial x}D\frac{\partial Z_b}{\partial x}+\frac{\partial }{\partial y}D\frac{\partial Z_b}{\partial y}$$

(5)

$$K=\left\{\begin{array}{c}\alpha U_{*}+(1-\alpha )W,\hspace{1em}{U}_{*}\ge W\\ W,\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{U}_{*}\le W\end{array}\right.; 0\le \alpha <1$$

(6)

$$D=\beta \hspace{0.17em}\tilde{S}h\hspace{0.17em}W$$

(7)

$$S_{*}={\alpha }_1\hspace{0.17em}\frac{\lambda \rho }{2{\rho }_s}\frac{{(\left|\overrightarrow{U}\right|-U_N)}^2}{gh}\left(\frac{0.13}{tg\phi }+0.01\frac{\left|\overrightarrow{U}\right|}{W}\right), \lambda =2g{n}^2{h}^{-1/3}$$

(8)

where $t$ is the time; $h$ is the flow depth; U and V are the components of the flow velocity along the X and Y axis, respectively; $\left|\overrightarrow{U}\right|=\sqrt{U^2+V^2}$; $S$ is the volumetric concentration of sediment particles in the flow; $S_{*}$ is the equilibrium volume concentration of particles (saturation concentration), taken according to the modified Bagnold formula; $K$ is the intensity coefficient of vertical sediment exchange between the bottom and the stream; $p$ is soil porosity (the ratio of the volume of pores to the volume of the entire soil with pores); ${\rho }_s,\rho$ is the densities of soil and water, respectively; φ is the angle of internal friction of the soil; $W$ is the hydraulic soil coarseness; $U_{*}$ is the dynamic speed; $\left|\overrightarrow{U}\right|$, $U_N$ is the module of the average vertical flow velocity and the non-shear velocity, respectively; $\lambda$ is the coefficient of hydraulic friction, calculated using the Manning formula; $n$ is the roughness coefficient.

For determining the non-shear velocity in the calculations, the formula provided by representatives of channel hydraulics [23,24] is utilized, which, with consideration of standard coefficient values, is expressed in two formally equivalent forms:

$$U_N=\lg \frac{8.8 h}{d_{90}}\sqrt{\frac{2}{3.5}\left[({\rho }_s-\rho )g{d}_{50}+1.25C_y^H\right]}$$

(9)

$$U_N=0.18\sqrt{\frac{2}{3.5\lambda }\left[({\rho }_s-\rho )g{d}_{50}+1.25C_y^H\right]}$$

(10)

where $C_y^H$ is soil adhesion in t/m²; $d_{50}$ is the average diameter of soil particles, $d_{90}$ is the 90% of soil particles' diameter.

In the calculations of the diffusion coefficient (equation 7), three variants of formulas for $\widetilde{ S}$ are employed:

a) Based on the total equilibrium concentration of transportable and suspended sediments:

$$\tilde{S}=S_{*}$$

(11)

b) Derived from bottom equilibrium concentration:

$$\tilde{S}=\hspace{0.17em}{\alpha }_1\frac{\lambda \rho }{2{\rho }_s}\frac{{(\left|\overrightarrow{U}\right|-U_N)}^2\cdot 0.13}{gh\hspace{0.17em} tg\phi }$$

(12)

c) Utilizing bottom concentration "without square"

$$\tilde{S}=\hspace{0.17em}{\alpha }_1\frac{\lambda \rho }{2{\rho }_s}\frac{(\left|\overrightarrow{U}\right|-U_N)\left|\overrightarrow{U}\right|\cdot 0.13}{gh\hspace{0.17em} tg\phi }$$

(13)

The initial conditions are defined for the initial bottom surface $Z(x,y,0)$, the corresponding instantaneous fields of the velocity $\overrightarrow{U}(x,y,0)$, the depth h(x,y,0), and the concentration S (x,y,0).

The boundary conditions for (9-13) are as follows:

−

On solid boundaries, the condition of no flow is specified.

−

For liquid boundaries, flow rates or water levels are typically specified.

−

Water flows into the computational domain through the boundaries of the computational domain, and the precipitation concentration is set at these boundaries.

−

Complex boundary conditions were also sometimes used. The complex boundary conditions could link costs with levels and non-reflective boundary conditions.

The solution for equations (4) and (5) concerning sediment concentration and bottom marks was carried out using the finite volume method on mixed triangular-quadrangular grids combined with the Saint-Venant equations.

The developed numerical scheme aligns with the scheme for the continuity equation of the liquid phase, which helps prevent the occurrence of so-called dipoles as sources and sinks of mass. It utilizes a directed difference type scheme to eliminate unphysical oscillations in the bottom topography, maintains the transportability property, and implements a difference analog of mass conservation for the solid phase.

Equations (4) and (5) represent a minimal approach for modeling processes of bottom deformation in uneven and non-stationary river flows. Disregarding equation (4) implies an assumption that the concentration in the flow is close to equilibrium and does not permit the specification of boundary conditions for concentrations. Eliminating terms from equation (5) that describe the diffusion of bottom marks makes it impossible to account for the processes of transformation (flattening) of underwater slopes. Moreover, to describe three distinct physical processes (uplift-sedimentation, longitudinal sediment transport, transverse diffusion of bottom marks), at least three empirical coefficients are required. In this model, the empirical coefficients are α, β, α₁. The choice of the closing relations for model (6)-(13) remains an open question and necessitates further research.

Calculations for identical conditions were conducted, assuming a value of 0.5 for the test problems. The type of closing relations (7)-(13) and parameters were varied during the calculation process.

3.2. Numerical model

The calculations were conducted using (4) and (5). A rectangular grid consisting of 1800 cells measuring 0.1 m × 0.2 m was constructed for a tray with 18.0 m × 2.0 m dimensions. At the entrance boundary corresponding to the first row of cells, it was assumed that the bottom was not eroded, which matched the entrance section of the tray reinforced with a cement crust in the laboratory experiment.

In the first step, the water flow rate at the inlet boundary was set to Q = 112 l/s, and flow was set as clarified (S = 0). The parameters of the numerical model, α_1 and β, were selected during the calculation process to achieve the best agreement in the average diameter of the flume between the calculated profile of the eroded slope and the experimental data. The results of the calculations are presented in Figure 3, Figure 4 and Figure 5.

Figure 3, Figure 4 and Figure 5 illustrate the different scenarios and parameters used in the calculations and their impact on the flattening or leveling of the slope. The results demonstrate the importance of selecting appropriate parameters for achieving the best agreement with experimental data.

The interaction between the flow and the moving bottom in the zone of a quarry, which typically occupies only part of the river's width, has been studied in detail. In real conditions, the interaction of the flow and the moving bottom in the quarry area is spatial, making the physical picture of the processes more complex and not fully describable in a one-dimensional setting.

The results of model studies conducted at the State Hydrological Institute, Moscow, Russia, in a hydraulic tray (2.0 m × 50.0 m) with a moving bottom were utilized to test the models. The bottom slope in the tray was 0.00125. At the flume inlet, the flow rate was set to Q = 40 l/s. During the establishment of the flow, a stable ridge relief was formed, which acquired a three-dimensional structure and scaly character. The average water depth increased from 6.65 cm to 8.86 cm, and the average flow speed decreased from 30.1 cm/s to 22.5 cm/s.

After achieving a quasi-uniform regime in the hydraulic flume, a riverbed quarry with a width of 0.7 m, which occupied 35% of the flume's width, was mined. The quarry was located in the axial part of the tray, with a depth of about 14 cm initially. The thickness of the sand layer at the bottom of the quarry was 3.0 cm. The sediment supply at the entrance boundary was 1.82 l/h (S = 0.02 kg/m³), with an average diameter of the initial soil being 0.32 mm.

The experiment in the work [8] described the process of bottom transformation, with background erosion deformations developing in areas above and below the quarry. Above the quarry, erosion began in the axial part of the flow, while below the quarry, the erosion process occurred more intensively. The highest bottom marks were observed in the edge areas near the upper ledge, which moved downstream.

For the numerical simulation of these experimental studies, a rectangular grid of 5250 cells was constructed, and the roughness coefficient was set to 0.027 to match steady-state depths and fluid velocities to the experimental data without a quarry. The calculations were then carried out in the presence of a quarry. The diffusion coefficient was determined from the bottom concentration (12) and the non-shear velocity using formula (9) with β = 15 and α_1 = 1.0. The comparison of results between calculations and experiments showed almost complete qualitative and quantitative agreement.

This extensive analysis and comparison demonstrate the model's ability to accurately represent and predict the behavior of the flow and bottom interactions in complex real-world scenarios.

4. Conclusions

This study has developed a comprehensive mathematical model to depict riverbed deformation behavior accurately. The conclusions of the study can be summarized as follows:

• The proposed mathematical model (Section 3.1) and numerical model (Section 3.2) for bottom deformations in uneven and unsteady river flows are suitable for calculating deformable channels, making it applicable to scenarios such as the easily eroded bed of the Amu Darya River. The models are characterized by their simplicity, with a minimum number of equations and empirical constants.
• A two-dimensional mathematical and numerical models of a deformable channel were developed and successfully verified.
• Numerical studies were conducted to investigate channel deformation, providing insights into the nature and intensity of these deformations. Comparison with experimental data demonstrated good agreement in the experimental results. However, it's worth noting that the optimal values of empirical coefficients varied between different experiments, highlighting the need for further work to establish more universal closing relations for the model.
• Verifying the model using field data is crucial to enhance the model's applicability and reliability. Additionally, the extension of the model to account for hydraulically heterogeneous soils is a promising avenue, as the concentration transfer equation (equation 4) allows for independent calculations of different particle fractions, assuming the total concentration of particles in the flow is not excessively high.