## Abstract

Time-dependent, unsteady flow has been studied in prismatic open channels with symmetric trapezoidal and triangular cross sections and small bottom slope. The St Venant equations without lateral inflow have been discretized in explicit as well as in implicit form and solved numerically, for unsteady, subcritical flow. The inflow hydrograph used can be applied for different flood events by adjusting its parameters accordingly. The results presented are derived from the explicit schemes Lax-Diffusive, MacCormack, Lambda as well as the implicit Preissmann scheme, and are compared to those from the Muskingum-Cunge method and the widely used commercial software HEC-RAS. The peak flow computed by the Lax-Diffusive scheme was reduced at the downstream end of the channel and the arrival time of the peak increased if compared to the other methods. The Muskingun-Cunge method forecasted the shortest peak flow arrival time at the downstream end cross section. Mass conservation computed from inflow and outflow hydrographs has been confirmed, since the maximum error did not exceed 2.60%. All codes were implemented in house using Matlab^{®}.

## INTRODUCTION

*t*and the longitudinal distance

*x*measured from the upstream end of the channel. The dynamic wave model includes the continuity (mass conservation) equation, which, if there is no lateral inflow, is written as (Chaudhry 2008): and the momentum equation is written as: where

*y*

*=*

*y(x,t)*is the unknown flow depth,

*V*

*=*

*V*(

*x*,

*t*) is the cross-sectional average velocity to be determined,

*D*is the hydraulic depth (wetted area to free surface width ratio),

*S*is the longitudinal bottom slope,

_{o}*S*

_{f}= n^{2}

*V*|

*V*|/

*R*

^{4/3}is the friction slope computed from Manning formula, accounting also for reverse flow conditions, and

*g*is the gravitational acceleration. Equations (1) and (2) were derived under the following assumptions: (i) the velocity distribution is uniform at a cross section, (ii) the pressure distribution is hydrostatic, (iii) the bottom slope is small, and (iv) the water is homogeneous and incompressible.

The above equations are nonlinear without analytical solutions, except for special cases where the linearized form of them can be applied. Earlier studies have dealt with the numerical solution of the equations either in the direct or in the characteristic form. The direct form of equations has been solved by finite difference and finite element methods, in explicit form (Liggett & Woolhiser 1967; Strelkoff 1970; Retsinis *et al.* 2018), or implicit form (Preissmann 1961; Amein 1968; Baltzer & Lai 1968; Amein & Fang 1970; Fread 1973; Amein & Chu 1975; Liggett & Cunge 1975; Delis *et al.* 2000; Krishnappan & Altinakar 2005; Retsinis *et al.* 2018).

The Muskingum-Cunge method can be classified as a hydraulic method based on the kinematic wave model with added diffusion terms, representing the real behavior of the flood wave propagation, i.e., attenuation of the peak flowrate, temporal translation of the peak as well as diffusion of the flood wave (Akan 2006). The rooting behavior of the Muskingum method is based on the kinematic wave but the attenuation of flood wave is due to the numerical diffusion. The analysis showed that the applied approximation modifies the kinematic wave model to the form of advection-diffusion equation that is not consistent (Szymkiewicz 2010). The resulting equation for estimating the hydrograph downstream is the same as the one produced by the Muskingum method (Akan 2006).

Equations (1) and (2) model the unsteady flow conditions in prismatic open channels where either added rainfall or lateral inflow can be neglected. The aim of this work is the numerical solution of the complete dynamic wave model under subcritical flow conditions in a long prismatic open channel of symmetric trapezoid and triangular cross sections. Three explicit numerical schemes will be used, namely, the Lax-Diffusive, the MacCormack, and the Lambda as well as the implicit four-point Preissmann scheme. The numerical results will be compared to those obtained by the Muskingum-Cunge method and the HEC-RAS (Warner *et al.* 2010) river modeling software.

## NUMERICAL METHODS

The governing equations may be written either in conservative or non-conservative form. The conservative form should be preferred since it better conserves the transported flow variables and simulates more accurately the flood wave celerity (Chaudhry 2008).

*A*

*=*

*A*(

*x*,

*t*) is the wetted area, and is the depth of the center of gravity of the wetted cross section measured from free surface. The St Venant equations are of hyperbolic type, having two real characteristic directions where the positive characteristic direction is given by

*λ*=

^{+}*V*

*+*

*c*, and the negative one by

*λ*=

^{−}*V−c*, where

*c*

*=*

*c*(

*x*,

*t*) is the celerity of the propagating wave. In hyperbolic problems, the number of characteristic curves inserted in the computational domain determines the required number of initial and boundary conditions for a mathematically well posed problem (Chaudhry 2008). In the present work, where subcritical flow is examined, the characteristic directions are

*λ*> 0 and

^{+}*λ*< 0, so that at the horizontal line crossing the initial time (

^{−}*t*

*=*0), we need two initial conditions for the flow depth

*y*and the cross-sectional average velocity

*V*, and at the vertical lines passing through each end of the channel (

*x*

*=*0,

*x*

*=*

*L*,

*L*is the total length of the channel), we only need one boundary condition since only one characteristic curve enters the computational domain. The above system of partial differential equations with the appropriate auxiliary conditions has to be solved numerically using the method of finite differences.

A schematic representation of the channel's geometry with both the trapezoid and triangular cross section is shown in Figure 1. The first step is the discretization of the physical domain. To do this, we divided the total length of the channel *L*, by a number of *n _{x}* nodes (including the end nodes) with index

*i*denoting the space location, and the total simulation time

*T*by

*n*nodes with index

_{t}*j*referring to the time level, all equally spaced for simplicity, resulting in a uniform grid with Δ

*x*and Δ

*t*being the space and temporal intervals, respectively, as shown in Figure 2. Next, the approximation of the partial derivatives and the variables met in Equation (3) was done with suitable finite difference expressions showing consistency, stability, and convergence for better overall accuracy of the method. Finally, we either computed directly all the unknown variables, or in the implicit schemes we solved the resulting algebraic system of finite difference equations with direct or iterative techniques, thus evaluating the flow variables. The different numerical schemes implemented for this work are presented next.

### Lax-Diffusive scheme

*et al.*2018): where for brevity we used

*f*for both variables

*y*and

*V*. Substituting the above into Equation (3) it yields:

From Equation (4) one can compute the values of *A*, *Q**=**VA*, at the *j* + 1 time level from the flow variables at time *j* level known either from the initial conditions or the previous time step calculations.

### MacCormack scheme

This is a two-step predictor-corrector scheme that is second-order accurate in space and time. In this algorithm, backward finite differences are used for the spatial derivatives during the predictor part, and forward finite differences for the corrector part. This sequence of differentiation may be reversed in modifications of the method. Thus, we discretize Equation (3) as follows (Chaudhry 2008).

*Predictor part:*where the superscript ‘*’ refers to the flow variables estimated at the predictor part. Substitution of the above expressions into Equation (3) and after some simplification it yields the following equation of finite differences:

The computed value of gives the variables and from which we determine the values of , . These values are then used in the corrector part to compute and .

*j*+ 1 time level is given by:

### Lambda scheme

*λ*-form. After some analysis we get (Chaudhry 2008): where superscripts ‘ + ’ and ‘ − ’ define the evaluation of the partial derivatives along the positive and negative characteristics, respectively. We discretize Equations (8) and (9) according to the sign of the characteristic direction using backward finite differences in the positive direction and forward finite differences in the negative one. The scheme is based on predictor and corrector parts. In the predictor part, the spatial derivative along the positive characteristic direction is discretized using three nodal values whereas in the corrector part, the spatial derivative along the negative characteristic direction is discretized using again three grid points with second-order accuracy in space and first-order in time. After all it yields the following (Chaudhry 2008).

*j*+ 1 time level are In the case of explicit algorithms, suitable one-sided finite differences have to be used for discretization of the equations at the boundary nodes according to the signal direction. Due to the explicit nature of the above schemes, the time interval is limited according to the Courant number, which must be less than or equal to unity at all grid points and time steps, indicating that the actual flood wave velocity |

*V*±

*c*| must be smaller than the computational velocity

*dx/dt*(Chaudhry 2008).

### Preissmann scheme

*j*and the unknown time level

*j*+ 1, weighted by factor

*α*, with 0 ≤

*α*≤ 1 as follows, which is second-order accurate in space and time and unconditionally stable if

*α*≥ 0.5 (Akan 2006):

Equation (16) has four unknowns, the , , , , where is the free surface elevation measured from a horizontal datum at grid location *i* and time level *j* + 1. Writing the two equations for each reach, we end up with 2·(*n _{x}* − 1) equations. However, we have 2·

*n*unknowns. The remaining two equations are provided from the boundary conditions, and then the system of nonlinear algebraic equations is solved simultaneously, providing the solution of unknown flow variables at

_{x}*j*+ 1 time stage.

In the present work, the iterative Newton-Raphson method was selected for the solution of the system (Akan 2006; Retsinis *et al.* 2018). We guess the values of the unknown flow variables , , , at each grid point and calculate the Jacobian matrix which has as elements the partial derivatives of Equation (16) and the boundary conditions with respect to , , , . The iterations stop once the corrections and for all grid points are within the limits we have set.

All the numerical codes discussed earlier have been implemented in house, using Matlab^{®} computational environment.

### Muskingum-Cunge method

*Q*

*=*

*Q*(

*x*,

*t*), given by the following equation (Szymkiewicz 2010): where

*C*

*=*

*C*(

*Q*) is the kinematic wave celerity and

*D*

*=*

*D*(

*Q*) is the coefficient of hydraulic diffusivity. Linearizing Equation (17) and discretizing it with finite differences, the final relationship for routing of the upstream hydrograph is given by the following equation (same as the Muskingum method (Akan 2006)): where

*I*is the inflow and

*Q*the outflow, with the subscripts 1 and 2 corresponding to the time levels n and

*n*+ 1, respectively, , , ,

*Q*is a reference discharge,

_{o}*S*is the longitudinal bottom slope of the channel,

_{o}*L*

_{1}is the length of the channel reach, while all variables with subscript ‘o’ correspond to those calculated from the reference discharge and

*m*

*=*5/3 (Akan 2006). The peak of the inflow hydrograph was selected as the reference discharge showing small sensitivity in the final results. To improve accuracy it is suggested that Δ

*t*must be smaller than one-fifth of the time interval from the beginning until the peak of the inflow hydrograph, and the length of the channel section for Muskingum-Cunge computations must be limited according to the following relationship (Akan 2006):

## RESULTS

*L*

*=*8,000 m long, with symmetric trapezoid cross section having a bottom width of 1 m, or triangular cross section, with side slopes 2H:1 V in both cases. The longitudinal bottom slope was

*S*= 0.001, small enough to ensure subcritical flow conditions, and Manning's roughness coefficient was 0.025, a typical value for open channels constructed by concrete. The initial condition at the beginning of time (

_{o}*t*

*=*0) was considered steady-uniform flow with discharge 10 m

^{3}/s. The upstream boundary condition was a known hydrograph with a mathematical expression (Akbari & Firoozi 2010): where

*Q*(

*x*

*=*0,

*t*) is the variable inflow with time

*t*at the upstream end of the channel,

*Q*

_{min}= 10 m

^{3}/s is the base flow,

*Q*

_{max}= 100 m

^{3}/s resulting in a maximum flow rate of 110 m

^{3}/s,

*t*

_{max}= 7,500 s (2.08 hr) is the time when peak flow occurs, and

*T*

*=*30,000 s (8.33 hr) is the total simulation time. The open channel was considered to be long enough to reach uniform depth so the friction slope equals the longitudinal bottom slope,

*S*=

_{f}*S*in any time stage at the downstream end of the channel. For all of the numerical methods used, the spatial and temporal discretization included

_{o}*n*= 9 and

_{x}*n*= 60,000 nodes, respectively, resulting in a uniform computational grid, also satisfying the Courant criterion in the case of explicit schemes.

_{t}In the Muskingum-Cunge method, the selected reference discharge was 100 m^{3}/s yielding a maximum reach length of 1,081 m calculated by Equation (19) for the trapezoid cross section, while *L*_{max} = 1,082 m for the triangular cross section. The spatial and temporal discretization was chosen as *L*_{1} = 1,000 m and Δ*t**=* 200 s, respectively.

For the trapezoid cross section, the temporal evolution of flow rate, flow depth, and mean velocity is presented for all methods at downstream locations *x**=* 0, *x**=**L*/4, *x**=**L*/2, *x**=* 3*L*/4, and *x**=**L*, measured from the upstream end of the channel, in Figures 3–16. Also, Figures 17–20 show the spatial variation of flow depth with distance, measured from the upstream end of the channel for the Lax-Diffusive, MacCormack, Lambda, and Preissmann methods, respectively, while Figure 21 displays the variation of the maximum flow depth for all times computed by all methods except for the Muskingum-Cunge one, along the length of the channel.

Regarding the triangular cross section, similar results for the flow rate, depth, and average cross-sectional velocity are presented at the downstream locations *x**=**L*/2, *x**=* 3*L*/4, *x**=**L*, measured from the upstream end of the channel in Figures 22–30. Figure 31 depicts the maximum value of flow depth for all times computed by all methods except for the Muskingum-Cunge one along the length of the channel.

For the trapezoid cross section, the flow rates computed from the MacCormack, Lambda, Preissmann, and HEC-RAS methods are almost identical, while Lax-Diffusive scheme underestimates the peak discharge. The Muskingum-Cunge method overestimates the peak discharge at all downstream locations. Both Lax-Diffusive and Muskingum-Cunge showed steeper increase in the discharge compared to the other methods. All methods except for the Lax-Diffusive one compute a small reduction of the peak discharge at the downstream end of the channel. The Lax-Diffusive computed peak is 103.2 m^{3}/s due to friction losses, while the arrival time of the peak discharge at the downstream end for the Lax-Diffusive, MacCormack, Lambda, Preissmann, Muskingum-Cunge, and HEC-RAS method was 3.84 hr, 3.37 hr, 3.39 hr, 3.40 hr, 3.11 hr, and 3.39 hr, respectively. It is obvious that the Lax-Diffusive scheme predicts the latest arrival time of maximum flow rate, while the Muskingun-Cunge method forecasts the earliest time. The results from MacCormack, Lambda, and Preissmann schemes as well as HEC-RAS compare well regarding the flow depth, while the results from Lax-Diffusive deviate slightly. Similar conclusions can be drawn for the average cross-sectional velocity where the results from Lax-Diffusive are different from those of all other methods, while this difference is reduced further downstream. The mass conservation equation using numerical integration with trapezoid rule for the computed hydrographs downstream was confirmed. More specifically, the maximum percentage error relative to the water mass entering the upstream cross-section for Lax-Diffusive, MacCormack, Lambda, and Preissmann schemes was 2.59, 1.41, 1.64, and 1.58, respectively.

Regarding the triangular cross section, similar conclusions are drawn for the flow rate, flow depth, and the average cross-sectional velocity. The arrival times of the peak discharge at the downstream end for the Lax-Diffusive, MacCormack, Lambda, Preissmann, Muskingum-Cunge, and HEC-RAS methods are 3.85 hr, 3.37 hr, 3.38 hr, 3.39 hr, 3.11 hr, and 3.40 hr, respectively. Again, the water mass conservation was confirmed while the maximum percentage error in the case of Lax-Diffusive, MacCormack, Lambda, and Preissmann scheme was 2.61, 1.42, 1.65, and 1.44, respectively, if compared to the water mass inserted in the upstream end of the channel.

## CONCLUSIONS

The flood wave propagation in long prismatic open channels of symmetric trapezoid and triangular cross sections under subcritical flow conditions has been studied. The results from four numerical methods, three explicit ones, namely, the Lax-Diffusive, MacCormack, and Lambda methods, as well as the implicit Preissmann scheme, are compared to those from the hydraulic Muskingum-Cunge method and the HEC-RAS commercial software program. All codes were implemented in Matlab environment. Due to the explicit character of Lax-Diffusive, MacCormack, and Lambda schemes, the maximum time interval was restricted according to the Courant criterion, whereas in contrast, the Preissmann scheme was unconditionally stable. The weighting factor selected for accuracy reasons was 0.5, the convergence criteria regarded the unknown discharge allowed error of 0.001 m^{3}/s, and the linear system of algebraic equations was solved by a direct method. The structure of the implicit code was improved resulting in a low number of iterations per time level, which led to accelerated computational times. The numerical results have been compared to those obtained from the commercial HEC-RAS software and the hydraulic/hydrologic Muskingum-Cunge method.

The largest error in mass conservation equation that appeared in the Lax-Diffusive method is probably due to the fact that it is a one-step algorithm, if compared to the two-step schemes. All the results from the numerical methods, the Muskingum-Cunge method and HEC-RAS software compare well.

The major innovation for this work was the comparison of different explicit or implicit numerical schemes with the Muskingum-Cunge hydraulic/hydrologic method and the HEC-RAS software. Furthermore, great effort was made for the coding of all numerical schemes in suitable programming language, revealing practical aspects and difficulties of code implementation useful for future works. Finally it proved that the explicit schemes produced satisfactory results if compared with the results of the more demanding implicit schemes.

To conclude, the numerical results compare well to those from use of the Muskingum-Cunge method and one-dimensional unsteady flow calculations software HEC-RAS, requiring more effort than the simpler Muskingum-Cunge method which can be used for the preliminary design of long prismatic open channels under unsteady flow conditions.