Third Order Reconstruction of the KP Scheme for Model of River Tinnelva

The Saint-Venant equation/Shallow Water Equation is used to simulate flow of river, flow of liquid in an open channel, tsunami etc. The Kurganov-Petrova (KP) scheme which was developed based on the local speed of discontinuity propagation, can be used to solve hyperbolic type partial differential equations (PDEs), hence can be used to solve the Saint-Venant equation. The KP scheme is semi discrete: PDEs are discretized in the spatial domain, resulting in a set of Ordinary Differential Equations (ODEs). In this study, the common 2 order KP scheme is extended into 3 order scheme while following the Weighted Essentially Non-Oscillatory (WENO) and Central WENO (CWENO) reconstruction steps. Both the 2 order and 3 order schemes have been used in simulation in order to check the suitability of the KP schemes to solve hyperbolic type PDEs. The simulation results indicated that the 3 order KP scheme shows some better stability compared to the 2 order scheme. Computational time for the 3 order KP scheme for variable step-length ode solvers in MATLAB is less compared to the computational time of the 2 order KP scheme. In addition, it was confirmed that the order of the time integrators essentially should be lower compared to the order of the spatial discretization. However, for computation of abrupt step changes, the 2 order KP scheme shows a more accurate solution.


Introduction
The Kurganov-Petrova (KP) scheme is a 2 nd order scheme, developed for solving hyperbolic Partial Differential Equations (PDEs) (Kurganov and Petrova, 2007).The Saint-Venant equation or commonly known as the Shallow Water Equation is a hyperbolic type PDE which is suitable to model flow of water in open channels, rivers, etc. (Fayssal et al., 2015).As a case study, flow of water in a reach of river Tinnelva between hydropower stations Årlifoss and downstream station Grønvollfoss, will be simulated using the Saint-Venant equation.Water height at Grønvollfoss dam due to changes in the volumetric flow out of the Årlifoss station will be computed in this study.
In this study, the 2 nd order KP scheme is being extended into a 3 rd order scheme by following the Weighted Essentially Non-Oscillatory (WENO) and Central WENO (CWENO) reconstruction procedures.Subsequently, both schemes will be used to discretize the Saint-Venant equation spatially.Different time integrators will be used to solve the resulting Ordinary Differential Equations (ODEs) from the spatial discretization: variable step-length solvers in MATLAB such as ode15s, ode23s, ode23t and ode45, as well as fixed step-length solvers such as the Euler method, the 2 nd order Runge Kutta (RK), and the 4 th order RK methods.The solution of this set of ODEs gives the final water height variation at Grønvollfoss dam.This paper is arranged as follows.A basic introduction to the Saint-Venant equation and the 2 nd order KP scheme will be given in Sections 2. In Section 3, the 3 rd order reconstruction steps will be discussed together with WENO and CWENO constructions.Flow of water in river Tinnelva is simulated in Section 4 for both the 2 nd order and 3 rd order KP schemes.Section 5 holds a discussion of results: simulated results are compared for both 2 nd and 3 rd order schemes in order to check improvements of the accuracy due to the increment of the order of the KP scheme.Numerical stability of the schemes is also discussed in this section.Finally, some conclusions are drawn in Section 6.
2 Governing equations with 2 nd order KP scheme The governing Saint-Venant equation is discussed in this section, together with the basis of 2 nd order KP scheme.

Mathematical model
The Saint-Venant equation/ Shallow Water Equation has been used to simulate flow of water in open channels, rivers, etc. for decades (Fayssal et al., 2015).Basic one-dimensional conservation laws such as momentum and mass balance provided the original basis for the Saint-Venant equation; subsequently it has been derived by integration of the full three-dimensional flow equations (Fayssal et al., 2015).
The Saint-Venant equation with source term can be posed as follows (Sharma, 2015) ∂U ∂t + ∂F (x, t, U ) ∂x = S(x, t, U ), where: . (4) Here, z is the water level above a datum, B is the bottom elevation from the datum, q is volumetric flow rate per unit width, w is the width of the river, n is Manning's roughness coefficient, and g is acceleration due to gravity.The S terms reflect source terms: including expressions of friction which give resistance against flow.In the 2 nd order KP scheme, properties of z and q are computed using the midpoint rule (Kurganov and Petrova, 2007).Hence, for the source term discretization, average values of the considered properties will be taken at each control volume (CV).

2.2
The 2 nd order KP scheme In the Finite Volume Method (FVM), the properties under consideration are averaged for each CVs (Versteeg and Malalasekera, 2007).During transients (shock, hydraulic jump, wave flow etc.), arbitrarily selected consecutive CVs have different average values of the fluxes.Hence low order reconstructions such as linear reconstruction, produce two different property values for a single interface of a CV, subsequently causing oscillations in the solution domain (Kurganov and Petrova, 2007).Such phenomenon is known as discontinuity or as Riemann problem (Kurganov and Tadmor, 1999).Due to the discontinuity, most of the numerical methods either failed to produce an accurate result or failed to achieve fast convergence.Hence, numerous efforts were done to develop a proper numerical scheme which could achieve proper/expected accuracy together with fast convergence.
In the development of the KP scheme, the precise local speed of discontinuity propagation is considered along with small CVs of variable size (Titarev and Toro, 2004).Hence the KP scheme is able to meet the demand of a smooth solution.In the numerical scheme development history, a staggered CV concept has been introduced to handle discontinuity at CV interfaces (Titarev and Toro, 2004).During a transient, the local velocities are usually different on each side of a CV interface.Therefore, altered staggering at both sides of the CV is reasonable.In the KP scheme, actual CV interfaces have been shifted to positive and negative direction, at a small time ( t) by considering the local velocity of discontinuity propagation.Then a new set of CVs have been created, here we refer them as "virtual CVs".Then, for each non-uniform CVs, a piecewise linear reconstruction is performed.Later on, the linearly reconstructed values are projected on the original uniform CVs while assuming the limit ( t → ∞) (Titarev and Toro, 2004).
In the KP scheme, properties are indexed by plus (+) and minus (-) with reference to the property flux.The local speed of discontinuity propagation has been calculated by considering the Jacobian matrix of the governing equations (Kurganov and Petrova, 2007).In order to achieve higher resolution and a well-balanced scheme, the Total Variation Diminishing (TVD) concept together with the flux limiter concept is used in the development of the 2 nd order KP scheme.The standard "minmod" limiter was used in the original development of the KP scheme.However, many alternative flux limiters can also be used (Titarev and Toro, 2004).The KP scheme does not need Riemann solvers: remedies introduced to manipulate the Riemann problem by German mathematician Bernhard Riemann, which consist of number of intermediate calculations to select proper value for a single interface.Such intermediate calculations consume considerably large amount of computational time.Elimination of the Riemann solvers in the KP scheme leads to less computational time (Kurganov and Tadmor, 1999).Numerical viscosity with the KP scheme is lower compared to the Nessyahu Tadmor (NT) scheme (Titarev and Toro, 2004).The 2 nd order KP scheme discretizes the Saint-Venant equation spatially, yielding a collection of ODEs in time.These ODEs (without the source term) can be written as follows (Sharma, 2015).
where ūj is the cell center average values, H j± 1 2 (t) are the central upwind numerical fluxes at the cell interfaces, defined as: are the one-sided local speeds of propagation and u ± j± 1 2 are property fluxes at indexed positions.

rd order reconstruction
Polynomial reconstruction plays a major role in the development of the KP scheme.In the 2 nd order KP scheme, constructed polynomials are 2 nd order accurate (Kurganov and Petrova, 2007).By following CWENO reconstruction: quadratic reconstruction, a 3 rd order interpolant could be obtained.This interpolant is a combination of two linear functions at the left and right boundaries of a CV and a parabolic interpolant at the central coordinates of the CV (Kurganov and Levy, 2000).If geometrical arrangement of a nodal group: here referred as stencils are smooth stencils, such interpolant show a 3 rd order accuracy (Kurganov and Levy, 2000).The WENO scheme on the other hand uses linear reconstructions.

WENO scheme
The WENO and the Essentially Non-Oscillatory (ENO) techniques were developed for constructing a higher order numerical scheme (Liu et al., 1994).The ENO scheme was developed by Harten et al. (1987).The WENO scheme was developed by Liu et al. (1994).Both schemes use adaptive stencils, thus automatically achieving a higher order of accuracy for the solution near the discontinuity (Zhang and Shu, 2009).The WENO scheme is based on decomposition of the local characteristic with flux splitting to avoid oscillations (Titarev and Toro, 2004).Both schemes approximate the solution by using nonlinear adaptive procedures.Hence, they avoid discontinuities in the interpolation procedure.These schemes are very common in shock and complicated smooth solution structures (Titarev and Toro, 2004;Zhang and Shu, 2009).
At transition conditions, the WENO scheme has a higher order of accuracy than the ENO in the same set of smooth coordinates (Kurganov and Levy, 2000).In the WENO scheme, instead of using the single candidate stencil, a linear combination of all candidate stencils are used (Shi et al., 2002).Calculation of the weights to each candidate stencil is a crucial factor for the success of the WENO reconstruction.Some advantages of the WENO scheme compared to ENO reconstruction can be summarized as follows: (Titarev and Toro, 2004;Kurganov and Levy, 2000;Zhang and Shu, 2009;Shi et al., 2002).
• Compared to ENO reconstruction, the WENO reconstruction has higher order of accuracy for the same set of stencils (Zhang and Shu, 2009).
• WENO reconstruction has 3 rd order of accuracy for piecewise linear reconstruction and 5 th order for piecewise quadratic reconstruction, while the ENO construction gives 2 nd order and 3 rd order accuracy for similar reconstructions (Kurganov and Levy, 2000).
Steps in the development of a 3 rd order WENO reconstruction, with a worked example, will be discussed in the sequel.

Steps in WENO reconstruction
A one-dimensional (1D) hyperbolic conservation law for a general system can be written as in Equation 8 (Titarev and Toro, 2004;Zhang and Shu, 2009;Shi et al., 2002;Shu, 2003).
Here U (x, t) is vector of known conservative variables and F (U ) is the physical flux vector.In the Finite Volume Method (FVM), average physical flux, F (U ) for the CV are calculated.Physical flux can be written as an integral form.For arbitrarily chosen CV, single property u(t) of the flux can be written as Based on Equation 9, property average for a CV arbitrarily located at spatial domain; j, can be written as: difference of fluxes at CV interfaces divided by its' spatial distance, x j [6].
Equation 10 can be evaluated by a suitable polynomial which achieves a desired degree of accuracy, hence the values of the polynomial at x j+ 1 2 and x j− 1 2 can be obtained.Values of reconstructed polynomial essentially should be matched with numerical values of x j+ 1 2 and x j− 1 2 .The polynomial reconstruction in general form for a selected position can be written as In the upwind discretization 1 scheme, the direction of the flux has to be considered (Versteeg and Malalasekera, 2007).According to the flow direction, single interfaces compose two values indexed as positive (+) and negative (-).Those values at x j+ 1 2 can be written as Equation 12.
Here f (u − , u + ) is referred to as a numerical flux.In Computational Fluid Dynamics (CFD), values of u − j+ 1 2 and u + j+ 1 2 can be obtained through piecewise reconstructions.

Steps of piecewise polynomial reconstruction
The property u(t) for a CV can be written as Equation 9. (Shu, 2003).In WENO scheme, three consecutive CVs will be considered for the evaluation of a value at a given position.Coordinates of the spatial domain can be illustrated as in Figure 1.
1 Difference between upwind and central schemes is defined by considering piecewise reconstruction.In an upwind scheme, reconstructed values are based on the middle point of the CVs while a central scheme is based on staggered average at the CV interfaces (Shu, 2003).In addition, the direction of the property flux is being considered (Versteeg and Malalasekera, 2007).Property u j+ 1 2 will be a combination of ūj−1 , ūj and ūj+1 in WENO scheme.For the considered stencils, polynomial p(x) can be obtained by integrating Equations 13 to 15 (Kurganov and Levy, 2000).
To find the desired point value, a combination of all three properties is weighted (Kurganov and Levy, 2000;Liu et al., 1994;Liu and Tadmor, 1998).Thus, (16) This averaging has 3 rd order of accuracy if the flux is smooth for the selected stencils.Using such approximation with fixed stencils leads to a higher order linear scheme, which will be oscillatoric in the presence of a shock condition (Godunov theorem) [13].Therefore, the method has to be normalized for a general system.Hence, five consecutive stencils have been chosen which produce three different combinations: here the x j−2 , x j−1 and x j combination is named 0, x j−1 , x j and x j+1 is named 1 and the combination of x j , x j+1 and x j+2 is named 2.Such stencil combination as in the WENO reconstruction can be illustrated as in Figure 2.These combinations are convex combinations (Arshed and Hoffmann, 2013).
Based on Figure 2, the computed approximation for the different sub stencils can be written as (Kurganov and Levy, 2000;Liu et al., 1994;Shi et al., 2002;Liu and Tadmor, 1998) If function u(x) is smooth in the three sub stencils, then three approximations u 0 in 3 rd order of accuracy.

Weighing technique of WENO scheme
In WENO scheme, computation of point value of u j+ 1 2 requires three intermediate equations: Equation 17 to 19.These three equations are written based on five consecutive CVs.Subsequently, value of each equation has to be weighed in order to get the final approximate solution to the desired point, u j+ 1 2 .A general form of the averaging can be expressed as The coefficients γ 0 , γ 1 and γ 2 are the linear weights.
In WENO scheme, a robust choice of such weights will be (Kurganov and Levy, 2000) If flux (F (U )) is smooth for all three stencils, then the reconstruction achieves 5 th order of accuracy.Then the choice of linear weight has to satisfy the condition γ k > 0 as long as sum of all weights equal to 1. 5 th order of accuracy for the smooth region is certain for the WENO scheme.However, at critical regions where wave/shock occur, accuracy of the developed WENO scheme is reduced to 3 rd order and in special case order of accuracy reduce even further, to at least 2 nd order (Arshed and Hoffmann, 2013).This possible accuracy drop can be restored by following nonlinear weighing system (Arshed and Hoffmann, 2013).
3.5 Nonlinear weights w 0 , w 1 and w 2 Discontinuity usually causes numerical dissipation.Two ways of minimizing the numerical dissipation can be performed with the WENO scheme: the upwind optimal stencil and nonlinear adaptation mechanism (Arshed and Hoffmann, 2013).Optimization techniques which manipulate optimal linear weights delays the dissipation in the upwind stencils.However, such optimizations techniques reduce the 5 th order of accuracy into 3 rd order of accuracy in the WENO scheme.
Henceforth a new strategy has to be used to maintain at least 3 rd order or higher order accuracy in the WENO scheme, while minimizing dissipation.Matching convex combinations of all possible candidate stencils by nonlinear weights could achieve higher order accuracy than 2 nd order: 2 nd order of accuracy being the worst case scenario (Arshed and Hoffmann, 2013).
With nonlinear weights Equation 20 can be written as (Kurganov and Levy, 2000) u These nonlinear weights are adapted for the stencils with discontinuity.These weighing techniques achieve optimal linear weights: where the nonlinear weights of w 0 , w 1 and w 2 are closer to the linear weights of γ 0 , γ 1 and γ 2 at smooth stencil coordinates (Arshed and Hoffmann, 2013).
Here O x 2 signifies that the remaining part is roughly proportional 2 to ( x) 2 when ( x) 2 → 0, then Equation ( 23) becomes, If property u(x) has a discontinuity, then the corresponding weight w k is considered to have a very small value.
The following expressions for the nonlinear weights can be shown to give a robust algorithm, (Kurganov and Levy, 2000) Here, = 10 −6 .This adaptive nature of the nonlinear weights is originally attributed to the local smoothness indicators (Arshed and Hoffmann, 2013).The three smoothness indicators, β in Equation 27 can be written as (Kurganov and Levy, 2000) 2 In general, if z(x) is a quantity that depends on the value of x, then as x → 0, the equation z(x) = CO (x p ) becomes |z(x)| ≤ Cx p (Sharma, 2015)

Steps in CWENO reconstruction
The WENO scheme will achieve 3 rd order in accuracy even for a linear reconstruction.In the 2 nd order KP scheme, the polynomial which represents CV average values is a linear function; the polynomial reconstruction is therefore a linear reconstruction.However, by following CWENO reconstruction, 3 rd order of accuracy can be achieved by merely changing the polynomial reconstruction into a quadratic polynomial reconstruction (Levy et al., 2000).In such quadratic polynomial reconstruction, interpolants are combinations of two linear functions at each plus and minus interfaces of the CV and one central convex combination (Levy et al., 2000).Simpson's rule is used in defining a 3 rd order polynomial in the CWENO reconstruction.
In CWENO reconstruction, the flux limiter has no influence.Hence, the flux limiters will be removed from the reconstruction step.The 2 nd order KP scheme is a central upwind scheme; however, the CWENO scheme is based on a central difference (CD) scheme (Kurganov and Levy, 2000).The CD scheme for a hyperbolic PDE can be written as follows, Equation 32 implies that the conserved variables change both in space and with time.Integrating in time from t n to t n+1 , the central scheme can be written as ūn+1 The value of u n j can be written as a polynomial which gives the desired order of accuracy, A general form of a quadratic polynomial which assures third order accuracy can be written as Equation 35 (Liu and Tadmor, 1998).Such polynomial function is nothing other than the Taylor series expansion of the given function.
Based on this quadratic polynomial, the 3 rd order reconstruction of the space coordinates can be written as Equation 36, ūn Subsequently, exact numerical fluxes in the time domain can be found by using Simpson's quadrature rule (Liu and Tadmor, 1998).Simpson's rule to solve an integral for a general case can be written as (Ujevic, 2007) Based on the general formula of Simpson's rule, time integration can be written as This method requires point values for intermediate time steps.Such point values can be approximated by Taylor series expansion or the RK method (Levy et al., 2000).The steps derived by the RK method can be re-written as Equation 39to 41 (Kurganov and Levy, 2000), Then based on all reconstruction ūn+1 For a smooth coordinate which does not have large gradients, polynomial which gives 3 rd order accuracy can be written as Equation 43.
For the non-oscillatory quadratic piecewise parabolic reconstruction, optimum polynomial becomes P opt,j (Levy et al., 2000).Then P opt,j can be represented as Equation 44.Optimum polynomial P opt,j is a parabola that interpolates the point values of ūn j−1 , ūn j and ūn j+1 .
Here l denotes CVs located at spatial domain j − 1 , j and j + 1.Based on quadratic polynomial reconstruction, P opt,j can be written as Here u n j is the reconstructed point value, u j is the point value of reconstructed gradient or slope, and u j is the point value of the reconstructed second order derivative (Liu and Tadmor, 1998).u n j , u j and u j can be written as, This reconstruction may cause oscillation at sharp gradients where discontinuity occurs (Levy et al., 2000).Therefore, WENO reconstruction is used to weigh each polynomial expression.In general, each polynomial with required weight can be written as Equation 49, Here, weights w j i are to be selected in order to assure the following conditions, i w j i = 1; w j i ≤ 0 (50) i refers to left L, center C and right R positions of a CV i.e. i = L, C, R. Then P L and P R become linear functions at the left and right boundaries, and P C becomes a quadratic parabolic expression.For the requirement of the conservation of the property, P R interpolates the CV adjacent to the right side of the stencils.Hence P R can be written as (Kurganov and Levy, 2000).
(51) Piecewise linear interpolation for the right side can be written as In the same way, the left side can be expressed as For the polynomials P R (x) and P L (x), weights can be taken as w R and w L respectively.Then the weight of the central polynomial P c (x) can be calculated using Equation 50.Then the whole expression with their weights can be written as (Kurganov and Levy, 2000;Levy et al., 2000) ) For 3 rd order of accuracy, weights w L and w R are chosen to be 1 4 .Hence, w C , which is the center point weight, becomes 0.5.Considering each weight, a parabolic expression can be written as More precisely P C (x) can be written as (Kurganov and Levy, 2000;Levy et al., 2000)

CWENO weight calculation procedure
For better stability and convergence nonlinear weighing techniques can be applied here.Henceforth actual weights can be calculated by the following steps.Then the weights are selected as (Kurganov and Levy, 2000;Levy et al., 2000).
Here α i can be written as, where = 10 −6 .The value of p will be chosen to provide the highest accuracy in the smooth areas and to assure non-oscillatory behavior at the discontinuities.
A typical empirical value for p is 2. Local smoothness indicators referred as β in Equation 58 can be written as a general form expressed in Equation 59.
Thus we have, 3.8 3 rd order KP scheme In the KP scheme, local speed of wave propagation can be written as follows based on the Jacobian of the flux function.
(63) Here λ refers to the eigenvalue of the flux function.Then local speed of wave propagation can be calculated as A quadratic polynomial can be written in general as, According to Kurganov and Levy (Kurganov and Levy, 2000), integrals for the CV over and Values of wn+1 j+ 1 2 and wn+1 j are average of property at virtual CVs at t = t n+1 .A third order polynomial of piecewise non oscillatory interpolant of CWENO reconstruction can be written as (Kurganov and Levy, 2000) wn+1 Then ūn+1 j can be written as Equation 71 (Kurganov and Levy, 2000).Here ūn+1 j represents the actual CV of the spatial domain.
For 1D problems, the semidiscrete approximation can be written by taking t → 0 as As t → 0 then the virtual CVs that were deliberately created to capture discontinuity have zero width.
When t → 0, the local speed a j+ 1 2 (t) (Kurganov and Levy, 2000) can be written as Equation 76, In Equation 66, A j , B j and C j can be written as, The final reconstructed equations can be written as 3.9 Source term discretization using Simpson's method The source term of the Saint-Venant equation can be written as (Sharma, 2015), According to Simpson's rule, properties have to be evaluated at x j− 1 2 , x j+ 1 2 and x j .Properties at x j denoted as S C are similar to Equation 84.Properties q, z and B are representing the average values of the CV.Then,

Simulation Study
A reach of river Tinnelva in Telemark, Norway is considered in the study.The focus is on the propagation of water wave, and the water level at the Grønvollfoss dam due to changes in the outflow of water from a power station located 5 km upstream.In order to simplify the complex river geometry, the width of the river  for the considered reach is assumed to be constant.Table 1 summarizes relevant quantities for the simulation study.
The bottom topography of the river is assumed to have three consecutive slopes.The assumed bottom slope is as in Figure 3.The Saint-Venant equation has been discretized by using both the 2 nd order and 3 rd order KP scheme.Then the resulting set of ODEs will be discretized in time by variable step-length solvers in MATLAB :ode23t, ode23s, ode45, ode15s and fixed step-length solvers :The Euler method, RK2 and RK4.

Results and Discussion
Simulated results for the 2 nd order and 3 rd order KP scheme and numerical stability of each solver with variable step-length and fixed step-length will be discussed.For easiness, simulated results and numerical stability will be discussed under two different headings.Hydraulic jump for both scheme will also be discussed under separate heading.Height (m) The Euler method-for the 2nd order and 3rd order KP scheme 2nd order Euler 3rd order Euler Figure 4: The Euler method for 2 nd and 3 rd order KP scheme ( t = 0.25s).

Simulation results
Solution with the Euler method for the 2 nd order and 3 rd order KP scheme are shown in Figure 4.According to Figure 4, the solutions seem to closely match each other.However, the 2 nd order KP scheme shows a small overshoot at the first wave peak.The 3 rd order scheme shows a smoother solution compared to the 2 nd order scheme.It was observed that in the 2 nd order KP scheme when the step length t increases up to 0.7s, the fixed step-length Euler method exhibits oscillatory nature (Dissanayake et al., 2016).However, for the 3 rd order scheme, a non-oscillatory solution is achieved until t is increased to 1.4s i.e. for t ≤ 1.4s.However increment of t should satisfy the CFL condition in order to achieve convergence LeVeque (1999); Griffiths et al. (2015).CFL condition can be written as, Here, C is a dimensionless number, u is the magnitude of the velocity and x is the length of CV.For Height (m) The 4th order Runge Kutta for the 2nd order and 3rd order KP scheme 2nd order RK4 3rd order RK4 Figure 6: The RK4 method for 2 nd and 3 rd order KP scheme ( t = 0.25s).
the upwind scheme, C max = 1 (Griffiths et al., 2015).This observation indicates that the choice of t together with CFL (Courant-Friedrichs-Lewy) condition has impact on solution stability.Results for fixed steplength solvers RK2 and RK4 are shown in Figure 5 and Figure 6.In both figures, the 2 nd order KP scheme shows a small overshoot at the first peak.However, both schemes produce more or less similar solutions.
The results for both the 2 nd order and the 4 th order schemes using the RK4 solver are plotted in Figure 6.The results obtained from all the fixed step-length solvers have been plotted in Figure 7 to observe how much they deviate from each other.A close inspection shows negligible deviations besides the overshoot occurring at the first peak of the solution.An exploded view of the overshoot for all fixed step-length solvers Figure 7: All fixed step-length solvers for 2 nd and 3 rd order KP scheme ( t = 0.25s).
with the 2 nd order scheme, is shown in Figure 8. Overshoot of the 2 nd order scheme at the first peak point can be clearly seen in the exploded view.This exploded view shows that all fixed step-length solvers produce very similar results.The variable step-length solvers also produce similar results.ode15s whose order can be varied from 1st order to 5th order (Shampine et al., 2003) shows small oscillatory behavior for the 3 rd order KP scheme (see Figure 9).Overshoot in the first peak with the 2 nd order KP scheme is also observed for the variable step-length solvers.Other variable step-length solvers: ode23s, ode23t, ode45 show similar behavior; other than the small overshoot in the solution with the 2 nd order scheme, the results appear to closely match each other.Figure 10 shows simulation results for ode23s solver.Both the ode23s plotted in Figure 10 and the ode23t plotted in Figure 11 produce similar results with negligible deviation in the 2 nd order scheme.Figure 11 shows the results of the ode23t solver for both the 3 rd order and the 2 nd order schemes.Variable step-length solver ode45 also produces similar results and there is no significant difference observed other than the small overshoot with the 2 nd order KP scheme.Simulation results have been plotted in Figure 12.The solution from all variable step-length solvers are plotted in Figure 13 from where it can be seen that all the variable step-length solvers produce very similar results.
Computation time needed for simulating 150 minutes of the river system is given in Table 2.According to the table, variable step-length solvers use less computation time for the 3 rd order KP scheme compared to the 2 nd order KP scheme.Less computation time might be a result of the elimination of the flux lim-  Figure 14: Numerical oscillation for all variable steplength solvers for 2 nd and 3 rd order KP scheme (Dissanayake et al., 2016).
iters in the polynomial reconstruction in the 3 rd order KP scheme.However, fixed step-length solvers on the other hand use more computational time for the 3 rd order KP scheme.
Table 3 shows minimum and maximum [min, max] time steps which are automatically taken by the variable step-length solvers in MATLAB.According to this table, variable step length solvers in the 3 rd order KP scheme have taken large time steps compared to the 2 nd order KP scheme.The solution smoothing ability of the 3 rd order KP scheme is higher compared to the 2 nd order KP scheme.Hence, 3 rd order scheme might eliminate unnecessary calculation of small oscillations, ultimately leading to larger time steps.For both the 2 nd order and the 3 rd order KP scheme, the steady state water level at the Grønvollfoss dam was computed.The results from these two schemes are very close to each other.For all solvers, the steady state water level was more or less the same at 17.09 m.

Numerical stability analysis
The simulation results with the 3 rd order KP schemes were checked for numerical stability.In the computation of the volumetric flow rate, oscillatory nature was observed for variable step-length solvers for the 3 rd order KP scheme.Similar behavior has been previously observed in the volumetric flow rate calculations with the 2 nd order KP scheme (Dissanayake et al., 2016).In addition, when the order of the time integrator is higher than the order of the spatial discretization, such oscillatory nature escalates.Figure 14 shows numerical oscillation in variable step-length solvers in MATLAB for the 2 nd order KP scheme.According to Figure 14,  the ode45 solver exhibits a more oscillatory solution.However with the 2 nd order KP scheme, for all the fixed step-length solvers: the Euler method, the RK2 and RK4, such oscillation is negligible (see Figure 15).
The 3 rd order KP scheme shows some similar behavior as the 2 nd order KP scheme.However, for the variable step-length solver ode45, it produces minor oscillation compared to the 2 nd order KP scheme.The solver ode15s with the 3 rd order KP scheme still exhibits marked oscillation.Comparing the solutions of both the 2 nd and the 3 rd order KP schemes, it seems appropriate to choose the order of the time integrators equal to or lower than the order of the spatial discretization.The solution of the 3 rd order KP scheme for all variable step-length solvers is plotted in the Figure 16.Numerical oscillation for all variable step-length solvers for the 3 rd order KP scheme can be observed in Figure 16.
With the fixed step-length solvers (the Euler method, the RK2 and the RK4), the solutions have negligible oscillations and are very close to each other as shown in Figure 17.However, the 3 rd order KP scheme shows some smoothing of the solution as compared to 2 nd order KP scheme.This observation indicates that the 3 rd order scheme is less suitable than 2 nd order scheme when the true solution has sudden (step) changes.

Hydraulic Jump computation
Hydraulic jump is a phenomenon which can be observed in open channel flow, river flow, etc. (Salame et al., 2015).A hydraulic jump condition occurs where the flow velocity is decreased abruptly (Hawary et al., 2016).With high speed flow suddenly reduced (as in water spill out from a dam into a still water bed or slow moving river), part of its kinetic energy will be converted into potential energy and subsequently cause a height increment in the liquid.Such a phenomenon highly depends on the flow velocity: with a change from super critical flow to sub critical flow, a hydraulic  Flow velocity of a channel may change from sub critical to critical at a steep change of the bottom slope.If the bottom becomes flat after this steep slope section, the flow changes from super critical flow into subcritical flow, hence hydraulic jump may be observed.In such a situation, a suitable numerical scheme for the Saint-Venant equation should be able to capture the hydraulic jump.In order to check the ability of capturing a hydraulic jump, both the 2 nd order KP scheme and the 3 rd order KP scheme were checked with a modified river geometry.The bottom of the river is assumed flat after a steep drop of the bed; here the third/last section of the river bed is assumed horizontal.Table 4 summarizes the modifications that were introduced to the river reach.All other quantities except the river bed geometry were kept unchanged.
For the 2 nd order KP scheme, simulation with RK2 as the time integrator captured a hydraulic jump as shown in Figure 18.A similar hydraulic jump should be observed with the 3 rd order KP scheme, unless the assumed smooth solution of the higher order method falsely removes the jump.In order to check the hydraulic jump with the modified parameters, simulation was done using the 3 rd order KP scheme by using RK2 as a time integrator.The result is plotted in Figure 19.
According to Figure 19, the 3 rd order scheme has failed to capture the hydraulic jump after the steep bottom section.Hence, the extra smoothing of the 3 rd order KP scheme sometimes has this undesirable side effect of failing to reveal the true characteristic of the system.On the other hand according to Figure 18, the 2 nd order KP scheme has successfully captured the hydraulic jump condition.
In order to observe the hydraulic jump, flow velocity has to be changed from supercritical to subcritical.Flow velocity changes can be studied by observing the change in the Froude number.The Froude number (Fr number) is a dimensionless number which explains speed-length ratio, and can be written as (xia Li et al., 2015), Froude number for the 2nd order KP scheme

Froude number
Figure 20: F r number for 2 nd order KP scheme.
Here µ 0 is the characteristic flow velocity, g 0 is the characteristic of the external field (gravitational constant) and l 0 is the characteristic length.F r number for the current scenario of the river reach can be written as, Here c is characteristic water wave propagation velocity, which can be written as, where g is gravitational constant and h is water wave height.F r number for the 2 nd order KP scheme is plotted in Figure 20.F r number increases over 1: F r ≥ 1 means flow change from subcritical to supercritical.According to Figure 20 for the for modified river geometry, the flow becomes super critical at the steep slope.Hence a hydraulic jump should occur.Such a hydraulic jump clearly can be seen with the 2 nd order KP scheme.
The computation of the Froude number with the 3 rd order KP scheme is shown in Figure 21.It can be clearly seen that the Froude number never becomes greater than 1 i.e. there is no change of flow characteristics (from supercritical to subcritical) and hence no hydraulic jump.According to Figure 20 and Figure 21, the 2 nd order KP scheme is preferable when the system has sudden step changes in its inputs: volumetric flow rate, step changes etc. Due to extra smoothing ability shown by the 3 rd order KP scheme, usage of such scheme is limited.

Conclusion
Based on the simulated results, the 3 rd order KP scheme is less oscillatory than the 2 nd order KP scheme.As the 3 rd order scheme uses a quadratic polynomial in the reconstruction, a smoother solution can be expected.Computational time for the 3 rd order KP scheme for all variable step-length solvers is smaller compared to the 2 nd order KP scheme.Computational time mainly depends on the number of calculations inside the algorithm.The 2 nd order KP scheme used a flux limiter, hence an extra set of calculations are associated with the 2 nd order KP scheme.In the 3 rd order scheme, the flux limiters have been eliminated in the reconstruction procedure.However, for fixed step-length solvers, a higher computational time was observed with the 3 rd order scheme.A high oscillatory solution for the volumetric flow rate was observed for the ode45 solver with the 2 nd order scheme.On the other hand, in the 3 rd order scheme the ode15s solver showed oscillatory nature.While closely considering this oscillatory nature, it can be seen that when the order of the spatial discretization increases, the order of the time integrators can also be increased in tandem without oscillation in the final solution.Based on this observation, it can be concluded that the order of the time integrator should be essentially of lower order compared with the order of the spatial discretization.This observation confirms the observations of Dissanayake et al. (2016).
On the other hand, even though the 3 rd order KP scheme appears to be more accurate and stable, this scheme produces too smooth solution with abrupt changes of the inputs.Hence, it can be concluded that the developed 3 rd order scheme is less suitable for systems with abrupt changes in the inputs.This result is supported by the lack of the 3 rd order scheme to exhibit a hydraulic jump in the solution.

Figure 13 :
Figure13: All variable step-length solvers for 2 nd and 3 rd order KP scheme.

Figure 17 :
Figure17: Numerical oscillation for all fixed steplength solvers for 3 rd order KP scheme.

Figure 21 :
Figure21: F r number for 3 rd order KP scheme.

Table 2 :
Computational time required for the 2 nd and 3 rd order KP scheme.Figure 9: ode15s for 2 nd and 3 rd order KP scheme.Figure 11: ode23t for 2 nd and 3 rd order KP scheme.Figure 12: ode45 for 2 nd and 3 rd order KP scheme.

Table 3 :
Computational time with variable step-length MATLAB solvers for the 2 nd and the 3 rd order KP scheme.