A semi-implicit approach for sediment transport models with gravitational effects

Applied Mathematics and Computation 421 (2022) 126938 Contents lists available at ScienceDirect Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc A semi-implicit approach for sediment transport models with gravitational effects J. Garres-Díaz a , ∗, E.D. Fernández-Nieto b , G. Narbona-Reina b , 1 a Dpto. Matemáticas. Edificio Albert Einstein, Universidad de Córdoba, Spain b Dpto. Matemática Aplicada I, Universidad de Sevilla, Spain a r t i c l e i n f o Article history: Received 29 April 2021 Revised 15 November 2021 Accepted 8 January 2022 Available online 24 January 2022 Keywords: Sediment bedload transport Shallow water exner model Gravitational effects Semi-implicit scheme a b s t r a c t In this work efficient semi-implicit methods for sediment bedload transport models with gravitational effects under subcritical regimes is proposed. Several families of models with gravitational effects are presented and rewritten under a general formulation that allows us to apply the semi-implicit method. In the numerical tests we focus on the application of a generalization of the Ashida-Michiue model, which includes the gradient of both the bedload and the fluid surface. Analytical steady states solutions (both lake at rest and non vanishing velocity) are deduced an approximated with the proposed scheme. In all the presented tests, the computational efforts are notably reduced thanks to the proposed method without losing the accuracy in the results. ©2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) 1. Introduction The sediment transport phenomena is a problem of interest in the environmental field. Besides the typical application in the morphology of a river, the sphere of action of the sedimentation process is wide and it has a deep impact on its environment. Due to the flow of the river some particles are swept along the current. This material is taken mainly from the riversides and also from the bottom, being deposited downstream when the flow becomes weaker, either again on the riversides or in its mouth. This effect change the morphology of the river environment sometimes affecting to crops or to natural protected areas but also to fluvial navigation and coastal zones with the consequent impact in the economy and environmental aspects. It is also important the effect of the sediment transport process in the building of hydraulic structures like dams or bridges. Either the long-term erosion process or for example in torrential rain events or during the thaw, some of the sediment of the bottom can be dragged and may affect to the stability of these structures. Normally most of these important phenomena are long-term processes so there is a special interest on the development of efficient predictive techniques in order to prevent hazards or to better manage the surroundings. Sediment transport occurs mainly through two phenomena: bedload and suspended load. Bedload entails the transport of sediment particles rolling or sliding on the bed and jumping into the flow and then resting on the bed again. Particles transported by suspension are supported by the surrounding fluid during a significant part of the current and may also be ∗Corresponding author. E-mail addresses: [email protected] (J. Garres-Díaz), [email protected] , [email protected] (E.D. Fernández-Nieto). 1 All authors were partially supported by the Spanish project RTI2018-096064-B-C22 and FEDER. https://doi.org/10.1016/j.amc.2022.126938 0 096-30 03/© 2022 The Author(s). Published by Elsevier Inc. This is an open access article under the CC BY-NC-ND license ( http://creativecommons.org/licenses/by-nc-nd/4.0/ ) J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 deposited. In this paper we focus on bedload sediment transport including gravitational effects that is an essential but not trivial task to take into account in the mathematical modeling. When the bedload phenomena is under study, we must consider both, the hydrodynamical component and the morphodynamical one that are coupled. The Saint Venant Exner system is commonly used to describe the bedload sediment transport and it is a system composed on the well known Saint Venant (or Shallow Water) system for the hydrodynamics and a continuity equation to update the morphodynamical part. This continuity equation is called Exner equation [16] , and it is defined in terms of the solid transport discharge that has to be prescribed. In the literature several models have been defined empirically to give this closure [1,17,25,33,35,40] . Even if these models are largely used they provide some inconvenient properties for the system as the lack of a dissipative energy or the loss of the mass conservation, other important limitation being the validation range for just nearly horizontal sediment beds. An alternative to these empirical models has been presented in [18] where a complete Saint Venant Exner type model is derived from the Navier-Stokes equations, this model encompassing the former disadvantages. In classical models, the solid transport discharge is defined in terms of the shear stress and the critical shear stress, or its normalized form called Shields parameter. Bedload models predict that there is no transport of sediment particles whenever the shear stress is smaller than the critical shear stress, see for instance [1,16,17,25,33,35,40] . The Shields parameter is defined as the ratio between the agitating and stabilizing forces on a deposited sediment particle. Lysne [32] showed experimentally for inclined sediment beds that the gravity is an important contributing action as agitating force, see also [20] . Nevertheless classical formulas consider that the shear stress is the unique agitating force. Therefore, gravitational forces play an essential role into the sediment transport, so they must be considered for the definition of an appropriate solid discharge for applications in general sloping beds. Several forms to include gravitational effects in the solid transport discharge are found in the literature. The simplest one is to include a second order derivative in the Exner equation, which is neglected when the slope of the bed is smaller than the tangent of the repose angle of the material (see Tassi et al. [39] ). Nevertheless, the form in which gravitational effects are usually included is by modifying either the definition of the shear stress or the critical shear stress. Although these two kind of modifications might coincide in some particular cases for 1D horizontal domains, they are not equal in general. Here we consider the inclusion of the gravitational effect in the solid transport discharge, which is defined in terms of the effective shear stress (see e.g. Fowler et al. [19] ). This is discussed in detail by Fernández-Nieto et al. [18] showing that the solid discharge naturally incorporates gravitational effects, which are included in the effective shear stress. Let us also mention that the inclusion of gravitational effects in the effective shear stress is analogous for scalar or vector systems, that is, 1D or 2D horizontal domains. Nevertheless, when the modification of the critical shear stress is considered, the inclusion of gravitational effects for 2D domains is not an easy extension of the 1D case. See Kovacs and Parker [31] , Seminara et al. [38] and Parker et al. [36] among others. From the numerical point of view, bedload transport models have been usually discretized by means of explicit schemes in collocated and staggered meshes (see [3,6,7,18,26,34] among many others). These problems are characterized by two different time scales: a small characteristic time for the hydrodynamical counterpart and a large characteristic time for the morphodynamical contribution. It makes the computational cost to be very high when explicit discretizations are used, even for simple models (e.g. Grass formula) without gravitational effects. In Bilanceri et al. [4] a comparison of explicit and implicit methods for bedload transport (without gravitational forces) is made, as function of the low/medium/high sedimentfluid interaction. They claim that the fully implicit scheme is too much expensive, and a linearization of the method based on the automatic differentiation software Tapenade [27] is made to solve the Grass model, showing that the computational cost is decreased, specially in the low interaction case. Just a few works are devoted to numerical approximation of bedload models with gravitational effects. In these models, the difficulty from a numerical point of view is the presence of a non-linear elliptic counterpart in the Exner equation. As a consequence, its explicit discretization is computationally expensive. Tassi et al. [39] used a Discontinuous Galerkin method to approximate a Grass model with a diffusion term depending on the free surface gradient. Later, Morales de Luna et al. [34] used a duality method based on the Bermúmez-Moreno algorithm to solve the morphodynamical component of a Meyer-Peter & Müller model with gravitational effects, which is computationally more expensive. Therefore, semi-implicit schemes are an interesting and promising alternative in the framework of bedload transport models. The method considered here is the one introduced by Casulli and co-worker for hydrostatic [8,10,11,13] and nonhydrostatic [9,12,14] free surface flows in z-coordinates and isopycnal cooordinates. In Bonaventura et al. [5] and Garres-Díaz and Bonaventura [23] this method was adapted to vertical multilayer discretization, showing its efficiency when it is applied to bedload transport with the simple Grass formula, and also for variable density flows. Such a semi-implicit approach was used in Garegnani et al. [21,22] , where an analysis of the coupled and decoupled approach of the Exner equation is made, introducing also a semi-implicit method [37] for the system with movable bed, although they did not consider gravitational effects. The main contribution of this work is a low-computational cost semi-implicit scheme for sediment transport problems solved through the Saint Venant Exner system introduced in [18] taking into consideration also gravitational effects, in long-time scale, i.e., slow processes. In particular, we propose a first order method ( Θ-method) and a more accurate method (IMEX-ARK2). It is based on a reformulation of the solid transport discharge by rewriting the sign function. Furthermore, the proposed method could be easily adapted to a wide range of families of bedload models with gravitational effects, which is 2 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Fig. 1. Computational domain and notation. also described in this work. Thus, to our knowledge, this is the first efficient numerical scheme for general bedload models with gravitational effects. Section 2 is devoted to present the model we use, as well as the reformulation of the solid transport discharge that we propose to apply the semi-implicit scheme to a wide family of models with gravitational effects. In Section 3 the semiimplicit method based on the Θ-method is detailed, and the numerical experiments are in Section 4 . Finally, we present some conclusions in Section 5 . In Appendix A the IMEX discretization of the system is presented. 2. Saint Venant Exner system with gravitational effects In this section we present the initial system for the hydrodynamical and morphodynamical counterparts. The solid transport discharge including gravitational effects is defined in Section 2.2 . Once the chosen model is exposed, the reformulation of the solid transport discharge is proposed in Section 2.3 , where we also analyze the steady solutions of the proposed model. Finally, we show in Section 2.4 how this reformulation could be easily adapted to a wide range of families of bedload models, depending on both the definition of the solid transport discharge and the effective shear stress. 2.1. Initial system The base system is the well-known Saint Venant Exner system, which is obtained from the coupling of the Saint Venant system (hydrodynamical component) and the Exner equation (morphodynamical component). This model is deduced from the non-dimensional Navier-Stokes equations for the hydrodynamical part together with the Reynolds equation for the evolution of the granular layer, a detailed derivation and analysis of this model can be seen in Fernández-Nieto et al. [18] . Considering a 1D incompressible fluid with constant density ρ∈ R , the system reads  ∂ t h + ∂ x q = 0 , ∂ t q + ∂ x q 2 /h + gh∂ x ( h + b ) = −F, ∂ t b + ∂ x Q b = 0 , (1) where (x, t) are the space and time variables, h (x, t) and u (x, t) the height and averaged horizontal velocity of the fluid, and q = hu is the discharge. The gravity acceleration is denoted by gand Fis a friction term between the fluid and the sediment layer b(x ) which will be defined later. This sediment is transported according to the solid transport discharge Q b that is defined by the chosen model for the bedload transport. In the deduction presented in [18] the thickness of the sediment layer (denoted by bhere and h 2 in [18] ) is subdivided in a lower static layer and a small upper movable layer whose height is h m (see Fig. 1 ). Simplified models are also presented there, which is the framework considered in this work. Assuming a quasi-uniform regime, the erosion rate equals the deposition rate and h m ca be estimated in terms of the rest of unknowns of the problem. Let us denote the shear stress by τ, which can be written, in the framework of depth-averaged models, as τ= ρgh Swith Sa friction term, whose more common definitions are through the Manning ( S = n 2 u | u | h −4 / 3 , n the Manning coefficient) or Darcy-Weisbach ( S = ξu | u | / 8 gh , ξthe Darcy-Weisbach coefficient) laws, depending on the averaged velocity u . Then, we can write τ ρ= C | u | u, with C = gn 2 h −1 / 3 (Manning law), C = ξ/ 8 (Darcy-Weisbach law) . (2) An empirical formula is used to define the solid transport discharge as a closure to the system, for instance Grass [25] , Meyer-Peter and Müller [33] , Ashida and Michiue [1] , among many others. For instance, when the classical Ashida 3 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 & Michiue’s model is chosen, it is written in nondimensional form as Q b Q = sgn (τ) 17 1 −ϕ ( θ−θc ) +  θ− θc , (3) where ϕis the porosity of the sediment bed, sgn (·) is the sign function, (·) + the positive part, and Q = d s  g ( 1 /r −1 ) d s the characteristic discharge that is defined in terms of the gravity, the density ratio r = ρ/ρs , with ρs and d s the density and the diameter of the sediment particles. Finally, θis the so-called Shields parameter, defined by θ= | τ| d 2 s g ( ρs −ρ) d 3 s , (4) with θc the critical Shield stress. A general formulation that includes a great number of classical models for the bedload solid transport discharge without gravitational effects may be written under the following compact form (see [18] ) Q b Q = sgn (τ) α1 1 −ϕ θβ1 ( θ−α2 θc ) β2 +  θ−α3  θc β3 , (5) with αi , βi i = 1,2,3 positive constants. Some important drawbacks of these classical models for the bedload transport are: (i) they do not take into account gravitational effects, since they are derived using as hypothesis that the sediment free surface is almost constant ∂ x b ≈0 ; (ii) they do not satisfy a dissipative energy balance; (iii) the mass of sediment may be not conserved. In this work we study gravitational effects in models satisfying a dissipative energy balance. We consider the case of quasi-uniform regimes, where the movable thickness of sediment is set in terms of the velocity of the fluid. In these cases the mass conservation will not be guaranteed but it can be easily solved as detailed in the following subsection. Models including gravitational effects are presented in following subsections. 2.2. Solid transport discharge with gravitational effects In this subsection, following [18] , we present the generalization of the Ashida-Michiue and Meyer-Peter & Müller models under the assumption of a quasi-uniform regime. The general definition of the solid transport discharge is Q b = h m v b  (1 /r −1) gd s , (6) where h m is the thickness of the movable bed and v b the averaged sediment velocity. This velocity is defined in terms of the effective shear stress ( τeff) as follows v b = sgn (τeff)(  θeff − θc ) + , (7) with θeff = | τeff| /ρ (1 /r −1) gd s , where τeff must be properly defined. Here, following [18] , we consider (see discussion in Section 2.4 ) τeff ρ= τ ρ−ϑgd s r ∂ x ( rh + b ) , with ϑ = θc tan δ, (8) being δthe angle of repose of the material. Frequent values for these constants are θc = 0 . 047 and ϑ = 0 . 1 , that is θc /ϑ ≈ tan 25 ◦(see Fredsœ in [20] ). Other possibilities are used in [15] and references therein. We can establish the relation between θand θeff by computing τfrom the previous equation, obtaining θeff = sgn (u ) θ−ϑ 1 −r ∂ x ( rh + b ) , where sgn (u ) is identified as the sgn (τ) and we used the definition of τin (2) . In order to set the definition of h m we will consider the case of a quasi-uniform regime, where h m is defined by a closed formula as a function of the erosion and deposition rates (see [18] and references therein for further details). A possibility is to define h m = K e d s K d ( 1 −ϕ ) ( θeff −θc ) + , (9) with K e , K d constant parameters related to the erosion and deposition effects respectively. Let us remark that Fernández- Luque and Van Beek [17] observed in the experiments this linear relation between the thickness of the movable bed and the shear stress. Most of classical formulas consider this relation to obtain the solid transport discharge. This linear relation was also observed in Bagnold [2] by investigating the momentum transference because of the sediment particles and fixed bed interaction (see also [15] ). Notice that the fact of using a closure formula for h m as above may be the mass conservation fail, 4 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 see [34] . Nevertheless, this problem is solved by defining h m as the minimum between the definition (9) and the thickness of the erodible sediment layer. In this work we consider that this layer is large enough, since a fixed bedrock has not been taking into account. By using the definition of h m (9) we obtain the following formula of the solid transport discharge Q b Q = sgn ( τeff) K e K d ( 1 −ϕ ) ( θeff −θc ) + (  θeff − θc ) + where Q = d s  (1 /r −1) gd s . (10) Note that this is a generalized Ashida & Michiue’s model (3) , where the ratio between erosion and deposition effects is set to K e /K d = 17 . Another possibility is to define h m as follows (see [38] ) h m = Kd s 1 −ϕ  (θeff −θc ) + (  θeff +  θc ) . (11) In this case we obtain that the solid transport discharge is Q b Q = sgn ( τeff) K ( 1 −ϕ ) ( θeff −θc ) 3 / 2 + . (12) Note that in this case we obtain a generalized Meyer-Peter & Müller model, where the constant parameter is set to K = 8 . Finally, in [18] it is also deduced the following definition of the friction term between the fluid and the sediment layer, which is proportional to the difference of velocities F =  gh m C r ∂ x (r h + b) + (1 −r) sgn (τeff) tan δif θeff > θc , τ/ρotherwise. (13) Remark 1. Notice that this definition of Fcoincides with the following definition of the friction between the fluid and the sediment layer: F = C| u −v | (u −v ) , where v = u − gh m /r |P | 1 / 2 sgn (P ) if θeff > θc , 0 otherwise , being v the velocity of the sediment layer and P = ∂ x (r h + b) + (1 −r) sgn (τeff) tan δ, for more details see [18] . In the next subsection the final model is reformulated in order to properly introduce a semi-implicit numerical scheme. 2.3. Reformulation of the bedload model with gravitational effects Hereinafter we consider the model defined by (1), (10) and (13) , corresponding to the generalization of the Ashida- Michiue Saint Venant Exner model including gravitational effects and a friction term proportional to the difference of velocities between the fluid and the sediment layer. Our goal in this work is to develop a semi-implicit scheme based on splitting the system as a stiff part (gravitational terms) with a high impact over the stability CFL condition and a non-stiff contribution. Looking for a more convenient writing of the system in that sense using (2), (8) and setting sgn (W ) = W/ | W | , the solid transport discharge (10) is reformulated as Q b = τeff ρ q b =  q b C | u | u −k 1 ∂ x ( b + h ) −k 2 ∂ x b (14a) with k 1 = ϑgd s , k 2 = ϑgd s (1 /r −1) (14b) positive constants, and  q b = ρQ | τeff| K e K d ( 1 −ϕ ) ( θeff −θc ) +  θeff − θc . (14c) For convenience to apply the numerical scheme, the final system is written in terms of the free surface level. Then, defining η= b + h and combining the first and third equations in (1) , the final system reads ⎧ ⎨ ⎩ ∂ t η+ ∂ x q +  q b Cu | u | − q b k 1 ∂ x η− q b k 2 ∂ x b = 0 , ∂ t q + ∂ x q 2 /h + g ( h + C h m ) ∂ x η= −g C h m (1 /r −1) ( ∂ x b + sgn (v b ) tan δ) −χτ/ρ, ∂ t b + ∂ x  q b Cu | u | − q b k 1 ∂ x η− q b k 2 ∂ x b = 0 , (15) where we identified sgn (v b ) as the sgn (τeff) and χ= 1 θc ≥θeff 0 otherwise . 5 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Regarding the stationary solutions of system (15) , it is difficult to obtain an explicit expression in general. However, it is interesting to firstly analyze the simple case concerning lake at rest solution that is a steady solution for the Saint Venant system. That is, those ones verifying b + h = η0 constant , q = 0 . (16) It is easy to check that they are also steady solutions of previous system if gravitational effects are not considered ( k 1 = k 2 = 0 ), since we recover τ= 0 and therefore θ= 0 and Q b = 0 . That means that solutions given by (16) are steady solutions independently of the sediment layer profile b(x ) . Obviously, this is a limitation of classical models and non-physical solution will be kept. On the contrary, this is not the case for the proposed model with gravitational effects. Actually, we have the following result. Proposition 1. Let η0 constant be the free surface level of a fluid at rest ( u = 0 ), with b(x ) the sediment layer and δthe angle of repose of the sediment. Then (η0 , 0 , b(x )) is a steady solution if and only if | ∂ x b | ≤tan δ. That is, a lake at rest solution is steady if and only if the sediment layer has a slope lower than the slope given by the angle of repose of the sediment. Proof. Let us start proving that the condition | ∂ x b | ≤tan δis equivalent to θeff ≤θc under the hypothesis ∂ x η= 0 , u = 0 . In this case θeff = k 2 | ∂ x b | g ( 1 /r −1 ) d s = ϑgd s (1 /r −1) g ( 1 /r −1 ) d s | ∂ x b | = θc tan δ| ∂ x b | ≤θc , and this inequality holds if and only if | ∂ x b | ≤tan δ. Using the condition θeff ≤θc , it is easy to check that Q b = h m = 0 (see (9) and (10) ), and therefore ∂ t η= ∂ t q = ∂ t b = 0 , i.e, (η0 , 0 , b(x )) is a steady solution. Assuming now that (η0 , 0 , b(x )) is a steady solution, we obtain ∂ x ( k 2  q b ∂ x b ) = 0 and h m ( ∂ x b + sgn (v b ) tan δ) = 0 . From the first condition, either ∂ x b = 0 and therefore τeff = 0 = θeff ≤θc and also h m = 0 , or  q b = 0 and therefore ( θeff −θc ) + = 0 , which ends the proof.  Let us now study steady solutions in the general case, neglecting the friction term χτ/ρin the momentum equation in (15) . Theorem 1. Let ( η(x ) , q 0 , b(x ) ) be the values of the free surface, discharge and sediment layer satisfying that Q b = 0 , with q 0 constant. Then, it is a steady solution of system (15) if and only if it satisfies −sgn ( τeff) β∂ x b ≤tan δ1 −sgn ( τeff) C | q 0 | q 0 gd s θc ( 1 /r −1 ) 1 ( η−b ) 2 , ∂ x η= α∂ x b, (17) with α= −q 2 0 g ( η−b ) 3 −q 2 0 and β= 1 + α (1 /r −1) . Proof. Let us start by noticing that the steady solutions of system (15) with no sediment transport are given by ∂ x q = 0 , ∂ x u 2 2 + gη= 0 , Q b = 0 . (18) The first condition implies that q (x ) = q 0 constant. From the second condition in (18) , and writing u as q 0 / (η−b) , we deduce that ∂ x η= α∂ x b, with α= −q 2 0 g ( η−b ) 3 −q 2 0 . Furthermore, since Q b = 0 , the only possibility is that θeff ≤θc holds, which leads to the inequality  tan δC | u | u gd s θc (1 /r −1) −β∂ x b ≤tan δ, where β= 1 + α (1 /r −1) . (19) Now, we can solve the inequality for the unknown ∂ x b, taking into account that the sign of the left hand side coincides with sgn ( τeff) , obtaining inequality (17) , which ends the proof.  Notice that the Proposition 1 is a particular case of the Theorem 1 where q 0 = 0 , obtaining α= 0 and β= 1 in (19) and η(x ) = η0 constant. When the limit case is considered in (17) , i.e. when the equality holds, the explicit expressions of ηand bare the solutions of the resulting nonlinear ODE system. Note that η, bin (17) can be also found by solving just the initial value 6 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 problem for the sediment, and updating the free surface value at each step using conservation of the energy (middle condition in (18) ), that is η= b + h where h is the greater root (for subcritical solutions) of the third-order degree polynomial 2 gh 3 + 2(gb −K 0 ) h 2 + q 2 0 = 0 , where K 0 = u 2 0 / 2 + gη0 . In the results presented in this subsection we have focused on the generalized Ashida-Michiue model for the sake of simplicity. However, many other Saint Venant Exner models including gravitational effects could be written in the same form as (14) - (15) , for instance the generalized Meyer-Peter & Müller model defined by (12) for K = 8 . In the next subsection other possible generalizations to include gravitational effects are discussed, which can be easily written under the form (14) - (15) . Consequently, the same numerical technique that is proposed in Section 3 can be applied for this family of models. 2.4. A general formulation of bedload sediment transport models with gravitational effects In this section we briefly present several approaches to include gravitational effects in classical models, and the relations between them. As we commented before, it is suitable for the development of the semi-implicit method presented in Section 3 if the model is written under the form (14) . Thus, all the discussed models below will be expressed in the same way. In classical models, with the purpose of including gravitational effects, we shall replace τwith τeff and θwith θeff in the solid transport discharge Q b . For example, when considering the family of classical models defined by the solid transport discharge (5) , the corresponding models with gravitational effects are defined by Q b Q = sgn (τeff) α1 1 −ϕ θβ1 eff ( θeff −α2 θc ) β2 +  θeff −α3  θc β3 , (20) being α1 , α2 , α3 , β1 , β2 and β3 non-negative constants depending on each model. In addition, note that the effective shear stress must be defined to write previous models under the form (14) . Actually, a different family of bedload models is obtained for each definition of τeff proposed in the literature. In the next lines we present some relevant definitions of τeff in order to show the differences and the similarities, but also to give a justification of our choice given in (8) . Let us begin with the effective shear stress introduced in Fowler et al. [19] that we denote by τF eff and which is given by τF eff ρ= τ ρ−λg d s (1 /r −1) ∂ x b, (21) where τis a quadratic law defined as in (2) and λ= 1 . Later, Morales de Luna et al. [34] proposed to define λ= ϑ, where ϑis defined in (8) , with the aim of recovering lake at rest steady solutions associated to the repose angle of the sediment. Whatever the value of λ, we find a family of models that can be reformulated defining Q b as in (14a) , and  q b , k 1 , k 2 as follows  q b = ρQ τF eff α1 1 −ϕ (θF eff) β1 θF eff −α2 θc β2 +  θF eff −α3  θc β3 , (22) with θF eff = | τF eff| /ρ (1 /r −1) gd s , k 1 = 0 and k 2 = λg d s (1 /r −1) . Note that the definition of k 2 is the same that in previous cases (14b) , but k 1 is neglected. This means that the free surface gradient has no influence on the solid transport discharge when τF eff is considered. Equivalently, τeff in (8) matches with τF eff if the free surface of the fluid is constant, ∂ x (b + h ) = 0 . Let us now go deeper in the definition of τeff to understand the expression (8) that we use in this work. In Fenández- Nieto et al. [18] the Navier-Stokes system is asymptotically analized to deduce a Saint Venant Exner system, leading to propose several models including gravitational effects and satisfying a dissipative energy balance. These gravitational effects are included through the considered effective shear stress. The law considered for the drag force between the fluid and granular layers determines both the model and the effective shear stress. In particular, two expressions for the effective shear stress are deduced corresponding to a linear or quadratic friction law, denoted by τL eff or τQ eff respectively, that we present next. In the linear friction case, it is defined by τL eff ρ= ϑd s h m τL ρ−ϑgd s r ∂ x ( rh + b ) , with τL ρ= Ch m u  g(1 /r −1) d s ϑd s , (23) or equivalently, τL eff ρ= Cu  g(1 /r −1) d s −ϑgd s r ∂ x ( rh + b ) . (24) 7 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Now the family of models is given by Q b = τL eff ρ q b =  q b Cu  g(1 /r −1) d s −k 1 ∂ x ( b + h ) −k 2 ∂ x b , (25) where  q b = ρQ τL eff α1 1 −ϕ (θL eff) β1 θL eff −α2 θc β2 +  θL eff −α3  θc β3 , (26) and with k 1 and k 2 as in (14b) for θL eff = | τL eff| /ρ (1 /r−1) gd s . In the quadratic friction case, the effective shear stress is (see [18] ) τQ eff ρ= gd s (1 /r −1) | | , (27) where = sgn (u )  | τQ | /ρ  (1 /r −1) gd s −  ϑ 1 −r ∂ x ( rh + b ) sgn ( ∂ x ( rh + b ) ) , with τQ defined in terms of τ(2) as τQ = h m ϑd s τ. (28) The effective shear stress τQ eff (27) gives other family of bedload models, where sgn (τQ eff) = sgn () and θQ eff = | | 2 , which can be reformulated under a similar expression to (14a) , concretely we obtain Q b =  q b u  ϑd s h m C −k 1 ∂ x ( b + h ) −k 2 ∂ x b , (29a) with  q b = Q u  ϑd s h m C −k 1 ∂ x ( b + h ) −k 2 ∂ x b  α1 1 −ϕ (θQ eff) β1 θQ eff −α2 θc β2 +  θQ eff −α3  θc β3 (29b) and k 1 =  ϑgd s r  | ∂ x (rh + b) | , k 2 =  ϑgd s r  | ∂ x (rh + b) | ( 1 /r −1 ) . (29c) Note that this model is much more complicate that any of the models defined by (22) or (25) . Firstly, because the values k 1 , k 2 are no more constant. Secondly, note that h m can be defined by (9) or (11) , in terms of θQ eff −θc + , and as a consequence θQ eff is implicitly defined. Let us now make a comparison of the three definitions (21), (23) and (27) given for the effective shear stress. When the effective shear stress of the linear model (24) is compared with the definition (21) given by Fowler [19] , two differences are observed. Firstly, the first term in the definition changes, it is quadratic in the velocity for the Fowler’s model and linear in (24) , which is consistent with the hypothesis of a linear drag force between fluid and sediment layer. Secondly, the gravitational terms are different, although they are equal when the free surface is constant and λ= ϑ. It can be easily seen by writing ϑgd s r ∂ x ( rh + b ) = ϑgd s (1 /r −1) ∂ x b + ϑgd s ∂ x ( b + h ) . Therefore, the definition of the gravitational effects introduced in [18] is more general, since it takes into account the weight of the upper fluid and the coupled effect of the variations of the free surface and sediment. Furthermore, when looking at the effective shear stress deduced for the quadratic model (27) we do not find the same definition than in the Fowler’s model, even if they both correspond to a quadratic friction law. Moreover the gravitational terms are also different, being much more complex in the quadratic model. The effective shear stress (8) used in this work may be interpreted as a linearized version of the quadratic effective shear stress, corresponding to the definition τL eff (23) when the friction term τL is replaced by τQ defined in (28) , thus obtaining our definition τeff ρ= τ ρ−ϑgd s r ∂ x ( rh + b ) . 8 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 As commented above, this model can be also seen as an enhanced Fowler’s model where the gravitational effects take into account the free surface gradient. Finally, the corresponding family of models for the proposed τeff defined in (8) can be reformulated easily keeping Q b as in (14a) , that is, Q b = τeff ρ q b =  q b C | u | u −k 1 ∂ x ( b + h ) −k 2 ∂ x b , but where  q b in (14c) is replaced with  q b = ρQ | τeff| α1 1 −ϕ θβ1 eff ( θeff −α2 θc ) β2 +  θeff −α3  θc β3 , (30) k 1 and k 2 taken the same values introduced in (14b) . 3. Semi-implicit approach We are interested in this work on the large-time scale, i.e., very slow processes where the characteristic time associated to the sediment is very large, and usually the fluid-sediment interaction is weak. This is the case for example of the movement of a dune in a river or a lake. In these situations it is common to have subcritical regimes. Furthermore, note that second order space derivatives of the free surface and the sediment layer appear in system (15) . Then, when discretizing it using an explicit scheme it leads to a very restrictive stability condition (CFL) and therefore to a huge computational cost because of the small time-steps. In this section, we develop an efficient semi-implicit numerical scheme to relax the CFL condition by removing the gravitational contributions, which can be applied to any model with gravitational effect written under formulation (14) - or (29) -and(15) . The spatial and time discretizations are described through the finite volume method. Let us start with the spatial discretization, where we consider a uniform mesh step x without loss of generality. Then, we subdivide the computational domain into control volumes denoted by V i = x i −1 / 2 , x i +1 / 2 with center x i = x i −1 / 2 + x i +1 / 2 / 2 , for i ∈ I , with # ( I ) = M. We consider a staggered mesh, that is, the discrete free surface and sediment variables ( ηi , b i ) are defined at the center of the control volumes ( x i ), whereas the discrete discharge values ( q i +1 / 2 ) are at the interfaces ( x i +1 / 2 ). This C-grid staggering has the advantage that the linear system resulting from the semi-implicit time discretization is more compact. We propose to discretize system (14), (15) using the semi-implicit scheme introduced in [8,10,11] based on the Θ- method. This strategy allows us to remove the celerity contribution to the CFL condition. Thus, defining h i = ηi −b i , u i = (q i −1 / 2 + q i +1 / 2 ) /h i and the Courant numbers C v el = t x max i ∈I | u i | , C cel = t x max i ∈I | u i | +  gh i , with tthe time step, the explicit method has a stability restriction that can be approximated by C cel < 1 , whereas the semi-implicit method relaxes this condition to C v el < 1 . Therefore, in subcritical regime ( | u | /  gh 1 ) this approach allows us to give larger time steps since the restrictive contribution is the gravitational one. The Θ-method (hereinafter Θ-method) reads w n +1 = w n + t Θf t n +1 , w n +1 + t ( 1 −Θ) f ( t n , w n ) , for an arbitrary ODE system w  = f(t, w ) , being Θ∈ [0 , 1] the implicitness parameter. Note that in the limit cases Θ= 0 and Θ= 1 , the explicit and implicit Euler methods are recovered. Although Θ-method is unconditionally stable for Θ∈ [1 / 2 , 1] it is usually choosen Θslightly larger than 1/2, in order to allow for some damping of the fastest linear modes and nonlinear effects. Since the goal of this work is not to propose high order semi-implicit schemes, and for the sake of clarity, we perform a discretization based on the Θ-method for Θ> 1 / 2 , which is first order accurate. However, this procedure can be adapted to more accurate time discretizations if needed via an IMEX-ARK (IMplicit EXplicit Additive Runge Kutta) methods (see [30] ), following the description made in [5,23] . We show in the numerical tests that the Θ-method produces good enough results, in fact we obtain order of accuracy 1.6. Here we describe the scheme corresponding to the Θ-method for the sake of brevity and simplicity in the reading of the paper. We add in Appendix A a discretization based on a formally second order IMEX method. In the numerical tests we show that in a test with regular initial conditions we achieve second order via this IMEX-ARK2 method (see Section 4.5 ). In the following we describe how the system (15) is discretized. The key point is considering all the convective terms in a explicit way, while terms involving the derivative of the free surface ( ∂ x η) and the sediment layer ( ∂ x b) are discretized using the semi-implicit method. Firstly, the discrete momentum equation is written as q n +1 i +1 / 2 = G n i +1 / 2 −t x g( h + Ch m ) n i +1 / 2 ηn +1 i +1 −ηn +1 i −t x g( 1 /r −1 ) ( Ch m ) n i +1 / 2 b n +1 i +1 −b n +1 i , (31) 9 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Fig. 6. Test 4.4 . Evolution of the sediment till t f = 46 days (1.5 months), with the Θ-method and t = 14 s ( C cel ≈8 . 1 ). Dashed red line is the initial condition, dotted green lines correspond to intermediate times and solid blue line is the solution at final time. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). Fig. 7. Test 4.4 . Evolution of the sediment velocity v b (7) (top figure) and the fluid velocity u (bottom figures) for a period of tide t ∈ [120 , 132] h. Green lines correspond to no bedload transport. This test involves a long-time simulation and consequently a big computational effort. It must be pointed out that we are using here a coarse mesh to reduce the computational cost, and that our goal is to show the speed-ups reached for the semi-implicit method with respect to the explicit one. The errors made by the Θ-method with respect to the explicit RK3 method are not significative for the current configurations, as we see in Table 7 . The Courant numbers and the speed-ups reached are shown in Table 8 . We see that the semi-implicit method is 15 times faster than the RK3 method without a significant loss of accuracy, which is an important result. As in previous test, we want to comment that we can go further in C cel , this violating the stability condition based on C v el . 4.5. Accuracy test We perform an order test to evaluate the accuracy of the proposed methods, both the first and the second order (see Appendix A ) semi-implicit approaches. To this end, we consider a test similar to Section 4.3 , but considering very regular initial conditions and without friction between layers, that is, setting C = 0 just in the momentum equation in system (15) . As material properties, we take the ones in Table 3 except for the Manning coefficient and the angle of repose. We choose for 16 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Table 7 Test 4.4 . Erosion process test: Relative errors with respect to the explicit RK3 with C cel = 0 . 1 , and Courant numbers reached at final time t f = 46 days by the semi-implicit method. t(s) C cel C v el Err η[ l 2 ] Err η[ l ∞ ] Err b [ l 2 ] Err b [ l ∞ ] Err u [ l 2 ] Err u [ l ∞ ] (×10 −5 ) (×10 −5 ) (×10 −4 ) (×10 −4 ) (×10 −4 ) (×10 −4 ) 5 2.9 0.35 0.3 0.5 0.4 0.5 1.8 1.2 10 5.8 0.70 0.7 0.9 0.9 0.9 3.6 2.5 15 8.7 1.0 1.0 1.4 1.4 1.5 5.4 3.7 20 11.6 1.4 1.3 1.9 1.8 1.8 7.2 4.9 25 14.4 1.7 1.6 2.3 2.2 2.4 8.9 6.1 Table 8 Test 4.4 . Erosion process test: Speed-ups and Courant numbers reached at t f = 46 days by the semi-implicit method. Method t(s) C cel Comp. time (s) Speed-up RK3 - 0.9 829.3 (13.8 min) 1 Θ-method 5 2.9 266.3 (4.4 min) 3.1 Θ-method 10 5.8 134.4 (2.2 min) 6.2 Θ-method 15 8.7 89.8 (1.5 min) 9.2 Θ-method 20 11.6 66.7 (1.1 min) 12.4 Θ-method 25 14.4 53.7 (0.89 min) 15.4 Fig. 8. Test 4.5 . Initial condition for the height (solid blue line) and the sediment layer (dot-dashed brown line). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). this test n = 0 . 03 and δ= 15 ◦in order to increase the sediment motion and properly evaluate the accuracy of the sediment approximation. We take x ∈ [ −20 0 m, 20 0 m ] , and the initial discharge and sediment layers given by q 0 (x ) = 10 m 2 /s, b 0 (x ) = 2 e x 8 / 10 14 m. In order to have regular initial condition, the initial height is computed as the analytical steady state of the subcritical flow resulting when considering a fixed bottom (no bedload transport), as shown in Fig. 8 . Finally, subcritical boundary conditions are considered, taking q ext = 10 m 2 /s upstream and ηext = 10 m downstream. We perform a time-space test order, starting from a mesh step δX = 4 m (100 nodes) and time step δT = 0 . 192 s , and taking as reference solution the one computed with 3200 nodes ( x = δX / 32 ) and time step 0.0 06 s ( t = δT / 32 ), that is, we keep the Courant number C cel = 0 . 52 . Tables 9 and 10 show the errors and order of accuracy for the semi-implicit methods at t = 50 and t = 3600 s. We see that second order accuracy is achieved for the IMEX-ARK2 at short and long times, while the Θ-method (with Θ= 0 . 55 ) exposed in Section 3 has order 1.6 approximately. Remark 3. We remark that in order to achieve second order accuracy of the IMEX scheme, a key point is the approximation of h i +1 / 2 . In Section 3 we have defined this value as the upwind value since it is a more stable choice in general situations. However, if this approximation is used in the IMEX-ARK2 with the second order upwind formula for the convective term (34) , just order 1.6 is recovered for all the variables. Then, we need to use a second order reconstruction to obtain second order accuracy (in this case it is enough to use h i +1 / 2 = (h i + h i +1 ) / 2 since we have a regular solution). This effect of the interpolation of h i +1 / 2 over the accuracy of the scheme in staggered meshes has been recently studied in [28] . 17 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Table 9 Test 4.5 . Errors and orders ( l 1 -norm) for the free surface, discharge and sediment for the IMEXARK2 and the Θ-method at t = 50 s , with δX = 4 m and δT = 0 . 192 s . ( x , t) Err ηOrder ηErr q Order q Err b Order b IMEX-ARK2 ( δX, δT ) 4.64 ×10 −2 0.00 1.98 ×10 −1 0.00 2.80 ×10 −2 0.00 ( δX/ 2 , δT / 2 ) 1.19 ×10 −2 1.97 6.47 ×10 −2 1.62 7.22 ×10 −3 1.96 ( δX/ 4 , δT / 4 ) 3.20 ×10 −3 1.89 1.80 ×10 −2 1.85 1.80 ×10 −3 2.00 ( δX/ 8 , δT / 8 ) 8.35 ×10 −4 1.94 4.90 ×10 −3 1.88 4.32 ×10 −4 2.06 ( δX/ 16 , δT / 16 ) 1.91 ×10 −4 2.13 1.22 ×10 −3 2.01 8.66 ×10 −5 2.32 Θ-method ( δX, δT ) 5.21 ×10 −1 - 2.33 - 1.36 ×10 −1 - ( δX/ 2 , δT / 2 ) 2.58 ×10 −1 1.01 1.26 0.89 6.73 ×10 −2 1.01 ( δX/ 4 , δT / 4 ) 1.23 ×10 −1 1.07 6.30 ×10 −1 1.00 3.16 ×10 −2 1.09 ( δX/ 8 , δT / 8 ) 5.38 ×10 −2 1.19 2.81 ×10 −1 1.17 1.36 ×10 −2 1.22 ( δX/ 16 , δT / 16 ) 1.82 ×10 −2 1.57 9.57 ×10 −2 1.55 4.53 ×10 −3 1.58 Table 10 Test 4.5 . Errors and orders ( l 1 -norm) for the free surface, discharge and sediment for the IMEXARK2 and the Θ-method at t = 3600 s , with δX = 4 m and δT = 0 . 192 s . ( x , t) Err ηOrder ηErr q Order q Err b Order b IMEX-ARK2 ( δX, δT ) 3.11 ×10 −2 - 5.56 ×10 −3 - 1.27 - ( δX/ 2 , δT / 2 ) 7.73 ×10 −3 2.01 1.44 ×10 −3 1.95 3.17 ×10 −1 2.00 ( δX/ 4 , δT / 4 ) 1.91 ×10 −3 2.02 3.70 ×10 −4 1.96 7.83 ×10 −2 2.02 ( δX/ 8 , δT / 8 ) 4.57 ×10 −4 2.07 9.93 ×10 −5 1.90 1.87 ×10 −2 2.07 ( δX/ 16 , δT / 16 ) 9.20 ×10 −5 2.31 2.58 ×10 −5 1.94 3.74 ×10 −3 2.32 Θ-method ( δX, δT ) 3.59 ×10 −1 - 3.55 ×10 −2 - 6.72 - ( δX/ 2 , δT / 2 ) 1.84 ×10 −1 0.97 1.78 ×10 −2 1.00 3.44 0.96 ( δX/ 4 , δT / 4 ) 8.88 ×10 −2 1.05 8.49 ×10 −3 1.07 1.66 1.05 ( δX/ 8 , δT / 8 ) 3.89 ×10 −2 1.19 3.70 ×10 −3 1.20 7.27 ×10 −1 1.20 ( δX/ 16 , δT / 16 ) 1.31 ×10 −2 1.57 1.25 ×10 −3 1.57 2.44 ×10 −1 1.57 5. Conclusions An efficient semi-implicit scheme for sediment transport models with gravitational effects under subcritical regimes has been proposed. For the sake of simplicity, here we have chosen the generalization of the Ashida & Michiue’s model, which includes gravitational effects through the definition of τeff. However, this method can be immediately adapted to other models, for both the solid transport discharge and the friction model, by redefining Q b in terms of  q b in (14) . Thus, the proposed approach can be adapted to a wide range of families of bedload transport models, as presented in Section 2.4 . We have also shown that the definition of the effective shear stress, τeff, considered in this paper can be seen as an improved formulation of the one proposed by Fowler et al. [19] . An efficient scheme based on the Θ-method has been proposed following [10] , where an iterative Gauss-Seidel algorithm is needed to solve the resulting linear system on (η, b) . This method is easily adapted to a second order method based on an IMEX approach. For the considered model an explicit expression for steady states verifying Q b = 0 . Gravitational terms play a key role on these steady states, since non-physical solutions are obtained if these gravitational terms are neglected. This behavior has been shown in the numerical tests, for both lake at rest and u  = 0 steady states. In particular, for solutions with u  = 0 the slope of the steady states is larger than the angle of repose of the sediment. These gravitational terms also determine the shape of a dune that is swept along by a flow, leading to more realistic (rounded) shapes of the advancing front, where the non-physical shock is corrected. Here we also have shown that a semianalytic approach, where ηand q are assumed as constant values, gives reasonable results in this case. An application to erosion processes have been also performed, where the time of simulation is very large. Finally, a test order is done, showing the expected orders for the proposed semi-implicit methods. In all the cases, we reduce the computational time of simulations with the proposed semi-implicit method while the accuracy is not notably degraded. Acknowledgements This research has been partially supported by the Spanish Government and FEDER through the research project RTI2018- 096064-B-C22. 18 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 Table A1 Butcher tableaux of the explicit ARK2 method. 0 0 2 ∓√ 2 2 ∓√ 2 0 1 1 −(3 + 2 √ 2 ) / 6 (3 + 2 √ 2 ) / 6 0 ±1 2 √ 2 ±1 2 √ 2 1 ∓1 √ 2 c l a lm b l Table A2 Butcher tableaux of the implicit ARK2 method. 0 0 2 ∓√ 2 1 ∓1 √ 2 1 ∓1 √ 2 1 ±1 2 √ 2 ±1 2 √ 2 1 ∓1 √ 2 ±1 2 √ 2 ±1 2 √ 2 1 ∓1 √ 2 c l  a lm b l Appendix A. A more accurate IMEX-ARK2 discretization In this appendix we describe a formally second order semi-implicit discretization based on the IMEX-ARK2 method. We would like to remark that it is difficult to reach second order in practical applications due to the freezing of some coefficient ( h n i , q n b ) in order to not solving a nonlinear system. Other reasons as the complexity of the solid transport discharge used (Ashida-Michiue’s model), where we have sign and positive part functions, could also affect the accuracy of the method. Nevertheless, in the case of regular solutions, avoiding the previous difficulties, the second order accuracy is achived, as has been shown in Section 4.5 . In order to apply the IMEX discretization, the ODE system must be written in additive form as the sum of a stiff (dealt implicitly) and non-stiff (explicitly discretized) contributions as dw dt = f s (t, w ) + f ns (t, w ) . Using the spatial semi-discretization of system (15) , these contributions are f k (t, w ) = f η k (t, w ) , f q k (t, w ) , f b k (t, w )  , k = s, ns, whose stiff components are basically the free surface and sediment gradients xf η s (t, w ) = −q i +1 / 2 −q b,i +1 / 2 x k 1 ( ηi +1 −ηi ) + k 2 ( b i +1 −b i )  −q i −1 / 2 + q b,i −1 / 2 x k 1 ( ηi −ηi −1 ) + k 2 ( b i −b i −1 ) , xf q s (t, w ) = −g ( h + Ch m ) i +1 / 2 ( ηi +1 −ηi ) + g ( Ch m ) i +1 / 2 1 r −1 ( b i +1 −b i ) , xf b s (t, w ) = q b,i +1 / 2 x k 1 ( ηi +1 −ηi ) + k 2 ( b i +1 −b i ) −q b,i −1 / 2 x k 1 ( ηi −ηi −1 ) + k 2 ( b i −b i −1 ) , and the non-stiff terms are convective contributions xf η ns (t, w ) = −q b,i +1 / 2 C i +1 / 2 | u i +1 / 2 | u i +1 / 2 −q b,i −1 / 2 C i −1 / 2 | u i −1 / 2 | u i −1 / 2 , f q ns (t, w ) = −∂ x q 2 /h i +1 / 2 + g ( 1 /r −1 ) ( Ch m ) i +1 / 2 sgn v b,i +1 / 2 tan δ + χi +1 / 2 C i +1 / 2 u i +1 / 2 u i +1 / 2 , xf b s (t, w ) = −q b,i +1 / 2 C i +1 / 2 | u i +1 / 2 | u i +1 / 2 −q b,i −1 / 2 C i −1 / 2 | u i −1 / 2 | u i −1 / 2 . The IMEX method consists of using a consistent combination of an explicit and an implicit method satisfying some coupling conditions. We consider the second order method proposed in [24] , which is defined by the coefficients a lm , c l , b l (explicit counterpart) and  a lm (implicit counterpart) given by the Butcher tableaux (see Tables A.11 and A.12 ). The implicit 19 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 part corresponds with the TR-BDF2 scheme [29] , which is L-stable, and the explicit part is stable under the CFL stability condition based on the velocity (instead of the celerity thanks to the semi-implicit approach). In general, a s-order IMEX-ARK method is written as follows. For l = 1 , . . . , s , the intermediate state w (l) is computed as w (l) = w n + t l−1  m =1 a lm f ns t n + c m t, w (m ) +  a lm f s t n + c m t, w (m )  + t  a ll f s t n + c l t, w (l) , (37) where note that all the terms in the right hand side are explicit contributions, except for the last term, which is multiplied by  a ll and is implicitly discretized. Then, the state w n +1 is w n +1 = w n + t s  l=1 b l f ns t n + c l t, w (l) + f s t n + c l t, w (l) . Let us detail the stages for the particular case of the IMEX-ARK2 that we use. For the first stage we trivially have w (1) = w n . For the second stage, q (2) i +1 / 2 = G (1) i +1 / 2 −t x g  a 22 ( h + Ch m ) (1) i +1 / 2 η(2) i +1 −η(2) i  −t x g  a 22 ( 1 /r −1 ) ( Ch m ) (1) i +1 / 2 b (2) i +1 −b (2) i , (38) with G (1) i +1 / 2 = q n i +1 / 2 + a 21 t f q ns t n + c 1 t , w (1) +  a 21 t f q s t n + c 1 t , w (1) . In addition, it is important to remark that to obtain a second order time-space discretization, the first order discretization of the convective term (32) in the non-stiff contribution for the discharge is replaced by the second order upwind formula x ∂ x q 2 /h i +1 / 2 = ⎧ ⎨ ⎩ (qu ) i −3 2 −4(qu ) i −1 2 + 3(qu ) i + 1 2 / 2 , if u i + 1 2 > 0 , −(qu ) i + 5 2 −4(qu ) i + 3 2 + 3(qu ) i + 1 2 / 2 , if u i + 1 2 < 0 . (39) Now, the discrete continuity and sediment equations are x η(2) i = x ηn i + t x a 21 f η ns t n + c 1 t , w (1) +  a 21 f η s t n + c 1 t , w (1)  −t  a 22  q (2) i +1 / 2 − q b (1) i +1 / 2 x k 1 η(2) i +1 −η(2) i + k 2 b (2) i +1 −b (2) i  −q (2) i −1 / 2 +  q b (1) i −1 / 2 x k 1 η(2) i −η(2) i −1 + k 2 b (2) i −b (2) i −1  (40) and x b (2) i = x b n i + t x a 21 f b ns t n + c 1 t , w (1) +  a 21 f b s t n + c 1 t , w (1)  + t x  a 22  q b (1) i +1 / 2 k 1 η(2) i +1 −η(2) i + k 2 b (2) i +1 −b (2) i  − q b (1) i −1 / 2 k 1 η(2) i −η(2) i −1 + k 2 b (2) i −b (2) i −1 . (41) As in the case of the Θ-method, the values of q (2) i +1 / 2 (33) are injected in (35) and then the new values of the free surface η(2) i and sediment b (2) i are found solving a linear system with the Gauss-Seidel algorithm described in Section 3 . Next, these values are used to compute q (2) i +1 / 2 by (33) . Note that in this case the coefficients of the linear system depend on the coefficient  a 22 instead of the parameter Θ. In the third step we write q ( 3 ) i +1 / 2 = G ( 2 ) i +1 / 2 −t x g  a 33 ( h + Ch m ) ( 2 ) i +1 / 2 η( 3 ) i +1 −η( 3 ) i −t x g  a 3 ( 1 /r −1 ) ( Ch m ) ( 2 ) i +1 / 2 b ( 3 ) i +1 −b ( 3 ) i , (42) with G (2) i +1 / 2 = q n i +1 / 2 + t a 31 f q ns t n + c 1 t, w (1) +  a 31 f q s t n + c 1 t, w (1)  +t a 32 f q ns t n + c 2 t, w (2) +  a 32 f q s t n + c 2 t, w (2)  20 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 collecting the explicit terms. For the free surface and the sediment we have x η(3) i = x ηn i + t x a 31 f η ns t n + c 1 t , w (1) +  a 31 f η s t n + c 1 t , w (1)  + t x a 32 f η ns t n + c 2 t , w (2) +  a 32 f η s t n + c 2 t , w (2)  −t  a 33  q (3) i +1 / 2 − q b (2) i +1 / 2 x k 1 η(3) i +1 −η(3) i + k 2 b (3) i +1 −b (3) i  −q (3) i −1 / 2 +  q b (2) i −1 / 2 x k 1 η(3) i −η(3) i −1 + k 2 b (3) i −b (3) i −1  and x b (3) i = x b n i + t x a 31 f b ns t n + c 1 t , w (1) +  a 31 f b s t n + c 1 t , w (1)  + t x a 32 f b ns t n + c 2 t , w (2) +  a 32 f b s t n + c 2 t , w (2)  + t x  a 33  q b (2) i +1 / 2 k 1 η(3) i +1 −η(3) i + k 2 b (3) i +1 −b (3) i  − q b (2) i −1 / 2 k 1 η(3) i −η(3) i −1 + k 2 b (3) i −b (3) i −1 . Again the values η(3) i and b (3) i are found solving a linear system and then q (3) i +1 / 2 are updated. Then, the new states are computed as ηn +1 i = ηn i + t 3  l=1 b l f η s (t n + c l t, w (l) ) + f η ns (t n + c l t, w (l) ) , q n +1 i +1 / 2 = q n i +1 / 2 + t 3  l=1 b l f q s (t n + c l t, w (l) ) + f q ns (t n + c l t, w (l) ) , b n +1 i = b n i + t 3  l=1 b l f b s (t n + c l t, w (l) ) + f b ns (t n + c l t, w (l) ) . Remark 4. Notice that, as in the case of the Θ-method, the height h , the coefficients Ch m and  q b are linearized in the stiff term f s (t n + c l t, w (l) ) in order to not solving global nonlinear systems, using in each stage the last available state w (l−1) to compute these stiff contributions. In subcritical regimes, which is when semi-implicit approaches are useful and relevant, this assumption is justified, and it has been widely used in environmental fluid dynamics (see e.g. [37] ). References [1] K. Ashida , M. Michiue , Study on hydraulic resistance and bed-load transport rate in alluvial streams, Proc. Jpn. Soc. Civil Eng. 1972 (206) (1972) 59–69 . [2] R.A. Bagnold , The flow of cohesionless grains in fluids, Royal Society of London Philosophical transactions. Series A. Mathematical and Physical Sciences, no. 964. Royal Society of London, 1956 . [3] F. Benkhaldoun , M. Seaïd , S. Sahmim , Mathematical development and verification of a finite volume model for morphodynamic flow applications, Adv. Appl. Math. Mech. 3 (4) (2011) 470–492 . [4] M. Bilanceri , F. Beux , I. Elmahi , H. Guillard , M.V. Salvetti , Linearised Implicit Time-Advancing Applied to Sediment Transport Simulations, Research Report RR-7492, INRIA, 2010 . [5] L. Bonaventura , E.D. Fernández-Nieto , J. Garres-Díaz , G. Narbona-Reina , Multilayer shallow water models with locally variable number of layers and semi-implicit time discretization, J. Comput. Phys. 364 (2018) 209–234 . [6] M.J.C. Díaz , E.D. Fernández-Nieto , A.M. Ferreiro , Sediment transport models in Shallow Water equations and numerical approach by high order finite volume methods, Comput. Fluids 37 (3) (2008) 299–316 . [7] M.J. Castro Díaz , E.D. Fernández-Nieto , A.M. Ferreiro , C. Parés , Two-dimensional sediment transport models in shallow water equations. A second order finite volume approach on unstructured meshes, Comput. Methods Appl. Mech. Eng. 198 (33–36) (2009) 2520–2538 . [8] V. Casulli , Numerical simulation of three-dimensional free surface flow in isopycnal coordinates, Int. J. Numer. MethodsFluids 25 (1997) 645–658 . [9] V. Casulli , A semi-implicit numerical method for the free-surface Navier-Stokes equations, Int. J. Numer. Methods Fluids 74 (8) (2013) 605–622 . [10] V. Casulli , E. Cattani , Stability, accuracy and efficiency of a semi-implicit method for three-dimensional shallow water flow, Comput. Math. Appl. 27 (4) (1994) 99–112 . [11] V. Casulli , R.T. Cheng , Semi-implicit finite difference methods for three-dimensional shallow water flow, Int. J. Numer. Methods Fluids 15 (6) (1992) 629–648 . [12] V. Casulli , G.S. Stelling , Numerical simulation of 3D quasi-hydrostatic, free-surface flows, J. Hydraul. Eng. 124 (7) (1998) 678–686 . [13] V. Casulli , R.A. Walters , An unstructured grid, three-dimensional model based on the shallow water equations, Int. J. Numer. MethodsFluids 32 (20 0 0) 331–348 . [14] V. Casulli , P. Zanolli , Semi-implicit numerical modelling of non-hydrostatic free-surface flows for environmental problems, Math. Comput. Model. 36 (2002) 1131–1149 . [15] F. Charru , Selection of the ripple length on a granular bed sheared by a liquid flow, Phys. Fluids 18 (2006) 121508 . [16] F. Exner , Über die Wechselwirkung zwischen Wasser und Geschiebe in Flüssen, Sitzungsber. d. Akad. d. Wiss. pt. IIa. Bd. 134 (1925) . [17] R. Fernandez-Luque , R. van Beek , Erosion and transport of bedload sediment, J. Hydraul. Res. 14 (2) (1976) 127–144 . [18] E.D. Fernández-Nieto , T.M. de Luna , G. Narbona-Reina , J.D. Zabsonré, Formal deduction of the Saint-Venant–Exner model including arbitrarily sloping sediment beds and associated energy, ESAIM Math. Model. Numer.Anal. 51 (1) (2017) 115–145 . [19] A.C. Fowler , N. Kopteva , C. Oakley , The formation of river channel, SIAM J. Appl. Math. 67 (2007) 1016–1040 . [20] J. Fredsœ, On the development of dunes in erodible channels, J. Fluid Mech. 64 (1) (1974) 1–16 . [21] G. Garegnani , G. Rosatti , L. Bonaventura , Free surface flows over mobile bed: mathematical analysis and numerical modeling of coupled and decoupled approaches, Commun. Appl. Ind.Math. 1 (3) (2011) . 21 J. Garres-Díaz, E.D. Fernández-Nieto and G. Narbona-Reina Applied Mathematics and Computation 421 (2022) 126938 [22] G. Garegnani , G. Rosatti , L. Bonaventura , On the range of validity of the Exner-based models for mobile-bed river flow simulations, J. Hydraulic Res. 51 (2013) 380–391 . [23] J. Garres-Díaz , L. Bonaventura , Flexible and efficient discretizations of multilayer models with variable density, Appl. Math. Comput. 402 (2021) 126097 . [24] F.X. Giraldo , J.F. Kelly , E.M. Constantinescu , Implicit-explicit formulations of a three-dimensional nonhydrostatic unified model of the atmosphere (NUMA), SIAM J. Sci. Comput. 35 (2013) . [25] A.J. Grass, Sediment transport by waves and currents, 1981. [26] P.H. Gunawan , R. Eymard , S.R. Pudjaprasetya , Staggered scheme for the Exner–Shallow Water equations, Comput. Geosci. 19 (6) (2015) 1197–1206 . [27] L. Hascoët , V. Pascual , TAPENADE 2.1 User’s Guide, Technical Report RT-0300, INRIA, 2004 . [28] R. Herbin, J.-C. Latché, Y. Nasseri, N. Therme, A consistent quasi-second order staggered scheme for the two-dimensional shallow water equations, 2021. [29] M.E. Hosea , L.F. Shampine , Analysis and implementation of TR-BDF2, Appl. Numer. Math. 20 (1996) 21–37 . [30] C.A. Kennedy , M.H. Carpenter , Additive Runge-Kutta schemes for convection-diffusion-reaction equations, Appl. Numer. Math. 44 (2003) 139–181 . [31] A. Kovacs , G. Parker , A new vectorial bedload formulation and its application to the time evolution of straight river channels, J. Fluid Mech. 267 (1994) 153–183 . [32] D.K. Lysne , Movement of sand in tunnels, J. Hydraulics Div. 95 (6) (1969) 1835–1846 . [33] E. Meyer-Peter , R. Müller , Formulas for bed-load transport, in: Rep. 2nd Meet. Int. Assoc. Hydraul. Struct. Res., Stockolm, 1948, pp. 39–64 . [34] T. Morales de Luna , M.J.C. Díaz , C.P. Madroñal , A duality method for sediment transport based on a modified Meyer-Peter & Müller model, J. Sci. Comput. 48 (1–3) (2010) 258–273 . [35] P. Nielsen , Coastal Bottom Boundary Layers and Sediment Transport, World Scientific, 1992 . [36] G. Parker , G. Seminara , L. Solari , Bed load at low shields stress on arbitrarily sloping beds: alternative entrainment formulation, Water Resour. Res. 39 (7) (2003) . [37] G. Rosatti , L. Bonaventura , A. Deponti , G. Garegnani , An accurate and efficient semi-implicit method for section-averaged free-surface flow modelling, Int. J. Numer. Methods Fluids 65 (2011) 448–473 . [38] G. Seminara , L. Solari , G. Parker , Bed load at low shields stress on arbitrarily sloping beds: failure of the bagnold hypothesis, Water Resour. Res. 38 (11) (2002) 31–1–31–16 . [39] P.A. Tassi , S. Rhebergen , C.A. Vionnet , O. Bokhove , A discontinuous Galerkin finite element model for river bed evolution under shallow flows, Comput. Methods Appl. Mech. Eng. 197 (33–40) (2008) 2930–2947 . [40] L.C.V. Rijn , Sediment transport, Part I: bed load transport, J. Hydraul. Eng. 110 (10) (1984) 1431–1456 . 22