Advances in computational morphodynamics using the International River Interface Cooperative (iRIC) software

Results from computational morphodynamics modeling of coupled flow–bed–sediment systems are described for 10 applications as a review of recent advances in the field. Each of these applications is drawn from solvers included in the public‐domain International River Interface Cooperative (iRIC) software package. For mesoscale river features such as bars, predictions of alternate and higher mode river bars are shown for flows with equilibrium sediment supply and for a single case of oversupplied sediment. For microscale bed features such as bedforms, computational results are shown for the development and evolution of two‐dimensional bedforms using a simple closure‐based two‐dimensional model, for two‐ and three‐dimensional ripples and dunes using a three‐dimensional large‐eddy simulation flow model coupled to a physics‐based particle transport model, and for the development of bed streaks using a three‐dimensional unsteady Reynolds‐averaged Navier–Stokes solver with a simple sediment‐transport treatment. Finally, macroscale or channel evolution treatments are used to examine the temporal development of meandering channels, a failure model for cantilevered banks, the effect of bank vegetation on channel width, the development of channel networks in tidal systems, and the evolution of bedrock channels. In all examples, computational morphodynamics results from iRIC solvers compare well to observations of natural bed morphology. For each of the three scales investigated here, brief suggestions for future work and potential research directions are offered. © 2019 The Authors Earth Surface Processes and Landforms Published by John Wiley & Sons Ltd


Introduction
In this article, we present a suite of research results in computational morphodynamics carried out by the International River Interface Cooperative (iRIC) group. Most of the work here centers on riverine morphodynamics and we apply the somewhat typical, although perhaps arbitrary (as discussed later), division of the field into mesoscale, microscale and macroscale features, meaning bars, bedforms, and channels. The order chosen is based on the historical evolution of the field of computational morphodynamics, as the first successes of this sort of endeavor were in the field of bar evolution, followed by bedform initiation and development, and only recently moving to channel-scale morphodynamics. This order also reflects increasing physical and computational complexity along with an associated demand for computer time, memory, and efficiency. Nevertheless, all the computations given here have been achieved on typical laptop/desktop computers and the iRIC solvers are suitable for practical application in real rivers using such equipment.
One of the hallmarks of the iRIC group is a strong commitment to make the computational morphodynamics tools we develop available to students and practitioners in as timely a manner as possible; the work presented here is no exception. Most of the results presented in this article can be reproduced using software in the iRIC public-domain suite of tools and tutorials available at www.i-ric.org (see Nelson, Shimizu, et al. [2016b] for a list of current solvers and capabilities); for those few results that use unpublished solvers or modifications to existing solvers, be assured that they will be made publicly available.

Mesoscale: Bars
So-called mesoscale bed features in rivers are generally river bars; these features include point bars as are commonly seen in meandering rivers, but also incorporate the rich gamut of alternate bars, mid-channel bars, confluence bars, braid bars, separation bars, and so forth that are found in almost all rivers. Their important distinction relative to other sedimentary features is that the basic mechanics of bars can generally be reasonably well explained with a shallow-water flow model, albeit often with some simple treatment of curvature-driven secondary flows. In this section, we introduce two examples of the computational prediction of bar behavior, the first of which assumes equilibrium sediment supply and the second of which briefly examines a small part of the rich behavior that stems from sediment oversupply or supply limitation. Notably, all the computations shown here can be accomplished in just a few minutes on a personal computer (PC) computer, as the computational demand is low for features that can be treated with vertically averaged approaches.

Sand bars in simple flows
Sand bars are one of the most fundamental morphological signatures of rivers (Figures 1a and 1b) and have received a great deal of attention from a wide range of research communities because of their important role in hydraulics, sediment transport and sequestration, and riparian ecology. Modeling of sand bars is a classical problem and may well be the most well-studied and understood research topic in river morphodynamics. Simple, idealized riverine sand bars have typically been classified by their origin as either free or forced bars. Free bars are an intrinsic response of the system produced by the interaction of flow, sediment transport, and morphological change. Early work on the genesis and mechanics of free bars was carried out by applying linear stability analysis (e.g. Callander, 1969) to understand both the physics of wavelength selection and the range of conditions for bed instability. However, forced bars are produced by external forces on the system such as planimetric channel geometry (e.g. Zolezzi et al., 2012), including curvature effects and width variations. Although separation of bars into free and forced categories is appealingly simple, more recent work suggests that these are endpoints of a spectrum of bar behavior, and hybrid bars displaying characteristics of both free and forced bars are common (Duró et al., 2016). Several review papers regarding sand bar dynamics have been published (e.g. Nelson, 1990;Tubino et al., 1999;Seminara, 2010), so readers can refer to these articles to understand the basic physics of the problem, perspectives, and future challenges. For schematic diagrams of bar types and regime graphs providing bar characteristics as a function of physical parameters, the reader is referred to the excellent discussion in Kleinhans and Van den Berg (2011). Here, we focus on recent progress and applications of numerical modeling to the creation and evolution of sand bars and, given the large body of work on point bars and other forced bars, we concentrate on spontaneously arising free bars in sand and gravel bed rivers.
The linear theories suggest that a shallow water (i.e. hydrostatic) flow model is a reasonable approximation for this problem and that changes in vertical flow structure are not necessary for the production of rivers bars; these analyses also suggest that consideration of several physical effects is required to correctly obtain the initial fastest growing wavelength of the bars (e.g. Nelson, 1990). This small-amplitude theoretical framework has been extended to understand finite amplitude bars (Colombini et al., 1987;Schielen et al., 1993). These theoretical studies suggest that a model of shallow water flow coupled with a bedload transport model including local slope effects is the minimal, physics-based model that can explain initiation and development of free bars. With this background and help of recent advances of computational power, fully non-linear numerical models have been proposed to directly simulate non-linear free bar dynamics (e.g. Shimizu and Itakura, 1989;Nelson, 1990;Federici and Seminara, 2003;Defina, 2003;Takebayashi and Egashira, 2008;Crosato et al., 2012;Siviglia et al., 2013;Kimura et al., 2015;Kimura, 2016a, Iwasaki, Nelson, et al., 2017a).
In this section, several examples of the numerical simulations of free bar morphodynamics are presented. First, a computational example demonstrates the formation and development of alternate bars (Figure 1a), the lowest mode bars appearing in channels. The computation is performed in a straight channel with a length of 100 m, constant channel width of 0.9 m, and initial slope of 0.005. The bed is composed of single grain-size sediment with diameter of 0.76 mm. A steady water discharge of 6.4 l/s and associated equilibrium bedload transport are imposed at the upstream end. The experimental result of this situation is demonstrated in Figure 1c, showing clear formation of alternate bar trains (Akahori et al., 2011). Figure 1c also shows that alternate bars formed at the early stage of the experiment (top panel of Figure 1c) are elongated during the development process (bottom panel of Figure 1c), as also documented by other experiments (e.g. Fujita and Muramoto, 1985). The linear stability analysis accurately predicts the initially formed sand bar wavelength, but fails to explain bar elongation (Nelson, 1990). Figure 1d shows computational results of alternate bar formation and development for the same conditions and at the same times as the experimental results shown in Figure 1c. This result clearly shows that the simple numerical model accurately reproduces observed free bar dynamics including the initial formation of alternate bars and bar development to the equilibrium state. This capability of numerical models for reproducing alternate bar morphodynamics are well known and several public software packages for river morphodynamics are able to capture these morphodynamic features: Delft3D (http://www.deltares. nl, Lesser et al., 2004), TELEMAC (http://www.opentelemac. org, Mendoza et al., 2017), andiRIC (http://www.i-ric.org, Nelson, Shimizu, et al., 2016b). The results shown in Figure 1d are computed using the solver Nays2DH within the iRIC publicdomain interface. We further extend this same model to a somewhat more interesting case by examining the formation and development of multiple row free bars. Multiple row bars (as seen in the photograph in Figure 1b) tend to develop in channels with relatively high width/depth ratios. In linear analyses these morphodynamic features are described by the order m of the lateral Fourier mode (Parker, 1976;Fredsøe, 1978;Colombini et al., 1987;Watanabe, 2007;Seminara, 2010;Pornpromin and Izumi, 2011). Here, so-called double row bars (m = 2) are numerically simulated. The computational settings are the same as the alternate bar simulation shown earlier except for the channel width/depth ratio and the discharge. The channel width and the discharge are doubled to 1.8 m and 12.8 l/s, respectively, for this case, but the initial water depth is the same (i.e. we use the same unit discharge in both calculations). The initiation and development of the bars are shown in Figure 1e. At the early stage of the simulation, double row bars form in the channel, but subsequently, the simulated double row bar changes to single-row alternate bars. Such reduction of the lateral mode of bars is also seen in laboratory experiments (e.g. Fujita and Muramoto, 1985). The numerical model captures this kind of morphological complexity, although this point needs further investigation (Watanabe, 2007;Pornpronmin and Izumi, 2011;Seminara, 2010).
These simple examples with well-sorted sediment demonstrate the basic success of the numerical modeling of free bars. The model has also been extended to more complex problems: unsteady flow conditions (e.g. Kobayashi et al., 2006), graded sediment (e.g. Takebayashi and Egashira, 2008) and many applications to bar responses in real rivers including characteristics of both free and forced bars with good results. Extending these approaches to the more realistic conditions of changing sediment supply introduces additional effects, as demonstrated in the next example.

Sand bars with sediment supply limitation
Temporal variations in sediment supply to a river can be related to natural disasters, land use, or upstream river management. These variations produce changes in channel sediment storage and channel patterns of topography, sediment transport rate, bar size and bed surface texture (Wilcock, 1998). Understanding the supply-driven morphological changes in a river with varying sediment inputs is critically important for maintaining or restoring the river. Despite this importance, most computational studies of river morphology assume equilibrium transport conditions at the upstream end of a real or idealized study reach. Equilibrium sediment supply is probably almost never met in river reaches except over timescales that are quite long relative to river adjustment and only when averaged over significant variations in the hydrograph, so morphological equilibrium in terms of bar amplitudes, spatial grain size distributions, and other physical river characteristics is likely to be rare at any given time. To better understand the role of disequilibrium sediment supply, Nays2DH was applied to understand the morphological processes in a braided river with varying upstream sediment supply. This model solves the two-dimensional depth-averaged flow equations and calculates sediment transport, sediment sorting, and geomorphic changes considering bank erosion and vegetation. This model has been widely applied in scientific research for rivers; readers are referred to Iwasaki et al. (2015) for more details.
To investigate the evolution processes of braided rivers with disequilibrium sediment supply, the model was applied to a straight, 3.0 m wide, 40 m long channel with erodible banks. The slope in the simulation was set to 0.01 and discharge was 0.012 m 3 /s. Mean diameter of sediment was 1.25 mm. The sediment-transport rate was calculated using Meyer-Peter and Müller's (1948) formula. Figure 2 shows the development of a braided river. As in the case of Figure 1e, multiple row bars developed and migrated downstream in the initial stage, as shown in Figure 2a. Multiple bars merged into large central bars, which strongly deflected the flow through topographic steering. In regions of strong flow convergence such as confluences, the bed was deeply scoured by the flow. The channel migrated laterally through erosion of the banks where high velocities occur in near-bank regions (Figures 2b and 2c). Over time, new channels are formed and old channel are filled with sediment. As time increased, the channel network is completely reformed (Figures 2d and 2e). The numerical model reproduced the essential features of braided rivers, such as the generation of new channels and abandonment of old channels, bifurcation and confluence of channels, and lateral migration of channels (Jang and Shimizu, 2005a). This numerical experiment (Run 1) was conducted as a baseline to investigate the development of bars and geomorphic changes in braided channels with uniform grain sizes with sediment supply variation. Thus, after this initial run during which the equilibrium transport rate was supplied, that supply was increased above the equilibrium value. Sediment supply was set to 120% of the equilibrium value for Run 2, 150% for Run 3, and 170% for Run 4. The computational time step for all the runs was 0.01 second; there were 80 grid points in the stream-wise direction and 240 in the transverse direction. The aspect ratio of the grid was 4.5. Alternate bars 1.0 cm in height and 5 m in wavelength were introduced as an initial perturbation to encourage the development of bars. Figure 3 illustrates simulation results of channel morphology with the sediment supply variation for each run. As the sediment input increases for Runs 1, 2, 3 and 4, the channels are progressively aggraded in the upstream part of the reach, which is due to the storage of excess supplied sediments in the channel. Even for Runs 3 and 4, almost all aggradation is in the upstream part of the reach, whereas the channel in the mid-reach and downstream reach shows little difference relative to the equilibrium supply result (Run 1). Aggradation magnitude varied throughout the length of the channel due to local variations in channel patterns and lower channel changes (Germanoski and Schumm, 1993). The bed aggradation affects channel morphology by increasing the number of channels and the Braiding Index, which is defined here as the average number of local elevation minima in a cross-section (typically averaged in the streamwise direction over at least a few channel widths). As sediment input increases, which causes aggradation and the development of individual sediment storage sites (or braid bars), the Braiding Index increases (Ashmore, 1991;Germanoski and Schumm, 1993).
The lower channel networks become more complex with increasing sediment supply ( Figure 4). The averaged Braiding Index for the entire channel is 2.03 for Run 1, 2.47 for Run 2, 2.83 for Run 3, and 2.92 for Run 4. The Braiding Index is 1.8 for Run 1, 2.5 for Run 2, 3.2 for Run 3, and 3.7 for Run 4 using only the upstream reach, where aggradation is most pronounced. However, the Braiding Index is 2.2 for Runs 1 and 3, 2.9 for Run 2, and 2.7 for Run 4 for in the middle reach, and there is no consistent streamwise variation at either the middle or downstream reach. Thus, the numerical experiment confirms the contention by Germanoski and Schumm (1993) that the Braiding Index increases with increased sediment input, which means that the channel is more complex with increasing sediment supply. The channel cross profile is affected by the variation of sediment supply. Braiding Index increases during aggradation, which causes a more complex profile in cross-sections.
The numerical model reproduces the features of braided rivers, such as the generation of new channels and abandonment  of old channels, bifurcation and confluence of channels, and lateral migration of channels. The model also confirms the more complex adjustment of the channel to the oversupply of sediment; a converse effect is found in the case of undersupply, as one might expect. Overall, this simple example shows how computational bar morphodynamics models can be used to examine more complex (realistic) river behaviors associated with oversupply or undersupply of sediment to a river reach.

Future work
Despite the successes of simple computational bar morphodynamics models, there are remaining questions about the effects of the relatively simple submodels used in the overall framework for computing bar initiation and finiteamplitude behavior, and various aspects of the simulation should be revisited to improve modeling of both free bars and other morphodynamic features. For example, the approximations employed (shallow water flow model, equilibrium bedload transport model, limited or non-existent secondary flow treatments, unproven gravitational effect treatments, etc.) present limitations. Several non-linear numerical simulations have indicated that simulated bar morphologies (especially multiple bars and time evolution of the morphodynamics) are somewhat sensitive to the choice of the submodels and their parameters (e.g. Schuurman et al., 2013;Mosselman and Le, 2016;Iwasaki, Shimizu, and Kimura, 2016a). Within this context, the local slope effect on bedload transport and secondary flow modeling have received the greatest attention (e.g. Iwasaki, Shimizu, and Kimura, 2016a;Baar et al., 2018), but this sensitivity may also be caused by the interaction or combined effects of the submodels and their parameters. More sophisticated models will be required for better understanding of free bar dynamics and the predictive behaviors of numerical models that have been used in the last few decades.
Most of the morphodynamics models in iRIC and similar models from other research groups incorporate mixed grain size sediment-transport models. Grain sorting patterns on bars and routing of mixed sizes over bars is a critically important part of morphologic adjustment that is not well understood. This is an area where our computational capability has run ahead of our observations, and better data sets on spatial grain-size sorting (both on and in the bed) are needed to test and improve mixed-size sediment modeling capabilities (e.g. Nelson et al., 2015aNelson et al., , 2015b. The development and use of more sophisticated models suggested earlier may only make minor improvements in the accuracy of the results for reproducing bar dimensions for simple situations. However, these attempts are likely to introduce a whole new layer of complexity associated with the mutual interactions of different types of bedforms. More sophisticated physics-based models will capture other morphological signatures in addition to free bars. One obvious area for research is the interaction between bars and bedforms (dunes and antidunes, both two-and three-dimensional). Colombini and Stocchino (2012) examine the interaction of three-dimensional dunes and free bars, and the coexistence of free bars and threedimensional antidunes at the linear level, but the non-linear behaviors and fully developed morphologies remain unknown. Application of high resolution hydrodynamic models and detailed discrete element sediment-transport models (e.g. Nabi et al., 2013a;Khosronejad et al., 2015;Schmeeckle, 2017) or advanced depth-integrated approaches (e.g. Uchida and Fukuoka, 2017;Yabe, 2016b, Iwasaki, Inoue, et al., 2017b) for reach-scale river morphodynamics are promising for future developments and new understanding.
Numerical models of sand bars have been well validated with results from controlled laboratory experiments. In addition to small-scale experiments, long-term sand bar morphodynamics have been recently investigated with the help of long-term field observation and aerial photographs (e.g. Adami et al., 2016). Moreover, some real-time measurements of the bar dynamics during a flood using non-contact measurement techniques (Nelson, Kinzel, et al., 2016a) are becoming available. These attempts will provide useful datasets for further validation of the numerical models in complicated real river systems. Not only comparison with the experimental and field data, but also theoretical analyses of sand bars will help us to interpret complicated non-linear behaviors of the numerical model in a simple manner (e.g. Vanzo et al., 2011;Siviglia et al., 2013;Adami et al., 2016). Diverse research methodologies incorporating careful field observations with state of the art computational methods will provide an increasingly comprehensive understanding of bar dynamics.
Finally, future work should greatly expand the nascent field of bar-vegetation and bar-habitat work. Although not really discussed in the work described earlier, this area of research is critically important for problem solving in real rivers. The interactions between biological parts of the system and physical ones are arguably more important than making improvements on the physical side of the problem, at least for real-world management issues in rivers.

Microscale: Bedforms
Although at first acquaintance the principal difference between mesoscale bars and microscale bedforms might appear to be size, such an understanding misses an important distinction that has a more interesting physical basis. Generally, or at least in their more common forms, bedforms arise spontaneously in sediment-transporting flows due to changes in vertical structure of the flow associated with spatial accelerations. Although the same physics can certainly play a role in bars, just as complex horizontal flow routing in the absence of vertical structure changes can play a role in bedforms, separation of the two along these physical lines has been a tremendous advantage for making progress over the past few decades. Thus, bedforms are separated from their bar counterparts by the need for a modeling approach that incorporates vertical structure changes and non-hydrostatic effects. Thus, a shallow-water (vertically averaged) flow model is inadequate; modeling microscale features generally requires either a two-dimensional streamwisevertical approach or a fully three-dimensional flow model.
We provide examples of both here. In both cases, computational effort is larger than that found for mesoscale features; typical computer run times for the examples shown below are from minutes for the two-dimensional case to hours for the three-dimensional case using typical personal computers.

Bedform development with a two-dimensional model
Understanding and quantifying the formation and evolution processes of microscale alluvial bedforms are of great importance not only from the fundamental perspective of river science, but also from the standpoint of real-world application and management of river-related issues. The evolution and transition of microscale bedforms, such as ripples, dunes, flatbed transition and antidunes, are important for flood stage prediction due to the importance of bedform-generated bed resistance and hysteresis effects, optimization of the measures for maintenance of river channels for inland navigability, and effective sediment passage during flushing operation in reservoirs and sediment bypass channels. Bedforms are a fundamentally important consideration in any computational treatment of river flow or sediment transport.
Substantial contributions have been made by our group and collaborators on developing computational models and conducting various studies to explore bedform evolution and transition processes including large-scale experiments and real-world applications. Some of the major studies include Nelson et al. (2005Nelson et al. ( , 2010, Giri and Shimizu (2006), Shimizu et al. (2009), Yamaguchi and Izumi (2005) Giri et al. (2011Giri et al. ( , 2015, Nabi et al. (2012Nabi et al. ( , 2013aNabi et al. ( , 2013b. We begin by describing a simple two-dimensional approach for predicting the initiation and growth of bedforms in sediment-transporting flows. The vertical two-dimensional model incorporates non-hydrostatic pressure, a non-linear turbulence closure, a free-surface flow condition, and both equilibrium and non-equilibrium sediment-transport formulations. This is the simplest approach that can simulate the bedform evolution process and bedform flow resistance in a physically correct manner. The solver has been described and verified in a number of previous studies . The predictions of the flow model have been compared to detailed laser-Doppler velocimetry (LDV) measurements over fixed beds . Results have also been compared to a variety of large-eddy simulation (LES) and Reynolds-averaged Navier-Stokes (RANS) computations Nelson et al., 2005); comparison of predicted bedform geometry with experimental values of wavelength and height show that this model performs well despite the relative simplicity of the approach. Importantly, this model is computationally much less demanding than three-dimensional models.
Application of equilibrium and non-equilibrium sedimenttransport models to bedform modeling was examined in the study by Yamaguchi et al. (2009). This study clarifies the significance of local bed slope and non-equilibrium effects on the bedload transport mechanism and implies that either mechanism is sufficient for reproducing bedform behavior. The work contrasts the equilibrium transport model with local bed slope effect against a non-equilibrium transport model in terms of their effects on bedform evolution. The formulation also provides a non-equilibrium bedload transport approach by incorporating a formulation for sediment saltation distance as a length scale that characterizes the microscale non-equilibrium transport process. This approach is compared to equilibrium bedload transport formula using the Meyer-Peter and Müller (1948) bedload equation including the local bed slope effect proposed by Fredsøe (1974).
The sediment step-length is composed of many smaller movements of the particles with intermittent stops and starts, which implies that the microscale processes within the steplength are not completely represented. Previous studies indicate that the saltation length, which is the actual hop distance of a sediment particle, is related to near-bed flow conditions. In one of our previous studies , the steplength relation was assumed to be a function of the bed shear stress to treat the particle movement scale as a direct function of the local flow conditions. This assumption allowed simulating the bedform transition process from dunes to flatbed under varying flow conditions. In the most recent approach, we have introduced the saltation length proposed by Sklar and Dietrich (2004) as the appropriate measure of near-bed microscale nonequilibrium transport instead of the overall step-length. We have also modified the definition of pick-up rate as sediment entrainment rate according to the saltation length, since formulations for the pick-up rate, proposed by various researchers, are developed using the total step-length rather than the saltation distance.
A large-scale experimental case has been simulated to assess the model performance. The experiment was carried out in an experimental facility built within the Tokachi River, Japan ( Figure 5). A comparison between simulated and measured (using an acoustic Doppler current profiler [ADCP]) vertical two-dimensional velocity is shown in Figure 6. Both measured and simulated results show the characteristic features of the flow field over bedforms, with very low velocity along the recirculation zone behind the crest and considerably higher velocity near the free surface over the bedform crest. However, the observed two-dimensional flow-field and the magnitude near the bed and the free surface are not entirely consistent with the computational result. Such discrepancies can be attributed to the complex flow field, the lack of temporal averaging, and the bedforms in the experiment, which were at least partly affected by the vertical piles on the left-side wall of the experimental channel. Since the numerical model is widthaveraged vertical two-dimensional, the approach cannot treat the flow complexity and dissipation induced by the side wall, which is inherently three-dimensional in nature. The simulated length of the flow separation zone behind the bedform is about 2 m, which is about five times the bedform height. This is consistent with the observed extent of flow separation behind bedforms in previous empirical studies.
The channel bed measured by a multi-beam echo-sounder is shown in Figure 7. The left bank with vertical piles appeared to induce a noticeable effect on the bedform feature near this bank, while there was almost no effect of the side wall near the right bank. The wavelength and height of the bedforms were found to vary between 10 and 12 m and 0.2-0.4 m respectively.
The model with an equilibrium transport model without bed slope correction could not simulate the bedform evolution. For the case using an equilibrium bedload-transport model with slope effect (Case 1), the computed bedforms reached a quasi-equilibrium state rather quickly without noticeable development of the wavelength (see Figure 8, upper left). For the simulated case using non-equilibrium bedload transport models (Cases 2 and 3), the bedforms (particularly the wave length) attained dynamic equilibrium only after about three hours of computation ( Figure 8, lower left). Generally, the gradual increase in wavelength during the evolution process shown in the measurements (right plot in Figure 8) appeared to be better replicated by the non-equilibrium model. However, more careful measurements, perhaps over a longer reach, would be necessary to draw a firm conclusion.
The computational result yielded a maximum bedform height of about 0.5 m for the fully developed bedforms in all simulation cases. Likewise, the mean bedform length was simulated to be 15 m and 12 m for Cases 2 and 3 respectively. However, in the experiment, the maximum wavelength and height was about 10 m and 0.4 m respectively. These discrepancies can be attributed to several factors, such as the possibility that the bedforms had yet to attain the equilibrium state in the experiment or that processes that cannot be adequately treated in a two-dimensional model were important, especially considering the weakly three-dimensional structure of the observed bedforms. For example, lateral dissipation induced by the wall disturbances and three-dimensional features in the experiments, which obviously could not be replicated by a vertical two-dimensional model, may have had some significant effect. Moreover, the wall disturbance may lead to an increase in the time required for bedform development due to the additional flow resistance not treated in the two-dimensional model. The periodic boundary condition used in the computation may also lead to early development of the bedforms in the numerical computation. The wavelength of the bedforms in the experiment was about two to three times the water depth, while experience and empirical bedform predictors suggest that the wavelength of the bedforms should be about four to seven times the water depth.
The water depth variation associated with the evolution of form drag exerted by the bedforms was almost identical for all the cases, showing about 3% to 7% difference compared to the experiment. This discrepancy is not very large and can be attributed to the graded sediment effect in the experimental channel with a maximum diameter of about 7 cm, as well as the potential influence of the vertical side wall with piles, which may lead to some additional factors that increased the flow resistance in the experiment.
The spatial distribution of instantaneous dimensionless bed shear stress, pick-up and deposition as well as the dimensionless bedload transport rates over a bedform for various stages of development showed that these parameters are noticeably different for the models with equilibrium transport with bed slope correction and the non-equilibrium transport (the results are not shown here). Although both mechanisms may be important, the issue of how to incorporate these mechanisms in a morphodynamic model to represent the correct physical behavior of bedform evolution in different stages remains ambiguous. The major differences in the two approaches are in how the wavelength and height evolve in a temporal sense. Therefore, further investigation is necessary to test the situations where there is greater variation in the bedform characteristics. For example, cases with a time-varying hydrograph in which the adjustment of wavelength and height during transient conditions is more pronounced would provide a better test. This issue is important in the context of recent advancement in numerical computation of bedform morphodynamics and deserves further attention with more analysis and careful observations. Bedform development with a three-dimensional large-eddy simulation model The development of the bedform model described in the previous section is limited by both the assumption of the twodimensionality and the associated inability of the model to treat the details of the turbulence field. Many real bedforms are three-dimensional and a more general view of the flow field requires a better understanding of the turbulence field than closure-type models can provide. These considerations both suggest that a more computationally intense, turbulenceresolving three-dimensional model is required for treating bedform formation and evolution. Thus, a large-eddy resolving flow solver was combined with a space and time resolved sediment-transport approach to investigate bedforms at a more detailed level.
This approach goes back to the most fundamental reasoning regarding the coupled flow-bed-sediment system. The bed morphology is the result of sediment transport through erosion and deposition. The sediment motion is governed by the turbulent flow, and the details of the turbulent flow are sensitive to the bed morphology. Together, these interactions form an inherently complex problem. Turbulence produces near-bed velocity and pressure fluctuations that give rise to fluctuations in the forces on and the resultant motion of sediment particles, as a result of which sediment may be picked up and transported by the fluid. However, the presence of sediment particles in the flow may generate or suppress turbulence, thus creating a dynamic feedback system. The change of bed also affects the flow motion and the turbulence structures, and it generates a strongly coupled mechanism, which is not yet fully solved in any general sense.
In the conventional approaches we have introduced so far in this article, sediment transport has mainly been described in terms of bulk fluxes related to mean flow properties. Some of these sediment transport formulae include a threshold of motion (the critical shear stress). However, using a more resolved flow treatment allows a temporally resolved approach wherein the transport of non-cohesive sediment particles on a loose sediment bed can now be simulated using deterministic approaches (Nabi et al., 2013a(Nabi et al., , 2013b. This approach is more advanced as the sediment transport and the bed deformation can be simulated at process level, in a physics-based way, which gives a better insight into the formation of microscale bed morphodynamics (ripples and dunes). Detailed simulation of flow, sediment transport and bed morphology is realized by employing LES for the turbulent flow (Nabi et al., 2012) and a model for sediment transport based on spherical rigid particles moving in a Langrangian framework (Nabi et al., 2013a). The bed deformation is the result of sediment pick-up and deposition at each finite area, which are the computational cells of the model (Nabi et al., 2013b).
The flow is solved on a multi-level Cartesian grid with local refinement to gain computational time efficiency. The grid is adaptive and can be refined close to the bed. The governing flow equations were discretized by a finite volume method, and the resulting sparse matrix is solved by a geometrical multigrid method. The flow and gravity forces exerted on sediment particles are determined and sediment pick-up, sediment transport in the water column, sediment transport by sliding, and sediment deposition are calculated based on first principles rather than empirical flux laws. The bed change at each cell is determined by mass balance between the picked-up and the deposited sediment particles.
This solver was applied to simulate the formation of ripples and dunes and tested for the capability to replicate existing flume measurements. Figure 9 shows the simulated bed change under constant discharge. The bed evolves from a totally flat plane (Figure 9a). Under the effect of turbulence fluctuations, the bed is disturbed (Figure 9b), and the disturbances on the bed migrate and interact to form ripple-like features (Figure 9c). The ripples in turn grow and interact and eventually form dunes (Figures 9d and 9e). Those dunes grow and ultimately reach a dynamic equilibrium (Figure 9f). Figure 10 shows the turbulence coherent structure ( Figure 10a) and the suspended sediment above the bed (Figure 10b). The horseshoe vortices near the bed are a factor that significantly enhances the sediment pick-up. The legs of the horseshoe structures are associated with high positive and negative vorticity in opposite directions, generating a strong vertical vorticity, which gives a significant upwelling to the flow and hence produces sediment entrainment. In some instances, the intensity of this vertical fluid upwelling is as strong as 30% of the bulk flow velocity (Grigoriadis et al., 2009). Figure 11 show the evolution of the bed (along the center of the channel) in time. The merging and splitting of dunes can be clearly observed. The dunes with smaller amplitude have higher celerity than those of larger ones, which leads to an amalgamation process as described by Venditti et al. (2005).
A natural river involves discharge variations that influence flow, sediment transport, and the resultant river bed topography. A constant discharge ultimately yields a statistically unique bedform pattern. As the discharge varies, the bedforms change with transition phases. For a variable discharge, the form of the bed is not unique due to a time lag between the discharge gradient and the bed deformation that leads to hysteresis effects. The model proposed by Nabi et al. (2012Nabi et al. ( , 2013aNabi et al. ( , 2013b was further extended to model the bed deformation under variable discharges. Figure 12 shows the applied discharge hydrograph and the results for the form drag for such a schematized hydrograph, with a flow period of 4000 seconds with uniform grain sizes of 500 μm. The resulting form drag diagram has the shape of an open loop, which is the well-known hysteresis effect associated with bedform disequilibrium (Simons et al., 1961). As the sediment size increases, the loop gets wider, i.e. the hysteresis effect gets stronger. For more details, we refer the interested readers to Nabi et al. (2012).
Under certain variable hydrographs the bedforms may disappear in the so-called flatbed transition. The occurrence of an upper plane-bed regime has been a topic of interest in several previous studies, both experimentally (e.g. Naqshband et al., 2014) and numerically (e.g. Shimizu et al. 2009). The model described here captures the transition to a flatbed in upper plane-bed regime . The starting point of the simulation was a flat bed under a relatively low discharge of 0.12 m 3 /s. Once the flow reached its statistically steady state, the sediment transport and morphologic evolution are turned on and then the discharge is linearly increased in a time span of 2000 seconds to a maximum value of 0.48 m 3 /s. This is associated with bulk velocities of 0.3 m/s and 1.19 m/s respectively, Reynolds numbers in the range of 120 000 to 476 000 and Froude numbers in the range of 0.15 to 0.6. Figure 13a shows the path of these simulations in the rising and falling flow stages in the phase diagram by Southard and Boguchwal (1990), and Figure 13b shows the implemented discharge hydrograph. The maximum discharge for these simulations lies in the upper plane bed regime. Figure 14 shows the generation of bedforms during the rising stage. Dunes are generated during the rising stage of the discharge (Figure 14a). Once the discharge reaches the maximum value (in the upper-plane regime) the flattening of the bed was observed in the simulation (Figure 14b). The flattening is associated with complex interactions between flow and sediment transport. At the transition to plane bed regime, sediment was eroded from the crest of the bedforms, and the sediment was mainly deposited behind the lee side (in the trough of the dunes). This is in contrast with the conditions for the dune regime in which the sediment goes mainly into suspension and after a certain travel time deposits on the stoss of the following dunes. After flattening of the bed (in the absence of dunes, Figure 14c), the sediment was primarily transported in suspension (not shown here). During the falling stage of discharge, the bedforms were regenerated in a realistic way (Figures 14d-14f).
While the intent of this model was originally to treat complex bedform dynamics, the approach is significantly more general. For example, the model has been extended to simulate bed morphology around piers by Nabi et al. (2013b) and Kim et al. (2014Kim et al. ( , 2015Kim et al. ( , 2017. Nabi et al. (2013b) simulated scouring around a single circular pier and compared the result with existing results from Khosronejad et al. (2012). Good agreement was found between the simulated results and the measured data. This study was extended to simulate local scour at two adjacent piers by Kim et al. (2014Kim et al. ( , 2017 and rigidemergent vegetation by Kim et al. (2015). As an example, Figure 15 shows the results from Kim et al. (2015) to investigate the effect of single pier, four piers and 16 piers (replicating vegetation density) on bed morphology. The boundaries in the streamwise and transverse direction are set as periodic, and the distance between the piers with diameter D are set to 10D (for one cylinder), 5D (four cylinders) and 2.5D (16 cylinders). Figure 15 shows the bedform for cases 10D, 5D and 2.5D, after their equilibrium state is reached. Note that the bedforms significantly differ when the distance between piers decreases. Those differences are associated with the interaction between the piers. As the piers get closer, the flow around piers is affected by the wake behind upstream piers. The sediment transport and the morphology around piers were primarily affected by two mechanisms, the acceleration of flow between piers and the effect of von Kármán vortex street from upstream piers. Figure 16 shows the simulated maximum scour profiles for the cases 10D, 5D, and 2.5D. The growth of the maximum scour depth is delayed as the cylinder interval increases (i.e. corresponding to a decrease in the vegetation density), while the maximum scour depths are similar.
Overall, the large-eddy approach combined with a physicsbased sediment model yields phenomenal results, but at a large computational cost. However, as parallel computing and cloud-based predictions become more common, this approach seems likely to revolutionize computational morphodynamics. There are, however, still applications for simpler models, including both two-and three-dimensional approaches that use advanced turbulence closures, as follows.
Streak development using a three-dimensional model with turbulence anisotropy Secondary currents are a classical topic in fluid engineering and can be classified into two types: the secondary current of the first kind (driven by an imbalance between hydrostatic pressure and convective accelerations in the vertical) and the second kind (driven by spatial gradients in turbulence intensities). Both give similar spiral structures to flows but the mechanisms are very different. Students are typically familiar with secondary flows of the first kind due to their important role in producing point bars in rivers, where curvature is the driving mechanism. The secondary current of the second kind is generally much weaker (and less well understood) than that of the first kind but is nevertheless important for river morphology because it is responsible for the generation of longitudinal streaks in the bed profile.
There is still controversy on the initiation of the secondary current of the second kind in open channels with large aspect ratios (= B/h > 5, where B is channel width and h is water depth), in which the effect of side walls is weak (Nezu et al., 1985). Some papers report that the secondary current in the center part of the smooth, shallow, rectangular channel with B/h > 5 disappears because the secondary current is initiated by the strong turbulence anisotropy near the side walls (Nezu et al., 1985). However, Blanckaert et al. (2010) reported the experimental result in an open channel flows with very large aspect ratio (B/h > 10), in which clear secondary flow cells were observed even at the center part of the channel. In this study, we attempt to reproduce the secondary current in a wide, shallow rectangular channel using a URANS (unsteady Reynolds-averaged Navier-Stokes) three-dimensional model  and show the capability of the solver for this phenomenon. We also try to simulate the generation of the streak structure for the case of a movable bed.
NaysCUBE is a fully non-hydrostatic three-dimensional solver for open channel flow based on the RANS equation ). The solver is available in the public domain within the iRIC river software package (www.i-ric.org). The governing equations are transformed using a generalized curvilinear coordinate formulation, such that any domain is transformed into a cube for computational solution. NaysCUBE also includes modules of bedload and suspended sediment transport and thus can simulate bed morphology induced by three-dimensional flow patterns. To simulate the secondary current of the second kind, it is essential to use a turbulence model that can consider the anisotropy of turbulence. NaysCUBE offers two turbulence models; the standard linear k-ε model and the second-order non-linear k-ε model. We use the latter approach in the results shown here. Figure 17a shows the computational result for the crosssectional velocity vectors and the color contour of the velocity magnitude; Figure 17b shows the observed vortex patterns in the experiment. In this computation, the second-order nonlinear k-ε model by Ali et al. (2007) is adopted as a turbulence model. The secondary currents are simulated well and the number of vortices matches the experimental result. It is known that the vortex tubes with a diameter of the flow depth align regularly in the span-wise direction except near the side walls in shallow, open channel flow. Near each side wall, two vortices are generated near the surface and near the bottom. The present simulated result reproduces this flow pattern and shows good agreement with experimental measurements (Suzuki et al., 2014).
One advantage of the numerical simulation is that unrealistic computational conditions can be considered, for instance, we can completely neglect the presence of side walls by setting a periodic boundary condition in the lateral direction. Figure 18 shows the computational result with the periodic boundary condition in the lateral direction. Clear secondary currents are generated and the scale of the diameter is also nearly the same as the water depth. The two vortices near the surface and near the bottom observed near the side walls in the previous case are not generated. As for the development process, it took more time (about 1000 seconds to reach the equilibrium state) starting from the initial conditions (parallel flow with no secondary flow), which is longer than the case with side walls (⋍ 300 seconds). This result indicated that the presence of the side walls is not the necessary condition to initiate the secondary current of the second kind but also that the presence of side walls accelerates the development of such flow features. Figure 19 shows the birds-eye view of the well-developed computational result with a movable bed. The hydraulic conditions are set for the case of the experiment by Blanckaert et al. (2010). These are the same conditions as used for the results in Figure 17 except for the inclusion of the movable bed. The color contour of the bed shows the bed elevation change from the initial flat bed. For the bedload flux model, we adopted the Ashida and Michiue (1972) model. For boundary conditions, a periodic boundary condition is used for the streamwise direction and side walls are included. A clear streak structure is generated at the bed. The erosion happens in the region of downward flow and the deposition happens in the upwelling region. The magnitude of vorticity in the streamwise direction becomes slightly smaller than for the case with a flat bed, probably because the wetted perimeter increases due to the wavy shape of the bed.
We made some considerations on the secondary current of the second kind in a wide, shallow, rectangular, straight, open channel. A three-dimensional solver NaysCUBE was used to simulate the flow. Computational results show that the secondary current is generated not only near the side walls but also at the center part of the channel, in a manner identical with the experimental results. The secondary current was also clearly simulated in case of the periodic boundary conditions in the lateral direction, in which the effect of side walls was completely neglected. When we changed the bed condition from a fixed flat bed to a movable sediment bed, a clear streak structure developed on that sediment bed in agreement with observations.

Future work
As noted earlier in the future work section on mesoscale bars, the interactions between bedforms and bars are important for improving computational predictions of bar behavior. The converse is also true; although most computational and experimental studies of bedforms have concentrated on situations without larger-scale non-uniformity, bedforms in rivers are frequently affected by the presence of bars, such as in the case of dunes migrating along the curved side slope of a point bar, where the interactions of bedforms with the point bar have a significant effect on the direction of sediment transport and sediment sorting Smith, 1983, 1984). Very little attention has been paid to the mechanisms of bedform adjustment in flows that are spatially non-uniform at larger scales. Understanding how wavelength, height, and shape adjust when bedforms migrate through regions of changing average depth and acceleration is not well understood, but accurate predictions of bedform fields in realistic flows require predictions of this nature. Although the work of Nabi (Nabi et al., 2012(Nabi et al., , 2013a(Nabi et al., , 2013b) is a step forward in predictive capability for riverine bedforms, such techniques still present significant and, in many cases, insurmountable computational demands for application to larger domains such as river channels. One of the immediate challenges is to work out ways to bring methods for bedform prediction into much larger-scale models. This is the key to prediction of the spatial and temporal variation of form drag, which is of central importance for improving river models in general.
For simpler approaches of bedform prediction such as the two-dimensional approach described earlier, better predictions almost certainly require a more detailed view of grain motion. Although the model we describe uses the simple saltation length predictor of Sklar and Dietrich (2004), recent experimental work shows that the distribution of particle hop lengths are more complicated (e.g. Roseberry et al., 2012, Fathel et al., 2015 than reflected in that simple approach. It seems likely that including a more realistic treatment of the sediment motion and the spatial lags associated with that motion can improve predictive capability for bedforms.
As in the case of bars, the linkages between bedforms and a wide suite of vegetation and riparian habitat issues deserve greater attention, especially with regard to the infiltration of finer materials into gravel spawning substrates in mixed sandgravel river systems. Prediction of bedform geometry in such systems is complex and poorly understood due to the factors of strong spatial gradients in bed roughness, non-equilibrium sediment supply, and interactions between channel and hyporheic flows. These same issues arise when considering bedforms in sand-gravel systems that have little topographic signature, but are identified primarily through spatial sorting patterns of sediment.

Macroscale: Channels
Up to this point, the examples we have considered have been primarily focused on features that evolve on the bed of the river, with little or no consideration of the planform shape of the channel itself. In this section, we present examples that specifically deal with planform changes in the channel form itself, meaning that the banks as well as the bed must now be treated in morphodynamic modeling. Generally, this requires thinking about longer timescales and the potential interaction of bed and bank processes. Computational demands for channel evolution models can be large for large domains, so judicious choices of timescale and spatial scale are necessary, as run times for channel evolution models shown here range from hours to a few days on typical PCs.

Meandering channels
The morphology of natural rivers is shaped by the interaction of flow, sediment transport, bed morphodynamics and bank erosion and deposition processes. In the past, many studies have tried to address the mechanisms that give rise to river planforms (Hasegawa, 1978(Hasegawa, , 1984Ikeda et al., 1981;Engelund, 1974), but most of those studies have been qualitative, with little directly predictive capability for the linked evolution of beds and banks in rivers.
Several meandering numerical models have employed the equation of bank migration developed by Hasegawa (1978Hasegawa ( , 1984 and Ikeda et al. (1981) to demonstrate the evolution of complex river planforms with channel cutoffs (Johannesson and Parker, 1989;Sun et al., 2001;Zolezzi and Seminara, 2001;Lanzoni and Seminara, 2006;Crosato, 2008;Frascati and Lanzoni, 2010). Several other numerical studies have made it possible to study bed topography under the assumption of non-cohesive bed material (Shimizu et al., 1996;Nagata et al., 2000). Also, some bank erosion models developed with the capability of simulating both cohesive and non-cohesive banks (Darby et al., 2002;Duan and Julien, 2005) have made it possible to understand the factors governing alluvial channels over relatively short timescales. However, in addition to these shortterm factors, there are other longer-term factors like vegetation growth that can affect a meandering channel. Importantly, although all these approaches have helped to understand and make limited predictions regarding the change in channel planform over time, none present a coherent approach linking both short-timescale bed evolution and longer-timescale bank effects to predict channel migration. In this section, a model is presented that can consider the effect of short-term factors as well as the key elements of longer-timescale processes.
As in the bar models described earlier, two-dimensional depth-averaged continuity and momentum equations are used in the flow model. Bed evolution is calculated using the Exner equation for sediment mass conservation. Both the critical Shields number and the critical shear velocity are computed using the formula of Iwagaki (1956). The near-bed velocity in the direction of the vertically averaged flow is computed as given by Engelund (1974). Bank migration is obtained by integrating the sediment continuity equation in the near-bank region in the cross-sectional direction (Parker et al., 2011;Asahi et al., 2013).
Although researchers tend to think about bank erosion when considering channel migration, the average preservation of width over the long term in migrating channels forces an equal emphasis on the poorly studied process of bank accretion. Land accretion affects the shape of rivers over long periods of time. In simple meandering channels, sediment deposition takes place on a point bar located at the convex part of the river bank. As time passes, this area becomes higher in elevation because of that deposition. As this area becomes higher, it is submerged progressively less often. Vegetation can grow in this area. Vegetation works to provide flow resistance (decreasing flow velocities) and catches fine sediment during floods, so the elevation continues to increase. With time, this area achieves the same height as the floodplain and is submerged only during floods (Figure 20, from Asahi et al., 2013). In estimating the amount of bank shift by land accretion, it is critical to determine the area that rarely submerges. Specification of this area requires specification of an appropriate discharge to ascertain the inundated area. However, it is not sufficient to simply assume all inundated area is subsequently accreted (i.e. becomes land) because the portion accreted depends on several complex processes, including details of the horizontal and vertical distribution of sediment size, the type or mix of vegetation, and even the local climate. The relationship between those factors and the land accretion phenomena is complicated and difficult to treat in a precise manner. To avoid these issues, we defined a new parameter f land which indicates the amalgamated result of the land accretion progress. If f land is equal to 1, then land accretion occurs in all the rarely submerged area. If f land is equal to 0, then no land accretion occurs. When the water depth goes below a specified minimum depth (h min ) and continues to stay below that value, the area is treated as dry. We compute the accreted portion of the inner bank as the multiple of f land with the width of dry area. As it is difficult to define f land in a process-based manner, we choose it empirically as a parameter. In this calculation, f land is set as 0.3 in accord with Asahi's study (Asahi, 2014). Although this is clearly a simplification, it reduces many complicated processes to a single parameter that varies between 0 and 1.
To treat channel cutoffs during the evolution process, the model domain searches for any crossing or meeting point of riverbanks. When an intersection is found, the part of the bank that is not part of that intersection is realigned (typically the outer bank in two bends) into a single bank, leaving a cutoff in the form of an oxbow lake. Cubic-spline interpolation is used to ensure that the new bank is smooth. After this path is determined, the coordinate system is remeshed, a new channel centerline is computed, and the bed elevation is interpolated onto the new grid coordinates. After this process is complete, channel evolution continues as the model recomputes hydraulic conditions using the new geometry (Asahi et al., 2013). Figure 21 shows migration of a channel with the calculation conditions specified as shown in Table I. The initial channel shape was a single sine-generated channel with a varying discharge hydrograph with just two representative values (essentially bankfull and base flow discharges). The calculations show a complex meandering pattern with streamwise and temporal variations in width. A similar meandering pattern can be observed in Omolon River, Siberia, as shown in Figure 22. As discharge variation plays an important role in formation of the meander shape, the effect of the amount of time to become land (T land ) including time for vegetation to grow, and the flood recurrence interval (return of flood condition is represented by T return ) is shown in Figure 23. The calculation conditions for the four cases shown in Figure 23 are shown in Table II. The calculation results show that with increase in the T land /T return ratio, the channel becomes wider, less sinuous and less prone to cut offs. This happens as the bank has more time for deposition and creation of floodplain, including the growth of vegetation, so the bank becomes less prone to erosion when the flood returns. The time for process of land accretion plays a vital role in formation of the channel shape.
Cantilever failure model for overhanging banks in curved channels The channel migration model described earlier treats bank erosion primarily in a parametric manner using user-specified timescales for the processes of accretion and erosion of banks. Improvement of this treatment requires a more physically based approach, but such a development is complicated by the presence of complex bank geometry including overhanging banks. To begin to address this problem and with an eye toward future development of more detailed methods for bank treatments, a model for failure of overhanging banks has been developed based on the concept of cantilever failures. Cantilever failures are the result of fluvial erosion in the lower part of cohesive riverbanks. As fluvial erosion progresses an overhanging block is developed in the upper part of the bank that eventually falls into the water. The dropped material (slump block) remains at the toe of the bank until it is removed by the fluvial processes of erosion. This process is repeated over and over; it has been suggested that this process provides a limit to bank erosion rates (Parker et al., 2011;Motta et al., 2014;Hackney et al., 2015). The numerical model presented here includes cantilever failure but also treats the removal process of the failed blocks, which has not been directly considered in previous models. The bank erosion model used for this work is developed and described by Patsinghasanee et al. (2018) based on an experimental case in a straight flume. The experiment considered a rectangular section with one cohesive bank that was a mixture of silt-clay with sand. The experiment and model considered steady uniform flow with scale and sidewall corrections. The physical mechanisms observed in the development of cantilever failures in the experiment were fluvial erosion, tension cracks, cantilever failure, slump block generation and decomposition, and sediment transport away from the banks (Patsinghasanee et al., 2018).
The model considered a triple grid approach, using a onedimensional coarse grid for the flow field, a one-dimensional fine grid for the bed deformation and sediment transport, and a two-dimensional fine grid for the cohesive banks (including the generation of the overhanging block and disintegration of slump blocks). More detailed information can be found in the publication of Patsinghasanee et al. (2018). To include the effect of curvature, the previous model was coupled with the flow field model of Garcia et al. (1994), which corrects the velocity profile for curvature effects. In addition, the effect of curvature-driven secondary flow was considered using the equation of Hasegawa (1984). The coupling process of the two models is described in the paper of Arnez et al. (2017). The numerical model was tested using both cohesive and non-cohesive banks.  Note: Discharge and timescales are given in Table II.

Y. SHIMIZU ET AL. 26
For non-cohesive banks, the numerical results were compared to the experimental results of Fukuoka et al. (1983), which were developed in an annular flume with a central angle of 305°with finite length. The equilibrium state was observed in the section located at 270°.
The simulation considered an infinitely long circular channel (wherein the equilibrium state should be observed) in order to compare the results with the experiments. The initial geometry of the channel corresponds to a trapezoidal section where the bank height is equal to 0.06 m, the inclination angle of the banks is 30°and the width is 0.05 m. The averaged grain size is 0.72 mm and the bed slope is 1/400. Figure 24 shows the results obtained with the numerical model at 60 and 180 minutes.
Three cases were simulated: Case 1 corresponding to the full coupled model; Case 2 ignores the effect of the correction in the velocity profile; Case 3 neglects the effect of secondary flow. The results showed that for all three cases there was an underestimation of the fluvial erosion in the outer bank compared to the experimental case. Nevertheless, the full coupled model showed the best agreement with the experimental results. For the cohesive case, we qualitatively analyzed a hypothetical circular channel with a radius of curvature of 1.5 m. The material of the bank was a mixture of sand and silt-clay (30%). The initial cross-section was rectangular with a height of 0.2 m and a width of 0.3 m. The mechanical properties of the bank material were similar to those used in the experiments of Patsinghasanee et al. (2018). The value of the discharge was constant and equal to 6.45 (l/s) and the bed slope was 1/1000. The results presented in Figure 25 show that the model was able to represent the different stages observed in cantilever failures. For the case of non-cohesive material, the discrepancies of the experimental and numerical results can be explained by the assumption that the virtual channel is infinitely long while the real channel in the experiment was limited in length. Nevertheless, the effect of the curvature and correction of the velocity profiles was shown by the asymmetric fluvial erosion in the results (towards the outer bank).
In the case of cohesive banks, the model qualitatively reproduced the process in the formation of overhanging blocks, their failure and the decomposition of slump blocks. The model presented here is powerful tool for analyzing cantilever failures. However, the model is limited by the fact that bedforms such as bars cannot be represented, so this approach needs to be coupled to a full morphodynamic simulation, including a more accurate flow field model such as a two-dimensional, depthaveraged flow model. With that coupling, application of this approach to real rivers should be straightforward.

Vegetation effects on channel morphology
Vegetation growing on mid-channel bars and floodplain plays an important role in channel evolution. With respect to hydraulic and environmental engineering, understanding the effects of physical parameters associated with vegetation on river morphodynamics is of great importance. Here, we describe the computational approach considering vegetation effects on both flow and morphology in a two-dimensional numerical model (Nays2DH within the iRIC interface).
Flow velocity and structures are affected by vegetation density, type, age, submerged and emergent conditions, etc. (i.e. Kouwen and Li, 1980;Nepf and Vivoni, 2000;Jang and Shimizu, 2007). To consider those vegetation effects on flow structures, we use the standard shallow-water conservation equation with an additional body force term to represent the drag of vegetation on the flow. Thus, vegetation effects are introduced as a vector drag force with component F vx /ρ and F vy / ρ as follows: where C D is the drag coefficient of vegetation, α s is the vegetation density per unit volume, n s is the number of vegetation stems, d s is the average diameter of stems or trunks, A s is the sampling grid width (Figure 26), h s is the minimum value of water depth or vegetation height to take into account that vegetation height varies with age. If vegetation height is greater than the calculated water depth, h s is equal to water depth of h. In the opposite case, h s is equal to vegetation height. Table II. Calculation condition for discharge and both T land and T return for the results shown in Figure 21 (Case 1) and Figure 23  Vegetation also has an influence on channel morphodynamics. Highly concentrated vegetation controls flow resistance, alters channel depth and width, changes the flow thalweg, alters suspended sediment accumulation, and reduces the number of active channels (i.e. Tsujimoto, 1999;Tal and Paola, 2010;Van Dijk, 2013;Akahori and Kasugai, 2014;Kyuka et al., 2017). In addition, vegetation strengthens banks by binding sediment particle by their roots (i.e. Ikeda and Izumi, 1990;Millar and Quick, 1993;Gran and Paola, 2001). These vegetation effects on bed morphodynamics are generally modeled as reduction of sediment transport rate associated with changes in the flow resistance and the effective angle of repose (i.e. Shimizu, 2005b, 2007;Iwasaki et al., 2015). As a simple model, bank vegetation effects are quantified to the angle of repose in the bank erosion model based on the concept of mass failure (Hasegawa, 1984) to simulate channel widening due to bank failure. Furthermore, dynamic vegetation processes including vegetation invasion, growth and removal are also important systems for channel morphodynamics. The authors have reported several ways for estimating the dynamic vegetation effects (i.e. Asahi et al., 2013;Uchida et al., 2015;Nagata et al., 2016;Kang et al., 2018), although this approach has not yet been incorporated into the iRIC Nays2DH model. Jang and Shimizu (2007) compared channel evolution without and with vegetation by means of the laboratory experiments and computational simulation as shown in Figures 27  and 28. The numerical model simulates well the scour hole where the main stream or thalweg is bordered by vegetation, and longitudinal deposition in the vegetated areas. In addition, the model showed reduction of the averaged channel width as the vegetation density increased, as illustrated by Figure 28.
The interactions between vegetation and channel evolution are important topics for both river engineering and environmental perspectives. The current numerical model accurately simulates the experimental situation. For further progress, it is important to develop more complex models that include natural river characteristics such as non-uniform material and dynamic vegetation effects over a range of vegetation types.

Tidal channel networks
Although a great deal of attention has been paid to the temporal evolution of river channels, tidal channels have received relatively little attention. As in the case of terrestrial landscapes and the associated channels, there is a large variety of landscapes and channel morphologies in tidal environments. The landscapes in tidal environments are affected by a larger number of natural processes (tide, wave, river flow, etc.), so that their dynamics and associated ecosystems are quite unique, but also vulnerable to environmental changes (e.g. climate changes, sea level rise, human impacts) (Tambroni and Seminara, 2006;Deheyn and Shaffer, 2007;Ganju et al., 2017). One of the most distinct morphological signatures in these zones is tidal channel networks (Figure 26a) (Allen, 2000;de Swart and Zimmerman, 2009). Tidal channel networks are commonly found in salt marshes and tidal flats, thus have major roles in controlling overall hydrodynamics, sediment and nutrient transport as well as habitat structures within the system (e.g. Wallace et al., 2005).
Development of reliable numerical models for better understanding and prediction of the morphodynamics of tidal channel networks and associated environmental feedbacks has been required to estimate and evaluate the impact of the environmental changes to the tidal systems (Fagherazzi et al., 2012;Coco et al., 2013;Zhou et al., 2017). In the last two decades, several numerical morphodynamic models have been proposed (e.g. Hibma et al., 2003;D'Alpaos et al., 2005;Marciano et al., 2005;Dissanayake et al., 2009;Tambroni et al., 2010;Xu et al., 2017). These models capture the fundamental properties of the channel networks (e.g. probability distribution of unchanneled length [Marani et al., 2003, see also for example D'Alpaos et al., 2005); however, there have been some difficulties in validating the model performance in terms of the development processes of a channel network because (1) field observations regarding initiation processes and long-term evolution of the network system are generally not available, and (2) controlled laboratory experiments of tidal landscapes had not been performed with exception of the experiments by Reynolds (1889Reynolds ( , 1890Reynolds ( , 1891. In fact, until recently there was a lack of knowledge concerning how to perform meaningful small-scale experiments or to transfer experimental observations to the real world in the context of morphodynamics of tidal channels. Because of the lack of prior work and the need for information, the last two decades have been a critical time for proposing novel laboratory experiments to understand the basic physics of tidal morphodynamics (Tambroni et al., 2005;Vlaswinkel and Cantelli, 2011;Stefanon et al., 2010Stefanon et al., , 2012Kleinhans et al., 2015). This recent work has greatly contributed information for  validating numerical models (e.g. Tambroni et al., 2010;Iwasaki et al., 2013;Zhou et al., 2014). Herein, we demonstrate the results of Iwasaki et al. (2013) as one of the examples of experimental and numerical modeling of the initiation and development processes of tidal channel networks. Iwasaki et al. (2013) experimentally modeled the initiation and development of a tidal channel network in a small flume (length of 1.8 m and width of 0.9 m) by downscaling the conditions of Notsuke marsh, Japan (prototype, see Figure 29a) using the dynamic similarity rule developed by Tambroni et al. (2005). The tidal fluctuation, which is the seaward boundary condition, is modeled as a sinusoidal water surface fluctuation with the amplitude of 7.5 mm and the period of two minutes. A light material, polyvinyl chloride powder, having a diameter of 0.12 mm and a density of 1480 kg/m 3 , was used as a bed material to obtain sufficient sediment mobility. Figure 29b shows the modeled equilibrium channel network obtained after 50 tidal cycles, exhibiting initial formation of some main channels from the seaward side and subsequent development of complex network features.
The observed morphodynamics of the channel network is reproduced by a numerical model. The model employs a depth-averaged shallow-water flow model with the use of an equilibrium bedload sediment-transport model (Ashida and Michiue, 1972) with a local slope effect (Hasegawa, 1989). Similar approaches have been utilized in other types of the numerical models for reproducing the morphodynamics of tidal channels (Tambroni et al., 2010(Tambroni et al., , 2017Zhou et al., 2014). A square computational grid with grid spacing of 2.5 mm by 2.5 mm is used in the computation. Figure 29c demonstrates the numerical result of the channel network. The simulated channel network is qualitatively in agreement with the experimental result in terms of the complexity of network shape. To further evaluate the model performance quantitatively, Iwasaki et al. (2013) used nearest edge distance, which is the shortest distance from any point in an unchanneled region to the nearest channel observed on a photograph (Edmonds et al., 2011), rather than unchanneled length (Marani et al., 2003) to quantify the network geometry. Figure 29d shows the probability density distribution of nearest edge distance in the experiment, indicating an exponential distribution in the same way to the unchanneled length. This feature is also recognized in the prototype (i.e. Notsuke marsh, Figure 29e), and the numerical model reasonably reproduces the statistical distribution of nearest edge distance observed in the experiment. The discrepancy between the simulation and experiment is seen in the range of large nearest edge distance, meaning that the simulation result has a larger probability for long nearest edge distance compared to the experiment. This indicates that the computational model has some difficulties in resolving smallscale channels, which is probably attributable to the computational grid resolution.
Many challenges remain for the study of tidal landscapes, particularly in the context of numerical modeling. Although the similarities and differences between terrestrial and tidal landscapes have been investigated and discussed largely (e.g. Marani et al., 2003;Finotello et al., 2018), the computational efforts for understanding the tidal landscapes still have not reached a level comparable to that found in the terrestrial counterpart. Several outstanding contributions on numerical tools have been found in the literature, such as that for sand bars in tidal channels (e.g. van der Wegen and Roelvink, 2008;Tambroni et al., 2017). We will need to extend the numerical model for other morphodynamic signatures in the tidal zones such as tidal meandering, gully erosion, and network development. To complement the brief summary of morphodyamic modeling of tidal landscapes we introduced here, the reader can also find future challenges in recently published review articles (e.g. Fagherazzi et al., 2012;Coco et al., 2013;Zhou et al., 2017).

Bedrock abrasion
Abrasion by saltating bedload particles is one of a few major mechanisms of bedrock channel erosion Dietrich, 2001, 2004;Finnegan et al., 2007;Whipple, 2007, 2010;Chatanantavet and Parker, 2009;Inoue et al., 2014). Sediment supply has two opposing effects on bedrock erosion: a 'tools effect' that increases as the number of bedload particles colliding with (and abrading) the bedrock increases, and a 'cover effect' that protects bedrock from abrasion as the amount of sediment deposited on the bedrock bed increases (Gilbert, 1877;Sklar and Dietrich, 1998). Therefore, understanding the spatial distribution of both alluvial cover and bedload transport is necessary for simulating river bed deformation in mixed bedrock-alluvial channels (Nelson et al., 2014). Figure 30 shows photographs of three bedrock rivers in Japan . The alluvial cover fraction decreases in order from the Shikaribetsu River to the Tachiya River and then to the Abashiri River. The number of longitudinal grooves observed in these rivers increases in the same order. These qualitative observations suggest that the number of longitudinal grooves is affected by the areal fraction of alluvial cover and the formation of alternate bars .
There have been a variety of modeling approaches applied to sediment-transporting bedrock channels. Hancock and Anderson (2002), Nelson and Seminara (2011), Finnegan and Dietrich (2011), and Limaye and Lamb (2014) conducted numerical models for predicting the bedrock valley formation; to determine valley type (narrow canyons or wide valleys with strath terraces), the ratio of vertical to lateral incision rate was numerically tested. Nelson et al. (2014) and Inoue, Parker, and Stark (2017a) proposed numerical models for simulating the spatial distribution of alluvial cover and bedrock abrasion in bedrock channel bends. The simulation results show that the formation of a point bar affects the shape of bedrock channel bends. Hodge et al. (2011) reproduced the formation of sediment cover in a straight channel by simulating the movement of individual bedload particles. However, these models did not reproduce multiple incisional grooves or free bars in mixed bedrock-alluvial rivers. Inoue et al. (2016) proposed a numerical model to predict the incisional/alluvial morphodynamics of a bedrock channel and investigated the relationship between bedrock erosional forms and sediment supply using the proposed model. Inoue et al. (2016) applied the two-dimensional depth-averaged model to simulate flow field. The alluvial layer thickness was calculated from a sediment continuity equation that includes both sediment deposited on the bedrock and sediment moving in the bedload layer (Luu et al., 2004;Inoue et al., 2014). The bedrock abrasion model proposed by Chatanantavet and Parker (2009) was used in two-dimensional form to calculate the bedrock erosion rate. The areal fraction of alluvial cover was predicted by a ratio between the alluvial layer thickness and the bedrock macro roughness in each calculation cell (Inoue et al., 2014;Zhang et al., 2015). In their simulations, the sediment supply to transport capacity ratio is large in the order of Run A (1.0), Run B (0.7) and Run C (0.2). Additional calculation conditions for the experiments are shown in Table III. The sediment transport capacity (0.002 m 2 /s) is estimated using the bedload transport formula of Fernandez Luque and van Beek (1976). In the model, small random perturbations are introduced on the initial planar bed along the section 0 km from the upstream end. For all cases considered here, the flow was constant at a bankfull value, so no variations of discharge were considered. Figure 31 shows the alluvial layer thickness and the bedrock erosion depth (reprinted from Inoue et al., 2016). The simulation results show that different erosional forms are formed depending on the characteristics of the sediment supply. When the sediment supply/capacity ratio is unity (Run A), alternate bars form on bedrock surface. These bars migrate downstream similar to alternate bars that form in a purely alluvial river. The bedrock is eroded as the pools (in which bedrock is exposed) migrate downstream. The alternate bars 'sandpaper' the bedrock surface to yield a form that is, on average, transversely symmetric, and is deeper toward the banks than along the center . The concurrent formation of alternate bars and a symmetric, incised bedrock surface results in the formation of a single meandering thalweg. The erosional bed morphology changes with decreasing sediment supply (Runs B and C). In these cases, the sediment supply is insufficient for the formation of alternate bars, so that the bedload concentrates in multiple strings due to small random perturbations set near the upstream end. This concentrated transport erodes multiple grooves in the bedrock surface. The number of multiple grooves increases and the spacing between grooves decreases with a decrease in the sediment supply.
In a natural bedrock river, turbulence-driven secondary flow (as described earlier) may affect the number of bedrock longitudinal grooves. This three-dimensional flow structure is not con- Exposed bedrock (P c = 0) Exposed bedrock (P c = 0) Exposed bedrock (P c = 0) sidered in the depth-averaged model proposed by Inoue et al. (2016). Experiments and three-dimensional simulations on groove spacing represent an exciting challenge for the future. Moreover, the model proposed by Inoue et al. (2016) assumes fixed, non-erodible banks. Fuller et al. (2016) showed that abrasion by bedload is an overriding factor for bank erosion. Migration of alternate bars may affect the lateral erosion rate of bedrock channel. In addition, some recent experiments and field studies have proposed wet-dry weathering increases the bedrock bank erodibility for saltation abrasion (Small et al., 2015;Inoue, Yamaguchi, and Nelson, 2017b). To analyze the morphodynamics of incised meandering valleys, future model development should incorporate these phenomena. Similarly, the differential effects of saltation impacts on the bed and banks appears to be a key ingredient to understanding many of the characteristics of bedrock channels, and future work will be directed at combining mechanics of bed and bank erosion to predict overall channel evolution for bedrock systems.

Future work
Researchers are only just beginning to scratch the surface of the fascinating problem of channel evolution. Even the best meander evolution models use highly parametric treatments of bank erosion; interactions with vegetation, bank cohesion effects on failure modes, variations in groundwater levels, and the degree of groundwater-surface water interactions are central players in this problem; however, coupling detailed physics at this level with large-scale channel evolution models is only just beginning. Even for the relatively well-studied case of meander evolution, there are subtle variations in planform channel geometry (meander shape) that are not well understood and need further investigation. Even more than in the case of bars and bedforms, channel morphodynamics is strongly coupled to riparian habitat and vegetation across a range of timescales. For example, the recent work of Kimura and Kitazono (2017) on driftwood movement clarifies the role of driftwood as a first-order effect in bar and bank morphodynamics. Observations also suggest that vegetation can have a first-order effect on the sediment-transport and erosion/deposition field that is not captured by simple drag models such as the one described earlier and should be treated with a more general view of the interaction of the vegetation with the turbulence field. Treating this and other similar problems will require fully coupled models of vegetation and the flow field. As this research field moves toward more realistic treatments of co-evolving morphodynamics and vegetation, groundwater models that incorporate physically correct connections to the channel flow will also be required. The iRIC group is working on these developments as well as trying to develop collaborative relationships with botanists and riparian ecologists to work together on better characterization of channel-vegetation couplings.
As in the case of both bars and bedforms, perhaps the greatest challenge is how to put the various morphodynamic scales together without creating a computational nightmare. We believe that the commitment of various scientists working on these problems to make codes and documentation available in the public domain is a key part of continued progress. As cloud computing and model coupling tools become more readily available and easier to use, community efforts to combine approaches at a variety of scales are likely to lead to future predictive treatments for real systems that are currently impractical or impossible to simulate.

Conclusions
The field of computational morphodynamics for natural flows is fascinating and covers a wide range of topics, as demonstrated from the limited examples mentioned earlier. Even though a good deal of progress has been made in terms of basic understanding and practical computational methodologies, the field is relatively young and has many future challenges. We believe that the continued development of scientific interaction among morphologists and a similar commitment to sharing computational code and interface tools, as exemplified by the iRIC project, are both fundamental parts of the continued development of this field. The examples shown here, although interesting and challenging in an academic sense, are also directly related to practical problems in rivers and other similar natural and manmade flows. The iRIC interface and solvers are freely available at www.i-ric.org, and we encourage students and practitioners to use these resources and encourage specialists in morphodynamics and related fields to make their own methods available within iRIC.