ENV5056 Numerical Modeling of

Flow and Contaminant Transport in Rivers

Numerical Solution of Flow

Equations

Asst. Prof. Dr. Orhan GÜNDÜZ

Saint-Venant Equations for River Flow

∂A ∂Q
+ −q = 0
∂t ∂x

∂Q ∂( βQ 2 / A)  ∂h 
+ + gA + Sf + Se  − βqυ x + Wf B = 0
∂t ∂x  ∂x 

Numerical techniques for the solution of expanded Saint-Venant

equations can be given as:

(i) method of characteristics

(ii) finite difference methods
(iii) finite element methods.

2
Of these methods, the finite element method is rarely used when flow is
approximated as one-dimensional such as in the case of Saint-Venant
equations. The other two methods have been commonly applied for the
numerical solution of one-dimensional unsteady flow since 1960s.

The finite difference methods can further be classified as explicit and

implicit techniques, each of which holds distinct numerical characteristics. A
major advantage of the implicit finite difference method over the method
of characteristic and the explicit finite difference technique is its inherent
stability without the requirement to satisfy the Courant condition, which
sets the criteria for the maximum allowable time step. This requirement to
satisfy Courant condition often makes the method of characteristics and
explicit techniques very inefficient in terms of the use of computer time.

Furthermore, certain implicit schemes such as the Preissmann scheme

(Preissmann, 1961) allow the use of variable time and spatial steps, which
make them extremely convenient for applications in routing of flood
hydrographs in river systems.

3
Considering these advantages, the implicit finite difference technique is
commonly used to solve the governing partial differential equations. In the
implicit formulation, all derivative terms and other parameters are approximated
by using the unknowns at the forward time line (j+1) as seen in the x-t grid
shown below. The implicit finite difference method advances the solution from
one time line to the next simultaneously for all points along the time line.
Time Line
Time, t Space Line

(i,j+1) (i+1,j+1)
j+1
Upstream ∆t Downstream

Boundary j Boundary
(i,j) (i+1,j)
Line Line

∆x
0
Distance, x
0 i i+1

The distance-time grid used to formulate the implicit finite difference scheme
4
Of the various implicit schemes that have been developed, the "weighted
four-point" scheme first used by Preissmann (1961) is very advantageous
since it can readily be used with unequal distance steps, which is particularly
important for natural waterways where channel characteristics are highly
variable even in short distances.

Similarly, the applicability of unequal time steps is another important

characteristic of the Preissmann scheme, particularly in the case of
hydrograph routing where floodwaters would generally rise relatively quickly
and recess gradually in time.

In this scheme, the four grid points from the (j)th and (j+1)th time lines are
used to approximate the terms in the differential equation. A weighing factor,
θ, is used in the approximation of all terms of the equation except for the time
derivatives in order to adjust the influence of the points (i) and (i+1).

5
The approximations of the derivatives and constant terms in the four-
point weighted difference scheme are given as follows for a dummy
parameter K.

∂K K ij++11 − K ij+1 K ij +1 − K i j
Time derivative =ψ + (1 − ψ )
j
∂t ∆t ∆t j

∂K K i j++11 − K i j +1 K ij+1 − K i j
Space derivative =θ + (1 − θ )
∂x ∆x i ∆x i

Other terms K i +j +11 + K i j +1 K i +j 1 + K i j

K =θ + (1 − θ )
2 2

where i and j are subscripts on x and t axis, respectively, ψ and θ are

weighing factors between 0 and 1, and ∆x and ∆t are incremental reach
length and time, respectively

6
It is possible to create different modifications of the implicit scheme for
various values of the weighing factors. Many researchers preferred using a ψ
value of 0.5 and approximated the time derivative at the center of grid
between (j)th and (j+1)th time lines On the other hand, different values of θ
were used depending on the particular application.

If a θ value of 0.5 is used, the corresponding scheme is called the "box-

scheme", which is a centered finite difference scheme. Similarly, a scheme
with a θ value of 1.0 is known as the "fully-implicit scheme" in space.

The weighted four-point implicit scheme is unconditionally stable for any

time step if the value of θ is selected between 0.5 and 1.0. In addition to
stability criteria, analysis on the influence of the weighing factor on the
accuracy of computations revealed that the accuracy decreases as θ departs
from 0.5 and approaches to 1.0. This effect became more pronounced as the
magnitude of the computational time step increased. Furthermore, analysis
also revealed that a θ value of between 0.55 and 0.6 provided unconditional
stability and good accuracy, which makes this scheme superior compared to
the explicit scheme that requires time steps of less than a critical value
determined by the Courant condition.

7
Finite Difference Formulation with Preissmann Scheme

The finite difference equations of the Saint-Venant equations are discretized

in the x-t plane using the approximations given above. Two algebraic
equations are obtained as a result of this approximation, representing the
partial differential equations of continuity and momentum.

In the approximations given below, some of the variables are defined at the
nodes of the grid (i.e., Q, h, and A), whereas the others are defined for the
reach (i.e., Sf, Se, sc, sm and ∆x).

The variables defined for the reach are denoted by the subscript (i+½) in
order to represent that these variables are defined half way between nodes
(i) and (i+1), or on the reach. It is also seen that the solution plane is
represented by a total of N nodes with node numbers starting from 1 and
running through N.

8
 Node Reach

Time, t

∆t
(i,j+1) (i+1,j+1)
j+1

∆x
j
Upstream (i,j) (i+1,j) Downstream

Boundary Boundary

i th

Grid

1 2 3 i i+1 N-2 N-1 N

Distance, x

The x-t solution plane used to generate the finite difference equations

9
The space derivatives in Saint-Venant equations are approximated as

∂Q Qi +j +11 − Qi j +1 Qi +j 1 − Qi j
=θ + (1 − θ )
∂x ∆xi +1/ 2 ∆xi +1/ 2

∂h hi j++11 − hi j +1 hij+1 − hij

=θ + (1 − θ )
∂x ∆x i +1 / 2 ∆x i +1 / 2

∂( β Q 2 / A ) ( βQ 2 / A) ij++11 − ( βQ 2 / A) ij +1 ( βQ 2 / A) ij+1 − ( βQ 2 / A) ij
=θ + (1 − θ )
∂x ∆x i +1 / 2 ∆x i +1 / 2

The time derivatives in Saint-Venant equations are approximated as

∂A Ai +j +11 − Ai +j 1 Ai j +1 − Ai j
= 0.5 + 0.5
∂t ∆t j ∆t j

∂Q Qi +j +11 − Qi +j 1 Qi j +1 − Qi j
= 0.5 + 0.5
∂t ∆t j ∆t j
10
Finally, the constant terms such as q, Sf, Se and A are approximated as

q = θq ij++11/ 2 + (1 − θ )q ij+1 / 2

j +1 j
S f = θS f i +1 / 2 + (1 − θ )S f i +1 / 2

j +1 j
S e = θS e i +1 / 2 + (1 − θ )S e i +1 / 2

Aij++11 + Aij +1 Aij+1 + Aij j +1 j

A =θ + (1 − θ ) = θ Ai + (1 − θ )A i
2 2

where the variables with subscripts (i+½) are defined for the reach
between nodes (i) and (i+1) and the variable with bar (−) represents an
average of the corresponding values at nodes (i) and (i+1).

11
The following explicit formulations are used to define the above variables

2
n i Qi Qi Qi Qi
(Sf )i +1 / 2 = 2 2 4/3
= 2
c1 Ai R i Ki

K ec i  Q  2 Q  
2
(S e )i +1 / 2 =   −   
2g∆x i +1 / 2  A  i +1  A  i 

Q i + Q i +1
Qi =
2

Ai + Ai +1
Ai =
2

Ai Ai
Ri = ≈
Pi Bi

B i + B i +1
Bi =
2

K i + K i +1
Ki =
2 12
 The Finite Difference Form of Continuity Equation
Substituting the necessary terms in the continuity equation yields:

Ai +j +11 − Ai +j 1 Ai j +1 − Ai j Qi +j +11 − Qi j +1 Qi +j 1 − Qi j
0.5 j
+ 0.5 j
+θ + (1 − θ ) − θ qij++1/1 2 − (1 − θ )qij+1/ 2 = 0
∆t ∆t ∆xi +1/ 2 ∆xi +1/ 2

If terms with common multipliers are rearranged, equation takes the

following form

Ai +j +11 − Ai +j 1 Ai j +1 − Ai j  Qi +j +11 − Qi j +1 j +1   Qi +j 1 − Qi j j 
0.5 + 0.5 + θ  − qi +1/ 2  + (1 − θ )  − qi +1/ 2  = 0
∆t j ∆t j  ∆ xi +1/ 2   ∆ x i +1/ 2 

One might obtain the final finite difference form of continuity

equation by multiplying this equation with ∆xi+1/2:

∆xi +1/ 2
θ Qi j++11 − Qi j +1 − ∆xi +1/ 2 qij++1/1 2  + (1 − θ ) Qi +j 1 − Qi j − ∆xi +1/ 2 qij+1/ 2  +  A
j  i +1
j +1
+ Ai
j +1
− A j
i +1 − A j

i  =0
2 ∆t
13
 The Finite Difference Form of Momentum Equation

Substituting the necessary terms in the momentum equation yields:

Qi +j +11 − Qi +j 1 Qi j +1 − Qi j
0.5 + 0.5 +
∆t j ∆t j
( β Q 2 / A)ij++11 − ( β Q 2 / A)ij +1 ( β Q 2 / A)ij+1 − ( β Q 2 / A)ij
θ + (1 − θ ) +
∆xi +1/ 2 ∆xi +1/ 2
 hi +j +11 − hi j +1 hi +j 1 − hi j 
(
g θA i
j +1
+ (1 − θ ) A θ

i
j
)
∆xi +1/ 2
+ (1 − θ )
∆xi +1/ 2
j +1 j j +1 j
+ θ S f i +1/ 2 + (1 − θ ) S f i +1/ 2 + θ Sei +1/ 2 + (1 − θ ) Sei +1/ 2  −

(θ (β qυ ) j +1
x i +1/ 2 + (1 − θ )( β qυ x )ij+1/ 2 ) = 0

14
 If terms with common multipliers are rearranged, equation takes the
following form

Qi +j +11 − Qi +j 1 Qi j +1 − Qi j
0.5 j
+ 0.5 j
+
∆t ∆t
 ( β Q 2 / A)ij++11 − ( β Q 2 / A)ij +1 j +1  h
j +1
i +1 − hi
j +1
j +1 j +1  j +1

θ + g Ai  + S f i +1/ 2 + Se i +1/ 2  − ( β qυ x )i +1/ 2  +
 ∆ xi +1/ 2  ∆ xi +1/ 2  
 ( β Q 2 / A)ij+1 − ( β Q 2 / A)ij j h
j
i +1 − hi
j
j j  j

(1 − θ )  + g Ai  + S f i +1/ 2 + Se i +1/ 2  − ( β qυ x )i +1/ 2  = 0
 ∆ xi +1/ 2  ∆ xi +1/ 2  

15
 One might obtain the final finite difference form of momentum
equation by multiplying this equation with ∆xi+1/2:

∆xi +1/ 2
j 
 Qi
j +1
+ 1 + Qi
j +1
− Qi
j
+ 1 − Qi
j
 +
2∆t

θ ( β Q 2 / A)ij++11 − ( β Q 2 / A)ij +1 + g Ai hi +j +11 − hi j +1 + ∆xi +1/ 2 S f i +1/ 2 + ∆xi +1/ 2 Se i +1/ 2 − ∆xi +1/ 2 ( β qυ x )ij++1/1 2  +
j +1

 ( j +1 j +1
) 

(1 − θ ) ( β Q 2 / A)ij+1 − ( β Q 2 / A)ij + g Ai hi +j 1 − hi j + ∆xi +1/ 2 S f i +1/ 2 + ∆xi +1/ 2 Se i +1/ 2 − ∆xi +1/ 2 ( β qυ x )ij+1/ 2  = 0
( )
j j j

 

16
The terms with subscript (j) in continuity and momentum equations are
known either from initial conditions or from the solution of Saint-Venant
equations at the previous time line.

Since cross sectional area (A) and channel top width (B) are functions of
water surface elevation (h), the only unknown terms in these equations are
discharge (Q) and water surface elevation (h) at the (j+1)th time line at nodes
(i) and (i+1).

Therefore, there are only four unknowns in these equations. All remaining
terms are either constants or are functions of these unknowns. The finite
difference forms of continuity and momentum equations are solved for each
grid shown before.

The resulting two algebraic equations obtained by the application of the

weighted four-point scheme are nonlinear and an iterative solution
technique is required.

17
As there are N-1 grids in a time line, a total of 2*(N-1) equations are
formed for one time line between the upstream and downstream
boundary.

There are two unknowns (Q and h) in each of the N nodes giving a total
of 2N unknown along each time line. The system of 2(N-1) equations
with 2N unknown require two additional equations to be determinate.

These two additional equations are supplied by the upstream and

downstream boundary conditions.

The resulting system of 2N non-linear equations with 2N unknowns is

commonly solved by the Newton-Raphson iterative technique to handle
the non-linearity in the equations.

18
The Boundary Conditions

The two additional equations required for the system to be determinate come
from the external boundaries at the upstream and downstream locales of the
reach. Two types of boundary conditions can be specified at the upstream
boundary:

1. Discharge Hydrograph
2. Stage (or Depth) Hydrograph

At the downstream boundary, it is possible to specify one of the five available

conditions given as:

1. Discharge Hydrograph
2. Stage (or Depth) Hydrograph
3. Single-Valued Rating Curve
4. Looped Rating Curve
5. Critical Flow Section

19
Upstream Boundary Condition

The upstream boundary conditions could be discharge or a stage hydrograph. A

discharge hydrograph is specified as

Q1j +1 − Q(t ) = 0

Similarly, a stage condition is specified as

h1j +1 − h(t ) = 0

where Q1j+1 and h1j+1 are the discharge and stage to be computed at the
upstream boundary node and Q(t) and h(t) are hydrographs inputted to the
model.

20
Downstream Boundary Condition

The downstream boundary condition can be a discharge or a stage hydrograph. It

is also possible to use a single or looped rating curve that define the relationship
between the stage (or depth) and discharge. Moreover, the downstream
boundary can also be a critical flow section. The corresponding equations for
these conditions are given as follows

Discharge Hydrograph

When a discharge hydrograph is used at the downstream boundary, one would

obtain the following equation

Q Nj +1 − Q (t ) = 0

where QNj+1 is the discharge to be computed at the downstream boundary

and Q(t) is the discharge hydrograph value.
21
Improper use of a discharge hydrograph at the downstream boundary can result
in gross errors. The imposed flow may exceed the capacity of the channel to
deliver water to that node. Moreover, if a discharge condition is also specified at
the upstream boundary, any error in the values of these flows is reflected in the
water levels which might sometimes cause the stream be partially or completely
dry. In this regard, a discharge condition at the downstream boundary is rarely
used in unsteady flow simulations

Stage Hydrograph

If a stage hydrograph is specified at the downstream boundary, one would obtain

hNj +1 − h(t ) = 0

where hNj+1 is the stage to be computed at the downstream boundary and h(t) is
the stage hydrograph value.

22
 Single-Valued Rating Curve

If a single-valued rating curve is used as the downstream boundary condition,

the boundary condition equation is written as

Q Nj +1 − Q' (t ) = 0

where QNj+1 is the discharge to be computed at the downstream boundary and

Q’(t) is the discharge value computed from the rating curve. If the rating curve
is expressed in a tabular, piecewise linear way, any intermediate discharge is
computed by linearly interpolating between the two data sets

Qk +1 − Qk j +1
Q' (t ) = Qk +
hk +1 − hk
hN − hk( )
where Qk, Qk+1, hk and hk+1 are consecutive tabular data sets of the rating
curve and hNj+1 is the stage at the downstream boundary.

23
 Looped Rating Curve

If a looped rating curve is used as the downstream boundary condition, the

boundary condition equation is again written as

Q Nj +1 − Q' (t ) = 0

where QNj+1 is the discharge to be computed at the downstream boundary

and Q’(t) is the discharge value computed from the rating curve. The rating
curve discharge value is computed using the Manning’s equation

j +1
 1.0 
Q' (t ) =  (
AR 2 / 3  S 1f / 2 ) j +1

 n N N

The loop is produced by using the friction slope rather than the channel
slope. The friction slope exceeds the bottom slope during the rising limb of
the hydrograph and falls below it during the recession limb.

24
The friction slope is computed using the full-dynamic wave form of the
momentum equation after neglecting the contributions of lateral inflow and
assuming that the momentum correction factor and sinuosity is equal to unity

1 ∂Q 1 ∂ (Q 2 / A) ∂h
Sf = − − −
gA ∂t gA ∂x ∂x

In its discretized form, this equation is written as:

Sf N
j +1
=−
1
  −
(
 Q Nj − Q Nj −1  1  Q 2 / A ) − (Q
j
N
2
/A )
N −1
j
 hj −hj 
− N N −1

gANj  ∆t  gA  ∆x   ∆x 


Since all the terms in this equation are known from the previous time lines, the
friction slope can easily be calculated and substituted in the Manning’s equation.
It should be noted that the loop rating equation allows the unsteady wave to pass
the downstream boundary with minimal disturbance by the boundary itself, and
hence, is extremely suitable when the boundary condition is located at an
arbitrary location in the waterway rather than at an actual control structure such
25
as a weir or a dam
 Critical Flow Section

If a critical flow section is selected as the downstream boundary condition,

the boundary condition equation is written as

Q Nj +1 − Qcr (t ) = 0
where QNj+1 is the discharge to be computed at the downstream boundary
and Qcr(t) is the critical discharge value computed according to the
following equation

g 3/ 2
Qcr = A
B
This condition is particularly suitable when the routing ends at an entrance
to a free fall or at a transition to a steep reach

26
 The Newton-Raphson Algorithm

The Newton-Raphson method is the most common iterative technique used

for the solution of a system of non-linear equations. It provides a very
efficient means of converging to a root given a sufficiently good initial
guess. Denoting the system of equations in vector form:

Fi ( x1 , x 2 , x 3 ,..., x 2 N ) = 0 i = 1,2,3,...,2 N

If x denotes the entire vector of unknown variables, xi, and F denotes the
entire vector of functions,

Fi, each of the functions Fi can be expanded in the
neighborhood of x as a Taylor series expansion.

2 N ∂F
Fi ( x + δ x) = Fi (x) + ∑



2
i
δ x j + O(δ x ) i = 1, 2,3,..., 2 N
j =1 ∂x j
27
The matrix of first partial derivatives appearing in the above equation is
called the Jacobian matrix:

∂Fi
J=
∂x j







2
F ( x + δ x ) = F ( x ) + J ⋅ δ x + O(δ x )

Neglecting the higher order terms and setting the left hand-side equal to
zero, one can obtain a set of linear equations that are solved for the
corrections δ x






J ⋅δ x = −F
This matrix equation is solved by a matrix solver such
as Gaussian
elimination or LU decomposition for the unknown, δ x , and an the
improved estimate of the solution is obtained by:

xnew = xold + δ x




28
The iteration process is continued until a predetermined convergence level
is achieved. It is important to assign different convergence criteria for Q and
h values. A generally accepted convergence criteria for h is 0.01 ft or 3 mm
(Fread, 1985). The convergence criteria for Q is computed according to the
following equality:

ε Q = ε h BV

where εQ and εh are convergence criteria for discharge and water surface
elevation, respectively; B is the top width of the cross section and V is the
flow velocity.

The convergence process depends on a good first estimate for the unknown
variables. A reasonably good estimate for the first time step is to use the
initial condition of discharge (Q) and water surface elevation (h). For any
other time step, the first estimates of the unknown variables can be obtained
by using the linearly extrapolated values from solutions at previous time
steps.

29
 Partial Derivatives of Finite Difference Equations
One of the discouraging issues about the Newton-Raphson technique is the
need to evaluate the partial derivative terms of the Jacobian matrix. When the
equations get complex, it might get very difficult to differentiate the difference
equations. However, it is this feature of the technique that provides fast
convergence to the root.

In this regard, the difference forms of continuity and momentum equations are
partially differentiated with respect to the unknown terms h and Q at the (j+1)th
time line for the reaches (i) and (i+1).

30
Partial Derivatives of the Continuity Equation:

The partial derivatives of the continuity (C) equation with respect to hi, hi+1,
Qi and Qi+1 (all at (j+1)th time line) are computed as follows:

∂C ∂C ∂Ai j +1 ∆xi +1/ 2

j +1
= j +1 j +1
= 
j  i
B j +1

∂hi ∂Ai ∂hi 2 ∆t

∂C ∂C ∂Ai j++11 ∆xi +1/ 2

j +1
= j +1 j +1
=  B j +1
j  i +1 

∂hi +1 ∂Ai +1 ∂hi +1 2 ∆t

∂C
j +1
= −θ
∂Q i

∂C
=θ
∂Q ij++11

31
 Partial Derivatives of the Momentum Equation:

The partial derivatives of the momentum (M) equation with respect to hi,
hi+1, Qi and Qi+1 (all at (j+1)th time line) are computed as follows:

 2 
j +1
j +1  ∂S
j +1
∂S
j +1  
 βQ B   − 1 + ∆x f i + 1 / 2 e i +1 / 2  
j +1 + g A i + ∆x +
∂M ∂M ∂Ai   A2   i +1 / 2
∂hi j +1 i +1 / 2
∂hi j +1  
j +1
= = θ  i   
∂hi ∂Aij +1 ∂hi j +1 

1
2
( j +1 j +1
gB i j +1 hij++11 − hij +1 + S f i +1 / 2 ∆x i +1 / 2 + S e i +1 / 2 ∆x i +1 / 2 )



where
j +1
∂S f i +1 / 2  2 ∆n j +1 5 B j +1 2 ∆B i 
j +1
i
S f i +1 / 2 
j +1 i
= − +
∂hi j +1  n ij +1 ∆hij +1 3 A ij +1 3B ij +1 ∆hi j +1 
 

j +1
∂S e i +1 / 2
=
(K )
ec i
j +1
(Q 2 ) ij +1 B ij +1
∂hi j +1 g∆x i +1 / 2 ( A 3 ) ij +1
32
   2 
j +1
j +1 
j +1
∂S f i +1 / 2 ∂S e i +1 / 2 
j +1 
−  βQ B   
j +1 + g A i 1 + ∆x i +1 / 2 + ∆x i +1 / 2 +
∂M ∂M ∂Ai +1   A2   j +1
∂hi +1 j +1 
∂hi +1  
= =θ   i +1  
∂hij++11 ∂Aij++11 ∂hi j++11 

1
2
gB i
j +1
+ 1 h i
j +1
+1 − h (
i
j +1
+ S f i
j +1
+ 1 / 2 ∆x i + 1 / 2 + S e i
j +1
+1 / 2 ∆x i +1 / 2 )



where

j +1
∂S f i +1 / 2  2 ∆n j +1 5 B j +1 2 ∆B i +1 
j +1
i
S f i +1 / 2 
j +1 i +1
= − +
∂hij++11  n ij +1 ∆hij++11 3 A ij +1 3B ij +1 ∆hi j++11 
 

j +1
∂S e i +1 / 2
=
(K )
ec i
j +1
(Q 2 ) ij++11 B ij++11
∂hi j++11 g∆x i +1 / 2 ( A 3 ) ij++11

33
 ∂M ∆xi +1/ 2   2 β Q  j +1 j +1 
j +1
∂S f i +1/ 2 ∂Se i +1/ 2  
j +1

= + θ −   + g Ai  ∆xi +1/ 2 + ∆xi +1/ 2 

∂Qi j +1
2∆t j
  A i  ∂Qi j +1
∂Qi  
j +1


where

j +1
∂S f i +1 / 2  1 ∆ n j +1 1 
i
S f i +1 / 2 
j +1
= +
∂Q i j +1  n ij +1 ∆Q j +1 Q j +1 
 i i 

j +1
∂S e i +1 / 2
=
( )
− K ec i
j +1
Q i j +1
∂Q i j +1 g∆x i +1 / 2 ( A 2 ) ij +1

34
 ∂M ∆x i +1 / 2  2 βQ  j +1 j +1 
j +1
∂S f i +1 / 2
j +1
∂S e i +1 / 2 
=
j +1
s m i +1 / 2 + θ   + g Ai  ∆x + ∆x i +1 / 2 
∂Q i j++11 2∆t j  A  i +1  i +1 / 2 ∂Q j +1 ∂Q j +1 
  i +1 i +1 

where

j +1
∂S f i +1 / 2  1 ∆ n j +1 1 
i
= S f i +1 / 2  j +1 + j +1 
j +1

∂Q ij++11  n i ∆Q j +1
Q i 
 i

j +1
∂S e i +1 / 2
=
( )
− K ec
j +1
i Q i j++11
∂Q i j++11 g∆x i +1 / 2 ( A 2 ) ij++11

35
Partial Derivatives of the Upstream Boundary Conditions:

If a discharge hydrograph is used as the upstream boundary condition, the

partial derivatives of Jacobian take the following form

∂UB
=1
∂Q1j +1

∂∂UB
UB
=0
∂h1j +1

However, if a stage hydrograph is used as the upstream boundary condition,

then the partial derivatives are computed as follows

∂UB
=0
∂Q1j +1

∂UB
=1
∂h1j +1
36
Partial Derivatives of the Downstream Boundary Conditions:

If a discharge hydrograph is used as the downstream boundary condition, the

partial derivatives take the following form

∂DB
=1
∂Q1j +1

∂DB
=0
∂h1j +1

and the following forms for stage hydrograph

∂DB
=0
∂Q1j +1

∂DB
=1
∂h1j +1

37
For the single-valued rating curve, the partial derivatives can be written as

∂DB
j +1
=1
∂Q1

∂DB  Qk +1 − Qk 
j +1
= − 
∂h1  hk +1 − hk 

If a looped rating curve is used, the partial derivatives are written as

j +1
∂DB  1 ∆n 
= 1 + Q '  
∂Q1j +1  n ∆Q N

j +1
∂DB  1 ∆n 5 B 2 ∆B 
= Q '  − + 
∂h1j +1  n ∆ h 3 A 3 B ∆h N

38
Finally, if a critical section is used at the downstream boundary, the partial
derivatives are written as

∂DB
=1
∂Q1j +1

j +1
∂DB  3 B 1 ∆B 
j +1
= Qcr  − + 
∂h1  2 A 2 B ∆h  N

39
 Solution of Finite Difference Equations
The system of non-linear equations can now be written in functional form in terms
of the unknowns Q and h at time level (j+1)

UB(Q1 , h1 ) = 0 This system of 2N non-linear

C1 (Q1 , h1 , Q 2 , h 2 ) = 0 equations with 2N unknowns is
M1 (Q1 , h1 , Q 2 , h 2 ) = 0 solved by Newton-Raphson
C 2 (Q 2 , h 2 , Q3 , h3 ) = 0 technique for each time step.
M 2 (Q 2 , h 2 , Q3 , h3 ) = 0
⋅ The computational procedure
⋅ requires a set of trial values to
C i (Q i , hi , Q i +1 , hi +1 ) = 0 initiate the iterative procedure.
M i (Q i , hi , Q i +1 , hi +1 ) = 0 These trial values of Q and h can
⋅ be the initial conditions (if j=1) or
⋅ the values from the previous
C N − 2 (QN − 2 , hN − 2 , QN −1 hN −1 ) = 0 time step (if j>1). Solving this set
M N − 2 (QN − 2 , hN − 2 , QN −1 hN −1 ) = 0 of equations results in 2N
C N −1 (QN −1 , hN −1 , QN , hN ) = 0
residuals. For kth iteration, these
M N −1 (QN −1 , hN −1 , QN , hN ) = 0
residuals can be expressed as 40
DB(QN , hN ) = 0
UB(Q1k , h1k ) = RUB k
C1 (Q1k , h1k , Q 2k , h 2k ) = RC1k
M1 (Q1k , h1k , Q 2k , h 2k ) = RM1k
C 2 (Q 2k , h 2k , Q3k , h3k ) = RC 2k
M 2 (Q 2k , h 2k , Q3k , h3k ) = RM 2k
⋅
⋅
C i (Q ik , hik , Q ik+1 , hik+1 ) = RC ik
M i (Q ik , hik , Q ik+1 , hik+1 ) = RM ik
⋅
⋅
C N − 2 (QNk − 2 , hNk − 2 , QNk −1 , hNk −1 ) = RC Nk − 2
M N − 2 (QNk − 2 , hNk − 2 , QNk −1 , hNk −1 ) = RM Nk − 2
C N −1 (QNk −1 , hNk −1 , QNk , hNk ) = RC Nk −1
M N −1 (QNk −1 , hNk −1 , QNk , hNk ) = RM Nk −1
DB(QNk , hNk ) = RDB k 41
 The solution is tracked by finding the values of the unknowns Q and h so
that the residuals given above are forced to zero or very close to zero.

∂UB ∂UB
δQ1 + δh1 = −RUB k
∂Q1 ∂h1
∂C1 ∂C ∂C1 ∂C
δQ1 + 1 δh1 + δQ 2 + 1 δh 2 = −RC1k
∂Q1 ∂h1 ∂Q 2 ∂h 2 Applying the Newton-
∂M1 ∂M1 ∂M1 ∂M1 Raphson technique by
δQ1 + δh1 + δQ 2 + δh 2 = −RM1k
∂Q1 ∂h1 ∂Q 2 ∂h 2 substituting in the
⋅
⋅
necessary parameters, the
∂C i ∂C i ∂C i ∂C i system of equations takes
δQ i + δh i + δQ i +1 + δhi +1 = −RC ik
∂Q i ∂hi ∂Q i +1 ∂hi +1 the form given here,
∂M i
δQ i +
∂M i
δh i +
∂M i
δQ i +1 +
∂M i
δhi +1 = −RM ik
which can easily be
∂Q i ∂hi ∂Q i +1 ∂hi +1 transformed into the
⋅ matrix form.
⋅
∂C N −1 ∂C N −1 ∂C N −1 ∂C N −1
δQN −1 + δhN −1 + δQ N + δhN = −RC Nk −1
∂QN −1 ∂hN −1 ∂QN ∂hN
∂M N −1 ∂M N −1 ∂M N −1 ∂M N −1
δQN −1 + δhN −1 + δQ N + δhN = −RM Nk −1
∂QN −1 ∂hN −1 ∂QN ∂hN
∂DBN ∂DBN
δQ N + δhN = −RDB k
∂QN ∂hN 42
∂UB ∂UB
0 0 0 ⋅ 0 0 0 0 δQ1 − RUB
∂Q1 ∂h1
∂C1 ∂C1 ∂C1 ∂C1
0 ⋅ 0 0 0 0 δh1 − RC1
∂Q1 ∂h1 ∂Q 2 ∂h 2
∂M1 ∂M1 ∂M1 ∂M1
0 ⋅ 0 0 0 0 δQ 2 − RM1
∂Q1 ∂h1 ∂Q 2 ∂h 2
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
∂C i ∂C i ∂C i ∂C i
0 0 ⋅ ⋅ 0 0 δQ i − RC i
∂Q i ∂hi ∂Q i +1 ∂hi +1
∂M i ∂M i ∂M i ∂M i
0 0 ⋅ ⋅ 0 0 δh i = − RM i
∂Q i ∂hi ∂Q i +1 ∂hi +1
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅ ⋅
∂C N −1 ∂C N −1 ∂C N −1 ∂C N −1
0 0 0 0 0 ⋅ δhN −1 − RC N −1
∂Q N −1 ∂hN −1 ∂QN ∂hN
∂M N −1 ∂M N −1 ∂M N −1 ∂M N −1
0 0 0 0 0 ⋅ δQ N − RM N −1
∂Q N −1 ∂hN −1 ∂QN ∂hN
∂DBN ∂DBN
0 0 0 0 0 ⋅ 0 0 δh N − RDB
∂QN ∂hN

43
A matrix elimination routine such as Gaussian elimination or LU decomposition
can be used to solve such a system given above. However, Fread (1971)
considered the sparse character of Jacobian coefficient matrix with a maximum
of 4 elements in a row and developed an efficient algorithm to solve such a
banded matrix. No matter which scheme is used to solve the matrix, the final
step is to compute the corrected values of the unknowns for the next iteration
according to the following formulae

Qik +1 = Qik + δQi

hik +1 = hik + δhi

After solving for the new values of iteration, the convergence criteria are
checked for all Q and h values in the vector as mentioned before

44
 Initialize all arrays
Flow Chart
Input global data

Input reach cross Read boundary conditions

section data and lateral inflows

Input lateral inflow Compute trial values of unknown

location and flow data variables (Q1k, h1k, …QNk, hNk) from
initial conditions or previous time steps

Compute initial conditions Solve partial derivatives in Jacobian matrix

using trial values of unknown variables
Advance to next time step:
j=j+1
Compute the residuals RUBk, RC1k,
RM1k, …, RCN-1k, RMN-1k, RDBk
Last NO
Time step Advance to next Solve system of equations for δQi and
iteration in Newton- δhi using a matrix solver
YES
Raphson method:
k=k+1 Compute values of Qik+1 and hik+1
Stop NO
and find the updated trial values
of unknown variables
Check convergence
for Q and h
YES

Write results 45
Reference
Gündüz, O. (2004). Coupled Flow and Contaminant Transport Modeling in
Large Watersheds. Ph.D. Dissertation, School of Civil and Environmental
Engineering, Georgia Institute of Technology, Atlanta, GA, USA, 467p

46