## Abstract

In this research, an explicit finite difference two steps scheme developed by Richtmyer is presented for simulation of 1D steady/unsteady flow and bed morphology in alluvial channels. Some laboratory experiments were performed in a mobile-bed flume for both of steady an unsteady flow to validate the proposed model. For improving the simulation, first, the role of sediment transport formulas, coupled–uncoupled approaches and simplification in the mass continuity equation are investigated. Comparisons between experimental data and model performance highlight the advantage of coupled method over the uncoupled method, domain role of sediment transport model in the bed morphological simulation and inaccuracy of aggradation processes due to simplifying the mixture continuity equation. Second, the importance of changing alluvial roughness is established by testing the model with calibrated and optimized friction factors. The inverse problem of estimating alluvial flow roughness is solved using a genetic algorithm (GA) optimization model coupled with alluvial flow model. The study demonstrates that the application of GA in the search for optimal values of roughness coefficients can significantly reduce computational errors and improve the computed water stage hydrographs. Role of the choice of the objective function is also mentioned.

## HIGHLIGHTS

An explicit finite difference two steps scheme developed by Richtmyer is presented for simulation of 1D steady/unsteady flow and bed morphology in alluvial channels.

Few laboratory experiments were performed in a mobile-bed flume for both of steady and unsteady flow to validate the proposed model.

The inverse problem of estimating alluvial flow roughness is solved using a genetic algorithm (GA).

### Graphical Abstract

## INTRODUCTION

Unsteady flow numerical simulation over mobile beds is theoretically challenging and practically important in the field of river engineering. So far, a large number of numerical river modeling methodologies of water-sediment-morphology fluvial systems have been developed (Van & Chua 2020). These methodologies have been applied as the primary methods in river engineering studies (Liao *et al.* 2014). The one-dimensional (1D) modeling of the free surface flow and bed evolution in alluvial channels are most often performed via applying the St. Venant equations (De St. Venant 1871) with the Exner equation (Graf 1971). In recent years, a great variety of numerical techniques have been applied in river engineering issues considering two basic aspects (Chen & Duan 2008). The first is associated with the numerical solution procedure of the governing equations system. Alluvial flows on erodible beds are different from streams on a fixed bed, because the stream may remove sediments from the bed or, on the contrary, sediments carried by the stream may deposit on the bed, which usually leads to river bed degradation or aggradation (Liang & Marche 2009; Daneshfaraz *et al.* 2021). It is supposed that the rate of bed morphological evolution has a lower order of magnitude compared to flow variations with low sediment concentrations (Estakhr *et al.* 2022). Therefore, the equation of water-sediment continuity is assumed to be identical for a fixed-bed clear water flow, regardless of the mobility of the channel bed, i.e. or the third term in Equation (1) is relinquished (Needham & Hey 1991; Sieben 1999; Sieben *et al.* 2009). As it will be described in the mathematical expressions in the next section, the water-sediment-morphology fluvial system is strongly coupled with the auxiliary closure relations. In available numerical approaches, these equations are often solved by asynchronous methods. In particular, at a given time step, the equations of mixture continuity and momentum are solved by supposing a small bed variation rate. In the next step, via the newly obtained flow variables, the sediment continuity equation is solved. Methods that contain asynchronous solutions are commonly referred to as uncoupled. Models involving the asynchronous solution are usually termed uncoupled. There are some-semi coupled methods that use iterative solutions of the bed and flow equations at a given time step (Palumbo *et al.* 2008). On the other hand, there is not a universal consensus about the best approach for morphological river modeling. The second aspect is the identification problem of alluvial channel friction factors. Mathematical modeling of open-channel flow systems involves wave propagation over channel reaches which is governed by the St. Venant equations, which describe the conservation of mass and momentum. The roughness coefficient used in the momentum equation is often considered as a lumped empirical coefficient that indicates the overall flow resistance of specific river channels. The literature on identifying and estimating parameters for unsteady open channel flow is sparse. Roux and Dartus investigated uncertainty in open-channel inverse problems using inundation observations (Roux & Dartus 2008). Sellier (2008) and Gessese *et al.* (2011) extended the inverse problem of reconstruction of river bed topography.

In the present research, the equations of 1D unsteady flow (i.e. St. Venant equations) and sediment continuity equation were solved numerically using the second-order accurate, explicit finite difference two steps scheme (Richtmyer & Morton 1967). The governing equations were solved in both coupled/uncoupled solution procedures. By evaluating the effect of bed load transport formulas and simplification of the continuity equation on simulation of alluvial flow, the model was validated with new laboratory experimental data. Then, the inverse problem of evaluating the open channel flow roughness was investigated using a genetic algorithm optimization model. Flow depths and discharges experimental datasets at several locations and time steps were used as inputs of the optimization model. In order to dynamically identify the roughness parameter the genetic algorithm and the unsteady alluvial flow model were coupled by time interval. Calibration was repeated with the selection of different objective functions.

### Governing equations

*t*= time;

*x*= stream wise coordinate;

*h*= flow depth; bed slope;

*u.h*= water discharge for unit width; bed load discharge for unit width;

*u*= cross-section- averaged stream wise velocity;

*z*= bed elevation;

*g*= gravitational acceleration;

*S*

_{f}*=*friction slope; and

*p*= porosity of the bed layer.

*n*denotes the Manning roughness coefficient, and stands for the hydraulic radius. The Equation (6) is used here for unit sediment discharge:where

*δ*= coefficients which depend on sediment characteristics. Equations (2)–(4) are a set of nonlinear hyperbolic relationships, and closed-form solutions can be only used for idealized cases. Accordingly, they may be solved using the numerical solutions. The modified Lax-Wendroff by Richtmyer and Morton scheme is a combination of the Lax-Friedrichs scheme and a midpoint leapfrog scheme, with each of the two steps being applied at half the time interval consecutively. The Richtmyer scheme has a second degree of accuracy both in space and in time. Figure 1 shows the solution grid for the two steps and boundary conditions. The system used in this study is a conservative one. By applying the above form in the system and considering the source conditions, the following steps can be expressed:

**Step 1:**

*et al.*(1981) can be applied. and parameters can be defined as below:

### Boundary conditions

*h*,

*u,*and

*z*values at nodes 1 and k + 1, one downstream and two upstream boundary conditions were applied. The inlet hydrograph of the flood was considered as the upstream boundary condition. Via the characteristics method, other dependent parameters which had not been determined by boundary conditions could be calculated. The downstream boundary condition was discharge-stage curves (Roushangar

*et al.*2010). The inclusion of the second boundary was not as straightforward as the former, and it was necessary to translate it into a relationship that the bed elevation could be computed at the upstream end. This was obtained by supposing a fictitious upstream node from node 1 and determining the sediment discharge at that point as . By applying the sediment continuity equation and the backward F.D. (Anderson

*et al.*2006) on the spatial differential part, the left side of Equation (20) will be obtained; so, the

*z*values at the unknown time level could be calculated:

**Stability**: From the stability point of view, the Richtmyer scheme should meet the Courant-Friedrichs-Lewy (CFL) requirement (Courant

*et al.*1928), which is given for a rectangular section as:

## MATERIALS AND METHODS

### Laboratory tests

### Aggradation under steady state flow and over loading

After establishing an equilibrium state, the stream flow which was capable to move the sediment was observed in the channel. Then, sediment was fed at upstream boundary, with a much higher rate than the flow transport capacity. The water supply rate was kept constant at 68 lit/sec, and the sediment supply at the flume entrance was 0.008 kg/s. The oversupply of sediment resulted in the sediment deposition and the aggradation of the stream bed near the flume entrance. The initial bed slope was S_{x} = 0.0005. The bed and water level profiles were recorded by traversing the bed and water level sensors along the flume length at various time intervals. The initial conditions were determined based on the water and bed level profiles measured at t = 0. The constant flow rate, bed level measured at the upstream, and water depth measured at downstream were considered as boundary conditions. The flume bed was covered with dunes during the experiments. The aggradation process occurred mainly in the first two mod channel length (with respect to the flume entrance). After passing this region, the movement of the dunes occurred, however, without significant change in the mean bed level. A better perspective of these features could be seen in Figure 2(b). The simulations of the proposed model were carried out to predict the average bed level and water surface profile, and good agreement was observed between the measurements and simulations.

### Degradation under steady state flow

_{x}= 0.0002. The sediment transport capacity of the flow allowed the immediate scouring after the false floor. The procedure of the simulations was generally the same as the modeling of the aggradation phenomenon described above. Again, the bed and water levels were monitored at different time intervals. Figure 3 shows the profiles of the water and bed levels after an elapsed time of 2 h. It can be seen from the figure that the simulations of the bed and water levels are predicted with higher accuracy.

### Aggradation and degradation under unsteady state flow

_{50}= 0.0011 m, m = 0.0168 and the constant sediment–supply rate of 0.02 kg/s at the flume entrance. Figure 5(a)–5(d) show the observed and simulated water and bed level profiles at different time intervals. The water and bed level profiles for t = 0 min were considered as initial conditions (see Figure 5(a)). Boundary conditions were the flow depth downstream and the sediment input and flow rate upstream.

### Asynchronous–synchronous solution procedure

In the coupled system of governing equations, selecting the appropriate numerical solution approach is so important. The water sediment-morphology fluvial system was strongly coupled, as mentioned in Equations (2)–(4) along with Equations (5) and (6). In existing numerical approaches, the mentioned equations are solved using an asynchronous procedure. In particular, at a given time step, at first, the mixture continuity and momentum equations are solved, supposing inconsiderable bed variation (or fixed bed). The sediment continuity equation is then solved through the newly achieved flow variables. Methods that include the asynchronous solution are commonly referred to as uncoupled (Cao & Carling 2002). The basis of the asynchronous solution strategy is the ‘fixed bed’ and ‘quasi-steady’ flow assumptions and often the flow is supposed to be steady when the channel bed evolution is investigated. On the other hand, the river bed is implicitly supposed to be fixed within a time step, whereas the flow over a mobile bed is of primary interest.

### Role of simplification (importance of )

Flow state . | Solution system . | |||||
---|---|---|---|---|---|---|

Uncoupled . | Coupled . | |||||

With () . | Without () . | |||||

%RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | |

Steady (aggradation) | 1.5 | 19 | 0.9 | 13 | 1.1 | 13.9 |

Steady (degradation) | 3.7 | 13 | 2.3 | 8 | 2.33 | 8.1 |

Unsteady (T = 1.5 h) | 3.9 | 18 | 4 | 18.4 | 5 | 20.2 |

Unsteady (T = 2 h) | 5.96 | 18.8 | 4.82 | 15.7 | 5.8 | 17.9 |

Unsteady (T = 4 h) | 6 | 30 | 2.55 | 23 | 3.4 | 25 |

Flow state . | Solution system . | |||||
---|---|---|---|---|---|---|

Uncoupled . | Coupled . | |||||

With () . | Without () . | |||||

%RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | |

Steady (aggradation) | 1.5 | 19 | 0.9 | 13 | 1.1 | 13.9 |

Steady (degradation) | 3.7 | 13 | 2.3 | 8 | 2.33 | 8.1 |

Unsteady (T = 1.5 h) | 3.9 | 18 | 4 | 18.4 | 5 | 20.2 |

Unsteady (T = 2 h) | 5.96 | 18.8 | 4.82 | 15.7 | 5.8 | 17.9 |

Unsteady (T = 4 h) | 6 | 30 | 2.55 | 23 | 3.4 | 25 |

Unsteady flow . | Solution system . | |||
---|---|---|---|---|

Uncoupled . | Coupled . | |||

%RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | |

T = 1.5 h | 2.3 | 15 | 2.1 | 14 |

T = 2 h | 3.4 | 12.2 | 2.5 | 10 |

T = 4 h | 2.6 | 21.8 | 1.33 | 18 |

Unsteady flow . | Solution system . | |||
---|---|---|---|---|

Uncoupled . | Coupled . | |||

%RAE_{1}
. | %RAE_{2}
. | %RAE_{1}
. | %RAE_{2}
. | |

T = 1.5 h | 2.3 | 15 | 2.1 | 14 |

T = 2 h | 3.4 | 12.2 | 2.5 | 10 |

T = 4 h | 2.6 | 21.8 | 1.33 | 18 |

### Bed load transport formulas

Bed load discharge is the main parameter in specifying the morphologic development of alluvial channel reaches. An essential task with the application of 1D numerical model is the use of an appropriate bed load discharge equation. Since most of the existing relationships depend on limited data and untested assumptions, the application of most relationships is limited to their development conditions. In this study, considering the validity of bed load transport formulas, Meyer-Peter & Müller (1948), Einstein (1950), and Ackers & White (1973) were applied for steady flow. The formulas mentioned above along with Song & Graf (1997) were used for unsteady flow.

### Inverse problem of channel friction factor optimization

*et al.*2009) the initial roughness coefficient was obtained 0.018. The range of

*n*variations was between 0.018 and 0.024. The roughness identification procedure is shown in Figure 6. Five objective functions (

*f*,

_{1}*f*,

_{2}*f*,

_{3}*f*,

_{4}*f*) were chosen. The objective functions were selected based on the comparison of observed and simulated values. If the value of the function is set higher than a prescribed tolerance value, this process continues iteratively through the correction calculation for the variable through the optimization algorithm. The objective function (Equation (24)) was the sum of squares of the difference between the observed and simulated values of flow discharge or depth and was applied to minimize the error criterion.where

_{5}*f*: error criterion or objective function, and : depth or discharge values at several stations.

### Optimization model: genetic algorithm (GA)

*et al.*1995). The most important property of GAs is their strong nature and the balance between efficiency and efficacy required to survive in different environments. It has been proved that GAs are successful in providing robust search in complex search spaces (Goldberg 2000). GA used here is only one of the variants of the method. For instance, the selection probabilities may be related directly to the function values instead of the ranking method used here. Further details on genetic algorithms can be found in Goldberg (Goldberg 2000) and Holland (Holland 1975). In this research, five objective functions (

*f*,

_{1}*f*,

_{2}*f*,

_{3}*f*,

_{4}*f*) were chosen for analysis:where and denote the observed flow depth and discharge values, respectively, and

_{5}*h*and

*Q*stand for the corresponding simulated values. Error percentages related to the optimized Manning coefficient were defined as .

*f*on Error percentages are presented in Figure 7.

_{1}–f_{5}## RESULTS AND DISCUSSION

There is an ongoing debate about using coupled or uncoupled approaches in alluvial flows. For most of the practical cases, the application of uncoupled methods may be justified, mainly due to different time scales in the motion of water and sediment, as well as the inclusion of the empirical relations with different degrees of approximations, which largely affect the accuracy of the outcomes. The coupled and uncoupled approaches were tested for steady and unsteady flow conditions. For unsteady flow, the uncoupled approach at t = 30 min and t = 90 min gave a close solution to the coupled method, however, at t = 120 min some oscillations were observed.

Strong oscillations of false bed and water level took place around t = 240 min. In the steady flow, at t = 120 min, it was found that the asynchronous solution procedure entailed significant inaccuracy for both aggradation and degradation processes. In the uncoupled procedure, the rapid changes of the bed level related to the sediment component were not introduced simultaneously to the flow component. Since the flow and sediment computations were done separately, the flow component received the new bed levels after *Δ*t. The obtained flow conditions (new water depths and velocities) related to the sudden variations of bed level significantly differed from the conditions extracted from the coupled method, which included a simultaneous solution of both water and sediment components at the previous time step. This deviation of the outputs was magnified through uncoupling of the components at the subsequent time steps, resulting in larger artificial oscillation. Conversely, in the coupled method, slight fluctuations of the bed level were gradually reduced due to simultaneous interference of the flow conditions.

In automatic optimization, a systematic iterative process is used to search a set of parameter values for which a chosen function assumes its extreme valves. The model implemented is an explicit finite difference scheme of Richtmyer and the process is automated using a genetic algorithm for estimating Manning's roughness coefficients using flow measurements data in a laboratory flume. The comparison between the water and bed level profiles (Figure 5(b)–5(d) and Table 2) monitored at various instants before/after optimization shows that the optimized model reduces computational errors, improves the simulation results, and performs satisfactorily in all cases (the evaluation herein is limited to subcritical flow conditions).

**)**are selected. Figure 9 presents a comparison between observed and computed stages before/after optimization. From the figures, it is seen that this model with optimized roughness can significantly improve the simulations of the water stage hydrograph. The computational effectiveness shown in Table 2 illustrates that the technique can reduce computational errors of stage simulation at each station.

In an unsteady flow, by optimization of the roughness coefficient, the form of the objective function can be elated. Both the optimization procedure and optimization criterion type have an effect on the duration and accuracy of computations. Also, if the selection process of these parameters has not been done properly, they may lead to incorrect values. In fact, the values, type, and quality of datasets play an important role in outcomes and often specify the selection of objective functions. The results of the performance of the objective functions, defined by *f _{1}–f_{5}* on error percentages (see Figure 7), show that the objective function significantly affects the quality of the calculations. The best outcomes are determined via functions based on relative errors (

*f*and

_{4}*f*). The impacts of the traditional

_{5}*f*and

_{1}*f*are relatively poor. The difference between

_{2}*f*and

_{1}*f*indicates that the objective function is particularly significant when both water stage and discharge are measured.

_{4}## CONCLUSIONS

In this research, an explicit finite difference two steps scheme developed by Richtmyer is presented for the simulation of 1D steady/unsteady flow and bed morphology in alluvial channels. Some new laboratory experiments were performed in a mobile-bed flume for both steady and unsteady flow to validate the proposed model. For improving the simulation, first, the role of the sediment transport formula, coupled/uncoupled approaches, and simplification in the mass continuity equation is investigated. Comparisons between experimental data and model performance highlight the advantage of the coupled method over the uncoupled method, the domain role of the sediment transport model in the bed morphological simulation, and the inaccuracy of aggradation processes due to simplifying the mixture continuity equation. Second, the importance of changing alluvial roughness is established by testing the model with calibrated and optimized friction factors. The inverse problem of estimating alluvial flow roughness is solved using a GA optimization model coupled with an alluvial flow model. The study demonstrates that the use of GA in the process of searching for optimal values of roughness coefficients can significantly reduce computational errors and improve the computed water stage hydrographs. The results of the performance of the objective functions show that the form of the objective function influences the quality of calculations and the best results are obtained by functions based on relative errors.

## FUNDING

The authors did not receive support from any organization for the submitted work.

## AUTHORS’ CONTRIBUTIONS

Kiyoumars Roushangar: Conceptualization, supervision, methodology, review & editing; Roghayeh Ghasempour: Project administration, investigation, data curation, methodology, writing; Hazi Mohammad Azamathulla: Conceptualization, methodology, review & editing.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

154(3), 207–219