Hydromorphological Modeling of Vertical Grain Sorting Process. Insights from the Günter Experiment 6 using TELEMAC & SISYPHE

I. INTRODUCTION

The main objective of river engineering is to predict river hydrodynamics and bed-level changes. The principal parameters of this field are discharge, cross-sectional area, water depth, roughness, and grain size. These parameters used in the hydrodynamic model predict the dynamics of rivers. Such river dynamics are generally characterized by roughness and grain size. To understand a complete hydromorphological process, the hydrodynamic models are coupled with the sediment flux and relative bed level changes .

Generally, sediment transport relation is described as the transport capacity of a river cross-section at a specific discharge. Hence, a well-calibrated hydromorphological model with particular geometry and bed grain size describes the river location and whether it is aggrading or degrading. A river at its origin has abundantly graveled size bed materials, and in the transport processes, the bedload is a dominant transport mode. As the river follows the valley's slope, the bed grain size decreases to clayey. At this location, suspended load is a dominant transport mode [2].

The study on river transport capacity started 70 years ago when Meyer-Peter and Mueller 1948 published a paper. They considered a uniform sediment model with one mean grain size [3].

However, the bed loads are widely graded, and Mueller already in 1943 noticed grain sorting effects. The Meyer-Peter and Mueller formula was still adopted as there was no other study on these effects [1].

Soon, the study on sediments and river engineering boomed with Harrison's study in

1950, which clearly showed the grain sorting effects. He described a static and mobile armor layer forming under different flow conditions. Both armors were the outcomes of bed level change. The static armor prevented the bed degradation, and there was a low sediment supply [4]. After Harrison, Gessler 1968 studied the static armor closely and developed an algorithm that could predict the armor layer's grain-size distribution. For this, he used the bed shear stress of the original bed material [5], and in 1971, he published a stability criterion for the armor layer [6]. After him, Günter 1971 performed a series of experiments and typically focused on stability criteria to define a threshold value for bed erosion [7], which raised questions on motion threshold and grain sorting processes. However, Einstein in 1950 had already elaborated the hiding and exposure function of mixed grain size [8], and Egiazaroff in 1971 formulated the hiding and exposure function of mixed grains on a hydrodynamic model [9], which simplified and answered the threshold of motion and transport rate. After that, Parker in 1990 and Sutherland in 1992 studied the hiding and exposure functions [10, 11]. Similarly, in 1998, Powell studied the patterns and process of sediment sorting in the graveled riverbed and reviewed the linkages between the pattern and sediment transport process [12]. In 2014, Jia and Zhong used the weighted implicit discretization method to simulate scouring and armoring in alluvial rivers [13]. Zhang et.al., 2020 used Laboratory method with two sets of flumes to analyze the armor layer and found out that statistics on bed armoring are cannot completely define erosions of the riverbed, however, recommended that other flow parameters must be simulated to completely understand the fluctuation in the bedload transport rate [17].

Like others, Gergely et. al., in 2017 used the 3D numerical flow model software called 'BASEGRAIN' coupled with combined application of two bedload transport (Wilcock and Crowe, Van Rijn) formulas and validated their model with the laboratory experiment. They argued that the numerical flow model coupled with combined application of bedload transport formulas improved the results of the numerical simulations . Similarly, Lee and Ahn in 2023 used HEC-RAS 1D to analyze the bed sorting process of the Geum River, Korea in two and three layers. They stated that their model performance on two layers of bed was reliable whereas for the higher layer the model was very unreliable and recommended sufficient investigations on the model .

Grain Sorting Process: the grain size distribution in the unconsolidated riverbed, which indicates the energy, rate, and duration of the deposition and transport process . Sediment transport theories state that critical shear stress for incipient motion depends on grain size, which applies only to uniform-sized materials. However, a riverbed contains a wide range of mixture classes where smaller particles hide under the larger particles, severely affecting transport processes. The studies on armoring tests suggest different grain mobilities. Günter, in his experiments, never used sediments as incoming flux. They instead started with conditions that resulted in erosion. This results in erosion upstream while maintaining a constant water level downstream. The result showed the sorted gain size, which implied that the ratio of smaller grains to larger grains was significantly high. So, a sorting process can be described as the variation of the coarsening of the bed surface.

Static Armor: Like in Günter's experiments, if sediments are not supplied at the inlet, a stable armor may form at the bed layer, which does not allow the erosion of the riverbed [5]. If the grain size is finer than the bed materials supplied, all the incoming materials will pass out without affecting the bed level. This does not allow the coarsening of the bed surface. Hence, a static armor layer formed on the surface of the bed layer prevents erosion and thus stops the bed level change.

Mobile Armor: A static armor layer cannot form when shear stress is very high. In such a case, the grain sorting process makes the top layer coarser, the bottom layer finer, and the layer is formed by the mobile armor layer . With this new layer, riverbeds can only be stabilized with the supply of sediments. Hence, a constant supply of sediment flux is required at the inlet.

Günter Experiments: Günter, in 1971, performed a series of experiments to define bed stability without supplying loads at the inlet and constant water at the outlet. The length of the rectangular flume was 40 meters, and the breath was 1 meter. He performed ten different experiments with varying discharge and bed slopes. As the slope was steep and there was no sediment flux at the inlet, the bed materials were eroded, and the bed level changed. Since constant water depth at the outlet was performed, no erosion occurred at the outlet, and most of the erosion took around the inlet.

When the slope was reduced, the bed shear stress was also reduced, and after many days of observations, the final stage was obtained as the threshold for the shear stress and the stability of the static armor layer. In his experiments, vertical grain shorting was observed. The sorted bed layers were described by the final coarser material at the top layer, which was originally finer where the erosion had flushed the smaller grain-sized material from the reach.

In this way, Günter used the physical model and measured the outcomes thoroughly. Meanwhile many researchers have tried to explain the vertical grain sorting process using numerical models. Some used 1D models and concluded that their model was not sufficient to define the vertical grain sorting process whereas some used 3D flow model coupled with combined bedload formulas which is usually time consuming.

So, this study tries to use a 2D numerical model based on the Günter experiment 6 to see how the parameters behave numerically. For that a Numerical Hydromorphological Model (NHM) was developed in TELEMAC coupled with SISYPHE. Then, a series of tests were conducted to find the effects of numerical parameters on transport processes. The model was then calibrated, and the armor layer was analyzed.

II. MATERIALS AND METHODS

Günter Experiment number 6 is formed in TELEMAC-2D coupled with SISYPHE. The length of the study domain is 40 meters, and breadth is 1 meter. The initial slope of the bed is 0.225% with 7-grain classes of bed materials. The discharge at the inlet is 48 lps, and the constant water level at the outlet is 0.0913 meters. The total simulation time is 691 hours [7].

The National Hydraulics and Environment Laboratory of the French Electricity Board developed TELEMAC-2D. It simulates surface flow in two dimensions using finite volume or finite element method in unstructured mesh. To simulate the flow and transport processes, a steering file containing the configuration of computation, a geometry file of geometrical mesh, a boundary conditions file describing the type of boundaries, an initial state file of computation, a bottom topography file describing the elevation of the riverbed, and the FORTRAN file which contains the specific programming are compiled in one folder.

Bluekenue is used to arrange all required files in a folder. Bluekenue is one of the many processors associated with TELEMAC. The experiment setup and mesh were generated in Bluekenue at a 0.35-meter nodal distance with some added hard points. Then, the slope was added to the mesh. A 2-D interpolation was used to interpolate the applied slope to each mesh node, and the default 'maximum distance' interpolator was used. The geometry file must be a SELAFIN object with a variable as the bottom slope so the TELEMAC recognizes it as the bathymetry for the project.

After this, the boundary conditions file is created in ASCII format. Bluekenue makes all boundaries solid, and to define an open boundary, discharge and water level data are required at the inlet and outlet, respectively. At the inlet, the 'open boundary with prescribed Q' and at the outlet, 'Open boundary with prescribed H' is used with the tracer code 'open boundary with prescribed tracer' applied to both ends. The boundary code for 'prescribed Q' is 455, and for 'prescribed H' is 544.

In addition, the number of grain classes, including their portion, were added in the FORTRAN file, and the simulation variables were described in the steering file. This file couples TELEMAC -2D with other files like FORTRAN, SISYPHE, etc. The other files, like boundary conditions files and geometry files, are also defined in the steering file. To generate the initial flow condition, the time step and the number of time steps were initially adopted 0.5 and 1000. Then the initial flow conditions were added to the initial condition file. From the experiment 6, the prescribed flow rate at the outlet and inlet of and , were adopted with the specified elevations at the outlet and inlet being 0.0813 and 0, respectively. Strickler's roughness coefficients are defined at bed and banks 50 and 100.

Where is Manning's roughness, and is the mean grain size.

TELEMAC-2D is then simulated on the Leibniz-Rechenzentrum der Bayerischen Akademie der Wissenschaften (LRZ) supercomputing computer at the Technical University of Munich, Germany. A short simulation was made to generate the initial condition file. Then, output SELAFIN files were animated on Bluekenue, and the stability of the model was checked. The stable output file is then used as an initial conditions file with the command "COMPUTATION CONTINUED = TRUE," previously applied initial conditions were commented out. Considering the stability of the numerical model, the time step of 0.5 seconds is used for the whole simulation period.

On the other hand, sediment data were defined on the SISYPHE file. The mean diameter of grain classes 1-7 sediments are 0.00102 m, 0.0020 m, 0.0031 m, 0.0041 m, 0.0052 m, 0.006 m, and 0.3 m respectively [10]. The active layer thickness is 0.06 m, and the number of bed layers is 3. Hunziker's bed load transport formula adapts the Meyer-Peter and Muller formula. The same friction coefficient is used on both steering files. The shield parameter for each grain class used is 0.047.

Hunziker's volumetric sediment transport per sediment class is:

This formula is only valid when . Is Shields parameter and is corrected Shields parameter. TELEMAC computes using the Critical Shields parameter, mean grain size of the surface layer, and the underlying layer. The formula also adopts a hiding function and is defined by . The volumetric sediment transport is computed for each grain class 'i'.

Finally, the SISYPHE file is coupled with the TELEMAC file for the first run. The desired simulation output can be achieved by changing the printout period. For the simulation, 12 hours as a printout period, which is 86400 timestep, is used. The output files are animated in Bluekenue, and the corresponding final water level, grain size, and slopes are compared to the original document of Günter's experiment. Shields coefficients and bed roughness are adjusted to calibrate the model, and the output data are compared with experimental data.

III. RESULTS

There are many parameters associated with hydromorphological modeling, and to find an exact relation between measured and simulated data, each parameter must be analyzed. Hunzikers' sediment transport model was initially analyzed at a separate Shield's coefficient, keeping bed roughness constant at Strickler's 50. The change in grain sorting and water depth is shown in Figure 1.

{"image_source":{"path":"images/e3977347fc1e3206c2eabbb8493e140062d4c5d22c9ee7e2bc6387fb1473fa47.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 1: This is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and Corresponding Water Depth at Constant Strickler's 50 and Shield's coefficients of 0.037, 0.04, and 0.043"}],"chart_footnote":[]}

After that, the reverse of the previous was analyzed. Shield's coefficient is kept constant at 0.040, and three analyses on Strickler's coefficients of 30, 40, and 50 were made. The plot of grain sorting and water depth is shown in Figure 2.

Observing the sensitivity of both Shields' and Strickler's coefficient,, both parameters are modified simultaneously to observe the effect on the river morphology. The impact on grain sorting and water depth is shown in Figure 3.

The thickness of the active layer was not completely defined in Günter's experiment, and as an initial guess, a thickness of 0.06 meters was adopted. Then, the different values of active layer thickness were analyzed at varying Shield's and Strickler's coefficients. The plot is shown in Figure 4.

{"image_source":{"path":"images/6991dae243469db5b298bfc091f5b1c0c721ec9733fff2ef67bbb0979e77eefa.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 2: This Is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and Corresponding Water Depth at Constant Shield's coefficient of 0.040 and Strickler's coefficientsat 30, 40, and 50"}],"chart_footnote":[]} {"image_source":{"path":"images/bae1a361210aa033d9c2224b641b5b321b9be05df83fe6ea73657e70b2cd3c2a.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 4: This Is a Plot of the Final Percentage of Grain Class in Different Thicknesses of the Active Layer of the Riverbed and the Corresponding Water Depth"}],"chart_footnote":[]}

{"image_source":{"path":"images/516519563c2ab1d8ff180e7dde919c7ef09b87d8986e83c7debe99c6a70fd660.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 3: This Is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and Corresponding Water Depth at Different Combinations of Shield's coefficient and Strickler's coefficients"}],"chart_footnote":[]}

Even after analyzing different parameters, some hidden parameters, like hiding function coefficients, could significantly impact sediment transport. So, a study on default value and one higher value was carried out, keeping other parameters constant. The plot of the analysis is shown in Figure 5.

{"image_source":{"path":"images/c7a881f9c6180216e2477cdd8f6c084410960deb6f9dc64af9614b5e8291f700.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 5: This Is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and Corresponding Water Depth at Two Different Hiding Function Constants of 6 and 5.5 Is a Default Value"}],"chart_footnote":[]}

Besides, there are many transport formulas associated with TELEMAC and SISYPHE. Various simulations were conducted to analyze them at constant Shield's and Strickler's coefficients. The result of the analysis on grain sorting and water level change is shown in Figure 6.

{"image_source":{"path":"images/54bf2a0cc843bd8d94debc086f762e77eb37a7f75b0be2adbd3516fa83f1f7c8.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 6: This Is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and Corresponding Water Depth at Different Bed Load Transport Formulae"}],"chart_footnote":[]} To achieve the water level as measured by Günter in his physical model, it is also necessary to analyze the model with different initial water levels. To do so, two water levels were analyzed, i.e., 8.13 cm and 7.13 cm. The results are plotted and shown in Figure 7.

Finding a perfect agreement between measured and simulated data is always tricky. The required water depth was achieved using Hunziker's bedload transport formula at Shield's coefficient of 0.047 and Strickler's roughness coefficient of 61 in thickness of active layer 2.5 mm. Still, the grain sorting did not match the physical measurements made by Günter. The evolution of bed grain, slope, and water depth is shown in Figure 8.

{"image_source":{"path":"images/4acc2eb0d50ead2c80d173e593db6c0959029d8d3bf5382af5ea1a783c3e7e39.jpg"},"content":"","chart_caption":[],"chart_footnote":[]} {"image_source":{"path":"images/ebbc0f404cf8a4980d1f1d96b013a68f6b20ff4827ac3393493dec6d071aa73a.jpg"},"content":"","chart_caption":[{"type":"text","content":"Figure 7: This Is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and the Corresponding Water Depth at Two Different Initial Water Level Conditions. 9.13 Cm Is the Last Water Level Measured Experimentally by Günter"},{"type":"text","content":"Figure 8: This is a Plot of the Final Percentage of Grain Class in the Active Layer of the Riverbed and the Corresponding Bed Slope and Water Depth at Different Time Steps. A Dotted Line With Levels at the Primary Axis Shows the Slope in Percentage. The Primary Axis Also Represents the Total Portion of Grain Classes."}],"chart_footnote":[]}

VI. DISCUSSION

Bed roughness is one of the sensitive parameters in river hydrodynamics, and Shield's coefficient is another sensitive parameter in riverbed transport processes. These two parameters are essential to developing a hydromorphological model. To understand their behavior on the grain sorting and water depth, both were analyzed separately and jointly.

At a higher Shield's coefficient, keeping the Strickler's coefficient constant at 50, a portion of the class 1 grains increases while the portion of the other classes decreases. This showed that with an increase in 'θcr,' the portion of the fine grain class increases, while the coarser grain class decreases in the active layer. The water level also decreases from 11.40 cm to 10.70 cm (see Figure 1), which implies that higher bed shear stress leads to increased sediment transport, resulting in a lower water level.

Similarly, at a higher Strickler's coefficient, keeping Shield's coefficient constant at 0.040, a portion of class 1 grain decreases, and a portion of the other classes increases. This showed that, with the increase of roughness, the portion of the finer grain class decreases, and the amount of coarser grain class increases, while the water depth also decreases with increased bed roughness. This implies that higher bed roughness affects sediment transport and reduces water depth (see Figure 2).

However, varying results were obtained when both coefficients were analyzed simultaneously (see Figure 3). This might be the combined effect of changes in bed shear stress and bed roughness, which leads to complex and non-linear interactions, resulting in outcomes that are not directly predictable based on individual variations in each coefficient.

The result above implied that, by increasing both the bed roughness and the Shield's coefficient, the development of numerical simulations could get closer to the measured results. Combinations of various increasing values of both parameters were analyzed. The plotted story (see Figure 3) showed that Strickler's coefficient between 60 - 70 and Shield's coefficient between 0.045 - 0.048 could give better results. After that, increasing values of both parameters showed fluctuation in grain sorting and water level.

The joint analysis on Shield's coefficient and Strickler's coefficient showed an excellent progression in the calibration of the model; however, the measured results were still far from the simulations' results. So, the thickness of the active layer was analyzed (see Figure 4). The plotted results showed that by decreasing the thickness of the active layer up to , the finer grain class and the water level decrease, and a portion of the coarser grain class increases. But if the coating is further decreased, the opposite phenomenon occurs where the water level and a bit of fine grain class rise rapidly. Figure 4 also showed that the thickness of the active layer is the most sensitive parameter in this study.

Besides, other parameters like the hiding function's constant could change the grain sorting process. For this, a default value was analyzed with a higher value (see Figure 5), and the result was almost the same. With higher values, the water level and portion of finer grain decreased but not significantly.

In addition to these, many available bed load transport formulae are analyzed separately (see Figure 6). The result showed that the portion of finer grain was higher, and the water level was lower in Hunziker than in others. The result also showed that Van Rijn, MPM, and Einstein-Brown formulas allocated the maximum portion for grain class 2 in active layers. It is also seen that Van Rijn, MPM, and Einstein-Brown recipes gave higher bits for coarser grain class, e.g., class 5 and 6, than Hunziker.

After that, the initial water level applied to the model was also analyzed and plotted (see Figure 7). The plot showed that the final water level also increased, and the finer grain class 1 portion increased with the increasing water level.

Finally, analyzing all the parameters, formulas, and initial conditions, two combinations were implemented. First, the Van Rijn formula with Shield's coefficient of 0.06, Strickler's coefficient of 117, 25 mm thick active layer, and 8.13 cm initial water level. Second, the Hunziker formula with Shield's coefficient of 0.047, Strickler's constant of 61, 2.55 mm thick active layer, and 8.13 cm initial water level. The latter was more appropriate, and the complete grain sorting process result, water level, and bed slope change were plotted at different time steps (see Figure 8).

In 12 hours of simulations, the portion of grain class 1 is more than 50%; the slope is 0.225%, and the water level is 8.13 cm. This is closer to the initial condition; this implied that nothing special was observed within the first 12 hours. After a day, the portion of grain class 1 significantly decreased to 38% while the portion of the other class increased. The water level also increased to 8.9 cm, but the slope decreased to 0.19%. The maximum bed shear stress was 1.89 N/m . After two days of simulations, a portion of grain class 1 further reduced to 29% while the water level was observed at level 9.1 cm, and the maximum bed shear stress increased to .

Similarly, at the end of day, the grain sorting process continued by decreasing the portion of grain class 1 and increasing others. The water level also increased to 9.14 cm and the slope to 0.3%. Then after the grain sorting and bed level change process retarded and can be seen on the plot at time step 144 hours. At the end of 288 hours, the change in grain sorting, water level, and bed slope is insignificant compared to the previous time step, and if compared to the last time step, the grain distributions, water level, and slope are constants. The maximum bed shear stress after 12 days was . This implied that after 288 hours or 12 days, the static armor layer was formed at the active layer's surface, protecting the grain from erosion. These 12 days can also be considered a threshold for creating a static armor layer.

However, the simulation results are still far from Günter's measurements, and the water level only comes closer to experimental data. Many parameters are associated with model calibration, and not every parameter is defined numerically in reality. A few were studied, and researchers have applied them as a numerical reference for numerical modeling. Those studies were a simplified prototype version, and many assumptions were made to find a relation. Those assumptions may not work for every scenario. So, other numerical parameters that could define the transport behavior better are required. An artificial neural network (ANN) can be used as an option. But it's always better to use two or more numerical models to see how the models represent reality. Also, there could be human errors in measuring the experimental model. A good model calibration is impossible due to the lack of numerical parameters and measurement errors. For this, updated parameters and precisely taken measurements will be required.

V. CONCLUSIONS

There are many formulas and parameters associated with them in numerical hydromorphological modeling. Choosing the best combination among them is one of the main tasks in the model calibrations. It was observed that local sensitivity analysis is not enough to calibrate a model, so a global sensitivity analysis was required. After performing the simulations out of best combinations, i.e., Hunziker formula with Shield's coefficient 0.047, Strickler's coefficient of 61, 2.55 mm thick active layer, and 8.13 cm initial water level gave a better output. It was seen that it took some time to form a static armor layer on the surface of the active layer, and in this model, it took 12 days. After the armor layer was created, the riverbed erosion completely stopped, and so did the grain sorting process, where the maximum bed shear stress was 1.98 N/m2. The final portion of grain class from 1 to 6 is 26%, 28%, 18%, 11%, 13%, and 4%, respectively; the bed slope is 0.33%, and the water depth is 9.15 cm. However, these results are still far from experimental results. So, using two or more numerical models for the analysis would be better. To find a good relation between them, updated parameters, their range, and precisely taken measurements will be required.

ACKNOWLEDGEMENT

I would like to acknowledge Dr. Markus Reisenbüchler for his invaluable guidance, Shubhechha Paudyal for her unwavering support, and the Technical University of Munich, Germany, for providing access to their resources.