Exact solution of the Linear Parabolic Approximation for flow-depth based diffusive flow routing

Accepted Manuscript Research papers Exact solution of the Linear Parabolic Approximation for flow-depth based diffusive flow routing Luigi Cimorelli, Luca Cozzolino, Andrea D'Aniello, Domenico Pianese PII: DOI: Reference:

S0022-1694(18)30448-7 https://doi.org/10.1016/j.jhydrol.2018.06.026 HYDROL 22875

To appear in:

Journal of Hydrology

Received Date: Revised Date: Accepted Date:

19 February 2018 4 June 2018 11 June 2018

Please cite this article as: Cimorelli, L., Cozzolino, L., D'Aniello, A., Pianese, D., Exact solution of the Linear Parabolic Approximation for flow-depth based diffusive flow routing, Journal of Hydrology (2018), doi: https:// doi.org/10.1016/j.jhydrol.2018.06.026

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

1

Exact solution of the Linear Parabolic Approximation for flow-depth based diffusive flow

2

routing

3

Luigi Cimorellia*, Luca Cozzolino b, Andrea D’Anielloa, Domenico Pianesea.

4 5 6

a

Department of Civil, Architectural and Environmental Engineering, Univ. of Naples Federico II.

7

b

Department of Engineering, Univ. of Naples Parthenope

8

*Corresponding author, mail to [email protected]

9 10

Keywords: Diffusive wave; Analytical solution; Backwater effects; Boundary conditions; Channel routing;

11

Simplified model.

12 13

Abstract

14 15

The Linearized Parabolic Approximation of the Saint Venant Equations is often used for flood

16

routing, but the corresponding analytic solutions usually neglect the downstream boundary

17

conditions. This is an issue because the presence of hydraulic structures at the end of a river reach,

18

or natural morphologic conditions, may influence the wave propagation. In order to take into

19

account realistic boundary conditions, namely a stage-hydrograph upstream and a stage-discharge

20

relationship downstream, a new set of exact solutions of the Linearized Parabolic Approximation of

21

the Saint Venant Equations with uniformly distributed lateral inflows is presented. This exact

22

solution is demonstrated by using it as a building block in a simplified flood routing model, whose

23

numerical results are compared with those supplied by laboratory experiments from literature.

24

Finally, the new solutions are used to analyze the range of validity of semi-infinite channel models.

25

The comparison shows that semi-infinite channel models are accurate when convective effects are

26

prevailing on diffusive effects, and the downstream boundary condition corresponds to uniform flow

27

conditions. In addition, the results show that semi-infinite channel models based on the knowledge 1

28

of the upstream stage-hydrograph can predict flow depths better than those making use of a flow

29

hydrograph, while being practically equivalent in predicting flow rates.

30 31

1 Introduction

32

The Saint Venant Equations (SVEs) are commonly used for the description of one-dimensional

33

unsteady flow in rivers or artificial channel networks (Cunge et al. 1980). The exact solution of the

34

SVEs is available for schematic cases only, and appropriate numerical methods, such as the Finite

35

Difference method (Vukovic and Sopta 2003, Xing and Shu 2006), the Finite Element method

36

(Szymkiewicz 1991), and the Finite Volume method (Garcia-Navarro and Vazquez-Cendon 2001,

37

Liang and Marche 2009, Cozzolino et al. 2012), can be used for real world applications. Despite the

38

availability of reliable numerical procedures, the application of the full one-dimensional Saint-

39

Venant numerical models may suffer from a number of drawbacks, such as the lack of cross-section

40

topographic data, obtained by traditional ground or LiDAR surveys (Hutton et al. 2012), the

41

incorrect application of internal and external boundary conditions (Cunge et al. 1980, Sobey 2001),

42

the arbitrariness in the choice of cross-sections location and orientation, and the uncertainties of the

43

friction coefficients (Hunter et al. 2007). These practical issues may overshadow the advantage of

44

an accurate physical representation. In addition, the Courant-Friedrichs-Lewy restriction of the time

45

step in explicit time-marching numerical models (Murillo et al. 2006) may advice against the use of

46

the SVEs in the case of forecasting applications, control problems, and design applications

47

involving Monte Carlo approaches (Neal et al. 2012). To overcome these drawbacks, simplified

48

models are considered an attractive alternative to fully dynamic models (Hunter et al. 2007, Bates et

49

al. 2010, Cheviron and Moussa 2016).

50

The Parabolic Approximation (PA) is a non-linear simplification of the SVEs, obtained by

51

neglecting inertial terms, which are often negligible in practical hydrologic applications (Chow et

52

al. 1988, Kazezyılmaz-Alhan 2012). The PA model has been used in numerous applications, as

53

river flow propagation (Cappelaere 1997), flood events simulation (Moussa and Bacquillon 2009), 2

54

drainage systems simulation (Xu and Schwanenberg 2016), flow propagation with significant lateral

55

contributions (Spada et al. 2017), optimal design of rural drainage systems (Cimorelli et al. 2013a).

56

Moreover, the applicability limits of the non-linear model have been actively studied (Ponce and

57

Simons 1977, Ponce et al. 1978, Fread 1985, Ferrick 1985, Dooge and Napiorkowski 1987, Chung

58

et al. 1993, Moussa and Bocquillon 2000, Tsai 2003, 2005, Fan and Li 2006, Prestininzi 2008,

59

Cimorelli et al. 2014).

60

In practical cases requiring massive and fast calculations, the Linearized Parabolic

61

Approximation (LPA), obtained by linearization of the non-linear PA model, may help to reduce the

62

computational burden. The analytic solutions of the LPA model have been used as a building block

63

in numerous applications, such as non-linear flow routing with simplified backwater modelling

64

(Todini and Bossi 1986, Cimorelli at al. 2013b), multilinear approaches (Becker and Kundzewicz

65

1987, Perumal et al. 2008, Perumal et al. 2010), flow routing modeling in semi-distributed rainfall-

66

runoff model (Todini 1996, Moussa et al. 2007), flow routing in non-uniform channels by means of

67

a cascade of diffusive linear models (Cimorelli et al. 2015), propagation of the uncertainty (Chang

68

and Yeh 2016a,b), river flow modelling in karstic areas (Charlier et al. 2015a,b), optimal

69

management of small agricultural catchments (Colin et al. 2011), influence of rainfall spatial

70

variability on rainfall–runoff modelling (Emmanuel et al. 2015). The main advantage of analytical

71

solutions is the unconditional stability of the associated numerical models, thus allowing for long

72

time steps, especially when accuracy is not an issue (Wang et al. 2014). Interestingly, the LPA is

73

related to other widely used simplified models such as the Muskingum-Cunge (MC) model (Cunge

74

1969, Weinmann and Laurenson 1979), that can be transformed into an equivalent LPA by

75

matching physical and numerical diffusion (Cunge 1969), even though it was originally derived

76

from the Kinematic Wave simplification of the SVEs.

77

In the literature, many efforts have been spent to develop simplified stage-hydrograph

78

routing models because upstream gauge information is usually available in the form of stage-

79

hydrographs. This is also convenient because upstream flow-hydrograph boundary conditions 3

80

require the evaluation of the rating curve, which is often lacking or inaccurate (Spada et al. 2017).

81

Following this approach, physically based stage-hydrograph models, such has the variable

82

parameter Muskingum stage-hydrograph, are now available (Perumal and Ranga Raju 1998a,b,

83

Perumal et al. 2007, 2009, 2010a,b). Unfortunately, the analytic solutions of the LPA with stage-

84

hydrograph imposed upstream do not take into account downstream boundary conditions (see the

85

review in Cimorelli et al. 2014), and the semi-infinite channel analytic solution by Hayami (1951) is

86

usually assumed in practical applications. The only exception is the model by Tingsachali and

87

Manandhar (1985), whose applicability is limited because it is based on the knowledge of the

88

downstream stage-hydrograph. As pointed out by numerous authors (Chung et al. 1993, Munier et

89

al. 2008, Cimorelli et al. 2014, 2015, Cozzolino et al. 2014), structures or channel geometric

90

singularities (weirs, sluice gates, bridge piles), which influence the stage-discharge relationship

91

downstream, may significantly affect the flow dynamics. This makes the semi-infinite channel

92

solution by Hayami (1951) inadequate in many cases.

93

Therefore, it is clear that the construction of the LPA analytic solution in finite-length

94

channels, with stage-hydrograph upstream boundary condition and given rating curve imposed

95

downstream, is an attractive task. Following the approach used in Cimorelli et al. (2014) for the

96

case of discharge-hydrograph upstream boundary conditions, this analytical solution is supplied in

97

the present work by applying the Laplace transform technique. The findings are compared with the

98

experimental results of a laboratory test, showing the suitability of the finite-length channel

99

analytical solution in simplified flow routing models.

100

Only the flow depth response of the semi-infinite channel LPA, with stage-hydrograph

101

imposed upstream with lateral inflow, is available in the literature (Hayami 1951). For this reason,

102

this solution is complemented in the present work by the corresponding discharge response to the

103

stage-hydrograph imposed upstream with lateral inflow. This allows the comparison with the new

104

finite-length channel solution, in order to evaluate the applicability limits of the Hayami (1951)

105

solution. The results of this comparison are of interest not only for the LPA-based flow routing 4

106

models, but also for those MC models (Perumal and Ranga Raju 1998a,b, Perumal et al. 2007,

107

2009, 2010a,b) where a stage-hydrograph is considered upstream.

108

The paper is organized as follows: in Section 2, the novel finite-length channel analytical

109

solutions are presented; in Section 3, the complete semi-infinite channel solution is derived for

110

upstream stage-hydrograph and lateral inflow; in Section 4, the finite-length and the semi-infinite

111

channel analytical solutions are compared by considering the results of a laboratory experiment

112

from the literature; in Section 5, the applicability limits of the semi-infinite channel solution is

113

accomplished by comparing its unit-step response with that of the finite-length channel; finally, the

114

conclusions are drawn in Section 6.

115 116

2 Finite-length channel analytical solution with stage-hydrograph imposed upstream

117 118

In the present section, the LPA model is presented, and the corresponding analytic solution in finite-

119

length channels, with stage-hydrograph upstream boundary condition and given rating curve

120

imposed downstream, is derived.

121 122

2.1 LPA governing equations

123

The Parabolic Approximation (PA) of the Saint Venant Equations for a prismatic channel can be

124

written as (Cunge et al. 1980)

125

126

h  Q B q   x t   h  S  S  0 0 f   x

,

(1)

127 128

where Q(x,t) is the discharge, h(x,t) is the flow depth, B(x,h) is the channel width, S 0 x  is the bed

129

slope, and S f Q, h, x  is the friction-slope. If Eq. (1) is linearized around a uniform state 5

130

characterized by flow depth hr, discharge Qr, and width Br, the Linearized Parabolic Approximation

131

(LPA) is obtained. With reference to a channel of length L, the dimensionless form of the LPA

132

model is (Cimorelli et al. 2014, 2015)

133

134

 Q* h*    q*   x* t * ,  *  h  r Q*  r h*  0 1 2   x*

(2)

135 136

where

137 138

t *  tC r L , x *  x L , h*  h' hr , Q*  Q' Qr , q *  qL Qr ,

(3)

139 140

are the dimensionless expressions of the time variable t, the space coordinate x, the first order flow

141

depth increment h’(x,t) = h(x,t) - hr , the first order discharge increment Q’(x,t) = Q(x,t) - Qr, and the

142

lateral inflow per unit of length q(x, t), respectively. The dimensionless coefficients appearing in

143

Eqs. (2) and (3) are defined as

144

145

Cr  

1 Br

 S f h   S   ,   hr Br C r , r1  L f  Q  S Q  Qr   f r

 Qr  S f    L , r2    r hr  h  r

(4)

146 147

where the partial derivatives of the slope friction Sf are calculated with reference to the flow depth

148

hr, the discharge Qr, and the reference abscissa xr. The celerity Cr is related to the hydraulic

149

diffusivity coefficient

150 151

Dr 

1

Br S f Q r

(5) 6

152 153

through the Pèclet number

154 155

Pe  Cr L Dr ,

(6)

156 157

which expresses the relative importance of flow advection and diffusion. It is easy to see from (4)

158

and (5) that Pe =  r1.

159 160

2.2 Finite-length channel general solution in the Laplace domain

161 162

The Laplace transform of Eq. (2) is

163

164

   Q    s h   h0*  q  * x ,    h  r Q  r h  0 2  x * 1

(7)

165 166

   where s is the complex Laplace variable, Qx* , s  , h x* , s  , and q x* , s , are the Laplace transforms

167

of Q* x* ,t * , h* x* ,t * , and q* x* ,t * , respectively, while h0* x *   h * x * ,0 is the initial condition of

168

the first order flow depth increment. Under the assumption that the lateral inflow is uniform in

169

space and the initial condition h0* x *  is null, the general solution of Eq. (7) is (Cimorelli et al.

170

2014)





 









171 172

 Q x* , s  *   *   x ,s h x , s 

   



 Q0, s  *    1 x , s qˆ s  ,  h 0, s 







(8)

173 7









174

    T T where Qx* , s  h x* , s  is the frequency domain solution vector, Q0, s  h 0, s  is the

175

  frequency domain upstream boundary condition vector, and q (s)  q ( x * , s) . In Eq. (8), the matrix

176 177

 x* , s   12 x* , s   x* , s    11 *  *  21x , s   22 x , s 

(9)

178 179

is the state-transition matrix, while x * , s   11 x * , s  11 x * , s  is a vector taking into account

180

the lateral inflow contribution. The expressions of  x* , s and x * , s  are deduced in Cimorelli et

181

al. (2014) and are reported in the Appendix A of the present work for the reader’s ease.

T





182

The Eq. (8) represents the Laplace domain general solution of Eq. (7), and it must be

183

complemented with the correct boundary conditions and lateral inflows to obtain the desired

184

particular solutions.

185 186

2.3 Particular solutions with upstream stage-hydrograph and downstream stage-discharge function

187 188

A stage-hydrograph h0, t  is assumed as upstream boundary condition, while the stage-discharge

189

relationship QL, t   f SG hL, t  is assumed downstream. In order to derive the analytical solution

190

of Eq. (7), the linearized dimensionless expressions of these boundary conditions are needed.

191

Upstream, the first order increment h' 0, t   h0, t   hr of the stage-hydrograph function is

192

considered, and the corresponding dimensionless form h* 0, t   h' 0, t  hr

193

obtained. Downstream, the stage-discharge relationship is linearized around the flow depth hr and is

194

then reduced in dimensionless form. The corresponding expression is

is immediately

195 196

 

Q* 1, t * 

 * * h 1, t  , 

(10) 8

197 198

where the dimensionless abscissa x* = 1 corresponds to x = L. The parameter   Br Cr  k , where

199

k  df SG hr  dh , quantifies how much the downstream boundary conditions deviates from the

200

steady uniform flow along the channel. The Laplace transform of Eq. (10) finally leads to the

201

Laplace domain downstream boundary condition

202 203

   Q 1, t *  h 1, t * .

 



 

(11)

204 205

The substitution of Eq. (11) into Eq. (8) for x* = 1 supplies:

206 207

  1, s    1 22 1, s  11 1, s    121 1, s     Q0, s    12 h 0 , s  q s  .  11 1, s    1 21 1, s   11 1, s    1 21 1, s 

(12)

208 209

The Eq. (12) represents the unique relationship between upstream flow depth and upstream

210

discharge that is compatible in a finite-length channel where Eq. (11) describes the downstream

211

boundary condition. The substitution of Eq. (12) into the first of Eq. (8) allows to eliminate Q0, s  ,

212

leading to the flow depth response of a finite-length channel to a stage-hydrograph imposed

213

upstream, with lateral inflow contribution:



214 215

   h x * , s   wu x * , s  h 0, s   wl x * , s  q s  .

(13)

216 217

In Eq. (13), the function wu x * , s  represents the flow-depth Laplace-domain response to an

218

upstream stage-hydrograph unit-pulse, while wl x * , s  represents the flow-depth Laplace-domain

219

response to a lateral inflow unit-pulse. These functions are defined as: 9

220

221





wu x , s  e *

Pe x* 2

 1  x*   1  x*    Pe 2    sinh  Z    Z cosh Z   2   2 , 1 1     Pe 2    sinh  Z    Z cosh Z 2  2 

(14)

222

 1  x*   1  x*  Pe 2    sinh  Z    Z cosh Z  1 e  2   2  wl x * , s   1 1 s s     Pe 2    sinh  Z    Z cosh Z 2  2 



223



Pe x* 2



Pe x * 1 2



 *  1   sinh x Z  2 Pe e 1  2   s 1  1  Pe 2    sinh  Z    Z cosh Z 2  2 

,

(15)

224 225

where Z  Pe 2  4Pe s .

226

In a similar manner, the substitution of Eq. (12) into the second of Eq. (8) leads to the

227

discharge response of a finite-length channel to a stage-hydrograph imposed upstream, with lateral

228

inflow contribution:

229 230

   Qx * , s   yu x * , s  h 0, s   yl x * , s  q s  .

(16)

231 232

In Eq. (16), yu x * , s  is the Laplace-domain discharge-response to an upstream stage-

233

hydrograph unit-pulse, while yl x * , s  is the response corresponding to a lateral inflow unit-pulse.

234

These functions are defined as:

235

10

236





y u x , s  e *

Pe x * 2

 1  x*   1  x*  Pe  2s sinh Z   Z cosh Z   2   2 , 1  1  Pe 2    sinh  Z    Z cosh Z 2  2 

(17)

237

  1  x*  Z   Z cosh Z  1 e  2   2  yl x * , s   s s 1  1  Pe 2    sinh  Z    Z cosh Z 2  2  . * *       x x * Pe x 1 Z   Pe sinh  Z   Z cosh    1e 2  2   2    s 1  1  Pe 2    sinh  Z    Z cosh Z 2  2 



238

Pe x * 2





Pe  2s sinh 1  x

*

(18)

239 240

The time domain expressions of Eqs. (14)-(15) and (17)-(18) are

241 242







  



(19)

  



(20)

  



(21)

  



(22)

wu x * , t *   c1,n x * exp  s n t * , n 0

243 244







wl x * , t *   c2,n x * exp  s n t * , n 0

245 246







yu x * , t *   c3,n x * exp  s n t * , n 0

247 248







yl x * , t *   c4,n x * exp  s n t * . n 0

249

11

250

where sn is the n-th pole of the corresponding Laplace-domain expression. The technique used to

251

obtain the coefficients ci ,n x *  (i = 1,..,4; n = 0, 1, 2…) and to calculate the poles is described in

252

Appendix B.

253

The unit-step response, frequently used to discretize inputs of any shape by means of

254

rectangles (Chow et al., 1988; Wang et al., 2014), is defined as the primitive of the unit pulse

255

response. Therefore, the Laplace Transform Inversion theorem (Ahlfors 1979) can be used to show

256

that the unit-step responses corresponding to the unit-pulse responses of Eqs. (19)-(22) are:

257

258

  

*  c e Pe    1e Pe x 1, n x Wu x , t   wu x , d   exp  s n t * ,  Pe sn e   1 n 0 0



*



*



t*



t*

*



*

*



*



*

*



(23)

259

260



Wl x , t   wl x , d  *

*

*

0

x*



e 

Pe x*



  P   1     c x  exp  s t ,

1



e

e   1 Pe

n 0

*

2,n

*

n

sn

(24)

261 *  c e P 3, n x  Yu x , t    y u x , d  P  exp  s n t *  , e    1 n 0 s n 0 t*

262

*

e

*

*

*

*

(25)

e

263

264

  

* 1  e Pe  Pe   1  c 4,n x Yl x , t   yl x , d  x   exp  s n t * . Pe sn Pe   e  1 n 0 0



*

*



t*



*

*



*

*







(26)

265 266

The Eqs. (23) and (25) represent the flow-depth and the discharge response, respectively, to

267

the stage-hydrograph unit-step input at the upstream end of the channel. Similarly, Eqs. (24) and

268

(26) represent the flow-depth and the discharge response, respectively, to the lateral-inflow unit-

269

step input when the upstream flow depth is given. The Eqs. (23)-(26) are valid under the condition

270

  2 . If   2 , the hyperbolic sines at the denominators of Eqs. (14)-(15) and (17)-(18) vanish,

271

and the time domain solutions are simplified (see Appendix B). 12

272 273

Remark 1. It is important to note that the Eqs. (24) and (26) differ from the Eqs. (38) and

274

(36) in Cimorelli et al. (2014), that represent the flow-depth and discharge response, respectively, to

275

the lateral-inflow unit-step input when the upstream discharge is given. The expressions of Eqs.

276

(23)-(26) are totally new.

277 278 279

3 Analytical solution in semi-infinite channel with stage-hydrograph imposed upstream

280 281

The flow-depth response of Eq. (7) is already known in the semi-infinite channel when a stage-

282

hydrograph is imposed upstream (Hayami 1951). In particular, the function

283

284



*

Wu , x , t

*



 x*  t * 1  erfc * 2  4t

 1 Pe x*  x*  t *  Pe   e erfc *  2  4t

 Pe  

(27)

285 286

represents the flow-depth response to the upstream stage-hydrograph unit-step input. The

287

corresponding discharge response is not available in the literature, and it is derived in the present-

288

section. First, the linearized parabolic rating curve is obtained from the second of Eq. (2):

289 290

Q *   h* 

 h* . Pe x*

(28)

291 292 293

If the position h * x * , t *   Wu , x * , t *  is made in Eq. (27), and the Eq. (27) is substituted in Eq. (28), the expression

294

13

295

Yu ,

x  t    x*  t *   *   x , t  erfc Pe   e 4t * 2 Pe t *  4t 



*

*



*



* 2

Pe

(29)

296 297

Is obtained after the position Yu , x * , t *   Q * x * , t *  . The Eq. (29) represents the discharge response

298

to the stage-hydrograph unit-step imposed upstream.

299 300

The substitution of the second of Eq. (7) into the first supplies the following flow-depth advection-diffusion equation (Hayami 1951)

301 302



1  2 h* h* h*  *  *  q* , 2 * Pe x x t

(30)

303 304

which formally coincides with the discharge advection-diffusion equation (ADE) by Cunge (1969).

305

For this reason, the solutions of Eq. (30) can be borrowed from the solutions of the discharge ADE

306

form. If the analytic derivation contained in Franchini and Todini (1988) is applied to Eq. (30), it

307

can be shown that the function

308

309





Wl , x * , t *  t * 

 x*  t * 1 * t  x * erfc * 2  4t





 1  x*  t * * Pe   t *  x * e Pe x erfc *  2  4t





 Pe  

(31)

310 311

is the flow-depth response to the uniform lateral-inflow unit-step input in the semi-infinite channel

312

for given flow-depth imposed upstream. If the position h * x * , t *   Wl , x * , t *  is made in Eq. (31),

313

the substitution in Eq. (28) supplies the expression

314

14



*

Yl , x , t 315



*



 *  t  x * e P x 2 Pe e

*

x  t 

* 2

 x*  t *   * 1  *      t   t  x  erfc P e   * 2 2 Pe   4t   x*  t *  erfc  Pe  *  4t 

t*    e Pe

*

4t

*

Pe

*

(32)

316 317

which represents the discharge response to the unit-step uniform lateral inflow in the semi-infinite

318

channel for given flow-depth imposed upstream.

319 320 321

Remark 2. The discharge responses of Eqs. (29) and (32) are novel and complement the flow-depth responses of Eqs. (27) and (31), which are already available in the literature.

322 323 324

4. Numeric application of the finite-length and semi-infinite channel analytical solutions

325 326

A simple linear flood routing model, based on the analytic findings of the preceding sections, is

327

here presented. The flood routing model is then used to reproduce the results of a laboratory test

328

case from the literature.

329 330

4.1 A simple linear flood routing model

331 332

In the case of uniform linear channel, a simplified version of the ILTDFR model (Cimorelli et al.

333

2015)

334

I t *  h * 0, t *  h 0, t *  hr hr be the dimensionless stage-hydrograph at the inlet of the channel.

335

This function is represented by the sequence of rectangular pulses

 

can



be

used

 



to

route

a

stage-hydrograph



336

15

imposed

upstream.

Let

337

1 Ik  * t

kt *

 I  d k  1, 2, ... ,

(33)

k 1t *

338 339

where t * is the dimensionless time interval. The flow-depth response of the channel at the station

340

x* and time t n*  nt * (n = 1, 2, …) is immediately calculated as

341 342





n

 







h * x * , t n*   I k W x * , n  k  1t *  W x * , n  k t * ,

(34)

k 1

343 344

where W x * , t *   Wu x * , t *  when the downstream boundary condition is taken into account, while

345

W x * , t *   Wu , x * , t *  whether the effect of the downstream boundary condition is neglected

346

(semi-infinite channel)

347 348

4.2 Numeric application

349 350

The linear model of Eq. (33) is applied to the experimental results by Rashid and Chaudrhy (1995).

351

The original experimental setup consisted of a flume, Lf = 21 m long, with uniform longitudinal bed

352

slope, compound rectangular cross-section (rectangular main channel with rectangular floodplains),

353

and Manning’s roughness coefficient dependent on the flow depth. An inclined sluice gate with bed

354

sill was present at the downstream end of the laboratory flume, and the corresponding stage-

355

discharge relation was expressed by

356 357

f SG  h   C  h  a  , m

(35)

358

16

359

where a is the sill height. In Eq. (35), the constants C = 9.35 and m=1.14 were experimentally

360

determined by Rashid and Chaudhry (1995). Nine gauging stations were located at different

361

distances from the flume inlet (Table 1), and the actual channel length used in the computations is L

362

= 18.6 m, because the upstream boundary condition is given by the flow-depth record at the

363

Gauging Station 1. The geometric characteristics of the flume and the initial conditions of the

364

experiments are reported in Rashid and Chaudrhy (1995), and their complete description is not

365

repeated here.

366 367

[Insert Table 1 about here]

368 369

First, the model of Eq. (33) is applied with W x * , t *   Wu x * , t *  (finite-length channel

370

model) to the Test 1, where the flow occupies the entire compound cross-section. The parameters

371

(Pe, α, and β), needed for the calculation of the finite-length channel response function Wu x * ,t * 

372

of Eq. (25), are evaluated with respect to the initial condition of Test 1, and no calibration is

373

performed. In particular, hr is computed as the average value of the initial backwater profile, Qr is

374

set equal to the initial discharge, and the reference conditions used for the downstream boundary

375

conditions are those corresponding to the reference discharge Qr. The results of the simulation with

376

time step t =10 s (thin continuous line) are compared in Figure 1 with the experimental results

377

(circles). From the figure, it is evident that the uniform linear channel model reproduces nicely the

378

flow-depth hydrograph at the Gauging Station 2 (Fig. 1a), which is close to the inlet, while the

379

essentials of the flow propagation phenomenon are captured at the Gauging Station 5 (Fig. 1b),

380

located at the middle of the channel. It is very important to observe that better results could be

381

obtained if a cascade of uniform diffusive channels were used, after a calibration of the parameters

382

(Cimorelli et. al. 2015).

17

383

In Figure 1, the results of a full Saint Venant finite-difference model by Rashid and

384

Chaudrhy (1995) are also reported. The comparison with the experimental results shows that the

385

parabolic linear model does not behave worse than the full Saint Venant equations. For this test

386

case, the discrepancy between the experimental results and the numerical results at the Gauging

387

Station 5 can be explained recalling that the cross-section is compound, with initial conditions

388

corresponding to more than 0.20 m of flow-depth in the main channel and less of 0.01 m of flow-

389

depth over the lateral floodplains.

390

The finite-length model is then applied to the Test 2, where the flow is entirely contained

391

into the main rectangular channel. The close inspection of Figure 2 shows that the uniform linear

392

channel model reproduces nicely the flow-depth hydrograph at both the Gauging Stations 2 (Fig.

393

2a) and 5 (Fig. 2b). The improvement of the numerical results can be explained recalling that now

394

the wetted area corresponds to a compact cross-section. In the same figure, the results of the

395

numerical model by Rashid and Chaudrhy (1995) are also reported. Again, the comparison shows

396

that the results supplied by the finite-length linear model are of comparable quality to those

397

supplied by the full Saint Venant model.

398 399

[Insert Figure 1 about here]

400

[Insert Figure 2 about here]

401 402 403

Tests

1

and

2

are

repeated

considering

the

semi-infinite

channel

model

404

W x * , t *   Wu , x * , t *  in Eq. (33). The results are represented in Figures 1 and 2, respectively, with

405

symbol. The comparison with the experimental results clearly shows that the semi-infinite channel

406

model behaves worse than the finite-length model.

407

18

408

5. Comparison between semi-infinite channel and finite-length channel analytical solutions

409 410

The findings of the previous section suggest that the finite-length channel analytical solutions are

411

more accurate than the semi-infinite channel analytical solutions in reproducing real flows.

412

Therefore, the applicability of the semi-infinite channel solutions is discussed by comparison with

413

the corresponding finite-length channel solutions. The analysis is accomplished considering the

414

unit-step responses, because these functions are used as a building block in practical numerical

415

methods where input hydrographs are represented by means of rectangular pulses sequences (Todini

416

and Bossi 1986, Chow et al. 1988, Cimorelli et al 2013b, 2015, 2016, Wang et al. 2014). Values of

417

Pe and  representative of realistic situations are considered, while the constant  = 1 is taken as

418

this coefficient is just a scaling factor. First, the dependence of the flow depth responses on Pe and

419

 is commented, and finally the Nash-Sutcliffe Efficiency Index is used for an objective

420

comparison.

421 422

5.1 Flow-depth response to a unit-step stage-hydrograph imposed upstream

423 424

The Eqs. (23) and (27) represent the flow-depth responses of finite-length and semi-infinite

425

channels, respectively, to a unit-step stage-hydrograph imposed upstream. The functions are plotted

426

in the plane (t*, h*) of Figure 3 for different values of the parameters Pe (0.5, 2.5 and 5) and  (0.5,

427

1, and 1.5), and at different locations x* (0.2, 0.4, 0.6, 0.8, 1). A continuous line is used for the

428

finite-length channel response, while a dashed line is used for the finite-length channel response.

429

First, the finite-length channel flow-depth response h*  Wu  x* , t *  of Eq. (23) is considered.

430

The value   0.5 (panels in the left column of Figure 3) is representative of a downstream

431

boundary condition where a free fall or a critical flow depth occurs at the downstream end.

19

432

Congruently with the physical meaning of the downstream boundary condition, the flow depth

433

decreases with x* at a given t* , for every value of Pe.

434

For   1 (panels in the central column of Figure 3), as time increases, the flow-depth tends

435

asymptotically to h*  1 everywhere. This behavior is consistent with   1 being representative of

436

a normal depth downstream boundary condition.

437

The value   1.5 (panels in the right column of Figure 3), is representative of an obstacle

438

downstream, such as a sluice gate or a weir. In this case, the Eq. (23) exhibits a noticeable

439

backwater effect propagating from downstream with increasing t*. Correspondingly, h*  Wu  x* , t * 

440

has a twisted behavior because the flow-depth decreases with x* for smaller values of t* (when the

441

wave entering upstream has not reached the downstream boundary condition), and increases with x *

442

for greater values of t* (because the wave is reflecting at the downstream boundary conditions).

443 444

[Insert Figure 3 about here]

445 446

The response h*  Wu ,  x* , t *  of Eq. (30) depends on Pe only, and not on . For this reason,

447

the corresponding plots in the left column of Figure 3 are equal to those of the central and the right

448

column, and no twisted behavior is exhibited. A closer inspection shows that lower values of Pe

449

cause a faster growth of h*  Wu ,  x* , t *  , but the asymptotic value of h* obtained for increasing t*

450

is not affected. In addition, h* decreases with x* at a given time t*.

451 452

4.2 Flow-depth response to a unit-step uniformly distributed lateral inflow

453 454

The functions h*  Wl  x* , t *  of Eq. (24) and h*  Wl  x* , t *  of Eq. (31) are plotted in the plane (t*,

455

h*) of Figure 4 with the same parameters Pe, , and x*, used for Figure 3. The finite-length channel

20

456

response to the lateral inflow is quite different from the response to the upstream stage-hydrograph.

457

When   0.5 (free fall condition at the downstream end), the time-asymptotic flow-depth value of

458

* * * h*  Wl x* , t * increases from x  0 to a maximum in x  0.8 , and then decreases from x  0.8

459

to x *  1 . This differs from the case of h*  Wu  x* , t *  in Figure 3, where the flow-depth is always

460

decreasing. When   1 (uniform flow at the downstream end), the time-asymptotic value of

461

h*  Wl x* , t *

462

value h* = 1 at every cross-section. Finally, when   1.5 (sluice gate or weir at the downstream

463

end), h*  Wl  x* , t *  increases regularly with x* at any instant t*, while h*  Wu  x* , t *  exhibits a

464

twisted behavior.









increases along the channel, while h*  Wu  x* , t *  exhibits the same asymptotic

465 466

[Insert Figure 4 about here]

467





468

The semi-infinite channel response h*  Wl , x* , t * can reproduce the results of Eq. (24)

469

only for high values of Pe, and close to the upstream channel end. When Pe  5 , the two solutions

470

behave similarly for x* from 0.2 to 0.6, and for each value of . When Pe  2.5 , the two solutions

471

exhibit similar behaviors for x* from 0.2 to 0.4, and for   1. Finally, the plots of Eqs. (24) and

472

(31) are very different at any position along the channel and for any value of  when Pe  0.5 .

473 474

4.3 Discharge response to a unit-step stage-hydrograph imposed upstream

475 476

The functions Q*  Yu  x* , t *  of Eq. (25) and Q*  Yu ,  x* , t *  of Eq. (29), are plotted in the plane

477

(t*, Q*) of Figure 5 with the same parameters of Figures 3 and 4. From Figure 5, it is evident that

478

the function Q*  Yu  x* , t *  exhibits a peak, originated at the channel inlet, that moves and decays 21

479

along the channel. Actually, when the inlet flow depth h* increases instantly from 0 to 1 (unit step

480

stage-hydrograph), the linearized looped rating curve of Eq. (28) prescribes that the derivative at the

481

right-hand side is a negative delta function. For this reason, the second term at the right-hand side of

482

Eq. (28) is a positive delta function that sums up with the first term (a unit step). For t > 0, this delta

483

function is convected and diffused along the channel, while the discharge entering upstream tends

484

to an asymptotic value that depends on h* and its space derivative. A closer inspection of Eq. (28)

485

shows that the inlet discharge peak weakens when Pe increases, and vanishes when Pe tends to

486

infinity, while the dependency on  of the time-asymptotic discharge value decreases with Pe.

487 488

[Insert figure 5 about here]

489 490

The function Q*  Yu ,  x* , t *  of Eq. (29) is inadequate to mimic the finite-channel model

491

when Pe = 0.5. Conversely, the solution is well reproduced almost everywhere along the channel,

492

with exclusion of x* = 1, when Pe  5 . Finally, the correspondence between Eqs. (25) and Eq. (29)

493

is strongly influenced by  when Pe  2.5 , but there is a reasonable correspondence for x *  0.6 .

494 495

4.4 Discharge response to unit-step lateral inflows.

496 497

A plot of the functions Q*  Yl  x* , t *  of Eq. (26) and Q*  Yl ,  x* , t *  of Eq. (32) is

498

reported in Figure 6 with the same parameters considered in the preceding subsections. The

499

inspection of the figure shows that the infinite-length and semi-infinite channel solutions exhibit

500

negative discharges in the vicinities of x* = 0, for low values of Pe, when uniform lateral inflow is

501

imposed This feature can be easily explained recalling that Eqs. (26) and (32) are found under the

502

assumption h*(0, t*) = 0 for t* > 0, which allows the separation of effects from the case of a stage-

503

hydrograph boundary condition imposed upstream. According to Eq. (28), this implies that the 22

504

discharge Q* along the channel is negative when Pe h *  h * x * . This is precisely what happens

505

for low values of Pe in the vicinities of x* = 0, where h* is small and h * x *  0 . Physically

506

speaking, the conditions in which these new solutions are derived can be representative of a water

507

body in the upstream cross section (a lake or the sea, a tank served by a pump) where the constant

508

level allows the lateral discharge injected along the channel to partly flow upstream. When a null or

509

positive discharge is expected at the upstream end, Eqs. (26) or (32) are not autonomous andmust be

510

complemented by a proper upstream stage-hydrograph. Interestingly, this also means that the

511

upstream stage-hydrograph cannot be assigned arbitrarily if lateral inflows are present. Therefore,

512

when the lateral inflow is the only input, the exact solution with discharge imposed upstream

513

(Cimorelli et al. 2014) should be preferred.

514 515

[Insert Figure 6 about here]

516 517

The Figure 6 shows that the function Q*  Yl ,  x* , t *  of Eq. (32) does not match the

518

function Q*  Yl  x* , t *  of Eq. (26) for Pe  0.5 . When Pe  2.5 , the two solutions behave

519

qualitatively in a similar manner, but there are significant differences for   0.5 and   1 , while

520

the differences are small for   1.5 . For Pe  5 the results given by the two solutions are similar

521

except for the location x*  1 .

522 523 524

4.5 Backwater effect analysis

525 526

To establish the range of validity of the semi-infinite channel solutions presented in Section 3, the

527

Efficiency Index (NS) by Nash and Sutcliffe (1970)

528 23

  f t   f t  NS  1    f t   f  I max

529

2

i 1 I max i 1

i

o

i

(36)

2

o

i

o

530 531

is used. In Eq. (36), f t i  and f o t i  are the approximating and the true function, respectively,

532

evaluated at the time instant ti (with i = 1, 2,..,Imax), while f o is the mean value of the true function.

533

The index NS ranges in (   , 1], where NS = 1 indicates a perfect match between the two models,

534

while NS  0 means that the prediction made by the approximating model has the same accuracy

535

of the time-averaged value of the true one. In this subsection, the values in x* = 1 of the semi-

536

infinite channel solutions are taken as the approximating models, while the step responses of Eqs.

537

(23)-(26) at the same abscissa are taken as the true models. All the computations are performed

538

considering Pe  0, 10 and   0, 3 . In the following, the interval NS > 0.7 is considered as

539

representative of a satisfactory matching (Moriasi et al., 2007).

540

The contour plot in Figure 7 represents the results of Eq. (36) when Eqs. (23) and (27) are

541

compared. It is evident from the figure that there is a region Pe ,    1, 10 0.6, 1.4 where

542

h*  Wu , x* , t * , (Eq. [27]) provides a good approximation of h*  Wu , x* , t *

543

practical point of view, this is significant because the interval 0.6    1.3 is wide enough to

544

represent the downstream boundary conditions usually encountered in river applications. This result

545

is consistent with the numerical application of Section 4: in experiment 1, Pe = 1.04 and  = 0.126,

546

while in experiment 2 Pe = 0.016 and  =0.0124, and in both cases the semi-infinite solution failed

547

to reproduce the numerical experiment.

548

The outcomes are quite different (Figure 8) when h*  Wl  x* , t *  (Eq. [24]) and h*  Wl ,  x* , t * 

549

(Eq. [31]) are compared, because in this case the region of validity of h*  Wl ,  x* , t *  is reduced to

550

Pe ,    2, 10 1.1, 1.5.







551 24



(Eq. [23]). From a

552

[Insert Figure 7 about here]

553

[Insert Figure 8 about here]

554 555

The contour plot of Eq. (36), computed using







Q*  Yu x* , t *



(Eq. [25]) and

556

Q*  Yu , x* , t *

557

boundary condition, the discharge response Q*  Yu ,  x* , t *  matches the finite-length channel

558

response Q*  Yu  x* , t *  better than the flow-depth response h*  Wu ,  x* , t *  matches the finite-

559

length model h*  Wu  x* , t *  , because the area of validity starts from a point Pe ,    0.5, 0.5 and

560

spreads up covering the whole range of considered  at Pe  2 . Actually, the propagation of the

561

discharge by means of semi-infinite channels model seems more robust than the propagation of the

562

flow depth, which is greatly influenced by the downstream boundary condition. A similar

563

conclusion can be deduced by looking at the contour plot of Figure 10, where Eq. (36) is calculated

564

considering Q*  Yl  x* , t *  (Eq. [26]) and Q*  Yl ,  x* , t *  (Eq. [32]).

(Eq. [29]), is depicted in Figure 9. Despite the use of the flow depth as upstream

565 566

[Insert Figure 9 about here]

567

[Insert Figure 10 about here]

568 569

5 Discussion

570 571

The exact solutions of Eq. (7) for the cases of finite-length and semi-infinite channels have been

572

systematically compared in Section 4, showing that the Péclet number Pe and the downstream

573

boundary condition coefficient  are the parameters controlling the suitability of the semi-infinite

574

channel solutions. The lower values of Pe correspond to short channels, or to cases where diffusivity

25

575

effects are not negligible with respect to the convection effects. Conversely, the higher values of Pe

576

correspond to long channels or to cases where the convection effects are predominant.

577

Congruently with the physical interpretation of Pe, the comparison between the finite-length

578

and the semi-infinite channel solutions shows that the downstream boundary conditions cannot be

579

neglected when Pe is low, or  is far from 1. Conversely, the use of the semi-infinite channel model

580

is reliable when Pe is high and  is close to 1. In addition, the analysis has shown that the semi-

581

infinite channel model with flow depth imposed upstream predicts Q*(x*, t*) better than h*(x*, t*),

582

because the flow depth is more directly influenced by the downstream boundary condition, but the

583

range of validity of h*(x*, t*) is satisfactory for large Pe. These results can be compared with the

584

findings contained in Cimorelli et al. (2014), where the range of validity of semi-infinite channel

585

solutions was considered when discharge was imposed upstream. For that case, it was concluded

586

that the range of validity of semi-infinite channel models for the prediction of h*(x*, t*) was very

587

narrow, while the range of validity for the prediction of Q*(x*, t*) was larger.

588

In conclusion, semi-infinite channel models that use a stage-hydrograph as upstream

589

boundary condition are able to predict the flow depths better than the semi-infinite channel models

590

where the upstream flow hydrograph is routed, while minor differences are found when the

591

prediction of the discharge is considered. For this reason, semi-infinite channel models with stage-

592

hydrograph imposed upstream should be used when an accurate prediction of the flow depth is

593

required.

594

These conclusions are not limited to the LPA model only, since the downstream boundary

595

conditions are neglected in many of the approximate flood routing models that are used in the

596

hydrology practice. One of the most employed of these models is the MC model, which is

597

equivalent to the Advection Diffusion Equation (ADE) form of the LPA in semi-infinite channels if

598

its parameters are properly chosen (Cunge et al., 1969; Tang et al., 1999; Wang et al., 2006). For

599

this reason, the applicability limits studied here also supply useful indications about the range of

26

600

applicability of the MC model with a stage-hydrograph boundary condition (Perumal and Ranga

601

Raju 1998a,b, Perumal et al. 2007, 2009, 2010a,b).

602 603 604

5 Conclusions

605 606

A novel set of analytical solutions for the linearized parabolic approximation (LPA) of the SVE has

607

been derived in this paper, considering a stage-hydrograph upstream boundary condition, uniform

608

lateral inflows, and a stage-discharge relationship as downstream boundary condition. These

609

solutions are motivated by the fact that downstream boundary conditions (control structures,

610

obstructions such as bridge piles, sills, jumps) can modify significantly the flow dynamics. For

611

further comparison, these analytical solutions have been complemented by the corresponding

612

analytical solutions in the semi-infinite channel.

613

The numerical results of a numerical application have suggested that a simple flood routing

614

model equipped with the finite-length channel exact solution behaves better than the semi-infinite

615

channel solution in reproducing the experimental results obtained in a laboratory channel. This

616

finding has suggested a systematic comparison of finite-length and semi-infinite channel analytic

617

solutions, in order to analyze the applicability of the routing models based on the semi-infinite

618

channel length assumptions. From the comparison, it is evident that the downstream boundary

619

conditions can be neglected only when the Péclet number is high, and the downstream rating curve

620

is not far from a normal depth one. Conversely, the semi-infinite channel models should not be used

621

if the Péclet number is low. In addition, the finite-length channel analytic solutions have shown that

622

the upstream stage-hydrograph and the lateral inflows concur in determining the discharge at the

623

upstream end of the channel. In particular, this means that the upstream stage-hydrograph and the

624

lateral inflows cannot be assigned autonomously if the discharge upstream is assumed positive.

27

625

Finally, when the semi-infinite channel assumption is made in flood routing models, the use of a

626

stage-hydrograph upstream should be preferred for the calculation of flow depths.

627

The authors hope that these results will encourage researchers to develop simplified models

628

based on the imposition of a stage-hydrograph upstream, rather than a flow-hydrograph, possibly

629

accounting for the downstream boundary conditions.

630 631 632

Acknowledgements

633 634

This research was partially funded by the University of Naples Parthenope through the funding

635

programs “Sostegno alla Ricerca individuale 2015-2017” and "Ricerca competitiva triennio 2016-

636

2018".

637 638

Appendix A

639

The components of the state transition matrix  x * , s  are defined as (Cimorelli et al. 2014)

640

641

 11x* , s  

e

Pe x* 2

  x*   x*     Z cosh Z  P sinh Z  ,  e   Z  2  2 

(A.1)

642

643

 12 x* , s   

Pe x* 2

2se Z

 x*  sinh  Z  , 2 

(A.2)

 x* sinh   2

(A.3)

644

645

 21 x * , s  

 2 Pe e

Pe x * 2

 Z

 Z  , 

646 28

647

 22 x* , s  

e

Pe x* 2

  x*   x*     Z cosh Z  P sinh Z  ,  e   Z  2  2 

(A.4)

648 649

where the complex variable Z is introduced in Sub-section 2.3.





The components of the complex vector  x * , s are defined as (Cimorelli et al. 2014)

650 651

652





11 x* , s 

Pe x* 2





x 1 e Pe  2 Pe s  sinh s sPe Z 2 2

*

Pe x* 2

 e  x*  Z   cosh Z  , s  2 

(A.5)

653

654

Pe x* 2

Pe x*

 x*  Pe 2  x*  1 e 21 x* , s   cosh Z   e sinh  Z  . s s 2  s Z 2 





(A.6)

655 656 657

Appendix B

658

In order to find the time domain expressions of Eqs. (15)-(16) and (18)-(19), the inverse Laplace

659

formula

660 661

f t  

1 f s e st ds 2i b

(B.1)

662 663

is used. In Eq. (B.1), f t  is a generic function in the time domain, while f s  is the corresponding

664

image in the Laplace complex domain. The integration contour b is a line parallel to the imaginary

665

axe, at the right of all the singularity of the function f s  . It can be proven that Eqs. (15)-(16) and

666

(18)-(19) are analytic in the complex domain, with the exception of a countable set of singular

667

points that are called poles sn (n = 0, 1, 2…). For this reason, the Eq. (B.1) can be computed by 29

668

means of the residue theorem, which states that the integral at the right hand of Eq. (B.1) is equal to

669

the summation of residues at the singularities of f s  times e st (Ahlfors, 1979):

670





671



f t    res e st f s ; s  s n .

(B.2)

n 0

672 673 674

Since the poles of the functions in Eqs. (15)-(16) and (18)-(19) are simple, the residues can be evaluated by means of

675 676





res e st f s ; s  s n  lim s  s n  f s e st .

(B.3)

s  sn

677

When   2 , the application of Eqs. (B.2)-(B.3) to Eqs. (15)-(16) and (18)-(19) leads to the

678 679

expression of the coefficients that appear in the sums of Eqs. (20)-(23):

680

681

 

c1,n x *  e

Pe x 2

*

 Zn

 1  x*  2

 sin

Z

2 n

  1  x*  Z n   Z n cos Z n    2  1      2 Pe cos Zn  2 

(B.4)



682

683

684

 

c 2,n x *  

e

Pe x * 2

 Zn s n

 1  x*  2

 sin 

  1  x* Z n   Z n cos   2

Z

2 n

P  e   x* Z n   2 Pe e 2  1  1 sin    2 1      2 Pe cos Zn  2 

(B.5)

685

30







 Z n  ,

686

e

 

c 3, n x *  

Pe x* 2

 Zn 



  1  x* Z n   Z n cos  2   2 1  2 Z n     2 Pe cos Zn  2 

Pe  2s n sin 1  x

*





 Z n  ,

(B.6)

687

 

c 4,n x *   688





e

Pe x * 1 2

e



Pe x * 2

 Zn s n

  1  x* Z n   Z n cos  2   2 1  2 Z n     2  Pe cos Zn  2 

Pe  2s n sin 1  x

*





x  Pe sin   1    Z n  2 s n 2

  x*  Z n   Z n cos Zn   2 1      2  Pe cos Zn  2  *

Z

n



  

 Z n  

.

(B.7)

689

In Eqs. (B.4)-(B.7),   Pe 2 1  1 and Z n  4 Pe sn  Pe 2 , while the poles sn are the roots

690 691

of the transcendental equation (see Appendix C):

692 693

1  1  Pe 2   sin Z n    Z n cos Zn   0 . 2  2 

(B.8)

694 695

The Eqs. (B.4)-(B.7) are null for n  0 when   0 , while preserve their expressions when

696

 2    0 . When   2 there is a pole in s   Pe 4 , and the coefficients corresponding to n = 0

697

are defined by

698

699

  e

Pe x* 2

 

Pe x* 1 2

c1,0 x

*

3x * , Pe

(B.9)

700

701

c 2, 0 x *  e





Pe

P 1     e 2 6x* e Pe  Pe   

,

(B.10)

31

702

703

Pe x* 2

  e

c3, 0 x

*





2  x*  1  , 2 Pe

3

(B.11)

704

705

 

c4,0 x *  3e



2  Px   1  e 2  x



Pe x* 1 2

Pe 2

*

Pe Pe   

*

 .

1 

(B.12)

706 707

When   2 there is a pole for 0  s   Pe 4 . In this case, the Eqs. (B.4)-(B.7) modify for

708

n  0 by exchanging the hyperbolic functions with their corresponding trigonometric functions, and

709

the minus sign is taken before Z n2 at the denominator. When   2 , the Eqs. (B4)-(B.7) must be replaced by:

710 711

712

  e

c1,n x

*

2n  1 sin x * 2n  1  

Pe x* 2

 

Pe

 

2

(B.13)

713

714

c 2,n

Pe x* 2

P  e  2n  1  sin x * 2n  1   n 1 x   (1) e 2    s n  Pe 2   

 

e

*

(B.14)

715

716

 

c3,n x * 

e

Pe x* 2

 2n  1  2n  1  

2 Pe

Pe

 2n  1   2n  1   cos x *    sin x *   2 2    

(B.15)

717

 

c 4,n x *  718



2 1

2e

Pe x* 2

 2n  1  2n  1

  2 2n  1  Pe 

n 1

2



Pe x 1 2 *



2

Pe

 2n  1   2n  1   cos x *    sin  x *   2 2    

  * 2n  1   2n  1      2n  1 cos x *    Pe sin  x 2 2  2n  1  Pe      2

e

2

2

32

.

(B.16)

719 720

Note that Eqs. (B.13)-(B.16) are identically null for n  0 .

721 722

Appendix C

723 724

Recalling the expression of Z defined in 3.1, the denominator of Eqs. (17)-(20) can be written as:

725 726

1 1 2  2 2  4 Pe s  Pe   4 Pe s  Pe cosh 4 Pe s  Pe  2  2 

 sinh

(C.1)

727 728

with  defined in Appendix B.

729

With the exception of few cases, the roots of Eq. (C.1) exist when the terms in the square root is

730

negative,

731 732

that



is,



when

s   Pe 4 .

Reminding

that





 

sinh  x  i sin x

and

 

2  x cosh  x  i x cos x (with i   1 ), and making the position z  4 Pe s  Pe , Eq. C.1

can be written as

733 734

1  iA sin z   z  2 

(C.2)

735 736

z where A  z 2   2 and  z   arctan  . Therefore, the poles correspond to the roots of the  

737

equation

738 739

1 z z  arctan   n 2  

(C.3)

740 33

741

with n  N .

742

Attention should be paid to the following situations:

743

1. when   0 , the poles of Eqs. (17)-(20) are given by the roots of Eq. (C3) with n  1 and the

744 745

residuals by Eqs. (B.4)-(B.7); 2. when   0 the poles of Eqs. (17)-(20) are analytical and corresponds to

746 747

sn  

2n  12  2  Pe 4 Pe

(C.4)

4

748 749 750 751 752 753

with n  1 , and the residuals are given by Eqs. (B.9)-(B.12); 3. when  2    0 , the poles of Eqs. (17)-(20) are given by the roots of Eq. (C.3) with n  0 and the residuals are still given by Eqs. (B.4)-(B.7); 4. when   2 there is a pole in sn   Pe 4 , while the other poles are still given by the roots of Eq. (C.3) with n  1 ;

754

5. when   2 , there is a pole corresponding to the root of Eq. (C.1), in the interval

755

s   Pe 4 ; 0, with residual given by modifying Eqs. (B.4)-(B.7) as specified in Appendix B,

756

while the other poles are still given by the roots of Eq. (C.3) with n  1 and residuals

757

corresponding to Eqs. (B.4)-(B.7).

758 759 760

References

761 762 763

Ahlfors LV. (1979) Complex analysis: an introduction to the theory of analytic functions of one complex variable. New York: McGraw-Hill.

34

764

Bates P.D., Horritt M.S., Fewtrell T.J. (2010) A simple inertial formulation fo the shallow

765

water equations for efficient two-dimensional flood inundation modelling, Journal of Hydrology

766

387, 33-45. DOI: 10.1016/j.jhydrol.2010.03.027.

767 768 769 770 771

Becker A, Kundzewicz ZW. (1987) Nonlinear flood routing with multilinear models. Water Resour Res; 23(6):1043–8. Cappelaere B. (1997) Accurate Diffusive Wave routing. ASCE Journal of Hydraulic Engineering 123(3), 174-181. DOI: 10.1061/(ASCE)0733-9429(1997)123:3(174). Chang C.M., Yeh H.D. (2016a) Probability density functions of the stream flow discharge in

772

linearized

773

10.1016/j.jhydrol.2016.10.033.

diffusion

wave

models.

Journal

of

Hydrology

543,

625-629.

DOI:

774

Chang C.M., Yeh H.D. (2016b) Stochastic modeling of variations in stream flow discharge

775

induced by random spatiotemporal fluctuations in lateral inflow rate. Stochastic Environmental

776

Research and Risk Assessment 30(6), 1635-1640. DOI: 10.1007/s00477-015-1170-x.

777

Charlier J. B., Moussa R., Bailly-Comte V., Desprats J.F., Ladouche B. (2015a). How Karst

778

Areas Amplify or Attenuate River Flood Peaks? A Response Using a Diffusive Wave Model with

779

Lateral Flows. In Hydrogeological and Environmental Investigations in Karst Systems (pp. 293-

780

301). Springer Berlin Heidelberg.

781

Charlier J. B., Moussa R., Bailly-Comte V., Danneville L., Desprats J.F., Ladouche B.,

782

Marchandise A. (2015b). Use of a flood-routing model to assess lateral flows in a karstic stream:

783

implications to the hydrogeological functioning of the Grands Causses area (Tarn River, Southern

784

France). Environmental Earth Sciences 74(12), 7605-7616. DOI: 10.1007/s12665-015-4704-0.

785

Cheviron B., Moussa R. (2016) Determinants of modelling choices for 1-D free-surface flow

786

and morphodynamics in hydrology and hydraulics: a review. Hydrology and Earth System Sciences

787

20, 1-32. DOI: 10.5194/hess-20-1-2016.

788

Chow VT, Maidment DR, Mays LW., (1988). Applied hydrology. McGraw-Hill.

35

789 790

Chung WH, Aldama AA, Smith JA., (1993). On the effects of downstream boundary conditions on diffusive flood routing. Adv Water Resour; 16:259–75.

791

Cimorelli L., Cozzolino L., Covelli C., Mucherino C., Palumbo A., Pianese D. (2013a)

792

Optimal design of rural drainage networks. ASCE Journal of Irrigation and Drainage Engineering

793

139(2), 137-144. DOI: 10.1061/(ASCE)IR.1943-4774.0000526.

794

Cimorelli L., Cozzolino L., Della Morte R., Pianese D. (2013b) An improved numerical

795

scheme for the approximate solution of the Parabolic Wave model. Journal of Hydroinformatics

796

15(3), 913-925. DOI: 10.2166/hydro.2013.130.

797

Cimorelli L., Cozzolino L., Della Morte R., Pianese D. (2014) Analytical solutions of the

798

linearized parabolic wave accounting for downstream boundary condition and uniform lateral

799

inflows. Advances in Water Resources 63, 57-76.

800

Cimorelli, L., Cozzolino, L., Della Morte, R., Pianese, D., & Singh, V. P. (2015) A new

801

frequency domain analytical solution of a cascade of diffusive channels for flood routing. Water

802

Resources Research, 51(4), 2393-2411.

803

Colin F., Guillome S., Tisseyre B. (2011) Small catchment agricultural management using

804

decision variables defined at catchment scale and a Fuzzy Rule-based system: a Mediterranean

805

vineyard case study. Water Resources Management 25(11), 2649-2668. DOI: 10.1007/s11269-011-

806

9831-0.

807

Cozzolino L., Della Morte R., Del Giudice G., Palumbo A., Pianese D. (2012) A well-

808

balanced spectral volume scheme with the wetting–drying property for the shallow-water equations.

809

Journal of Hydroinformatics, 14(3), 745-760.

810

Cozzolino, L., Cimorelli, L., Covelli, C., Della Morte, R., & Pianese, D. (2014). Boundary

811

conditions in finite volume schemes for the solution of shallow-water equations: The non-

812

submerged broad-crested weir. Journal of Hydroinformatics, 16(6), 1235-1249.

813 814

Cunge J.A, Holly F.M., Verwey A., (1980). Practical aspects of computational river hydraulics. Pitman, Boston. 36

815 816 817 818

Cunge, J.A. (1969). On the subject of a flood propagation computation method (Muskingum Method), IAHR Journal of Hydraulic Research 7(2), 205-230. Dooge J.C.I., Napiorkowski J.J. (1987). Applicability of diffusion analogy in flood routing. Acta Geophysica Polonica 35(1), 66-75.

819

Emmanuel I., Andrieu H., Leblois E., Janey N., Payrastre O. (2015) Influence of rainfall

820

spatial variability on rainfall–runoff modelling: Benefit of a simulation approach?. Journal of

821

Hydrology 531(Part 2), 337-348. DOI: 10.1016/j.jhydrol.2015.04.058.

822 823 824 825

Fan P, Li JC. (2006). Diffusive wave solutions for open channel flows with uniform and concentrated lateral inflow. Adv Water Resour;29(7):1000–19. Ferrick M.G. (1985) Analysis of river wave types, Water Resources Research 21(2), 209220. DOI: 10.1029/WR021i002p00209.

826

Todini, E., & Bossi, A. (1986). PAB (Parabolic and Backwater) an unconditionally stable

827

flood routing scheme particularly suited for real time forecasting and control. Journal of Hydraulic

828

Research, 24(5), 405-424.

829 830

Franchini, M., & Todini, E. (1988). PABL: A parabolic and backwater scheme with lateral inflow and outflow. In Fifth IAHR international symposium on stochastic hydraulics, Birmingham.

831

Fread D.L. (1985) Applicability criteria for kinematic and diffusion routing models,

832

Laboratory of Hydrology, National Weather Service, NOAA, U.S. Dept. of Commerce, Silver

833

Spring, Md.

834

Garcia-Navarro P., Vazquez-Cendon M.E. (2001) On numerical treatment of the source

835

terms in the shallow water equations, Computer & Fluids 29(8), 951-979. DOI: 10.1016/S0045-

836

7930(99)00038-9.

837

Hayami S. (1951). On the propagation of flood waves. Disaster Prev Res Inst Bull; 1:1–16.

838

Hunter N.M., Bates P.D., Horritt M.S., Wilson M.D. (2007) Simple spatially-distributed

839

models for predicting flood inundation: A review. Geomorphology 90(3-4), 208-225. DOI:

840

10.1016/j.geomorph.2006.10.021. 37

Kazezyılmaz-Alhan, C.M., (2012). An improved solution for diffusion waves to overland

841 842

flow. Appl. Math. Model. 36 (9), 4165–4172.

843

Liang Q., Marche F. (2009) Numerical resolution of well-balanced shallow water equations

844

with

845

10.1016/j.advwatres.2009.02.010.

complex

source

terms,

Advances

in

Water

Resources

32(6),

873-884.

DOI:

846

Moriasi, D. N., Arnold, J. G., Van Liew, M. W., Bingner, R. L., Harmel, R. D., & Veith, T.

847

L. (2007). Model evaluation guidelines for systematic quantification of accuracy in watershed

848

simulations. Transactions of the ASABE, 50(3), 885-900.

849 850 851 852

Moussa R. (1996) Analytical Hayami solution for the diffusive wave flood routing problem with lateral inflow. Hydrol Process; 10:1209–27. Moussa, R., Bocquillon, C. (2009). On the use of the diffusive wave for modelling extreme flood events with overbank flow in the floodplain. Journal of Hydrology, 374(1), 116-135.

853

Moussa R., Chahinian N., Bocquillon C. (2007) Distributed hydrological modelling of a

854

Mediterranean mountainous catchment – Model construction and multi-site validation. Journal of

855

Hydrology 337(1-2), 35-51. DOI: 10.1016/j.jhydrol.2007.01.028.

856 857

Munier S, Litrico X, Belaud G, Malaterre PO. (2008). Distributed approximation of openchannel flow routing accounting for backwater effects. Adv Water Resour;31(12):1590–602.

858

Murillo J., Garcìa-Navarro P., Brufau P., Burguete J. (2006) Extension of an explicit finite

859

volume method to large time steps (CFL>1): application to shallow water flows. International

860

Journal of Numerical Methods in Fluids 50(1), 63-102. DOI: 10.1002/fld.1036.

861 862

Nash, J. E. and J. V. Sutcliffe (1970). River flow forecasting through conceptual models part I - A discussion of principles, Journal of Hydrology.

863

Neal J., Villanueva I., wright N., Willis T., Fewtrell T., Bates P. (2012) How much physical

864

complexity is needed to model flood inundation?, Hydrological processes 26(15), 2264-2282. DOI:

865

10.1102/hyp.8339.

38

866 867 868 869

Perumal M, Ranga Raju KG (1998a). Variable parameter stage-hydrograph routing method: II. Evaluation. J Hydrol Eng ASCE 3(2):115–121 Perumal M, Ranga Raju KG (1998b). Variable parameter stage-hydrograph routing method: I. Theory. J Hydrol Eng ASCE 3(2):109–114

870

Perumal M, Moramarco T, Sahoo B, Barbetta S (2007). A methodology for discharge

871

estimation and rating curve development at ungauged river sites. Water Resour Res W02412

872

43(2):1–22. doi:10.1029/2005WR004609.

873

Perumal M., Moramarco T., Sahoo B., Barbetta S. (2008) Multilinear Diffusion Analogy

874

Model for Real-Time Streamflow Routing. 9th International Congress on Environmental Modelling

875

and Software. 3. https://scholarsarchive.byu.edu/iemssconference/2008/all/3.

876

Perumal M, Moramarco T, Sahoo B, Barbetta S (2010a). On the practical applicability of the

877

VPMS routing method for rating curve development at ungauged river sites. Water Resources

878

Research 46(3), W03522. doi: 10.1029/2009WR008103.

879

Perumal M., Moramarco T., Sahoo B., Barbetta S., Melone F. (2010b) Multilinear Diffusion

880

Analogy Model for Stage Hydrograph Routing. 9th International Congress on Environmental

881

Modelling and Software. 508. https://scholarsarchive.byu.edu/iemssconference/2010/all/508.

882

Perumal M, Sahoo B, Moramarco T, Barbetta S (2009). Multilinear Muskingum method for

883

stage-hydrograph routing in compound channels. J Hydrol Eng ASCE 14(7):663–670.

884

doi:10.1061/(ASCE)HE.1943-5584.0000029.

885 886 887 888

Ponce V.M., Simons D.B. (1977) Shallow wave propagation in open channel flow, ASCE Journal of the Hydraulic Division 103(HY12), 1461-1476. Ponce V.M., Li R.-M., Simons D.B. (1978) Applicability of kinematic and diffusion models, ASCE Journal of the Hydraulic Division 104(HY3), 353-360.

889

Prestininzi P. (2008) Suitability of the diffusive model for dam break simulation:

890

Application to a CADAM experiment, Journal of Hydrology 361(1-2), 172-185. DOI:

891

10.1016/j.jhydrol.2008.07.050. 39

892 893

Rashid, R. M., & Chaudhry, M. H. (1995). Flood routing in channels with flood plains. Journal of Hydrology, 171(1-2), 75-91. DOI: doi.org/10.1016/0022-1694(95)02693-J

894

Sobey R.J. (2001) Evaluation of numerical models of flood and tide propagation in

895

channels, ASCE Journal of Hydraulic Engineering 127(10), 805-824. DOI: 10.1061/(ASCE)0733-

896

9429(2001)127:10(805).

897

Spada, E., Sinagra, M., Tucciarelli, T., Barbetta, S., Moramarco, T., & Corato, G. (2017).

898

Assessment of river flow with significant lateral inflow through reverse routing modeling.

899

Hydrological Processes.

900

Szymkiewicz R. (1991) Finite-element method for the solution of the Saint Venant

901

equations in an open channel network, Journal of Hydrology 122(1-4), 275-287. DOI:

902

10.1016/0022-1694(91)90182-H.

903 904

Tang, X., Knight, D. W., & Samuels, P. G. (1999). Variable parameter Muskingum-Cunge method for flood routing in a compound channel. Journal of hydraulic research, 37(5), 591-614.

905 906

Tingsanchali T, Manandhar SK. (1985). Analytical diffusion model for flood routing. J Hydraul Eng;111:435–54.

907

Todini E, Bossi A. (1986). PAB (Parabolic and Backwater), an unconditionally stable flood

908

routing scheme particularly suited for real time forecasting and control. J Hydraul Res; 24(5):405–

909

24.

910

Todini E (1996) The Arno rainfall-runoff model. J Hydrol 175:339–382

911

Tsai, C.W. (2003). Applicability of kinematic, noninertia, and quasi-steady dynamic wave

912

models to unsteady flow routing. Journal of Hydraulic Engineering, 129(8), 613-627.

913

Vukovic S., Sopta L. (2003) Upwind schemes with exact conservation property for one-

914

dimensional open channel flow equations, SIAM Journal on Scientific Computing 24(5), 1630-

915

1649. DOI: 10.1137/S1064827501392211.

916 917

Wang, G. T., Yao, C., Okoren, C., & Chen, S. (2006). 4-Point FDF of Muskingum method based on the complete St Venant equations. Journal of hydrology, 324(1), 339-349. 40

918

Wang, L., Wu, J. Q., Elliot, W. J., Fiedler, F. R., and Lapin, S. (2014). Linear diffusion-

919

wave channel routing using a discrete Hayami convolution method. Journal of Hydrology, 509,

920

282-294.

921 922

Weinmann P.E., Laurenson E.M. (1979). Approximate flood routing methods: A review, ASCE Journal of the Hydraulic Division 105(HY12), 1521–1526.

923

Xing Y., Shu C.-W. (2006) High order finite difference WENO schemes with the exact

924

conservation property for the shallow water equations, Journal of Computational Physics 208(1),

925

206-227. DOI: 10.1016/j.jcp.2005.02.006.

926

Xu M., Schwanenberg D. (2017) Sequential and simultaneous model predictive control of a

927

drainage canal network using an implicit Diffusive Wave model. ASCE Journal of Irrigation and

928

Drainage Engineering 143(3), B4016003. DOI: 10.1061/(ASCE)IR.1943-4774.0001082.

929 930 931

41

932

List of Tables

933 934

Table1. Experimental Flume characteristics.

935 936

List of figures

937

Figure 1. Test 1 by Rashid and Chaudrhy (1995): Gauging Station 2 (a), Gauging Station 5 (b).

938

Figure 2. Test 2 by Rashid and Chaudrhy (1995): Gauging Station 2 (a), Gauging Station 5 (b).

939

Figure 3. Flow depth responses to an upstream unitary step stage-hydrograph: Eq. (23), continuous

940

line, and Eq. (27), dashed line.

941

Figure 4. Flow depth responses to a uniform lateral inflow: Eq. (24), continuous line, and Eq. (31),

942

dashed line.

943

Figure 5. Flow rate responses to an upstream unitary step stage-hydrograph: Eq. (27), continuous

944

line, and Eq. (33), dashed line.

945

Figure 6. Flow rate responses to a uniform lateral inflow: Eq. (28), continuous line, and Eq. (34),

946

dashed line.

947

Figure 7. Nash-Sutcliffe contour plot of flow depth responses to an upstream unitary step stage-

948

hydrograph.

949

Figure 8. Nash-Sutcliffe contour plot of flow depth responses to uniform lateral inflow.

950

Figure 9. Nash-Sutcliffe contour plot of flow depth responses to an upstream unitary step stage-

951

hydrograph.

952

Figure 10. Nash-Sutcliffe contour plot of flow depth responses to a uniform lateral inflow.

953 954 955

42

956

Gauge Station distance from inlet [m]

X [m]

Boundary conditions Upstream boundary

1

1.22

0

2 3 4 5 6 7 8 9

3.35 6.09 7.61 11.3 15.2 16.8 18.9 19.8

2.13 4.87 6.39 10.1 14 15.5 17.7 18.6

Inclined gate sill

957 958

Table1. Experimental Flume characteristics.

959

43

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Figure 6

Figure 7

Figure 8

Figure 9

Figure 10

Highlights

960 961 962 963 964

 Diffusive wave analytical solutions requiring upstream stage-hydrograph (flow depth based) are provided;

965

 Downstream boundary conditions are accounted for;

966

 Flow depth based simplified models requiring upstream stage-hydrograph predict the flow

967 968 969

depth better than those requiring an upstream flow hydrograph;  Flow depth-based solutions for lateral inflows are not suitable for flood propagation when upstream boundary conditions are not used;

970

44