Verification of Saint-Venant equations solution based on the lax diffusive method for flow routing in natural channels

Hydrodynamic models, based on the Saint-Venant equations, represent the transient flow in water systems, and simulate flow routing over time and space. Several numerical solutions have been applied to these expressions, but often differences are found in results among distinct procedures, considering that all numerical approaches have limitations. Natural channels are characterized by a dynamic flow behavior, with multiple uses of water and changes in morphology, which occur naturally and by human influence. In this context, this study seeks to complement the understanding of different assumptions and approximations for solution of the Saint-Venant equations, verifying results of the Lax diffusive scheme and comparing with HEC-RAS model. The case study consists of a river that covers an urban area (River Iguaçu, located in Curitiba, Paraná). Comparing the observed and simulated hydrographs, it was possible to assess that how the equations are solved, which determines the type of boundary conditions to be used, that may contribute to differences in simulated flows. Additionally, the calibration strategy and the hypothesis of trapezoidal cross-sections produced positive results. The development of this research should also encourage the progress in solving the hydrodynamic model trough the Lax diffusive scheme, particularly concerning computational efficiency. Explicit schemes, although more sensitive to the time step, have the advantage of simplicity to implement, which could be beneficial in the study of complex systems such as urban rivers.


INTRODUCTION
Comprehension and accurate prediction of river flow behavior have fundamental relevance in problems involving water resources, such as flood simulations and water quality studies.In natural channels, a wave propagates under unsteady state, due to processes related to dynamics of rainfall-discharge, reservoirs and sluices operation, modifications of channel geometry, basin drainage patterns or effluent releases and water withdrawals.In addition, land use, especially variable in urban areas, cause major changes in hydrological processes, resulting in significant impacts on rivers discharges (WANG et al., 2014).
In such a context, it is common to apply the Saint-Venant equations to represent one-dimensional problems in channels, such as discharge routing (SERRANO, 2016), flood prediction (ALEKSEEVSKII et al., 2014), dam ruptures (PENG, 2012), and surface and subsurface runoff (HUGHES; LANGEVIN; WHITE, 2015).
The Saint-Venant equations (hydrodynamic model), formed by the expressions for mass and momentum conservation, describe the transient flow, gradually varied, in a channel with irregular cross sections and lateral contribution.Its main advantage is a physical and temporal representation of the longitudinal flow closer to reality, besides the flexibility in simulating modifications in the studied system.
These expressions constitute a group of hyperbolic and non-linear partial differential equations and, therefore, have no analytical solution except for simplified cases.One of the most applied numerical approximation methods for solving the Saint-Venant equations is the characteristics method (DELPHI, 2011;MENDOZA et al., 2011;LOBEIRO, 2012), whose procedure offers the possibility to define boundary conditions more precisely (HASHEMI et al., 2007).
For problems with complex geometry, in two or three dimensions, the finite element or volume methods have also been applied successfully for the hydrodynamic model (ZARMEHI et al., 2011;PENG, 2012;FABIANI;OTA, 2013).These schemes have a relatively complex mathematical basis, and have the same problems common to other numerical approximations, such as approximation errors and conditional stability Finite difference methods, explicit or implicit, are the most commonly procedures for one-dimensional problems.Such schemes are based in the principle of transforming differential equations into algebraic expressions, in which derivatives are converted into finite differences (CHAUDHRY, 1979;GUNDUZ, 2004).
In general, they have first or second order of accuracy, are simple for computational implementation and generate results quickly (LI; JACKSON, 2007).
The usual explicit schemes are: Foward-Time-Centered-Space (FTCS), MacCormack, QUICK, Leap-Frog and Lax.Common implicit representations are Backward-Time-Centered-Space (BTCS), Crank-Nicolson and Preissman.The explicit method is less complex for implementation, but its stability is conditioned to the simulation time interval.For Paiva, Collischonn and Bravo (2011), hydrodynamic modeling through these schemes have become interesting because of the easiness of programming parallelization.The algorithm practicality is an advantage when considering nonlinear problems and simulation of modifications in complex systems, such as urban streams.
In the implicit procedure, the definition of temporal discretization is less restrictive than that of the explicit one, although some studies indicate loss of results quality when increasing the Courant number (GAJDOS; MANDELKERN, 1998;HASHEMI et al., 2007).In this type of method, the solution involves a system of equations, which sometimes causes an increase in the total simulation time due to the size of the matrices (KALITA; SARMA, 2012).
In the discussion presented by Wang et al. (2013), the authors draw attention to the fact that there are often significant differences of results between models due to the different hypothesis and solution methods tested.Since every modeling has simplifications and limitations, in order to achieve correct interpretations one must pay attention to the level of detail required by the problem studied, and understand the implications of the assumptions considered.
In this research verifications for the Saint-Venant equations through the Lax diffusive scheme are conducted, in addition to the study presented by Ferreira, Fernandes and Kaviski (2016), which obtained solutions of the hydrodynamic model associated to water quality (SIHQUAL model).For the same simulation conditions, are shown results generated with the HEC-RAS software (USACE, 2010), a model well established in the literature (PAIVA; COLLISCHONN; BRAVO, 2011;FAN et al., 2012;HAMEED;ALI, 2013).Therefore, the main objective of the study is to evaluate difference of results between distinct procedures for solving the hydrodynamic model, verified for a natural channel that drains an urban area (Iguaçu river), where the flow dynamics is affected by an intense anthropic activity and multiple uses of water.

Hydrodynamic model: governing equations
Considering the wide number of variables that characterize flow in rivers, and the complex channel geometry, some assumptions are made in order to apply the conservation principles represented by Saint-Venant (LIGGETT, 1975): i. unidimensional flow and incompressible fluid; ii. uniform velocity in each cross section, varying over the longitudinal direction; iii. vertical accelerations not considered -hydrostatic pressure distribution; iv. average slope of the channel bottom is sufficiently small to approximate the sine by the tangent of the inclination angle; v. friction losses are not significantly different from those in the steady flow; therefore, the Manning equation or similar may be used; vi. invariable bed channel (erosion or sediment depositions neglected); vii. longitudinal channel axis represented by a rectilinear reach with low bed slope.

Ferreira et al.
Through these simplifications, the Saint-Venant equations may be written as (LIGGETT, 1975): where B represents the cross section top width (m) -which varies with flow depth y (m) -, U is the longitudinal velocity of the flow (m/s), A is the cross section area (m 2 ), q is the lateral contribution per unit of channel length (m 3 /s.m),g is the acceleration of gravity (m/s 2 ), v L is the input velocity of the lateral contribution in the longitudinal direction (m/s), S 0 is the channel bottom slope (m/m) and S f refers to the friction slope (m/m).
In the HEC-RAS model, the Saint-Venant equations are solved in terms of discharge Q (USACE, 2010): where the term ∂z/∂x is the water surface slope.
The mass conservation equation represents the hydrographs damping effects, which occur due to variation of storage capacity.The expression for momentum conservation, in turn, considers the balance of forces acting on the mass of water, which includes gravity, friction, pressure and inertia of the flow.In the latter are expressed the translation effects.
The lateral contribution q represents the inputs and outputs of water.Point contributions include releases (domestic and industrial wastewater or tributaries) and withdraws.Examples of diffuse sources are precipitation and evaporation at surface water, infiltration into the soil and runoff.
Although Manning and Chézy's equations were developed for uniform and steady flow, it is accepted that they are well suited for resistance estimations in open channels with unsteady regime (CHOW, 1959;LIGGETT, 1975).Therefore, the Manning equation can be used to calculate friction slope (STEPIEN, 1984): ( ) where n represents the Manning roughness coefficient and R H the hydraulic radius (m).
According to USACE (1993), other equations may be used to evaluate the term S f , such as Einstein (1950), Simons andSentürk (1976), andASCE (1975).However, they are avoided due to the presence of sediment-related parameters, and the need for iterative solutions.
The Manning roughness coefficient represents the resistance due to friction in the channel, and several factors interfere in its determination, such as: flow event (drought or flood), presence of vegetation and obstructions (bridges and gates, for example), bed material, irregularity of the cross sections, and river alignment (presence of meanders or rectified stretches) (CHOW, 1959;ARCEMENT;SCHNEIDER, 1984).
Some studies point out that the correct characterization of channel bed (such as roughness and slope profile) is the main factor influencing the results of hydrodynamic simulations (SALEH et al., 2013;BENJANKAR;TONINA;MCKEAN, 2015).

Software Hec-Ras
The Hydrologic Engineering Center's -River Analysis System (HEC-RAS) model is a program of public domain, developed by the U.S. Army Corps of Engineers.In this research is used the version 4.1.0,that solves the equations for one-dimensional flow under steady and unsteady state, with modules for sediment transport and water quality (USACE, 2010).
For the Saint-Venant equations solution, the solution is based on the four-point implicit numerical scheme, also called box scheme.Through this method, for each river reach is generated a system of equations, and the simultaneous solution allow that information from the entire river impact the solution at any point.
Consequently, the simulation time interval may be significantly larger than in explicit schemes.Particularities of the method and the discretization of the Saint-Venant equations are presented in USACE (2010).

METHOD The river system input data
The case study is performed for approximately 84 km of Iguaçu River, located in the Upper Iguaçu basin (Paraná), using five known sections (IG2 -upstream, IG3, IG4, IG5, IG6 -downstream), presented in Figure 1.This region concentrates near 30% of Paraná's population (PARANÁ, 2013).Details of the study area and a description of the monitoring points are provided in Ferreira, Fernandes and Kaviski (2016).
In order to solve the hydrodynamic model, the following information are required: geometry of cross sections and distance between them, lateral contribution and Manning coefficient, besides initial and boundary conditions.
The lateral contribution was estimated by the difference between the hydrographs measured at the known points, divided by the distance between them.In consequence, the total balance between water outlets and inlets in the main channel is evaluated.
It was considered that the lateral contribution is evenly distributed along the length of each reach.Considering the lack of information, it was also admitted that the lateral flow entering or leaving occurs in the same main river velocity, a hypothesis also assumed by Steinstrasser (2005) in a similar study.
The initial conditions required for solving Equations 1 and 2 are depth, velocity and channel geometry for time t = 0 along the longitudinal direction.The boundary conditions are these same properties in upstream and downstream sections, known throughout the simulated period.Thus, to obtain the solutions are needed data from cross-sections, hydrographs and rating curves equations of the monitoring points.
The cross sections at the control points were approximated by the trapezoidal shape, assuming they are symmetrical.Consequently, the required information are only bottom width and side slope.In was also considered that these values do not vary over the time interval analyzed.
The observations of elevation and flow measurements were carried out from 1973 to 2013 at IG2, IG3, IG5 and IG6, and from 1999 to 2012 at IG4.However, the available data do not consider the amplitude of maximum flows reached and, therefore, extrapolation of the discharge curves (i.e.rating-curves) was necessary.The procedure used for extrapolation was the logarithmic method, in which the association between flow and depth is approximated by the expression (JACCON; CUDO, 1989): where y is the elevation for the discharge Q, y 0 represents the elevation for a null flow, and a and r are constants.Since the elevations are taken according to an arbitrary referential (fixed and materialized on the surface), y 0 is an unknown parameter, which can be understood as the reading of the water level in the ruler for which the flow is zero (PORTO et al., 2001).
The logarithmic method is widely used because of its simplicity, and is suitable for extrapolations of the upper part of rating curves (PEREIRA; SANTOS; FILL, 2003).The procedure is based on the equations of uniform regime, assuming a regular cross section.In this research it is considered that this approach represents the Iguaçu river conditions sufficiently well.
In order to determine the parameters a, r and y 0 , an iterative calculation is required, through the method of least squares.Table 1 shows the rating curves equations and the determination coefficients (R 2 ) between data observed and adjusted with the expression obtained in each section.

RESULTS AND DISCUSSIONS
The calibration using the Lax method, presented in Ferreira, Fernandes and Kaviski (2016), was performed using data from the year 2010.The verification, on the other hand, is applied in this study for the year 2011.
Manning roughness coefficient was used as calibration parameter, adjusted for each control section.The simulation starts with a characteristic value for the study area, and by a trial and error procedure and with comparisons between measured and simulated data, the required values are obtained.The criterion applied in this calibration was the Nash-Sutcliffe coefficient, E ns (NASH;SUTCLIFFE, 1970): where Q obs is the observed discharge, Q sim represents the discharge simulated with the model and Q med is the average of the flows observed in the period t = 1,2,...N (N is the total number of data).

Explicit scheme solution
The Lax diffusive scheme, presented by the mathematician Peter Lax in 1954, is an explicit finite difference method, simple to apply and with first order accuracy in time and second order in space (TANNEHILL et al., 1997).This method has been shown to be stable and convergent for the solution of the Saint-Venant equations in several studies (STEINSTRASSER, 2005;AKBARI;FIROOZI, 2010;KALITA;SARMA, 2012;PATOWARY;SARMA, 2013).
Complementary information regarding the simulations using the Lax diffusive method is presented in Ferreira, Fernandes and Kaviski (2016).
Since is a procedure with explicit representation, the stability of the generated solutions depends on the Courant condition, given by: where the variable c represents the celerity.
The spatial (Δx) and temporal intervals (Δt) defined were 500 m and 50 s, respectively, yielding mean C r values along the channel ranging from 0.29 to 0.76.
The calibrated Manning coefficients for each point are presented in Table 2, that also shows the Nash-Sutcliffe coefficients obtained for the simulations.The simulated and observed discharges at the five monitoring points for the years 2010 (calibration) and 2011 (verification) are compared in Figures 2 and 3.
The results show that the model simulates adequately the flow oscillations in the period considered, with Nash-Sutcliffe coefficients varying from 0.802 at point IG5 to 0.967 at point IG3.For the verification period ( 2011), E ns values ranged between 0.914 and 0.959.
The period of drought in 2010 extends from August to November, while 2011 shows the driest period between the months of April and July.This confirms that the performed calibration is able to reproduce variations from one year to another.
In the analysis of the hydrographs it is still possible to notice that some flow peaks were overestimated by the model, mainly in the points IG4 and IG5.Upstream of this area there are changes in the hydraulic characteristics, once the river begin to have meanders.Thus, the resulting velocity variations may not have been well represented.Another relevant question is that the adjusted rating curve equation may not be suitable for extrapolation of maximum flows, especially at points IG4 and IG5.

Solution through HEC-RAS
The input data in the HEC-RAS program (cross-sections, Manning coefficients and lateral contribution) were the same as those adopted in the numerical solution using the Lax method.The hydrograph in section IG2 was applied as upstream boundary condition.At the last downstream section (IG6), a series of daily values of depth was obtained using the observed flows and the rating curve equation (shown in Table 1).The initial conditions for the problem were defined as the observed discharges for the first day of 2010 and 2011 in each monitoring section.
In addition, cross sections were interpolated every 500 m (spatial discretization), while the Δt interval used was 60 seconds.Thus, are defined simulation conditions close to the solution using the Lax diffusive scheme.
Figures 4 and 5 show the simulated hydrographs for 2010 and 2011, and Table 3 presents the Nash-Sutcliffe coefficients obtained with these data.The average Courant numbers generated were between 0.61 and 1.40.
In the simulation with the HEC-RAS model there were no overestimations of the maximum flow observed in the numerical solution implemented, as previously described.
In the solution obtained with the Lax method, the discharge is generated from the simulated depths and the equations of the rating curves of each monitoring point; in the procedure assumed by the HEC-RAS, on the other hand, the flow is directly simulated.

Processing time
For the original simulation conditions (intervals Δt = 50 s and Δx = 500 m in the explicit solution, and Δt = 60 s and Δx = 500 m for the HEC-RAS solution), the processing time for Lax solution was approximately 23 min, against 8.30 min in HEC-RAS.This difference occurs mainly because the explicit method has been implemented in Matlab®, while the HEC-RAS is written in FORTRAN language, which processing is usually more efficient.In addition, the data output in the HEC-RAS was defined in daily intervals; when defining smaller periods, the total time increases considerably, generating large output files.

Cross sections
During the simulations a test was conducted with natural cross-sections, instead of trapezoidal sections in the HEC-RAS model (with a 30 min Δt interval).In Table 4, which shows the E ns coefficients obtained in this experiment, it is observed that this change of cross sections generated slightly different results, demonstrating that the trapezoidal geometry approximation is reasonable for the studied case.

Temporal discretization
Simulations with different time intervals are summarized in Table 4.The sensitivity of the explicit scheme to Δt intervals is elevated, and the tests generated unstable solutions.The HEC-RAS solution, on the other hand, proved to be flexible when choosing time discretization, generating results with Nash-Sutcliffe coefficients close to 1.0.

CONCLUSIONS
The solution of the Saint-Venant equations is a research topic that has calling the attention of researchers through time and solved applications in several distinct problems in open channel flows.This work allowed to assess some differences between different numerical schemes for resolution of these expressions, complementing the study by Ferreira, Fernandes   and Kaviski (2016), who presented a numerical solution a water quality hydrodynamic model.The diffusive Lax method, although a simplified computational procedure, was shown to be significantly more restrictive and sensitive to the change in the Δt interval.However, considering the advantages of this type of scheme, it is understood that its application is relevant and should be studied and improved.
The tests performed attested that the explicit method solution presented must be rewritten in a compiled language, such as FORTRAN, to optimize the simulation time.In addition, it was verified that the trapezoidal section hypothesis, which reduces the complexity of the algorithm implementation and the simulation time, is compatible for the case studied.
The comparison of the simulation and verification results between SIHQUAL and HEC-RAS showed that the calibration performed by the roughness coefficient is adequate, with very similar results between the different procedures.
This research also allowed to assess the effect of different solution forms of the Saint-Venant equations, mainly regarding the use of the rating curve as boundary condition.Due to the overestimation of the maximum flows observed only in the solution by the Lax diffusive method, there is an indication that the rating curves for monitoring points IG4 and IG5 are not suitable for extrapolations.
In addition, Domeneghetti, Castellarin and Brath (2012) point out that the uncertainties regarding the rating curve can induce unrealistic estimates for the roughness coefficient.Considering these aspects, defining a suitable relation for the dynamic discharge-elevation in a water body is still a task that offers challenges for the accurate prediction of levels and flows.

Ferreira
et al.Therefore, the potential differences of results are evidenced by the adoption of different boundary conditions in the solution of the Saint-Venant equations.

Table 1 .
Fitted equations for the rating curves at the monitoring stations and R 2 obtained.

Table 3 .
Nash-Sutcliffe (E ns ) coefficients for the HEC-RAS simulation.Points