A coupled 1-D/2-D model for simulating river sediment transport and bed evolution

For calculating flood hazard in complex situations, the 2-D models are more and more used, replacing the 1-D models because they permit a detailed view of the flooding patterns to be obtained. However, a 2-D calculation on large areas remains time-consuming. Thus, one alternative consists in coupling 1-D and 2-D models in order to limit 2-D modelling to the relevant areas. Often, the 1-D model is dedicated to the areas such as the river main channel in which the flow velocity direction can be fixed to a hydraulic axis all over the flood duration; conversely, the 2-D model is limited to the areas such as floodplains in which the flow pattern evolves with time. However, during floods, sediment transport may cause a change of the shape of the bed with consequences on the flowing conditions. Except strong variations of hydraulic parameters, the finer sediments follow the flow and settle where the flow slows down strongly. Then, for river morphology changes, one can focus on the coarser sediments that interact with the main channel bed during flood events while studying finer sediments can allow following the long-term evolution of the floodplain morphology (Juez et al. 2019). Consequently, many authors such as Hou et al. (2018) simultaneously introduce separate equations for suspended load and bed load. However, the problem is simplified here below considering only one sediment class of which the features evolve in space and time.

For detailed modelling of one flood, the processes of deposition, erosion and transport of the sediments should be simulated and the bed parameters should be updated along the simulation if one wishes to estimate the consequences of sediment transport on flow parameters. In order to cope with this objective of detailed modelling, models coupling shallow water equations with more or less complex sediment transport equations are still being developed in one (e.g., Franzini & Soares Frazao 2018) or two dimensions (e.g., Di Cristo et al. 2016). For operational use of such 2-D models, the issue of high computational time is still acute (Garcia-Navarro et al. 2019). Then, field studies are often carried out using a set of models but limiting 2-D hydrosedimentary model to restricted areas (e.g., Pinho et al. 2020).

For coupled models, in order to avoid increasing computational time, it is recommended to use shared memory instead of files for exchanging data. A standardization of the inputs and outputs such as the one proposed by OpenMI software (www.openmi.org) is likely to simplify the access to internal variables of other software without being obliged to know the structure of both codes. Here below, the problem is simplified in coupling two codes written in the same language (Fortran) and compiling the three parts (1-D code, 2-D code and coupling code) at the same time as one single code. However, one main difficulty remains in defining the exchange terms from the variables of one model or the other one with specific care to the time and location of the exchanges. For a coupled 1-D/2-D model, Chen et al. (2012) discuss the point for field application.

In order to present this new coupled model, the contents of the paper are as follows. First, the equations and the numerical schemes of the 1-D and 2-D hydro-sedimentary codes are described. Then, the method for coupling is detailed. In order to illustrate the ability of the new coupled model, a comparison with a full 2-D model is carried out on two field cases. The Agly river case shows the lateral coupling of a 1-D main channel with a wide 2-D floodplain protected by levees set as boundary between the two models. The Rhine river case shows the upstream/downstream coupling within a main channel.

The coupled model was built from the software RubarBE that solves 1-D de Saint Venant equations coupled to Exner equation (El Kadi Abderrezzak et al. 2008) and the software Rubar 20 TS that solves 2-D shallow water equations coupled to an advection-diffusion equation of the sediment concentration (Paquier 2010). In both codes, the single sediment class is characterized by two parameters: the median diameter D₅₀ and the standard deviation of the diameters σ, both parameters evolving in time and space; σ can be calculated by or D₈₄/D₅₀ assuming that the grain-size distribution is lognormal (n % of the mass of sediments have a diameter less than D_n). This latter description of the sediments is adapted to the river in which the transported sediments are homogeneous in terms of interaction with the flow and the bed.

In RubarBE and Rubar 20 TS, the system of equations is solved using an explicit second order Godunov-type numerical scheme, which makes possible that part of any model is temporarily dry and which permits integrating hydraulic jump as an ordinary calculation point. Locally if flow becomes complex or cannot be represented by shallow water equations, specific equations describing the features of hydraulic structures (e.g., weir-type equation) can replace or be added to the standard calculation using shallow water equations. Because the equations are quite similar, the terms exchanged during coupling are the fluxes of mass (water and sediment) and of momentum through the interfaces between the 1-D and the 2-D models. As both solvers are explicit and use an intermediate time at mid time step, the calculation of the fluxes can be performed at that intermediate time taking into account the estimated values of the variables of both models. Then, one can concentrate on the main difficulty that stands in defining a calculation method of these fluxes starting from the variables of both models at locations and times close to the location and time of the boundary problem. Before detailing this latter point, the solvers of the two codes are presented.

These latter flow equations are supplemented by equations defining the behaviour of the sediments.

The lag distance L_s characterizes the distance for sediment transport to reach its saturation rate for a flow condition (Bell & Sutherland 1983) but because L_s is related to the dimensions of sediment movement, bed forms and channel geometry, its definition is difficult. Its value can be set fixed or calculated using empirical equations such as the one proposed by Wu & Wang (2007). However, to be efficient in a 1-D model, this lag distance should be higher than the space step and the width of the active channel (in which sediment transport is concentrated).

The sediment model domain is defined at each cross section by three types of layers: the upper layer refers to the flow region in which sediment is routed downstream with the water flow; the second layer is the active layer, which contains the sediment particles available for transport; either, the volume of the active layer is defined at each time step by Q_s^capΔx/V, with V as the flow velocity and Δx as the cell length or the thickness of the sediment layer is defined from D₈₄; beneath the active layer is the riverbed substrate defined by several sediment layers described at each point of the cross section. For each layer and sub-layer, the sediment size distribution is represented by the two parameters D₅₀ (median diameter) and σ (standard deviation). Exchanges between layers are calculated according to sediment balance (Camenen et al. 2018). During this latter calculation, when two sediment masses (from any layer) are mixed, D₅₀ and σ of the mixture are averages calculated using empirical equations (Camenen et al. 2018). When a mass of sediments is shared between two parts, one can either choose the hypothesis that all the sediments have the same size or consider that there is a sorting process. In this latter case, the first part, which is finer, is transported by the flow, whereas the second one, which is coarser, settles down according to empirical equations that introduce adaptation lengths for median diameter and standard deviation that are generally one or two orders of magnitude higher than L_s (Camenen et al. 2018).

Second member G includes three kinds of terms:

• 1.

  the bottom slope term that is computed as a difference of pressures at constant water elevation;

• 2.

  the friction term that is computed in an implicit way;

• 3.

  the other terms treated in an explicit way.

Because the numerical scheme is explicit, the computational time step should be limited by numerical stability conditions. For the flow routing, the usual Courant–Friedrichs–Levy (CFL) condition is that the Courant number is limited to 1 but, in case of rapid evolution of the geometry, the time step should be reduced.

In a vertically averaged 2-D representation of the flow and the bed, the question of distribution of the erosion and deposition is no longer a problem if the processes can be reduced to what occurs in a vertical column. The behaviour of the sediments can be described in the same way as in a 1-D model, i.e., focusing on the median diameter and considering additional modelling to represent the evolution of the additional parameters such as the standard deviation and the critical shear stresses. When sediment masses are mixed or shared, the same empirical equations (Camenen et al. 2018) as in 1-D are used.

E is calculated at each node and then averaged over the cell (or directly calculated at the cell centre and then distributed between the nodes) in order to obtain sediment input or output. The bed elevation and the features of the sediment layers defined at each node evolve considering E at every time step.

Rubar 20TS uses an explicit finite volume scheme (El Kadi Abderrezzak et al. 2009) on a calculation mesh consisting of quadrilaterals and triangles (with the condition that any edge is the boundary of no more than one cell on each side) in order to adapt to a natural topography of a river. The second order Godunov-type scheme for the flow equations is similar to the one of RubarBE. The flux through one edge is generally computed solving a 1-D Riemann problem perpendicularly to the edge, but this latter calculation can be replaced by the equation of a hydraulic structure similarly to what is proposed as internal boundary conditions for levees (Echeverribar et al. 2019) or for bridges (Dazzi et al. 2020). For Equation (12), only first order is used: the flux is calculated using the hydraulic and sediment variables at time t_n; then, the bottom elevation (Equation (13)) is directly updated at time t_n+1 from time t_n. Generally, hydraulic structures are assumed to transfer the sediment discharge unchanged; however, if the structure represents a breach in an embankment, sediment can be produced inside the hydraulic structure using similar erosion equations as here above (Paquier & El Kadi Abderrezzak 2019).

Then, the numerical scheme includes six steps:

• 1.

  Computing the slope of each one of the four variables h (or z water elevation), hu, hv and C in every cell on x and y axis by the method of the least squares and applying limitations of slopes.

• 2.
  Computing the values of at the intermediate time in the middle of the edge of cell i

  

  (18)

  in which f₁ (respectively f₂) are the fluxes on x (respectively y) axis, S the second member of Equations (9)–(11), the slope of W along x axis, index L (respectively R) means left (respectively right) side of the edge.
• 3.
  At the intermediate time t_n+1/2, solving a 1-D Riemann problem (by an approximate method such as Roe-type linearization) in the direction normal to the edge (Equation (19) similar to Equations (9) + (10) + (11) on the x axis because these last equations do not vary through a rotation) in order to estimate the fluxes through the edges.

  

  (19)

  This step is eventually replaced by the use of an equation typical from a hydraulic structure that will provide the flow discharge from which the momentum fluxes are calculated.
• 4.

  The sediment discharge through an edge is obtained multiplying the flow discharge by the upstream concentration (obtained at step 1). In case of a breach represented by a hydraulic structure, the sediment discharge produced by the breaching process is added. The sediment flux ϕ through one edge is calculated as the sum of this latter sediment discharge and the diffusive term calculated at t_n (from hν multiplied by the gradient of the concentration estimated at the middle of the edge).

• 5.
  At every node k, the bottom elevation Z_k and the sediment layers are updated starting from the bed evolution rate E_k (i.e., Equation (20)) and using the same rules for the sediment parameters as in the 1-D model.

  

  (20)
• 6.
  At the intermediate time t_n+1/2, Equation (21) is applied in order to obtain the final value of the 3-D variable integrating second member on the surface of the cell and adding the corresponding contribution for the second member ⁠:

  

  (21)

  in which f₁ represents flux along x axis, summing on j cells which have a common edge with cell M_i (of area ⁠) and ɛ_ij equals 1 or −1 according to the orientation of edge m_ij (length ⁠) common to the cells M_i and M_j. The second member S includes: the gravity or slope terms or that are treated as fluxes in such a way as a horizontal water surface remains strictly horizontal, the bed friction terms assessed at the centre of the cell but calculated using an implicit way and the diffusion terms also treated as fluxes through the edges after calculation of the slopes on the cell by the method of the least squares. Similarly, the sediment concentration is calculated using Equation (22):

  

  (22)

  and the corresponding sediment parameters are calculated using the same rules for mixing and sharing as in the 1-D model.

For flow equations, the stability condition is a Courant number below 1. As for the 1-D model, a rapid evolution of the bed topography requires a reduced time step.

Owing to the similarity of the two numerical schemes, the following principles for coupling are proposed:

• 1.

  There is no overlap between the 1-D and the 2-D domains, which means that the boundary between the two domains is constituted of lines and/or areas where the exchange terms are calculated.

• 2.

  The exchange terms are the fluxes (water and sediment) through the boundaries between the areas dedicated to each model. These latter fluxes are the same ones as in the 2-D model (flow discharge, (water) momentum flux in the normal and transverse (relatively to the boundary) directions and sediment discharge); they can be characterized starting from the six variables of the 2-D model: h, u, v, C, d₅₀ and σ. The calculation of these fluxes should be performed taking into account the variables of both models. If the time steps are the same ones in both models, the variables of both models are known at the same time and an explicit scheme can be built calculating the fluxes at the intermediate time t_n+1/2. Otherwise, a time interpolation is necessary. Then, the remaining difficulty stands in defining a calculation method of these fluxes, which should take into account the physical processes at the boundary. In order to keep the more detailed approach and to take into account eventual discontinuities at the boundary of the models, the use of the Riemann problem of the 2-D model allows calculation of the flow fluxes at the 2-D edge and providing the results to the 1-D model (Paquier & Bazin 2013).

The approach is similar to the one proposed by Morales-Hernandez et al. (2015) for solute transport; however, the calculation of the concentration flux is slightly different. While Morales-Hernandez et al. (2015) proposed to upwind the concentration and then obtain the solute flux multiplying the flow flux by the value of the concentration at the centre of the 1-D or 2-D cell, the user of the present coupled model can take into account a gradient of concentration in the 2-D cell. The calculation method of the gradient of concentration (and of the flow variables) depend on the type of boundary of the 2-D model; if the boundary is of ‘wall’ type, the gradient is set to zero and the same scheme as the one of Morales-Hernandez et al. (2015) is applied.

As in the 2-D code, the Riemann problem can be replaced by hydraulic structures’ equations if it is more suitable to represent the physical processes.

• 3.

  In order to avoid complicated interpolations, a constraint is added on the building of the 1-D and 2-D meshes: the boundary on the 1-D side is constituted by one or several full 2-D edges (same condition as in Morales-Hernandez et al. (2015)). A difficulty remains in obtaining the values of the six 2-D variables on the 1-D side of one 2-D edge because only five 1-D variables are known (at a location that is not exactly the one of the edge): mainly, the two horizontal components of the flow velocity are known in the 2-D model while only the velocity intensity is known in the 1-D model. Consequently, an assumption is required to define both the intensity and the direction of the velocity on the 1-D side of the boundary edge and this influences the model results (Paquier & Sigrist 1997). More generally, this latter difficulty adds some constraints on the location of the boundary between the 1-D and 2-D domains if one wishes to keep hydraulic consistency.

In order to detail this latter point, nine possibilities for coupling are identified (Figure 1).

Cases 1 and 2: 1-D model upstream, 2-D model downstream, either solving Riemann problem (case 1) or using hydraulic structure (case 2). The boundary is a cross section for the 1-D model. The 1-D flow discharge is supposed to be perpendicular to the edges and is distributed in a similar way as in the 2-D model at the previous time step (or with a uniform velocity if no flow in the 2-D model). The water elevation downstream the 1-D model is the average of the water elevations on the edges of the 2-D model (or uniformly the 1-D water elevation if the flow is supercritical). The sediment flux entering the 2-D model is defined using the flow discharge at every edge as described here above and the concentration, diameter and standard deviation on the 1-D side (thus guessing that these latter variables are uniformly distributed along the boundary).

Cases 3 and 4: 2-D model upstream, 1-D model downstream, either solving Riemann problem (case 3) or using hydraulic structure (case 4). The boundary is a cross section for the 1-D model. The 1-D flow discharge is the sum of the 2-D discharges through the 2-D edges; the water elevation at the 2-D boundary is equal to the one calculated at the upstream cross section of the 1-D model. The 1-D flow discharge is the sum of the 2-D discharges through the 2-D edges. On each 2-D edge, a sediment concentration is associated with the 2-D flow discharge; these resulting sediment discharges are added to obtain one sediment discharge with an average diameter and average standard deviation, which becomes the input of the 1-D model.

Cases 5 and 6: 1-D model at left side, 2-D model on the right bank of 1-D model, either solving Riemann problem (case 5) or using hydraulic structure (case 6). On the 1-D side, the water elevation at the boundary is equal to the one calculated at the centre of the 1-D cell (between two cross sections). One assumes that the 1-D velocity at the boundary is parallel to the boundary on the 1-D side and has the same intensity as the value in the 1-D cell. This latter assumption should be questioned: generally, it leads to a too high velocity intensity; conversely, in case of strong lateral flow between the two models, the transfer is underestimated because the velocity in the direction perpendicular to the boundary cannot be neglected. The sum of the flow discharges through the 2-D edges constitutes a lateral input for the 1-D model. Depending on the flow direction at the boundary, the sediment is entering or leaving the 1-D model so that either the 1-D cell sediment features are used or an average of the 2-D sediment features at the 2-D side of the boundary edges. Note that there are two problems: is actually the 1-D concentration at the interface equal to the average concentration in the 1-D cell? If the flow is going out only on some edges and entering the other ones, are the average concentrations representative of the sediment flux in one way or the other way?

Cases 7 and 8: 1-D model at right side, 2-D model on the left bank of 1-D model, either solving Riemann problem or using hydraulic structure. The calculation method is similar to the previous case.

Case 9: Hydraulic structure defining another kind of coupling (e.g., 2-D model above 1-D model if the 1-D model represents an underground drainage network; the case of superposing 1-D and 2-D with common physical boundary as in Fernandez-Nieto et al. (2010) is excluded). For the 1-D model, the water elevation and the velocity are the ones at the 1-D cell centre and the direction of the velocity is defined as perpendicular to the line middle of the two cross sections defining the 1-D cell. If the flow enters the 1-D model, an average of the sediment features at the 2-D edges is used; if the flow enters the 2-D model, the sediment features of the 1-D model cell are used. The same two questions as for the previous case should find answers: what occurs if the flow on the 2-D edges is not in the same direction? Is the 1-D concentration at the interface equal to the average 1-D cell concentration? For this latter question, if the 1-D model is below the 2-D model and sediment is coarse (bed load transport), the answer is certainly that the 1-D concentration at the interface is much lower except if the structure is starting close to the bottom of the 1-D cross section. Generally, this case occurs for urban flood modelling and is not discussed here below; calculating the exchanged flow is still an active research question (Bazin et al. 2014; Rubinato et al. 2018).

For all the cases, the values of the variables on the 2-D side at the intermediate time are defined at the edge as usually in the 2-D model. They only depend on the way the boundary is set in the 2-D model because the type of boundary determines the calculation method of the slope in the adjacent cell. Finally, in the coupled model tested here below, the coupling formulation can be summed up as follows:

• 1.

  The flow fluxes exchanged between 1-D and 2-D models are calculated on each 2-D edge (included in the boundary) using the same methods as in the 2-D model, either 1-D Riemann problem or equation of a hydraulic structure. The concentration, median diameter and standard deviation are simply upwinded.

• 2.

  On the 2-D side of the boundary, the values of the six variables are calculated using the same method as elsewhere in the 2-D model.

• 3.

  On the 1-D side of the boundary, the water elevation, the concentration, the median diameter and the standard deviation are directly equalled to the values at the 1-D cell. In case of lateral (cases 5–8) or vertical (case 9) coupling, the two components of the velocity are calculated using the hypothesis of the same velocity intensity as at the 1-D cell centre and same direction as the boundary (lateral coupling) or as 1-D axis (vertical coupling). In the case of upstream/downstream coupling (cases 1–4), the 1-D flow velocity is assumed to be in the direction of the 1-D axis at the boundary; the distribution of the normal (respectively tangential) flow discharge between the edges is equal to the one of the normal (respectively tangential) flow discharge on the 2-D side except if no water on the 2-D side for all the edges (then, uniform velocity on the 1-D side).

The calculation algorithm is based on the assumption that one or several 2-D time steps correspond to the 1-D time step, although here below, the simplified assumption of equality of the time steps is used. In this latter case, the time step is equal in the two models so that there is only one intermediate time for which the sediment variables at the interface are calculated. Then, during one time step, calculation passes easily from 1-D to 2-D and from 2-D to 1-D as follows:

• 1.

  The 1-D model is calculated up to the intermediate time t_n+1/2 (steps 1–4 of the numerical scheme).

• 2.

  The 2-D model is fully calculated: at the edges that form the boundary with the 1-D model, the missing values (1-D side) are calculated as explained here above using the 1-D predictions at intermediate time (water elevation, velocity, concentration (ratio of sediment discharge over flow discharge), median diameter and standard deviation). Thus, the steps 3 and 4 of the 2-D numerical scheme can be performed as if the boundary edges were inside the 2-D model. They provide the fluxes for the 2-D model but also the ones of the 1-D model. These latter ones are obtained adding the fluxes of the various 2-D edges that share the neighbourhood of the same 1-D cell.

• 3.

  The time step is finished for the 1-D model (steps 5 and 6 of the 1-D numerical scheme).

As the methods used in the coupled model are the same ones as in the 1-D and 2-D models, the stability of the coupling scheme is a priori similar and a Courant number below 1 should be fixed in both the 1-D and the 2-D models that are linked by coupling. Then, the coupled model should save computational time reducing the number of cells and increasing the time step if small cells can be avoided. However, here below, for the comparison between the full 2-D model and the coupled model, a common time step is fixed and is set short enough in order that all the calculations are stable. In such a case of fixed time step, the computational time is partly proportional to the number of wetted cells and partly proportional to the number of wetted nodes.

The overtopping of Agly levees during the November 1999 flood in the southern part of France (Paquier et al. 2002) led to the flooding of a wide coastal plain, quite flat, in which flow was constrained by many embankments. The peak discharge reached about 2,000 m³/s (return period of about 100 years), which was higher than the conveyance of the main channel between the levees in the lower reach (about 1,500 m³/s). Although the river overflowed to the floodplain upstream from this reach, the levees were overtopped in several locations and particularly upstream from the bridge of the coastal road, which was blocked by wood and other vegetation deposits. Overtopping was one of the causes of flooding of large areas but only with a few centimetres of water. At one location, the overtopping induced dike breaching because of the levee material heterogeneity. Although it occurred at the start of the falling phase of the discharge hydrograph, the breach developed (about 40 m wide and 3 m high on average) and released a 2 m high wave that damaged a sewage plant in Saint Laurent de la Salanque and flooded the village.

A 2-D model of the whole floodplain (18,568 cells) was built to simulate the 1999 flood event. The cell size varies from about 5 m on the levees (eight points to describe one cross section of the main channel) to 200 m in the floodplain. The Strickler coefficients set from the land use vary from 7 m^1/3/s in the built-up areas to 40 m^1/3/s in the main channel. The water inputs are limited to the discharge hydrograph at the upstream end of the Agly main channel (accompanied by sediment inputs at the sediment transport capacity) and to the rainfall over the whole area. The downstream boundary condition is the sea level for the sea at east and the ponds at north. The sediment median diameter is set to 2 mm (σ= 10) for the main channel in which the finer sediments are regularly washed out while they constitute the alluvial deposits of the floodplain (d₅₀ = 0.2 mm and σ= 10). The sorting option is activated with adaptation lengths of 1,000 m for diameter and standard deviation instead of 200 m for the lag distance for the sediment transport capacity calculated using the equation from Meyer-Peter & Müller (1948) with the standard value of 0.047 for the dimensionless critical shear stress. At the breach location, the breach model used (representation of the type ‘hydraulic structure’) assumes that the levee crest goes down by about 0.7 m (at 6 m asl) at 9 a.m., actual time of the beginning of breaching.

Then, a coupled model was built keeping the 2-D model of the floodplain and describing the main channel using a 1-D model with eight-point cross sections, the points being exactly the same ones as the vertices of the initial 2-D model. All the sediment parameters used in the initial 2-D model are transferred unchanged in the coupled model. Note that coupling occurs along the two banks of the main channel at the crest of the levees. Generally, coupling cases are 5 or 7 but, at the breach location, coupling case 8 is used with the specific hydraulic structure corresponding to a simplified breach model; this hydraulic structure involves three 1-D cells and three 2-D cells in the calculation of the flow and sediment discharge through the breach.

The peak water elevation along the main channel during the 1999 flood (Figure 2) is quite well calculated by the coupled model while the 2-D model provides too low results, particularly along the reach between the distances 40,000 and 45,000. This difference provides a too small flooded area along the same reach for the 2-D model (Figures 3 and 4), particularly downstream the breach located at a distance of about 42,300 m (left bank at x= 54,500). One reason is the general erosion of the main channel bed with the 2-D model while the coupled model (1-D in the main channel) provides a varied bottom evolution (erosion followed by deposition, etc.). This varied evolution agrees more with the measured evolution although the latter one (which refers to a period of more than 20 years) is much higher and includes bank sliding in the outside of the curves (Figure 5). In the floodplain, the bottom evolution is quite similar in the two models with some erosion and deposition in the flood pathways. However, the processes are stronger upstream for the 2-D model and stronger downstream for the 1-D model, which agrees with the overflow to the floodplain (higher upstream in 2-D and higher downstream in the coupled model). At the location of the deposits or erosions in the floodplain, the diameter of the floodplain surface is similar in the two models and intermediate between the coarse sediment of the main channel (2 mm) and the fine sediment of the floodplain (0.2 mm) because of successive phases of deposition and erosion. In both models, the sorting process along the main channel is quite limited (about 1% of decrease of the median diameter).

Finally, the process of transfers from the main channel to the floodplain is very similar in the two models, which validates the use of the coupled model in such a lateral coupling including the case of hydraulic structure at the boundary. For the breach, the coupled model that provides a peak flow discharge of about 140 m³/s (instead of 80 m³/s for the 2-D model) seems to provide better results because the water elevation in the main channel is better predicted. However, if one looks at the whole flood, the coupled model provides more wrong wet cells than the 2-D model as shown in Table 1 (percentage shown is calculated using the sum of the number of dry cells inside the observed flooded area and the number of wetted cells outside this latter area divided by the number of cells of the observed flooded area). Nevertheless, this result depends strongly on the calibration of the models and there is a slight bias because the calibration of the flow model was performed on the 2-D model (without sediment transport). If the two models are run with the same time step, Table 1 shows also that the computational time is reduced by 20%, corresponding to the reduction of wetted cells from 2-D model to coupled model.

Issued from a study of the efficiency of introducing sediment to restore the river bed, a 5 km-long reach of the Rhine river was used to compare 1-D and 2-D modelling approaches (Paquier et al. 2013). The originality of this case is the building of an artificial bar close to the right bank of the main channel with a sediment diameter smaller than the one of the surrounding sediments at the surface of the armoured bed.

In October 2010, the artificial bar was built from 182.450 km to 183.070 km (y-coordinate between 850 and 1,350 in Figure 6), about 8 km downstream of the Kembs dam. Horizontal layers were successively displayed with a narrowing width in order to form a bar with an average side slope of about 40°. The total volume of the artificial bar reached slightly less than 22,000 m³ and the bar width about 15 m (i.e., less than 10% of the river main channel width). The sediment used for building this bar was excavated from the neighbouring right floodplain in which polders were created to mitigate floods. For the artificial bar, the final median sediment diameter (D₅₀) was about 21 mm and σ was estimated to 3.7, which illustrates the heterogeneity of the fluvial deposits in the floodplain. This resulting sediment of the bar was very different from that of the armoured bed for which the median diameter is about 100 mm (σ = 1.8) but quite close to the subsurface sediment of the main channel that has a median diameter of about 35 mm. The armoured bed does not move or only slightly moves for flow discharge below 2,000 m³/s, which usually allows simplification of the main transport process (both in the field and in the models) to a simple transport of the bar sediments over a fixed bed.

At the beginning of December 2010, high flows coming from the Rhine upper basin (Switzerland) exceeded the capacity of the derivation canal and a flood entered into the Rest Rhine and eroded the artificial bar. The flood lasted 5 days with two peaks at more than 1,000 m³/s. During the flood, about 80% of the bar volume was eroded, this process being followed by a deposition of the eroded sediments in the main channel near the bar or at a location immediately (less than 100 m) downstream from the bar.

The topography of the riverbed is defined by a set of cross sections of the main channel between the levees surveyed at the beginning of 2010 with an average distance between cross sections of about 200 m. For the 2 km-long reach including the bar, some additional surveying allows definition of the topography of the bar in October 2010 immediately after its building and in February 2011 (after the December 2010 flood that eroded the bar). Interpolation between the levelled cross sections and within these sections leads to a full description of the topography through 273 cross sections of 45 points (corresponding to an average distance of 5 m) separated by a distance of about 20 m. These latter cross sections are directly used in the 1-D model while, for the 2-D model, the corresponding points are used to define a quadrilateral mesh. The coupled model includes the 2-D mesh for the 1.1 km long central part of the reach (2,376 cells out of 11,968 cells) containing the bar and the 1-D mesh for the upstream and downstream parts (220 cross sections).

In order to make the comparison between results easier, the three models of the 2010 flood use the same reference parameters (although, for each model, this set of parameters is not the one that permits the best calibration to be obtained): a Strickler coefficient equal to 28 m^1/3/s, the (Meyer-Peter & Müller 1948) equation for the sediment transport capacity with a non-dimensional shear stress τ_c* equal to 0.047, a lag distance of 1 m (equivalent to no space lag), a corrected critical shear stress using (Ikeda 1982) an equation with a stability angle of 30° and the grain sorting option is not applied. The 1-D model distributes erosion according to Equation (6) while deposition is proportional to a power function of ⁠, both using an exponent m= 1.5 and the Merged Perpendicular Method to calculate τ_j.

The 1-D model provides an erosion of the whole bar while the 2-D model provides only an erosion of about 25% of the bar volume: only the upstream end of the bar is eroded (Figure 6). The coupled model provides nearly the same results as the 2-D model because no sediment comes from upstream and the coupling boundaries are very slightly influencing the flow at the bar neighbourhood. The 2-D (and thus coupled) model provides a too weak erosion of the artificial bar but erosion and deposition are better described so that the quality of the sediment modelling is close to the one of the 1-D model as shown in Table 2. In Table 2, the percentage of nodes with wrong evolution refers to the number of nodes eroded (if actually there was deposition) or with deposit (if actually erosion) divided by the total number of nodes for which there was elevation change.

In all the models, the first flood peak starts the erosion mainly over 300 m and the eroded sediments settle in the immediately downstream sections that are the sections of the downstream half of the bar; then, the second flood peak erodes (partly or totally, depending on the model) the rest of the bar and the previous deposits.

At the upstream section S1 (Figure 7), all the models erode the bar totally or nearly totally, which is more than actually. At the downstream section S2 (Figure 8), the 1-D model erodes totally the bar because the deposits from upstream were deposited in the main channel and then evacuated easily while in the 2-D calculation, most parts of the deposits were stored near the bar and cannot move easily during the second peak of the flood. Further downstream, there is a strong difference between the 2-D model that is disseminating the deposits at specific spots on both banks depending on the flow pattern and the coupled model for which passing to the 1-D model modifies the deposition pattern with a spreading of the deposits over half the width of the main channel.

In order to estimate how large is the disturbance created by the boundary between models at the upstream boundary of the 2-D model, a calculation was performed with the same flow conditions but guessing that the bottom was made of the sediments of median diameter equal to 35 mm, which can be easily eroded during the flood. In such a case, at the upstream end of the 2-D reach of the coupled model, the upstream sediment discharge was shared between the 2-D edges guessing the concentration is equal. The result is a large deposit at the upstream end of the 2-D reach with a much thicker layer where the velocity is lower. In order to make consistent the values on both sides of the boundary, one should reduce the 1-D sediment transport capacity (that it is slightly larger than the 2-D one although it is calculated in the same way) and provide another distribution of sediment between the various edges of the 2-D model (a uniform concentration in the cross section is an assumption that creates strong disturbance).

Although the 1-D model uses a distribution of bed shear stress across the section (Merged Perpendicular Method), it leads to higher sediment transport capacity than the 2-D model because the 1-D model uses the mean velocity in the cross section instead of various velocities across the river width. At the upstream/downstream boundary between a 1-D reach and a 2-D reach, this difference in the total (i.e., for the whole cross section) sediment transport capacity creates some disturbance. However, the main disturbance comes from the distribution of the concentration across the section.

The coupling of the 1-D model RubarBE with the 2-D model Rubar 20TS provides a 1-D/2-D coupled hydrosedimentary model adapted to simulate various field cases. The use of explicit numerical schemes permits treatment of the cases in which the boundary becomes dry or wet without stability problems. In the 1-D and 2-D components of the coupled model, the same representation of the sediments by a mean diameter and a standard deviation is used. Camenen et al. (2018) showed that this representation could be efficient for such complex questions as sediment sorting. Then, the extension of the coupling of 1-D and 2-D flow models to sediment transport is reduced to defining three variables at the interface of the models: concentration, mean diameter and standard deviation. The simple assumption of a uniform distribution of the sediment concentration at the 1-D/2-D boundary of one 1-D cell is tested on two field cases: one with lateral coupling (Agly river) and one with upstream/downstream coupling (Rhine river).

The comparisons between full 2-D modelling and coupled modelling on the two field cases show very similar results, which justifies the use of the coupled model. The main difficulty for the coupled model comes from the distribution of a 1-D flux between various edges of the 2-D model. One can expect to cope with this latter difficulty by adapting the interpolations at the 1-D/2-D interfaces of the coupled model and correcting the calculations of the sediment transport capacity in order to reduce the difference between the 1-D and 2-D models.

More generally, as for the coupling of hydrodynamic models, if getting parameters for the 2-D side of the boundary is easy, it is not so simple for the 1-D side, in which additional assumptions are necessary. Although the presented coupled model is stable and allows computational time to be served, further improvements are necessary to obtain suitable variables at the boundary on which coupling creates exchanges. In order to avoid discrepancy at the boundary, the distribution of the variables on the 1-D side should be adapted to the 2-D representation of the physical processes.

The PROFAS B+ program for the scholarship (Algerian-French cooperation) funded the stay of first author at INRAE.

Data cannot be made publicly available; readers should contact the corresponding author for details.