Morphodynamic Model Suitable for River Flow and Wave-Current Interaction

Morphodynamic models are a great support in water environment management and decision-making, as well as in integrated coastal zone planning. In the present paper, a 2DH model is presented, able to deal with both currents, waves and their mutual interaction. The model is briefly presented and the results of its application to some benchmark tests are discussed.


Introduction
Morphodynamic models are numerical models able to couple the hydrodynamic issue to the sediment transport analysis and to the morphological evolution that represents the changes in bottom height due to sediment transport. The models are a valuable tool in hydraulic and morphological studies, resulting particularly useful in comprehending and predicting sediment movement and the following morphological response in rivers, estuaries or coastal environments. For this reason, their use is becoming more and more common in support of water environment management and decision-making, as well as in integrated coastal zone planning. Two main approaches for the analysis of the sediment transport problem are known in literature: equilibrium and non-equilibrium approaches. The former assumes that sediment transport corresponds everywhere to the actual sediment transport capacity, while the latter studies the distribution of sediment concentration, solving a proper advection-diffusion equation. With particular regards to the modelling of coastal areas, the main contributions to the morphological evolution are due to currents, waves and their mutual interaction. Moreover, the morphological evolution in lower river courses can be particularly affected by the secondary flows, which take place in presence of curvilinear flows, typical for example of bends and meanders. A model, suitable to be applied to this context, should consider all these aspects.
In literature, several different morphodynamic models have been presented and discussed. In particular, the most useful ones in coastal applications are bidimensional shallow water models.
CMS (Coastal Modelling System), developed by the U.S. Army Corps of Engineers [1], [2] and MIKE, developed by the Danish Hydraulic Institute [3] are among the most popular commercial coastal morphodynamic models. Moreover, Open-telemac [4], as an open-source modelling suite, is also widely used. Among other models, recently presented in literature, are those by Ferrarin et al. [5] or Carniello et al. [6] or Begnudelli et al. [7].
Despite the large number of available models, some of them neglect the interaction between waves and current or they consider it only in a limited way. Others disregard secondary flow effects. Only a few models consider both aspects. In the last few years, 3D modelling has also been profoundly developed; however, 3D-models are still too heavy from a computational point of view to be applied to complex and wide domains.
With a view to developing a morphodynamic model able to be applied to coastal areas like river mouths and lagoons, a finite volume 2DH morphodynamic model is presented in the present paper as a first step, which can consider the effect of the regular wave field in both hydrodynamics and morphodynamics. Moreover, secondary flow effects on bed load are also considered.
The model is based on an accurate shock-capturing and C-property preserving scheme, also able to describe shock phenomena properly [8]. Finally, the present model has been developed in order to reach a balance between accuracy and simplicity: in fact, the model is suited for the study of complex 2D domains, with a minimum number of physically based parameters to be considered in sediment transport problems. In this work, a few different theories have been implemented and applied to some benchmark tests, in order to carry out a comparison between them, which permits one to focus on the most appropriate approach in each single case study.
In particular, the analysis focuses on morphodynamic evolution under the combined effects of current and waves, propagating in the same direction, and that under curvilinear flow.

Numerical model
The model presented in this paper is an extension of that adopted in Bosa et al. [8] in the case of current and wave motion acting simultaneously. The granular sediment transport as usual is split into bed and suspended load, where the former is dealt with an equilibrium approach and the latter with a nonequilibrium approach. Thus, a depth-averaged advection-diffusion equation associated with the concentration of the suspended load is added to the hydrodynamic equations to study the distribution and the propagation of sediments in the flow.
The resulting 2DH shallow water equations system is then: with: , , , ; = 0, , , ; , , , ; = 0,2 , , ; ( , , , In these equations, (t, x, y) are temporal and horizontal spatial coordinates and U variable vector, with h the water depth, (U, V) the mean velocity over the depth in x-and y-directions and C depthaveraged concentration of granular suspended sediment. (F, G) and (F t , G t ) are vectors of advective and turbulent fluxes, with g the gravity acceleration and  t horizontal eddy viscosity coefficient evaluated by means of a Smagorinsky approach. S is the source term vector, with  the water density and z b the bottom height; ( bx ,  by ), ( sx ,  sy ) and (f rsx , f rsy ) are the components along x and y of the mean bed shear stress  b , the surface wind stress  s and the forcing due to wave radiation stress f rs . E -D represents the erosion or deposition rate of the suspended load.
Mean and maximum bed shear stress due to the combined action of current and waves are evaluated through Soulsby's 13-parameters regression formula as functions of current shear stress  c and maximum wave shear stress  w [9]. In the present paper, the regression formula has been employed to represent the following models: Bijker, Grant and Madsen, Fredsøe, Huynh-Thanh and Temperville and Myrhaug and Slaattelid. Moreover, an additional set of parameters called DATA2 introduced by Soulsby [10] has also been implemented, which derives from fitting the regression formula on a wide experimental data set. In this context, current and wave bottom shear stress have been evaluated through Manning and Jonsson formulae respectively: k s being the Gauckler Strickler coefficient, f w the wave friction factor and U b the maximum bottom velocity amplitude, computed using the linear wave theory. The mean bed shear stress is supposed to be responsible for current velocity and hence involved in the momentum and advection-diffusion equations. The maximum bed shear stress determines the threshold of motion and the entrainment of sediments and therefore it has been used in the evaluation of E -D.
The bed level change is evaluated by means of the sediment continuity equation: n being sediment porosity and q b bed load transport, using an equilibrium approach and the formulae proposed by Van Rijn [11], Soulsby [10] and Meyer -Peter and Müller [12]. A transverse and a longitudinal component are added to the bed load to take into account transverse and longitudinal bed slope effect [13]. The source-sink term of transport and sediment continuity equation E -D is evaluated as: , (7) w s being the sediment fall velocity evaluated as in Soulsby [10]. Assuming as the reference height the upper end of the bed load layer, called a, the equilibrium concentration at this height is c a and the actual concentration at the same height is c 0 . c a is evaluated as in van Rijn [11], while c 0 is calculated from depth average sediment concentration C assuming a particular vertical concentration profile. In the present paper, van Rijn [11] and Soulsby [10] profiles have been implemented. The former considers, for the mixing coefficient, a parabolic-constant vertical distribution for currents ( c ) and constant-linearconstant vertical distribution for waves ( w ). The resulting mixing coefficient is due to their combination: . (8) The relationship between C and c 0 is based on an exponential concentration profile employing a depth-averaged mixing coefficient; thus, c 0 can be evaluated from C as: Soulsby considers a linear mixing coefficient vertical distribution, with a change in the vertical slope over the wave boundary thickness . The resulting vertical concentration has a power law profile and c 0 can be evaluated from C as: , being with u *m and u *max the friction velocities evaluated through mean and maximum bed shear stress, respectively. As a comparison, the equilibrium approach for suspended load has also been implemented. In this case, the advection-diffusion equation associated to the concentration of the suspended load is not used, the bed load transport q b in eq. (6) is substituted by the sum of bed and suspended load transport (q b + q s ) and the souce-sink term E -D vanishes. q s is evaluated using the van Rijn [14] formula.
Equations (1) and (6) are integrated over time with a Strang splitting scheme and in space with a finite volume method, based on HLLC Riemann solver associated to a hydrostatic variable reconstruction that assure the scheme is well balanced also in wet and dry conditions [15], [16], [17]. The turbulent fluxes and the bed load at the intercell are evaluated by means of a centred finite difference scheme. The resulting scheme is first order accurate both in space and time [18].

Secondary flow for morphological evolution
The present model has been developed to be applied to transitional environments in the future like river mouths, where fluvial and maritime phenomena coexist. Thus, typical river aspects should also be considered. For example, river bends cause an imbalance of the centripetal force acting on the fluid, which generates an outwards motion near the surface and an inwards motion at the bottom perpendicularly to the main flow, well-known as secondary flow. Thus, bends are characterized by a resulting helical flow, that has an important influence on bed topography [19]. In the present paper, the helical flow has been considered for the change in bed shear stress direction, which has an important influence on bed load. The algorithm here implemented is based on the Rozovskii theory [19], which admits a deflection of the bed shear stress from the main flow in presence of secondary flow of an angle  s :  being a parameter depending on the Gauckler Strickler coefficient and R s the radius of curvature of flow stream line. Secondary flow thus causes a further bed load component q bs transverse to the flow direction, which is then added to q b : tan , In regions where streamline curvature changes in space and time, secondary flow will adapt gradually [20]; thus, a transport equation for  s is considered: Equation (13) is finally solved using the same finite volume technique adopted for the concentration advection-diffusion equation.

Validation tests
The numerical model has been tested on laboratory experience, performed by van Rijn [21] and Yen and Lee [22] in different hydrodynamic conditions. A comparative analysis among bed load formulae and suspended load approaches is conducted, to check their reliability.

Trench migration in a flume
The trench migration test has been carried out by van Rijn [21] and it consists in the measurement of the bottom evolution of a trench perpendicular to the flow under steady current with and without wave motion (figure 1). In the present paper, two configurations were chosen among those originally presented by van Rijn; their geometric and hydrodynamic conditions are summed up in table 1 and sediments and morphodynamic characteristics are summarized in table 2. The former is characterized by current alone (T1 configuration) and the latter by current superimposed on a regular wave, propagating in the same direction (T2 configuration).

Trench migration under steady current (T1).
In the current alone test, in order to guarantee experimental hydrodynamic conditions, the inflow, k s and longitudinal bed slope were set to 0.1 m 3 /s, 50 m 1/3 /s and 0.4 ‰ respectively. Different computational methods for bed and suspended loads were considered and compared. Suspended load was studied analysing both equilibrium and non-equilibrium approaches, considering both van Rijn and Soulsby vertical concentration profiles. With regard to the bed load, Van Rijn, Soulsby and Meyer-Peter and Müller formulae were considered. As a sediment inflow boundary condition, an equilibrium suspended load was imposed, and the outflow was a cinematic condition. As mean grain diameter 160 μm was imposed and the reference height a was 0.0250 m, coherent with ripples height. Bottom porosity was set to 0.3.
In figure 2 the results in terms of bottom elevation after 7.5 h and 15 h are shown using the different above mentioned approaches for suspended load and van Rijn formula for bed load, together with van Rijn experimental data. The experiments show that the upstream slope of the trench advances while the trench bottom gradually fills and, instead, the downstream slope is eroded, beginning from the upper side. It can be seen that the equilibrium approach tends to overestimate trench migration, while the filling process is actually underestimated. In case of a non-equilibrium approach, the slopes migration is more gradual and the entire trench is progressively smoothed and filled. The trench shape is then in better agreement with experimental data. However, when the van Rijn vertical concentration profile is used, the infilling process is generally quicker. Soulsby vertical concentration profile seems to be generally in better agreement with measured data. After that, the behaviour of different bed load formulas was investigated, keeping a non-equilibrium approach for suspended load, with power-law Soulsby concentration distribution. The results obtained are compared in figure 3. Van Rijn and Meyer Peter and Müller formulas result in better agreement with experimental data, being both filling process and upstream slope position more coherent. However, the Van Rijn formula appears to describe slightly better the upstream trench slope. The Soulsby formula gives less accurate results, with a less filled trench and a less migrated upstream slope. Thus, the Van Rijn bed load formula seems to be more appropriate for this test.

Trench migration under current and regular waves, propagating in the same direction (T2).
The second trench test considers the superposition of a regular wave (wave period 1.5 s and wave height 0.08 m) onto a steady current, propagating in the same direction. The inflow, Gauckler Strickler coefficient and longitudinal bed slope were set to 0.014 m 3 /s, 40 m 1/3 /s and 0.018 ‰ respectively, in order to guarantee experimental hydrodynamic conditions. The sediment transport was characterized by a mean grain diameter of 100 μm, a bottom porosity of 0.4 and a reference height a of 0.01 m, coherent with ripples height. As in the previous test, an equilibrium suspended load was imposed at the inflow boundary and a cinematic condition at the outflow. Due to particular hydrodynamic conditions, bed load played a very small role, nevertheless it was computed through the van Rijn formula. The suspended load was evaluated using both equilibrium and non-equilibrium approaches, considering both van Rijn Distance, x (m) and Soulsby vertical concentration profiles. As a first approach DATA2 formulas to evaluate mean and maximum bed shear stress due to the combined action of current and waves were used. A uniform wave field had been imposed into the whole channel; thus, the simulation did not take into account the effect of the radiation stress forcing. Considering the bed level after 10 h (figure 4), a non-equilibrium approach gave similar results, whatever vertical concentration profile is adopted. The migration of upstream trench slope results here less evident than in experimental data. This fact is probably caused by the almost total absence of bed load in the numerical simulation. In the experiments the presence of a non-uniform granulometry probably gave rise to a bed load transport, which is generally responsible for trench migration. This could reveal that the mean sediment diameter is arguably not representative for characteristic bed load diameter. Nevertheless, the filling process is adequately simulated by the model: the profile is close to experimental data in the central part of the trench and at the downstream slope. Downstream to the trench, the flume in the numerical model is slightly more eroded, but the profile is anyway parallel to experimental data and coherent with them. The equilibrium approach fails again in the description of the trench evolution, the trench migration being broadly overestimated and the filling much lower than in experimental measurements. To appreciate the effect of different formulas for the evaluation of mean and maximum bed shear stress due to the combination of wave and current, the same test was carried out employing the Soulsby non-equilibrium approach for suspended load and varying the theory used to calculate mean and maximum bed shear stress. In every configuration, the equilibrium inflow concentration was imposed. The trench profile after 10 h is depicted in figure 5 using different mathematical theories and experimental data. Bijker model results as less accurate than others, with regard to the filling process. Also in the Fredsøe, Myrhaug and Slaattelid and Huyngh-Thahn and Temperville models the simulated trench appears to be slightly less filled than experimental data. Conversely, the Grant and Madsen model overestimates the process, being the trench bottom slightly higher. With regard to the migration process, all models underestimate the upstream trench feed, while the lowest part of the trench is slightly ahead compared to the experimental results. The downstream channel appears slightly more eroded in all configurations, but coherent with experimental data. From the present comparison, DATA2 formula seems to better describe the physical process.

Channel bend evolution under unsteady flow
Yen and Lee [22] studied the effect of an unsteady flow on the morphological evolution of a laboratory channel bend, where secondary flow and transverse slope play an important role. The channel geometry is depicted in figure 8. The longitudinal slope was 0.002, the mean radius was 4 m and the channel width was 1 m. A 20 cm layer of sand with a mean grain diameter of 1 mm was placed on the bed before each experiment began. The initial flow was set at 0.02 m 3 /s, which corresponds to the incipient motion condition for a grain diameter of 1 mm. Yen and Lee analysed five experiments, with different triangular inflow hydrographs; for their features, please refer to the original paper.
In the numerical tests performed, the same geometry and hydrodynamic conditions of the original Yen and Lee experiments were adopted. In order to guarantee such conditions, the Gauckler Strickler coefficient was calibrated at 57.47 m 1/3 s -1 . At the inflow, the hydrograph of the corresponding experiment was imposed and, as an outflow boundary condition, a cinematic condition was set. Channel lateral walls were modelled with wall boundary conditions. Due to particular hydrodynamic and sediment conditions, the suspended load was not relevant and hence it was neglected in this test. As a first approach, the van Rijn formula for bed load was applied.
Due to the bed load deviation under the combined effect of secondary flow and consequent transverse slope, deposition is concentrated at the inner wall, with maximum values in the first part of the bend. Similarly, scour phenomena are concentrated at the outer wall, where bed material is gradually eroded and transported into the inner areas. Maximum scour is located in the final part of the bend. At the end of each simulation, the bend is characterized by a fairly constant transverse slope.
In figure 6 the results of run1 are shown as the ratio between erosion/deposition and the initial water depth at the beginning of the test and compared with experimental data. In the following, this adimensional variable is mentioned as z b '. Maximum erosion, maximum deposition and the general morphology obtained by numerical simulation are fairly much in agreement with the experiments. Yen and Lee observed that maximum deposition occurred near the inner bank at a 75° section from the beginning of the bend, while maximum scour took place near the outlet bank at a 165° section from the bend beginning (section1 and 2 respectively in figure 6).
In figure 7, the comparison between experimental data and numerical result at these sections is illustrated for each run. In section 1, the transverse slope appears to be less sensitive to the hydrograph in the numerical model than in the physical experiment. In particular, the transverse slope seems to be slightly underestimated in run1. However, the mean transverse slope seems to be in good agreement with experimental data. In section 2, deposition in run1, run2 and run3 is well described by the numerical model, although erosion is underestimated. On the contrary, in run4 and run5, deposition is underestimated by the model, while the erosion is well predicted in run 4 and somewhat overestimated in run 5. It should be considered that Yen and Lee used a sand mixture with a significant variance, which probably has a relevant effect on the bend morphology. In order to obtain a comparison of different bed load formulas for the description of bend morphodynamics, some more simulations have been considered using Soulsby and Meyer-Peter and Müller formulas. In particular, two runs are here reported, the first one (run1) and the last one (run5), considered to be representative of the two extreme configurations. The results are depicted in figure 8, in terms of z b ' on the two considered sections. Van Rijn and Meyer-Peter and Müller formulas gave similar results: in run1, scour is slightly underestimated while in run5 both scour and deposition are more accentuated. However, in both sections the transverse slope is in fairly good agreement with

Conclusions
With a view to developing a model useful in coastal and riverine environments, a morphodynamic model is presented in the present paper and applied to some validation tests. The work focuses on two main features of the model: the capability to consider the effects of wave field in hydrodynamic and morphodynamic problems on the one hand and the secondary flow effect on sediment transport on the other.
Wave field influences both hydrodynamics, through wave-current bed shear stress and radiation stress forcing, and sediment transport, using bed and suspended load theories developed for wave and current combination.
Secondary flow influences bed load magnitude and direction, causing deposition in the inner side of bend and erosion in the outer side.
A trench migration test was conducted with and without a regular wave, propagating in the same direction of the current. The model's results show a better agreement with experimental data using a non-equilibrium approach for the suspended load. Moreover, while in the current and wave configuration the Soulsby and van Rijn vertical concentration profiles gave similar results, the current alone configuration was better described by Soulsby profile. With regard to bed load, the van Rijn formula seems to have better results.
Finally, in the case of a wave parallel to the current, the comparison of different bed shear stress theories evidences that the Soulsby DATA2 formula, based on a regression of experimental data, gave better results than the others.
The channel bend evolution test under unsteady conditions describes a typical situation of the lowest reach of a river. In this test, sediment transport occurs mainly close to the bed and, for these reasons, the influence of a helical flow due to the main-flow curvature plays an important role. Van Rijn, Soulsby and Meyer-Peter and Müller bed load formulas were considered and compared.
Numerical results show that the model can adequately describe the general morphological evolution. Moreover, Van Rijn and Meyer-Peter and Müller formulas seem to better fit experimental data, while the Soulsby formula obtains a less accurate solution, in particular in correspondence to a maximum scour section.
The next step to developing a morphological model suitable to study coastal and riverine environments will be the coupling of the developed model with a spectral model.