Influence of Erodible Beds on Shallow Water Hydrodynamics during Flood Events

Flooding has become the most common environmental hazard, causing casualties and severe economic losses. Mathematical models are a useful tool for flood control, and current computational resources let us simulate flood events with two-dimensional (2D) approaches. An open question is whether bed erosion must be accounted for when it comes to simulating flood events. In this paper we answer this question through numerical simulations using the 2D depth-averaged shallow-water equations. We analyze the effect of mobile beds on the flow patterns during flood events. We focus on channel confluences where water flow and sediment mobilization have a marked 2D behavior. We validate our numerical simulations with laboratory experiments of erodible beds with satisfactory results. Moreover, our sensitivity analysis indicates that the bed roughness model has a great influence on the simulated erosion and deposition patterns. We simulate the sediment transport and its influence on the water flow in a real river confluence during flood events. Our simulations show that the erosion and deposition processes play an important role on the water depth and flow velocity patterns. Accounting for the mobile bed leads to smoother water depth and velocity fields, as abrupt fields for the non-erodible model emerge from the irregular bed topography. Our study highlights the importance of accounting for erosion in the simulation of flood events, and the impact on the water depth and velocity fields.


Introduction
The numerical simulation of flood events to model fluvial hydraulics and flood inundation areas has become a key tool to assess flood risks, effectively manage and control flood events, and evaluate practices on flood prevention, protection and mitigation-the so-called flood control. Flooding is the most common environmental risk [1], causing heavy life and economic losses. Moreover, climate change is expected to increase the intensity and frequency of extreme rainfall events that drive floods [2,3].
One of the most widely used approaches to model fluvial hydraulics at a river scale is based on the numerical solution of the one-dimensional (1D) Saint Venant equations [4]. This approach requires low computational resources, allowing us to simulate long river length. Nevertheless, the simplifying assumptions involve several drawbacks: (1) the flow is one-dimensional, (2) velocity is uniform across the cross-section, and (3) the water level is horizontal across the cross-section. Despite the simplifying assumptions of these approaches, in-channel flows are satisfactorily simulated by 1D models. However, out-of-bank flows or singular elements such as channel junctions or bridge piers incorporating two-and three-dimensional (3D) effects may not be satisfactorily described by 1D approaches [5].
Currently, the development of computational capacity has allowed us to simulate fluvial hydraulics with approaches based on the two-dimensional (2D) Shallow-Water Equations. These approaches have emerged as a robust and efficient tool for the simulation of flow on both natural and man-made environments [6]. Moreover, sophisticated 2D hydrodynamic codes have recently become available, as for instance, Delf3D-FLOW [7] or Iber [8][9][10]. The approach relies on the assumption of depth-averaged velocity and, in most cases, quasi-hydrostatic flow pressure, which involves neglecting vertical acceleration. Although this assumption has been identified as a source of error in flows in non-prismatic channels [11], over curved beds in 1D [12] and 2D simulations [13,14], or channel configurations with vertical water circulation [15], the approach is suitable for modeling the hydrodynamics of flood events.
While 2D Shallow-Water models may naturally incorporate well resolved bed topography on flow hydrodynamics, most numerical studies of flood events so far have neglected changes in bed topography on the water flow dynamics during flood events. Much effort has been focused on the long-term bed erosion patterns and the evolution of river morphologies over time scales of years or decades, such as the formation of bars in alluvial channels [16], the bedload transport in river confluences [17,18], the riverbank migration and island dynamics [19], or at a basin scale [20]. However, erosion during flood events may significantly modify bed topography [21,22], even leading to the formation of canyons by a single flood event [23]. The hydrodynamics, sediment transport processes and river morphology are strongly coupled in rivers and channels with erodible beds and cannot be studied in isolation [24]. Mathematical models should then interlink these three processes through fully coupled approaches.
A paradigmatic case of bed erosion on both natural and man-made environments are the river and channel confluences. Junctions of two open channels or rivers are common configurations present in many hydraulic structures such as urban channels [25], tidal channels [26], tidal river confluences [27], irrigation channels or river confluences. The morphological progress of river bathymetry is significantly influenced by turbulent flow structures associated with dunes [28], junctions of channels with different depths [29] or discordant channels [30][31][32], indicating that these structures may play a key role in the upstream erosion patterns under subcritical flows. Flow dynamics at channel junctions have been characterized with laboratory or field experiments, aimed at the characterization of the flow-separation zone [33] and the effect of curved profiles [34], the contaminant transport with lab experiments [35] or field studies [36,37], the development computational approaches [38,39], or providing benchmark experimental data sets [40]. Erosion patterns at junctions have also been experimentally characterized in channels with [41,42] or without [43][44][45] sediment supply.
Flow dynamics at junctions are characterized by six regions [43], as shown in Figure 1; of particular interest among them are the separation zone in the main channel intermediately downstream of the junction [33], and the shear plane between the two flows in the main channel [46,47]. In the separation zone the water depth depresses and a recirculating flow appears. The shear plane is characterized by high turbulence levels that generate quasi-two-dimensional vortices with near-vertical axes [48,49].  Steady-state bed bathymetry at junctions presents two salient morphological features [50]: a scour-hole on the downstream corner of the junction, and a bar on the left bank closed to the downstream side of the confluence, denoted as a bank-attached lateral bar. Scour hole depth varies non-linearly with the confluence angle [45,51,52]. The bar appears due to the lower water depth and recirculating flow in the separation zone, as well as the flow detaches from the left sidewall downstream of the junction. Moreover a helicoidal flow emerges in the main channel downstream of the junction [50], which is responsible for the movement of water in the vertical direction and the concentration of the finest sediments in the bar [45].
In this study we assess the effect of mobile beds on the flow patterns during flood events. Our case study is a river confluence in a semiarid basin located in southeast Spain. Our main objective is to determine whether erosion should be accounted for or not to simulate fluvial hydraulics during flood events. To validate our numerical model, we compare our simulated results with laboratory experiments of sediment transport at channel confluences. In particular, we focus on two laboratory setups: (1) a set of experiments without mobile beds to validate the hydraulic model, and (2) an experiment of dynamics of bed morphology to validate the erosion model. We test those predictions against the experiments and identify the sources of model discrepancy. Lastly, we take advantage of the mathematical model to evaluate the effect of mesh size and bed properties on the erosion patterns.

Materials and Methods
We use the freely available Iber code [8,53]. Iber is a 2D hydraulic model that solves the hydrodynamic equations coupled to the sediment transport equations, among others. In the following sections, we describe the hydrodynamic equations and the erosion model.

Two-Dimensional Depth-Integrated Shallow-Water Equations (SWE)
Iber solves the 2D SWE using a finite volume scheme. This scheme is able to handle unstructured meshes, irregular topographies, friction losses, and wet-dry fronts [53]. The 2D SWE are derived from the Navier-Stokes equations by assuming quasi-hydrostatic flow and incompressibility of water. The 2D mass conservation equation in a Cartesian coordinate system reads: where h is the water depth, t is time, (U x , U y ) is the depth-averaged velocity in the x and y-directions, and S q are the sources or sinks. The linear momentum conservation equation in a Cartesian coordinate system is given by: and ∂hU y ∂t where g is the gravity, Z b is the channel-bottom elevation, τ b is the bed friction, ρ is the density of the water, and ν t is the eddy viscosity. We neglect wind surface tractions and Coriolis force. τ b is computed through the Manning formulation as: where j = x, y, and n is the Manning's coefficient. The last two terms in Equations (3) and (4) are the turbulent stresses. They are determined from the so-called kturbulence model of Rastogi and Rodi [54]. The model consists of two differential transport equations for the turbulence kinetic energy, k, and the rate of its dissipation, . The turbulence kinetic energy is given by: ∂k ∂t and the rate of its dissipation by:

Erosion Model
Iber erosion and sediment transport module solves the non-cohesive sediment transport equations with uniform granulometries, in a non-stationary regime. The hydrodynamics, sediment transport processes, and river morphology interplay at mobile-bed confluences [55,56] and cannot be studied in isolation [24].
Iber software couples these three processes. The sediment transport process is simulated through the 2D Exner equation [57][58][59]. The equation provides the bed elevation evolution in response to the erosion and sedimentation and is given in the form of bedload and suspended load. It reads: where Φ is the porosity, q b,j is the bedload discharge along the j direction, ω s is the fall velocity of the suspended sediments, E s is the dimensionless factor of the inlet of suspended sediment, and C b is the concentration of suspended sediment on the bed. Several formulations are available in Iber to compute the bedload solid flow rate due to the bedload transport. We have conducted our simulations with the van Rijn equation [60], given by: where q * b is the dimensionless bedload solid volumetric flow rate, D * is the dimensionless diameter of solid particles, τ * is the dimensionless shear stress acting on the particles, and τ * crit is the dimensionless critical shear stress.
The dimensionless bedload solid volumetric flow rate, q * b , reads: being q * b the bedload solid volumetric flow rate, s is the ratio between the densities of the particles and water, s = ρ s /ρ, and D is the characteristic diameter of sediments, usually taken as the median diameter D = D 50 .
The dimensionless diameter of solid particles, D * , reads: where ν is the kinematic viscosity of water. The dimensionless shear stress, τ * , has the form: where τ b is the bed friction shear stress, given by Equation (4). The solid bedload flow rate equation is combined with the Shields-stress criterion for incipient motion [61] and a bed slope correction for bedload transport. Moreover, Iber simulates the suspended sediment transport with the depth-averaged convection-diffusion equation for the sediment concentration. The model accounts for the turbulent diffusion and the equilibrium suspended concentration. The latter, together with the sedimentation velocity, are computed with the van Rijn formulation [62]. In the following, we do not consider the suspended sediment transport since it is negligible compared with the bedload sediment transport.

Solution of the System of PDEs
We solve the system of Partial Differential Equations (PDEs) with a second-order numerical scheme based on the Monotonic Upwind Scheme for Conservation Laws (MUSCL) scheme [63]. High-resolution finite volume methods [64][65][66][67][68] are suitable to minimize numerical dissipation and capture complex flow patterns [69].

Validation of the Hydraulic Model
We test our hydraulic model by comparing the numerical outputs with laboratory experiments reported on the literature [70]. Although the Iber code has been previously validated with experiments of open-channel contractions with several geometries and flow regimes [71], here we compare the outputs of the model with a set of lab experiments of open-channel junctions [70]. The experiment setup is a rectangular channel confluence of 0.30 m wide, as shown in Figure 2a. The channel was made of acrylic plates and metal sheets. The slope in both channels was 0.14%. Water was introduced in the system through both inlets located upstream in the main and secondary channels. Water level was measured at five points. Here, we simulate the experiments under subcritical flow. The confluence angle θ between the main and secondary channels θ is 30 • and 60 • . The simulated experiments are listed on Table 1, as well as the flow rates and water levels measured at five points.
Laboratory experiments have characterized the flow behavior at the junction of rectangular channels [33,34,[38][39][40]. Briefly, flow dynamic is characterized by six regions [43], of which the separation zone in the main channel intermediately downstream of the junction [33], and the shear plane between the two flows in the main channel are of particular interest [46,47].
Our model is able to reproduce the main flow features at channel junctions, as shown in Figure 2b,c. The momentum of the water in the secondary channel makes the flow detach from the sidewall as water gets in the main channel, and water depth decrease in the vicinity of the wall. Our model replicates this effect, as shown in Figure 2b where a depression in water depth appears at the downstream intersection of both channels and extends along the downstream direction. Another side effect of the intersection is the increase in water depth upstream the junction in both channels. Our model reproduces this effect, Figure 2b.
The interaction between both channel streams generates a shear plane that divides the main and lateral flow. We identify this plane with our numerical results by plotting the x-velocity gradient, given by ∂U x /∂y, in Figure 2c. The gradient is high along the intersection of both water streams in the main channel. As expected, the gradient is especially high downstream of the junction, in the separation zone. However, since our model is 2D depth-averaged, it is unable to reproduce the helicoidal flow that appears downstream of the junction in the main channel [50]. The helicoidal flow is responsible for the movement of water in the vertical direction and may influence the erosion patterns.
We compare the output of our numerical model with the results of the laboratory experiment listed ih,n Table 1 in Figure 2d,e for the channel junction with θ = 30 • . The agreement between the five experiments and the simulations is remarkable. Such agreement is also satisfactory for the channel junction with θ = 60 • , depicted in Figure 2f,g. Nevertheless, the numerical results at point P5 slightly deviate from the lab experiments for this channel configuration. However, the agreement between experiments and numerical simulations at point P5 is very satisfactory for the channel junction with θ = 60 • .
(a) Experimental setup Inlet

Validation of the Erosion Model
Our next validation experiment was conducted by Yuan et al. [45]. The authors studied the dynamics of bed morphology in channel confluences through laboratory experiments under subcritical flow conditions. The setup consisted of a perpendicular confluence of rectangular concordant-bed open-channels, the dimensions of which are shown in Figure 3a. The inlet flow rates of the simulated experiment are Q m = 9 l/s and Q t = 6 l/s. The bed sediment was poorly sorted with D 50 = 0.9 mm, and D 90 = 2.5 mm, and the bed was initially covered with a 5.5 cm-thick layer of sediment. In our numerical simulations, we got the best agreement between experiments and numerical simulations with the following properties of the sediments: density of grains ρ = 2500 kg/m 3 , friction angle φ = 25 • , and Manning's coefficient n = 0.018 s·m −1/3 . The steady-state bed elevation contour is shown in Figure 3b, we remind the initial bed elevation was 5.5 cm.
Our 2D numerical model is suitable to replicate this experiment. The eddy development process in the shear layer takes place in a quasi-horizontal plane [46] in concordant-bed confluences, and vortices are then quasi-two-dimensional with near-vertical axes [48,49]. The two-dimensional nature of the shear layer may disappear in discordant-bed confluences [30,31].
We depict the steady-state bed elevation contour plot computed with our mathematical model in Figure 3c. The profile presents two morphological features [50]: a scour-hole at the downstream corner of the confluence, and a sediment bar on the left bank immediately downstream of the confluence. The mathematical model is able to reproduce both features, as shown in Figure 3c.
The scour hole depth is similar in both models, experimental and numerical, with a maximum depth of about 2 cm. The experiment is highly sensitive to the depth, since it increases with the confluence angle [45,51] in a non-linear fashion [52]. However, the extension of the hole slightly differs in both models. The experiment shows the hole extends up to the central axis of the main channel, whereas the numerical model provides a hole that extends up to the left bank. Under certain conditions, experiments also show scour holes that expand up to the bank [50].
The sediment bar on the left bank downstream confluence, also denoted as a bank-attached lateral bar, is also reproduced in both the experiment and the numerical model. The sediment bar appears due to the lower water depth and recirculating flow in the separation zone, as well as from the flow detached from the left sidewall. The experiment results in a higher and more concentrated bar than the dimensions computed with the mathematical model. Since our mathematical model is depth-averaged, it is unable to reproduce the helicoidal flow that emerges downstream of the junction in the main channel [50]. The helicoidal flow is responsible for the movement of water in the vertical direction and may concentrate the sediments resulting in a higher and narrower bar. Moreover, the van Rijn model for bedload transport [60] characterizes the sediments with the median diameter (D 50 ) and disregards the variance of the diameter distribution. Bank-attached lateral bars are made of the finest sediments [45], which our model is unable to consider because the variance of the sediment diameter is ignored. The evolution of the bed elevation, water head, water velocity, and stream lines provided by the numerical model are depicted in Figure 4. Initially, the bed elevation is 5.5 cm, as shown in Figure 4a, which provides an almost constant water head everywhere except in the separation zone, i.e., the downstream corner of the confluence (Figure 4d). Under this situation, the velocity contour plot shows a clear shear layer in the separation zone, as shown in Figure 4g, and the highest velocities are given downstream of the confluence in the main channel. The streamline plot shows how water recirculates inside the circulation zone, as shown in Figure 4i. Moreover, this zone extends up to 1.6 m downstream of the junction.
After 3600 s, the bed topography has changed significantly. The scour hole depth is up to 3 cm and the extension is almost the definitive, and slight sedimentation appears downstream of the hole (Figure 4b). The water depth increases in the area of the scour hole due to the erosion, but increases downstream of the hole because of sedimentation (Figure 4e). The water velocity is also affected by the bed elevation evolution. The velocity decreases in the scour hole area thanks to the erosion, but also water slows down downstream of the hole due to the contraction of the recirculation zone that increases the effective channel width (Figure 4h). The contraction of the recirculation zone is well represented in the streamline plot (Figure 4j).
The steady-state is given at time 10,800 s. The bed elevation has evolved up to create a scour hole of about 4 cm deep that extends up to the left bank and a bar in the left bank downstream of the hole (Figure 4c). The water depth increases in the scour hole due to erosion, as shown in Figure 4f, and water slows down due to the larger depth (Figure 4i). The recirculation zone contracts slightly compared with the previously analyzed time step (Figure 4l).

Sensitivity Analysis
We take advantage of our mathematical model to study the erosion patterns as a function of the mesh computational size and the bed properties-friction angle of sediments and roughness. In the following subsection we present and discuss the results of a mesh refinement study and the effect of the friction angle of the sediments and roughness on the erosion patterns. The properties of the sediments and the bed are: D 50 = 0.9 mm, n = 0.018 s·m −1/3 , φ = 25 • , ρ = 2500 kg/m 3 , unless otherwise indicated.

Mesh Refinement Study
We test the grid convergence of our model results with a mesh refinement study. We simulate the erosion experiment depicted in Figure 3 using four mesh sizes, ∆δ = (0.5, 1, 2, 4) cm, in the erodible areas of the channel. The erosion patterns and two bed profiles at the x-coordinates 0 and 20 cm (see reference system in Figure 3a) are included in Figure 5.
The overall erosion patterns are grid independent, but quantitatively we observe smoother topographies (less erosion) in coarser grids. Note that the finest grid already has a grid size that is comparable with the characteristic size of the sediment particles (around 1 mm), so refining below the sub-millimeter scale may raise questions about the validity of the continuum hypothesis. For fixed hydro-mechanical properties the steady state bed elevation converges for four increasing mesh sizes, from ∆δ = 4 cm to ∆δ = 0.5 cm (Figure 5a-d). The deepest scour hole and highest downstream bar are computed with the finest mesh: as the mesh is refined, the scour hole depth increases.
The difference between the computed erosion depths decreases as the grid is refined in both profiles, as shown in Figure 5e,f, as it remains relatively small for the finest grids, ∆δ = 1 cm and ∆δ = 0.5 cm. Nevertheless, the scour hole depth is quite short for the coarsest mesh, 1 cm deep, in contrast with the estimated depth with the finest mesh, 4 cm deep. Based on these convergence results, we adopt a mesh size of 0.5 cm.

Bed Properties: Friction Angle of Sediments
The model assumes that the sediments can be described as a cohesion-less granular material. We study the effect of the internal friction angle, φ, of sediments on the erosion pattern through four mathematical simulations of the experiment depicted in Figure 3. We adopt four angles, from φ = 20 • up to φ = 35 • and run the mathematical simulations. The erosion patterns and two bed profiles at the x-coordinates 0 and 20 cm (see reference system in Figure 3a) are included in Figure 6.
The friction angle plays an important role in the erosion pattern: the higher friction angle, the steeper scour hole and bar. The sediment with the lowest angle, φ = 20 • , provides a flatter pattern, as shown in Figure 6a, than the sediment with the highest angle, φ = 35 • and Figure 6d, for both the scour hole and the bar. The deepest point of the scour hole is given near the downstream corner of the junction close to the left bank, and ranges from almost 5 cm for the sediment with the highest friction angle to 3.5 cm for the sand with the lowest φ (Figure 6e,f). Sediments with high friction angles allow for steeper steady state holes and, consequently, deeper ones. Nevertheless, the erosion in the right bank of the main channel, y-coordinate equal to 0 cm, is deeper as the friction angle decreases. Indeed, the low friction angle of the sand hampers a deep scour hole in the downstream corner of the junction and sediment is conveyed by sliding from the right bank to the left (Figure 6e). A similar behavior is given in the profile at x-coordinate equal to 20 cm, but in the left bank where y = 40 cm (Figure 6f). Since the sediment with the lowest angle does not allow for a steep hole, the erosion conveys up to the left bank.

Bed Properties: Effective Roughness
We simulate the roughness of the bed channel with the Manning's coefficient, n. We simulate four erosion experiments with four coefficients with values ranging from 0.015 to 0.025 s·m −1/3 . The erosion patterns and two bed profiles at the x-coordinates 0 and 20 cm (see reference system in Figure 3a) are depicted in Figure 7. Roughness has a significant influence on the erosion patterns. The larger the Manning's coefficient, the deeper and larger scour hole, and the lower and smaller downstream bar (Figure 7a-d).
The maximum depth of the scour hole varies from 3.5 cm for the experiment with the lowest friction to 5 cm for the simulation with the highest n. The shear stress due to friction between the flow and the bed is proportional to the roughness and increase as n does. Consequently, the scour hole is larger in both sections of the main channel as the n-coefficient increases (Figure 7e,f).
The bar in the main channel downstream of the junction is affected only slightly by roughness. As the n-coefficient increases, the area of the bar decreases and concentrates in both banks of the channel (Figure 7a compared to Figure 7d).
The effective roughness can be estimated through formulae as function of the particle size. Two examples are the Muller formula [72], n = D 1/6 90 /26, or the Strickler equation [73], n = D 1/6 50 /21.1. Both formulae provides n-values of 0.0142 and 0.0147, respectively, for the sediment size used in the experiments. We got the best agreement between experiments and numerical simulations for n = 0.018 s·m −1/3 , which is higher than the previous values. The discrepancy may arise from the grain-size distribution of the sediments, since the formulae account for the whole distribution with only one variable, D 90 or D 50 .

Results and Discussion. Model Application to Field-Scale Channel Junction
To test the implications of sediment erosion and transport for river hydraulics during flash floods, we simulate the erosion on a real river confluence during a flood event. The simulations allow us to study the interplay between bed erosion and water flow patterns and to discern whether erosion has to be taken into account in eroded beds during flood events. Our case study is a channelized junction in southeast Spain. The confluence is located within a semiarid basin with major soil erosion problems situated between Sierra Nevada and Sierra de Los Filabres mountain ranges in Almería, Spain [74,75].

Case Study
Our case study is the confluence of the Nacimiento and Andarax rivers. Both streams are located in Almería province, Spain, between Sierra Nevada and Sierra de Los Filabres mountain ranges. The Nacimiento river is a left tributary of the Andarax river, and both meet at Alhabia, Almería province, Spain. The junction is channelized and the bed of the river is sand, whereas banks are protected with vertical concrete walls. The topography of the simulated stretch is depicted in Figure 8a. Since the river typically has a very small natural discharge, we know the actual bathymetry, which is typically unknown in most practical scenarios. The width of the Andarax river before the confluence is 25 m, and after it is 40 m, and the width of the Nacimiento river is 30 m.  The Nacimiento river is 40 km in length at the junction. The catchment area is 600 km 2 , the terrain elevation at the junction is about 263 m and the maximum elevation is about 1100 m. The Almarax river is 33 km in length at the junction. The basin area is 288 km 2 , and the maximum elevation is approximately 1160 m. The median diameter of the sediments at the junction is D 50 = 1 mm, the porosity is 0.4, the density is 2650 kg/m 3 , and the friction angle is 30 • . The roughness of the rivers, measured with the Manning's coefficient, is n = 0.025 s·m −1/3 . We include an elevation map of our simulation domain in Figure 8a.
We estimate the annual maximum flood quantile for the desired return period through the map of maximum discharges in Spain [76,77]. The map divides Spain into a set of homogeneous regions and provides the most likely regression model to estimate the annual maximum flood quantiles for the desired return period. The map was built using flood series recorded at gauged basins [76]. We compute the flood quantiles downstream of the junction for the return periods, T, of 2, 10, and 25 years. These quantiles are 8.97, 65.84, and 173.06 m 3 /s, respectively. We consider that the flood event lasts 5 h with a parabolic shape, and assign the total flow rate to each basin proportional to its area. The simulated flood events in both basins for the selected return periods are plotted in Figure 8b. Since the overall water flow regime in our domain is subcritical, as boundary conditions we impose the total discharge at the upstream boundary, and critical flow (free discharge) at the downstream end. The outlet of the model is far enough from the junction to guarantee that water reaches uniform flow downstream the confluence. The imposed flow rates are plotted in Figure 8b.
We take advantage of our model and study the effect of the discharge ratio on the erosion patterns. In this sense, we study two situations for the total flood quantile of 10-year return period, i.e., 65.84 m 3 /s: (1) the flow distributes 80% in the Andarax river and 20% in the Nacimiento river, which gives a discharge ratio of 0.2; and (2) the flow distributes 20% in the Andarax river and 80% in the Nacimiento river, which gives a discharge ratio of 0.8. The evolution of the flood events, assuming 5-h duration and parabolic shape, are depicted in Figure 8c.

Effect of Bed Erosion on Hydrodynamics During Flood
We study the interplay between the bed erosion on the flow pattern and vice-versa. We simulate the 25-year return period flood event, as shown in Figure 8b, in our computational domain (Figure 8a). The erosion at the end of the flood event is represented in Figure 9a, and we include the water depth and velocity contour plots at time step 2.5 h in Figure 9b,c, respectively. We also depict the water depth and velocity contour plots at time step 2.5 h in Figure 9d,e, respectively, for a simulation with non-erodible bed.
The initial bed topography is irregular, as represented in Figure 8a. The cross sections of all the channels are proximately U-shape in the central part and flat in the banks. The irregular initial bed shape results in an irregular erosion pattern, as shown in Figure 9a. The junction between both rivers is curved to avoid the scour hole that appears in rectangular confluences.
The curved shape produces erosion in the left or inlet bank and sedimentation in the right or outer bank. Erosion experiments and numerical models of sediment transport in curved channels have demonstrated that erosion occurs in the opposite way: erosion in the outer bank and sedimentation in the inlet bank [78,79]. However, the influence of the lateral flow is not considered in these studies, as well as the bed topography was initially flat. The flow from the Andarax river pushes the Nacimiento flow and detaches it from the outer wall. This effect, together with the bed topography produces erosion on the inner bank. Moreover, the widening of the Andarax river at the confluence promotes the sedimentation at the junction, as two bars appear: the first one in the upstream corner of the junction and the second one in the central area of the Andarax river just downstream of the upstream corner.
The flood event slightly erodes the right bank downstream of the confluence, whereas several sediment ridges emerged in the mid-channel area as well as a bar in the left bank (Figure 9a). Sedimentation patterns downstream of the junction have been observed in laboratory experiments [45], and they are caused by the decrease in velocity since the shear layer vanishes [46].   Erosion has a great influence on the water depth and flow velocity patterns. The water depth pattern computed with a bed-mobile model at a given time step is smooth, as shown in Figure 9b, whereas the output provided by a model without bed erosion is very irregular, with areas with high depths close to other ones with low depths, as shown in Figure 9c. The abrupt pattern for the non-erodible model emerges from the irregular bed topography. Since the bed topography in the bed mobile model can evolve according to the flow conditions, the water depth pattern is smoother.
Water velocity patterns are also greatly influenced by erosion. The water velocity pattern computed with the bed-mobile model at a given time step is much smoother than the one provided by the fixed bed approach, as shown in Figure 9d,e. As in the water depth pattern, the initial abrupt bed topography causes the water velocity pattern to also be abrupt.

Effect of Flooding on Erosion Patterns
We take advantage of our model to study the effect of the flood event on the erosion pattern. We analyze the effect of the flood return period, and the flood discharge ratio-the distribution of the total flow between the two rivers-on the erosion patterns.
We simulate three flood events with return period of 5, 10, and 25 years, respectively. The flood return period influences the erosion or sedimentation intensity, but not the erosion pattern, which remains constant in all of them ( Figure 10). As the return period increases, so does the intensity of the erosion and sedimentation. We also study the effect of the flood discharge ratio, or total flow distribution between both rivers, on the erosion patterns. We simulate three new scenarios for the 10-year flood event: (1) the total flow distributes 80% in the Andarax river and 20% in the Nacimiento river, which results in the pattern included in Figure 11a, (2) the flow distributes 40% in the Andarax river and 60% in the Nacimiento river, which has the outcome of the pattern depicted in Figure 11b, and (3) the flow distributes 20% in the Andarax river and 80% in the Nacimiento river with the pattern depicted in Figure 11c. (c) Discharge flow ratio: 20% Andarax river and 80% Nacimiento river Erosion-Sedimentation (m) 1.00 0.33 -0.33 -1.00 Figure 11. Effect of the discharge ratio on the erosion patterns. Here we plot the erosion patterns under the 10-year flood event. In panel (a) the total flow distributes 80% in the Andarax river and 20% in the Nacimiento river, in panel (b) the distribution is 40% in the Andarax river and 60% in the Nacimiento river, and in panel (c) it is 20% in the Andarax river and 80% in the Nacimiento river.
The comparison between the erosion patterns for the three flood discharge ratios allow us to reach the following conclusions. As the flow rate in the Andarax river increases, the erosion in both banks of the river upstream the junction also does. In the same way, the erosion in the banks of the Nacimiento river upstream the junction rises up as the flow rate in the rives does.
The sedimentation patterns downstream of the junction are nearly the same for the three flood discharge ratios, since the flow rate in this stretch is identical in all the simulations. Small differences emerge in the tail of the bar, which may arise from the total volume of sediment eroded that is different in every situation.
The erosion in the right bank downstream of the junction increases as the flow in the Nacimiento river does. The Andarax river flow displaces the Nacimiento river flow at the junction, pushing the main stream away from the right bank downstream of the confluence. The larger the Nacimiento river flow rate, the lower the Andarax river flow and the more water hits the right bank that promotes erosion.
Two bars emerge from the upstream corner of the junction. The first one is parallel to the Andarax river and extends in some cases upstream along the Nacimiento river. This bar shrinks as the Nacimiento river flow increases. The second bar is perpendicular to the Andarax river and in some cases changes its direction and becomes parallel to the flow. Opposite to the other bar, this one contracts as the Nacimiento river flow increases.

Conclusions
Flooding has become the most frequent environmental hazard, causing significant losses of life and severe economic damage. Mathematical models are useful tools for flood risk forecasting and control, and the current computational resources allow us to simulate flood events with two-dimensional (2D) approaches. An open question is whether erosion should be accounted for or not when it comes to high-fidelity simulation of flood events.
In this paper we have addressed the question by assessing the effect of mobile beds on the flow patterns during flood events through numerical simulations. We have simulated the erosion in a river confluence located in a semiarid basin of southeast Spain. We have also validated our simulations by comparing our numerical results with laboratory experiments of channel confluences reported in the literature. Our 2D simulation is able to reproduce the main features of the erosion patterns in a rectangular channel confluence: a scour-hole on the downstream corner of the confluence, and a bar on the left bank immediately downstream of the junction. However, laboratory experiments have shown that a helicoidal flow emerges downstream of the junction in the main channel and concentrates the sediments. Since our mathematical model is 2D depth-averaged, it is unable to reproduce the helicoidal flow and discordance between experiments and simulations arises.
We have conducted a mesh refinement study that has shown the well-posedness of the problem. Moreover, our sensitivity analysis has concluded that the bed roughness has a great influence on the erosion patterns: the higher the effective bed roughness, the deeper and larger the scour hole, and the lower and smaller the downstream bar. The internal friction angle of the granular river bed sediments plays an important role in the erosion pattern: the larger the friction angle, the steeper the scour hole and bar, which in turn controls the overall hydrodynamics.
We have simulated the erosion process in a channelized river confluence during flood events. The case study is located in a semiarid basin in southern Spain. Erosion has a great influence on the water depth and flow velocity patterns. The water depth pattern computed with a bed-mobile model at a given time step is smooth, whereas the depth provided by a model without bed erosion is very abrupt. Water velocity patterns are also greatly affected by erosion: the pattern provided by a bed-erodible model is smoother than the one computed with the fixed bed approach.
We have analyzed the effect of the flood return period on the erosion pattern. The return period influences the erosion or sedimentation intensity, but not the erosion pattern, which remains constant for all the simulated periods. We have studied the distribution of a given flow rate between both rivers, i.e., the discharge ratio. The increase in flow rate in one river also increases the erosion upstream from the junction in that river. Moreover, the formation of bars in the junction is quite influenced by the discharge ratio. Nevertheless, the erosion pattern downstream of the junction is quite independent of the discharge ratio.
Our study elucidates the interplay between shallow water hydrodynamics and sediment transport during flood events. Erosion plays a key role in these processes, leading to smoother water profiles and overall smaller flow velocities. Current computational resources can deal with 2D numerical simulations of river flow that incorporate sediment transport. Hydrodynamic models in conjunction with sediment transport arise as a useful tool for engineers, practitioners, and stakeholders to evaluate flood control problems.