Comparison of numerical schemes of river flood routing with an inertial approximation of the Saint Venant equations

The one-dimensional flow routing inertial model, formulated as an explicit solution, has advantages over other explicit models used in hydrological models that simplify the Saint-Venant equations. The main advantage is a simple formulation with good results. However, the inertial model is restricted to a small time step to avoid numerical instability. This paper proposes six numerical schemes that modify the one-dimensional inertial model in order to increase the numerical stability of the solution. The proposed numerical schemes were compared to the original scheme in four situations of river’s slope (normal, low, high and very high) and in two situations where the river is subject to downstream effects (dam backwater and tides). The results are discussed in terms of stability, peak flow, processing time, volume conservation error and RMSE (Root Mean Square Error). In general, the schemes showed improvement relative to each type of application. In particular, the numerical scheme here called Prog Q(k+1)xQ(k+1) stood out presenting advantages with greater numerical stability in relation to the original scheme. However, this scheme was not successful in the tide simulation situation. In addition, it was observed that the inclusion of the hydraulic radius calculation without simplification in the numerical schemes improved the results without increasing the computational time.


INTRODUCTION
River's flow and the propagation of flood waves can be represented by hydrodynamic models that solve numerically the Saint-Venant equations.These are formed by the continuity equation (Equation 1), which represents the mass conservation, and by the dynamic equation (Equation 2), which represents the forces acting on the flow (CUNGE; HOLLY; VERWEY, 1980;TUCCI, 2005).The two equations are presented considering a one-dimensional flow. . (1) where, Q is the flow, h is the water depth, t is the time, B is the cross-section width, x is the longitudinal distance, g is the gravitational acceleration, A is the cross-section area, S 0 is the bed slope, and S f is the friction slope.
The left side of the dynamic equation (Equation 2) refers to the local inertial (1) and advective inertial terms (2).The right side represents the external and internal forces acting on the fluid: pressure (3), weight (4) and friction (5).
The Saint-Venant equations can be simplified by eliminating some terms of the dynamic equation, which result in benefits for the models such as simpler formulation, more accessible programming codes, and gain in computational efficiency.For example, Kinematic Wave, Muskingum and Muskingum-Cunge models (CUNGE, 1969;FREAD, 1993;CHANSON, 2004) consider only the terms of weight (4) and friction (5).These terms have a greater contribution to the flow acceleration than other (NÉELZ;PENDER, 2009).The Muskingum-Cunge model is often used as a flow propagation module in rainfall-runoff models (TUCCI, 2005;BEVEN, 2011).However, the traditional linear Muskingum-Cunge method does not adequately represent the flow in some situations where the flood wave celerity is not constant, or in situations with the backwater effect (PONCE, 1989).
Another simplification, considered slightly more complete than Kinematic Wave, Muskingum and Muskingum-Cunge models, is the result of using only the terms of pressure (3), weight (4) and friction (5), called non-inertial (YEN;TSAI, 2001) and diffusion wave models (YEN;TSAI, 2001;CHANSON, 2004).The diffusion wave model has the advantage of representing more adequately the attenuation of a flood wave propagating downstream.However, it has limitations in comparison with the complete Saint-Venant equations, as in representation of downstream effects (MONTERO et al., 2013).
The inertial model, or local inertial, disregards only the advective inertial term (2).This model, popularized by Bates, Horritt and Fewtrell (2010), was widely applied to represent two-dimensional flow in flood areas (NEAL et al., 2011;ALMEIDA et al., 2012;ALMEIDA;BATES, 2013).Montero et al. (2013) demonstrated that the inertial model has advantages over the diffusion wave model when compared to complete solutions of the Saint-Venant equation.Fan et al. (2014) tested the inertial formulation for representing one-dimensional flow in rivers.The authors showed their applicability in scenarios with high and low rivers' slope and subjected to downstream effects, such as backwater and tide.The results were comparable to the complete Saint-Venant equations.However, it was also observed that the model presents numerical instability for values of the Courant number much smaller than 1 (the stability condition requires that t gh x ∆ ∆ must be smaller than 1).
Furthermore, Monteiro et al. (2015) tested the inertial solution for the simulation of waves caused by the closure of floodgates, and also obtaining results comparable to the complete solutions of the Saint-Venant equations.
The main advantage of the inertial model is the easy implementation with simple code and good results in relation to the solution of the Saint-Venant equations.As a result, the model has been used in hydrological models and flood models by several authors (DOTTORI; TODINI, 2011;NEAL;SCHUMANN;BATES, 2012;ALMEIDA;BATES, 2013;COULTHARD et al., 2013;YAMAZAKI;ALMEIDA;BATES, 2013;SAMPSON et al., 2015).In Brazil, the authors Pontes et al. (2015Pontes et al. ( , 2017) ) presented coupling of the one-dimensional inertial model in a large-scale hydrological model (MGB-IPH).The authors pointed out that with this model it is possible to perform site simulations, as in floodplains, more adequately than when using other simplified methods.They showed that the scheme has robustness in coupling with hydrological models and, in addition, the parallel processing is easy.
As commented, several studies point to the inertial model as an alternative to the complete solution of the Saint-Venant equations and to other simplified models representing the one-dimensional flow.However, alternatives to improve the model's numerical stability may still be explored.Therefore, this paper proposes and tests some changes in the inertial formulation in order to improve the numerical stability.The numerical scheme of the one-dimensional inertial model, as applied by Fan et al. (2014), was altered in relation to the approximation of the numerical derivative, in relation to the numerical method and also in relation to the dynamic equation's friction term (term 5 in Equation 2).Six new schemes were proposed and tested considering the hydraulic radius without simplification ( ) and six numerical schemes with hydraulic radius simplification ( R h ≈ ).The performance of the schemes was compared with original scheme, presented by Bates, Horritt and Fewtrell (2010), and with the complete solution of the Saint Venant equations (HEC-RAS hydrodynamic model -USACE, 2016), i.e. the HEC-RAS model was considered the "true solution" in tests.This model solves the equations by an implicit finite difference scheme described in USACE (2016).For the sake of isonomy with the knowledge already generated and published on the one-dimensional inertial model application, scenarios identical to those used by Fan et al. (2014) were considered.The scenarios were: four situations of river's slope (normal, low, high and very high) and two situations where the river is subject to downstream effect (backwater and tide).

ONE-DIMENSIONAL INERTIAL MODEL
The dynamic equation (Equation 2) in the one-dimensional inertial model is rewritten (Equation 3) considering a rectangular river cross-section, representing the pressure and weight forces in terms of water surface level (y) and estimating the friction force using the Manning equation, where y is the water surface level, n is the Manning roughness coefficient, R is the hydraulic radius.
The hydraulic radius is approximated by the depth ( R h ≈ ), considering that the rivers cross sections have a width much larger than the depth.The derivatives of the dynamic equation are approximated by a numerical scheme of progressive finite differences in space and time, which results in Equation 4, as presented in Fan et al. (2014).
Where i and k are the indices in space and time, respectively.
x ∆ is the length of the channel section (m). The term is the depth in the contour of the channel stretch ( x ∆ ) and can be estimated by the difference of the maximum value of the water level ( y ) and the maximum value of the bottom level ( z ) in the center of the channel stretch, according to Equation 5.

[ ]
; ; From the explicit solution of Equation 4, the flow is estimated at each time step according to Equation 6. Bates, Horritt and Fewtrell (2010) considered the multiplication of the flow in time step k by the flow in time step k 1 + in friction term to increase the stability of the numerical scheme.
Based on the flow values ( , the continuity equation (Equation 1), discretized by a spatial and time progressive scheme (Equation 7), is applied to estimate the depth at all reaches i at the end of the time interval ( k 1 + ).From this value, the water level at the end of the time interval is computed by Equation 8.
Equation 6 requires the initial values of the depth ( 0 h ) and the water level ( 0 y ) in all reaches.For this, the same initial flow is assumed in all reaches and the depth is calculated by considering the continuous and uniform flow, that is, the Manning equation with simplification of the hydraulic radius ( R h ≈ ).In addition, boundary conditions are required in the first and last channel stretch.The typical upstream condition is a flow hydrograph and the downstream condition can be a rating curve or a normal depth.
The numerical scheme is restricted to the size of time step (Δt) and channel stretch (Δx) to avoid numerical instability.Thus, the choice of Δt and Δx must respect the Courant-Friedrichs-Levy condition (Courant number -C), given by Equation 9.
The depth (h) used in this calculation corresponds to the largest depth found in the river reach.Although this condition requires a value less than 1, Bates, Horritt and Fewtrell (2010) observed that the value must be less than 0.7 to ensure the numerical scheme stability in two-dimensional simulations.
After the initial development of the model, Almeida et al. (2012) proposed a modification in the numerical scheme to improve stability in low friction situations, such as in urban areas.In the new scheme, a weighting (θ) was introduced in the flow derivative in relation to the time.Thus, the flow value at time k and reach + is weighted considering the flow at time k and in the reach − , and + , according to Equation 10.
The weighting (θ) varies from 0 to 1. Using Equation 10 to approximate the numerical derivative, the flow in each time step in the new scheme is estimated according to Equation 11.
According to Almeida et al. (2012), the inclusion of the weighting allowed more stability in simulations performed by the authors.

MODIFICATIONS PROPOSED TO THE NUMERICAL SCHEME OF THE INERTIAL MODEL
In the inertial model formulation proposed by Bates, Horritt and Fewtrell (2010), the flow at time step k 1 + (Equation 6) is solved by the Euler method, i.e. a first order method where the first derivative provides a direct estimate of the slope for entire interval (Δt).In this method, the global truncation error is proportional to the time step size (CHAPRA; CANALE, 2012).The second-order Runge-Kutta method (RK2), equivalent to the Heun method without iteration, and the fourth-order Runge-Kutta method (RK4) are more accurate because the errors are of the order of 2 t ∆ and 4 t ∆ , respectively (CHAPRA; CANALE, 2012).Thus, the reduction of the step size decreases the error at a faster rate than in the Euler method (CHAPRA;CANALE, 2012).Another way to improve the accuracy is to use a centered derivative (Cent) instead of the progressive derivative (Prog).This derivative allows a second-order approximation, with an error of the order of 2 t ∆ .In this paper, the modifications of the numerical schemes were done in the derivative approximation (Prog and Cent), in the numerical method (Euler, RK2 ou RK4) and also in the friction term (Equation 2), specifically in the time index ( k ) of the flow The properties of the proposed numerical schemes, such as convergence, consistency and stability, were not evaluated.These analyzes indicate how close the numerical scheme is to the analytical solution, and if the errors are not amplified (POPESCU, 2014).Thus, the numerical schemes were only tested and evaluated in different numerical applications.In each test, 16 simulations were run, each one with the methods listed in Table 1 and described in detail below ("R" means that the hydraulic radius was not simplified).

Scheme Q(k+1) x Q(k-1):
The first modification made in the original numerical scheme was to approximate the flow by a time-centered derivative, which considers the flow in the time steps k+1 and k-1.Thus, the proposed scheme also considered the flow in these time steps . Rearranging, the flow k 1 estimated directly by Equation 12.This scheme will be identified by the abbreviation: In this numerical scheme, the friction term was calculated by the flow at time step k+1, k 1 k 1 . Equation 13shows the scheme for a progressive derivative, identified as Prog Q(k+1) x Q(k+1), and Equation 14 for a centered derivative, identified as Table 1.Summary of proposed methods and characteristics.
Note that both equations are second-degree polynomial equations: 2 ax bx c 0 + + =.Therefore, the flow can be calculated directly by the Bhaskara's formula, without the need for an iterative method for resolution.The desired estimate is always the largest root given by: ) is always negative, as shown below.
It is known that the roots sum is given by = and a 0 > ( , , , ∆ are always positive), we have Hence 2 x 0 < , otherwise the sum of the two roots would be positive.Thus, the root 2 x does not have physical importance as an always negative flow.On the other hand, the root 1 x can be positive or negative depending on the term c.As .x is positive if c 0 < , and negative if c 0 > .In fact, if the flow is negative, as occurs when the river is subject to the downstream effect, the term c will be positive and the flow at time step k+1, + , is negative.The scheme is restricted to a new stability condition, in addition to the Courant-Friedrichs-Levy condition, where the discriminant of the Bhaskara's formula must be positive ( 2 b 4ac 0 − > ) so that the root is not complex.Thus, the condition is given by . When the flow is positive ( 1x 0 > ) it was shown that a 0 > and c 0 < , satisfying the condition.However, when the flow is negative ( 1 x 0 < ), the term c is also negative and this condition is mainly dependent on the time step size ( t ∆ ), the channel stretch size ( x ∆ ) and the level variation between reaches ( ).

Scheme Q(k) x Q(k):
Finally, the original numerical scheme was modified by considering the flow at time step k . This scheme allows to calculate the flow directly (Euler) and by the second-and fourth-order Runge-Kutta methods.The flow in time step k+1, calculated by the Euler method, is given by Equation 15 and it is identified as Prog Q(k)xQ(k). .
In the second-order Runge-Kutta method, the flow determined by the Euler method is an intermediate result, where the slope estimate is calculated at the initial point (predictor equation).In the next step, a new estimate is calculated with the predictor result, i.e., the slope estimate is made at the end point.Then, the mean of the two derivatives is calculated to obtain an improved estimate of the slope over the whole range (broker equation).Therefore, Equation 15 is rearranged to define the derivative function ( ) 16), the predictor step is given by Equation 17(equivalent to Equation 15) and the corrector step by Equation 18.This numerical scheme is identified by the Multiple slope estimates are developed in the fourth-order Runge-Kutta method to obtain an improved mean slope (Equations 19 to 22).Each K value represents a slope and the final equation is a weighted average (Equation 23).This scheme is identified by the abbreviation Prog RK4 Q(k)xQ(k).Note that K 1 is the slope at the interval beginning and K 2 is the slope at the interval final, equivalent to the second-order Runge-Kutta method.

EXPERIMENTS
The described numerical schemes were evaluated in different applications to compare the performance in relation to the same schemes adapted by Fan et al. (2014) of the formulation proposed by Bates, Horritt and Fewtrell (2010) and the formulation proposed by Almeida et al. (2012).The schemes were implemented in MATLAB and the simulations were run on eight cores of four Quad-core 3.4 GHz Intel Xeon processors with 8 MB cache memory.The schemes were evaluated for stability (maximum Courant number for which the method remains stable), peak flow at the river's last reach (Qpeak), volume conservation error (VolErr), and the flow's mean square error (RSME), calculated based on the results of simulations performed in the HEC-RAS hydrodynamic model.In addition, each proposed scheme was solved with and without the hydraulic radius simplification, that is, R h ≈ and ( ) The applications were equivalent to the scenarios tested by Fan et al. (2014).The first test considered a real reach of the São Francisco River, located between Três Marias UHE and the Pirapora city (MG), without lateral contributions.Tests 2, 3 and 4 considered the same reach but with low, high and very high slopes, respectively.In tests 5 and 6, downstream effects were considered: backwater and tide.The hydrograph input was the same as used by Fan et al. (2014).The simulated period in all tests was 150 hours.Table 2 presents a summary of the data used in the six tests.

Test 1
This test considered a real reach of the São Francisco River between the Três Marias UHE and the Pirapora city (MG) with approximate values of 135 km in length, 300 m in width (rectangular cross-section), river's slope of 29.5 cm/km and Manning's roughness coefficient of 0.030.The spatial discretization (Δx) was 2 km and time discretization (Δt) was varied between 1 and 4 minutes, assuming C values equal to 0.26, 0.51, 0.77 and 1.03 at the time of greater depth.A constant slope was considered as a downstream contour condition.
The six formulations of the numerical models described (Equations 12,13,14,15,18 and 23), the formulation proposed by Bates, Horritt and Fewtrell (2010) (Equation 6) and the formulation proposed by Almeida et al. (2012) with a weighting factor equal to θ = 0.9 (Equation 11), adopted based on the authors' results, were tested with and without the hydraulic radius simplification.
All stable schemes presented values similar to the HEC-RAS model since the highest RMSE was 51.48 m 3 .s - for the Prog Q(k+1)xQ(k+1) scheme.Figure 1 shows only the region of the hydrograph peak flow in the river's last reach for the original scheme, the Cent Q(k+1)xQ(k-1), Prog Q(k+1)xQ(k+1), Prog RK2 Q(k)xQ(k) schemes and the HEC-RAS model, with and without the hydraulic radius simplification.Figure 2 shows the results of the parameters evaluated by varying the time step.In the first column, the results correspond to simulations with the hydraulic radius simplification, and in the second column without.All values are shown in Table A1 of Appendix A.
The hydraulic radius calculation without simplification did not increase the processing time of the schemes and their use reduced RMSE.As an example, the scheme of Bates, Horritt and Fewtrell (2010) presented a 60% lower RMSE when the hydraulic radius was calculated.Despite this, the volume conservation error increased from 0% to 0.05%, but still remained negligible.This parameter (ErrVol) did not vary with the time step nor with the numerical scheme.It was also observed that the calculation of hydraulic radius did not influence the numerical schemes' stability.
The Almeida et al. (2012) scheme presented no advantages compared to the Bates, Horritt and Fewtrell (2010) scheme in terms of stability (see Table A1 of Appendix A).These schemes were unstable with a time step of 3 min (C = 0.77).All schemes with centered derivative presented instability, with stable results only for the time step of 1 minute (C = 0.26).In addition, these schemes did not improve the RMSE and the volume conservation error, as was expected for a second-order approximation scheme with respect to a first-order scheme.

Test 2
This test considered the same conditions of the river reach used in test 1 but it used a river's low slope of 5 cm/km.The time step was varied in 1, 2, 3 and 3.5 minutes, corresponding to C values of 0.29, 0.60, 1.90 and 1.04, respectively, at the time of greater depth.As in test 1, the numerical schemes' performance was compared.Figure 3 shows the peak flow region in the last reach for the original scheme, the Prog Q(k+1)xQ(k+1), Prog RK2 Q(k)xQ(k) schemes and the HEC-RAS model.Figure 4 shows the performance of all schemes evaluated and Table A2 ( Appendix A) shows the values of Figure 4.
The Prog Q(k+1)xQ(k+1) scheme was the most stable, with stability up to the time step of 3.5 min (C = 1.04), and with an increase in processing time of 20% relative to the original scheme.In addition, it was the scheme that presented the highest RMSE values.
Schemes with centered derivative presented instability at all time steps adopted.The other schemes, such as Bates, Horritt and Fewtrell (2010) and Almeida et al. (2012), presented equal stability up to a time step of 3 min (C = 0.90).
The Prog Q(k)xQ(k) scheme presented the lowest RMSE with practically the same processing time as the original model.The Prog RK2 Q(k)xQ(k) and Prog RK4 Q(k)xQ(k) schemes presented RMSE values closer to the ones in Bates, Horritt and Fewtrell (2010), but the processing time was, on average, almost twice.
The volume conservation error did not changed with the time step but increased from 0.57% to 0.75% when the calculation of the hydraulic radius was included in the numerical scheme.This inclusion of this calculation did not increase the processing time, as observed in Figure 4, and led to the results closer to that of the HEC-RAS model, reducing RMSE by 64% compared to the Bates, Horritt and Fewtrell (2010) scheme.

Test 3
In this test, a high slope of 300 cm/km was considered.The time step was varied between 1 and 4 minutes corresponding to C values of 0.18, 0.37, 0.55, 0.74 at the time of greatest depth.Figure 5 shows a section of the hydrograph at the river's last reach for the original scheme, the Cent Q( k+1    The inclusion of the hydraulic radius calculation in numerical schemes reduced the RMSE by 20% in the original scheme and by 30% in the Prog Q(k+1)xQ(k+1) scheme.Moreover, the increase in the volume conservation error from 0 to 0.01% is negligible considering this calculation.

Test 4
The same river's reach of test 1 with a very high slope of 10 m/km was considered in this test.This slope changes the flow regime to the supercritical regime, where the original inertial model does not show good stability results (FAN et al., 2014).The time step was varied with values of 80, 120, 160 and 200 seconds, corresponding to the values of C at the time of greatest depth of 0.20, 0.31, 0.41, and 0.51, respectively.Figures 7 and 8     conservation error, processing time and RMSE.In addition, the processing time of these schemes was the same as the original scheme with a slightly higher RMSE, 31.24m 3 /s compared to 30.19 m 3 /s.
As in the other tests, the inclusion of the hydraulic radius calculation in the numerical schemes did not increase the processing time, reduced the RMSE from 30.19 to 25.75 in the model of Bates, Horritt and Fewtrell (2010) while the volume conservation error increased from 0 to 0.01% in all stable numerical schemes, but may be considered negligible.

Test 5
This application was identical to test 1 but considered a constant water level downstream and a constant initial level from the half of the reach as boundary conditions.These conditions represent a dam operating at constant water level in simplified form.Other models of flood propagation, such as Muskingum-Cunge, are invalidated in the presence of dam and backwater (FAN et al., 2014).The time step was varied with values of 1, 2, 2.2 and 2.5 minutes, corresponding to the maximum C values of 0.44, 0.88, 0.97 and > 1 at the moment of greatest depth (higher values computed during the simulation).Figure 9 shows a hydrograph of the last reach for the original, Prog Q(k+1)xQ( k+1) and Prog RK2 Q(k)xQ(k) schemes.Figure 10 shows the main results and Table A5 (Appendix A) the detailed results.
All numerical schemes, with the exception of the centered derivative schemes, showed stable results until the adopted time step of 132 sec (C = 0.97).The processing times of the Prog Q( k+1)xQ( k+1) and Prog Q(k)xQ(k) schemes presented an average increase of 25% in comparison to the original scheme.On the other hand, the Prog RK2 Q(k)xQ(k) and Prog RK4 Q(k)xQ(k) schemes took twice as long as the original scheme.
In this application, the volume conservation error did not change with the inclusion of the hydraulic radius calculation.In addition, the peak flow approached the HEC-RAS model's peak flow and the RMSE decreases.

Test 6
In this test, we considered a situation equal to test 2, low slope, but with different contour conditions.A constant flow of 500 m 3 /s was adopted as the upstream condition and a variable

CONCLUSIONS
This paper compared the performance of new numerical schemes with inertial approximation of the Saint-Venant equations.The differences between the numerical schemes in terms of level as the downstream condition.The water level downstream (y) was given by the equation , which represents a periodic level variation between 4m and 8m with a period of 12 hours.This downstream condition -tidal variation -causes a wave that propagates upstream and which can reverse the flow direction in the river's last reach during some moments along the tide cycle.The time step was varied between 1, 2, 3.5 and 4 minutes with maximum values of C equal to 0.27, 0.54, 0.94 and 1.07 (highest values computed during the simulation), respectively.Figure 11 shows a region of the hydrograph peak in the last reach for the original scheme and for the Prog RK2 Q(k)xQ(k) scheme.Figure 12 shows the main results and Table A6 (Appendix A) the detailed results.
Unlike the other applications, in this test the Prog Q(k+1)xQ(k+1) scheme was not stable in any adopted time step.As shown, the scheme requires that the discriminant of the Bhaskara's formula be positive and this criterion was not fulfilled even with a very small time step.
The most stable scheme was Prog RK2 Q(k)xQ(k) with C equal to 1.07 (time step of 4 minutes), followed by the other schemes, with C = 0.94 (3.5 min).Despite the improvement in stability, the Prog RK2 Q(k)xQ(k) scheme took, on average, twice the original scheme time and presented higher RMSE.
The volume conservation error and RMSE were equal in all stable schemes, approximately 0.97%, and had a slight increase when the calculation of the hydraulic radius was considered.Prog RK2 Q(k)xQ(k), did not show stability in these tests (3 and 4).For this scheme, the stability was equivalent to original model in tests 2 and 5 and greater in tests 1 and 6.The disadvantage of the schemes with second and fourth order approximations is the longer processing time.
The parabola scheme, whose solution is a second degree equation (Prog Q(k+1)xQ( k+1)), presented advantages in relation to the Bates, Horritt and Fewtrell (2010) scheme in most applications.In tests 1, 2, 3 and 4 the proposed model was the most stable, presenting good results with C values of 1.03, 1.04, 0.74 and 0.41, respectively, compared to the maximum values of 0.77, 0.90, 0.37 and 0.20 for the original scheme.In test 5 (dam and backwater), the stability was equal in both schemes (C = 0.97).However, in test 6 (tide), the parabola scheme was not stable at any time step.The proposed scheme is more unstable in negative flow situations.
The parabola scheme's processing time was equal to or slightly higher than the original model.This difference is in the tenths of a second scale and may be considered negligible.Moreover, in applications that do not occur negative flow, the greater stability makes this scheme faster because it allows to adopt a greater time step.As an example, in the supercritical regime application (test 4) the proposed scheme spent 0.58 sec (time step 160 sec) and the original scheme spent 1.12 sec (time step 80 sec).Regarding RMSE, despite the higher values in the parabola scheme in tests 2 (low slope) and 5 (dam and backwater), the difference with the original model is very small.The largest difference found was 2.3 m 3 /s (RMSE of 22.52 m 3 /s compared to 24.82 m 3 /s) in test 5.
Therefore, among all tests performed, we believe that the important information raised in the research is that the greater stability of the parabola scheme makes this scheme more attractive than the original model for use as flow propagation modulus in hydrological models.However this use is limited  stability, peak flow, volume conservation error and RMSE were investigated.All stable schemes presented similar results to the hydrodynamic model HEC-RAS, and can be used with good accuracy in comparion to the model that considers the complete equations of Sant-Venant.The advantages of the proposed schemes are ease of implementation, robustness in the coupling with hydrodynamic models and the possibility of parallelization.However, more comparison of the inertial model with implicit and explicit numerical schemes is recommended for future studies.
With respect to performance measures, Table 3 presents the best scheme(s) in each metric and in each test.The difference of the peak flow and the volume conservation error between schemes can be considered negligible.Furthermore, although the original scheme is faster, the difference in processing time is very small in relation to the Prog Q(k)xQ(k) and the Prog Q(k+1)xQ(k+1) schemes.
It was observed for all stable numerical schemes that the peak flow approaches the HEC-RAS model's peak flow and that the RMSE decreases when the hydraulic radius calculation is considered without simplification.In this way, the inertial model formulation is more similar to what is done in the HEC-RAS, which does not simplify the hydraulic radius calculation.In addition, the processing time did not increase with this calculation and the volume conservation error showed a negligible increase.
In general, the scheme proposed by Almeida et al. (2012) presented no advantages in relation to the Bates, Horritt and Fewtrell (2010) scheme.In all tests, both models presented instability at the same time step and both RMSE and peak flow were equivalent.The processing time was equal or greater in the Almeida et al. (2012) scheme.
Schemes with centered numerical derivatives showed instability for low C values in tests 1, 2, 5 and 6.Only in the applications of test 3 (high slope) and test 4 (very high slope) the Cent Q(k+1)xQ(k+1) scheme showed an advantage over the others.On the other hand, the second order approximation scheme, in some cases where there are negative flows, which reduces one of the main applications of models using the inertial model.This disadvantage could easily be overcome by using an adaptive mixed method which recognizes the possibility of a flow reversal and adapts temporarily to another scheme and then returns to the parabola scheme as soon as the negative flow disappears.
Finally, it is important to emphasize that, according to the results, the inclusion of the hydraulic radius calculation without simplification in the numerical schemes improved the results without increasing the computational time.That is, the inclusion of this consideration approximates the results of physical reality without additional performance costs.
APPENDIX A. SUPPLEMENTARY DATA.

Figure 2 .
Figure 2. Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a river's slope equal to 29.5 cm/km.

Figure 3 .
Figure 3. Region of the hydrograph peak flow in the river's last reach for the Bates, Horritt and Fewtrell (2010), Prog Q(k+1)xQ(k+1) and Prog RK2 Q(k)xQ(k) schemes with and without the hydraulic radius simplification considering a river's slope equal to 5 cm/km (the red line is subscribed by the green line).
schemes and HEC-RAS.Figures 6 shows the main results and TableA3(Appendix A) the detailed results.As in the previous tests, the Prog Q(k+1)xQ(k+1) scheme was the most stable, reaching the value of C equal to 0.74 compared to 0.37 in theBates, Horritt and Fewtrell (2010) andAlmeida et al. (2012) schemes.The increase in processing time compared to the original scheme was 20%.The Cent Q(k+1)xQ(k+1) scheme also presented better performance than the original model, with stability up to the present the main results and Table A4 (Appendix A) the detailed results.Only the numerical schemes Prog Q(k+1)xQ(k+1), Cent Q(k+1)xQ(k+1),Bates, Horritt and Fewtrell (2010) andAlmeida et al. (2012) presented stable results.The last two were stable only in the time step adopted in the simulation of 80 sec, corresponding to C = 0.20.The other two schemes were stable up to 160 sec (C = 0.41).Despite the improvement, where the Courant number increased from 0.20 to 0.41, the value is still small compared to other models of flood propagation where C can reach 1.The Prog Q(k+1)xQ(k+1) and Cent Q(k+1)xQ(k+1) schemes showed practically identical performances with near peak flow, volume

Figure 5 .
Figure 5. Region of the hydrograph peak flow in the river's last reach for the Bates, Horritt and Fewtrell (2010), Cent Q(k+1)xQ(k-1) and Prog Q(k+1)xQ(k+1) schemes with and without the hydraulic radius simplification considering a river's slope equal to 300 cm/km.

Figure 6 .
Figure 6.Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a river's slope equal to 300 cm/km.

Figure 4 .
Figure 4. Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a river's slope equal to 5 cm/km.

Figure 7 .
Figure 7. Region of the hydrograph peak flow in the river's last reach for the Bates, Horritt and Fewtrell (2010) and Prog Q(k+1)xQ(k+1) schemes with and without the hydraulic radius simplification considering a river's slope equal to 10 cm/km.

Figure 8 .
Figure 8. Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a river's slope equal to 10 m/km.

Figure 9 .
Figure 9. Region of the hydrograph peak flow in the river's last reach for the Bates, Horritt and Fewtrell (2010), Prog Q(k+1)xQ(k+1) and Prog RK2 Q(k)xQ(k) schemes with and without the hydraulic radius simplification considering a boundary condition downstream of dam and backwater (the red line is subscribed by the green line).

Figure 11 .
Figure 11.Region of the hydrograph peak flow in the river's last reach for the Bates, Horritt and Fewtrell (2010) and Prog RK2 Q(k)xQ(k) schemes with and without the hydraulic radius simplification considering a boundary condition with tide variation (the red line is subscribed by the green line).

Figure 12 .
Figure 12.Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a boundary condition with tide variation.

Figure 10 .
Figure 10.Peak flow, volume conservation error, processing time and RMSE in relation to the time step for different formulations of numerical schemes considering a boundary condition downstream of dam and backwater.
Results of test 2 (low slope) for different numerical schemes.
Results of test 3 (high slope) for different numerical schemes.
Results of test 5 (backwater) for different numerical schemes.Results of test 6 (tide) for different numerical schemes.

Table 2 .
Summary of the parameters used in tests.
*Flow velocity calculated by the initial flow, width and initial depth (Manning's Eq. with simplification of the hydraulic radius).

Table 3 .
The best numerical scheme(s) in each test in relation to metrics.

Table A1 .
Results of test 1 (real reach) for different numerical schemes.

Table A4 .
Results of test 4 (very high slope) for different numerical schemes.