Xem mẫu
- CHAPTER 8
Non-linear Deterministic Models
8.1 NON-LINEARITY IN HYDROLOGY
If we examine the basic physical equations governing the various hydrologic
processes, we find that these equations (and hence the processes they represent)
are non-linear. Consequently, we face the distinct possibility that all of the approaches of
linear analysis discussed in Chapters 4, 5, 6 and 7 may be irrelevant to real hydrologic
problems, save as a prelude to the development of non-linear methods. Accordingly, in the
Non-linear
present chapter we take up this question of non-linearity and ask ourselves whether we
methods
can determine under what circumstances the effects of non-linearity will be most marked
and also whether we can adapt the methods of linear analysis described in previous
chapters to the non-linear case. While knowledge of linear methods of analysis is
valuable in such an examination, we must avoid the tendency to carry over into non-linear
analysis certain preconceptions, which are valid only for the linear case.
The basic equations for the one-dimensional analysis of unsteady flow in open
channels are the continuity equation and the equation for the conservation of linear
momentum. The continuity equation can be written as:
Q A
(8.1)
r ( x, t )
x t
where Q is the discharge, A the area of flow, and r(x, t) the rate of lateral inflow. The
above equation is a linear one and consequently poses no difficulties for us in this
regard. The second equation used in the one-dimensional analysis of unsteady free-
surface flow is that based on the conservation of linear momentum, which reads
y u u 1 u u
(8.2)
S0 S f r ( x, t )
x g x g t gy
where y is the depth of flow, u is the mean velocity, S0 is the bottom slope and Sf is
the friction slope. This dynamic equation is highly nonlinear. Consequently, it is not
possible to obtain closed-form solutions for problems governed by equations (8.1)
and (8.2). The extent of the non-linearity can be appreciated if we examine the
special case of discharge in an infinitely wide channel with Chezy friction, in which
case the continuity equation takes the form
q y
(8.3)
r ( x, t )
x t
where q = uy is the where q = u y is the discharge per unit width; and the momentum
equation takes the form
u2
y u u 1 u u
(8.4a)
S0 2 r ( x, t )
x g x g t C y gy
which appears to be non-linear in only three of its six terms. However, if we multiply
through by g y, f ive of the six terms of the equation are seen to be non-linear. If, in
- 143 -
- addition, we express u in terms of q and y, which are the dependent variables in the
linear continuity equation, we obtain
y q q g
( gy 3 q 2 ) 2qy y 2 S0 gy 3 2 q 2 (8.4b)
x x t C
in which every term is seen to be highly non-linear (see Appendix D).
On the basis of the above equations we would expect such processes as flood
Unsteady flow with a
routing, which is a case of unsteady flow with a free surface, to be characterised by
free surface
highly non-linear behaviour. However, practically all the classical methods of flood
routing commonly used in applied hydrology are linear methods. In contrast most of
the methods used in applied hydrology to analyse overland flow (which is another
case of unsteady free surface flow) are non-linear in character.
The basic equations for sub-surface flow are also non-linear in form. For the
case of one-dimensional unsteady vertical flow in the unsaturated zone, the basic
equation (often known as Richards equation) was developed in Section 7.1 above
and given as equation (7.11) on page 129. This equation reads, in its diffusivity form,
as
c c
(8.5)
D (c) z z [ K (c)] t
z
where c is the moisture content, K(c) the hydraulic conductivity of the unsaturated
soil at a moisture content c, and D(c) the hydraulic diffusivity of the unsaturated soil
Sub-surface flow at a moisture content c. Since the hydraulic con- ductivity is usually a non-linear
function of c and the hydraulic diffusivity varies with c, equation (8.5) above is clearly
a non-linear equation and so represents a non-linear process.
The equation for horizontal flow in the saturated zone was also derived in
Section 7.1 as equation (7.6) on page 128 and is usually known as the Boussinesq
Boussinesq equation
equation:
h h
(8.6)
K h r ( x, t ) f
x x t
where h is the height of the water table above a horizontal impervious layer, K is the
saturated hydraulic conductivity, f the drainable porosity of the soil, and r(x, t) is the
rate of recharge at the water table. This equation is clearly non-linear, because of
the first term on the left-hand side of the equation, and also because the usual
assumption that f may be taken as constant, is open to serious doubt.
There are so many uncertainties in the derivation of unit hydrographs that
reliable and significant data on the existence of non-linearity in surface response is
not readily available. However, there have been some interesting results which have
Non-linearity of been published in the literature and which need to be taken into account in any
catchment response attempt to evaluate the non-linearity of catchment response. Minshall (1960) derived
unit hydrographs for a small catchment of 27 acres in Illinois for five storms whose
average intensity varied from 0.95 inches per hour to 4.75 inches per hour. If this small
catchment acted in a linear fashion, the unit hydrographs should have been essentially
the same for each of the five storms.
Actually, as shown in Figure 8.1 the peak of the unit hydrograph showed a more
than threefold variation, being higher for the greater rainfall intensity. The time to
- 144 -
- peak also showed a threefold variation, being smaller for the larger rainfall
intensities. Minshall's results are a clear indication of non-linear behaviour.
Amorocho and Orlob (1961) and Amorocho and Brandsetter (1971) published
data for a very small la basin, in which the artificial rainfall was carefully controlled
and the runoff accurately measured (see Figure 8.2). The test basin consisted of a
thin layer of gravel placed over an impervious surface. The results for varying rates
of input showed a clearly non-linear response as indicated by the set of three
experimental results, in which the cumulative outflows are not proportional to the
corresponding cumulative inflows (see Figure 8.3). Both Minshall's and Arriorocho's
data will be discussed later in Section 8.5, which deals with the concept of spatially
uniform non-linearity.
- 145 -
- Ishihara and Takasao (1963), in a paper on the applicability of unit hydrograph
methods, showed results for the Yura river basin at Ono, which is a river basin of
346 square kilometres. Figure 8.4 presents the relationship, which they obtained,
between the mean rainfall intensity and the time of rise between the start of the
equivalent mean rainfall and the peak of the flow hydrograph. They interpreted the
events for a mean rainfall intensity above 10 millimetres per hour as representing
essentially surface runoff, and the events for a mean rainfall intensity of less than 8
millimetres per hour as representing essentially sub-surface runoff. It is noteworthy,
that their results, as presented in Figure 8.4, show remarkably little variation in the
- 146 -
- time to peak for rainfall intensities from 4 millimetres per hour to 18 millimetres per
hour, but show a distinct variation for smaller rainfall intensities. The results from
Ishihara and Takasao also indicate a linear relationship between peak runoff and
mean intensity of rainfall at higher values of rainfall intensity. These two results
taken together would indicate that for conditions similar to those in the Yura river
catchment, the unit hydrograph approach might be reliable for high intensities, but
not for low ones.
Pilgrim (1966) measured the time of particle travel in a catchment area of 96 square
Time of particle travel
miles by means of radioactive tracers. His results showing the relationship between
time of travel and level of discharge are presented in Figure 8.5. As in the case of the
Japanese result, we see here an essential constancy of time of travel at higher rates of
discharge and a tendency for the time of travel to be inversely proportional to the discharge
at lower discharges.
We conclude from this brief summary, both from the basic equation of physical
hydrology and from experimental data, that there are sufficient indications of non-linearity to
justify an investigation of the extent to which non-linearity affects the techniques commonly
used in applied hydrology.
There are many possible approaches to the analysis of non-linear processes and
systems. One approach is to analyse each input-output event as if it were linear and
then to examine the effects of the level of input on the results obtained. Linearisation
can be applied to all three basic approaches used in hydrology: black-box analysis,
conceptual models, or solution of the basic equations. The Linearisation approach is
discussed in Section 8.3 below.
- 147 -
- A second line of approach to non-linear systems is to accept the non-linearity and to
attempt to develop a non-linear method of black-box analysis, which would be a
Non-linear method of
generalisation of linear black-box analysis as discussed in Chapter 4. This method is
black-box analysis
summarised in Section 8.4. Still another approach would be to attempt to find simple non-linear
conceptual models, which would simulate the operation of non-linear systems with the
same degree of accuracy as achieved by the simple linear conceptual models presented in
Chapter 5. This is the subject of Sections 8.2 and 8.5. A final approach would be to
accept the full complexity of the complete non-linear equation and to seek solutions by
numerical methods. This last approach is outside the scope of this book.
8.2 THE PROBLEM OF OVERLAND FLOW
Overland flow is an interesting example of a hydrologic process, which
appears to require a non-linear method of solution. It would appear that because it
occurs early in the runoff cycle, the inherent non-linearity of the process is not
Inherent non-
dampened out in any way as appears to occur to some extent in the question of catchment
linearity
runoff. A physical picture of overland flow is shown in Figures 8.6 and 8.7 together with a few
of the classical experimental results of Izzard (1946).
For the two-dimensional problem of lateral inflow the equation of continuity is written as
q y
(8.7)
r ( x, t )
x t
where q is the rate of overland flow per unit width, y is the depth of overland flow and r is the
rate of lateral inflow per unit area. The equation for the conservation of linear momentum is
written as (from equation 8.4a)
y u u 1 u g
(8.8)
S0 S f 2 r ( x , t )
x g x g t gy
- 148 -
- where u is the velocity of overland flow, S0 is the slope of the plane and Sf is the
friction slope.
The classical problem of over flow is the particular case where the lateral inflow is
uniform along the plane and takes the form of a unit step function. There are several parts to
the complete solution of this problem.
Firstly. there is the steady-state problem of determining the water surface profile
when the outflow at the downstream end of the plane increases sufficiently to balance the
inflow over the surface of the plane. Secondly, there is the problem of determining
the rising hydrograph of outflow before this equilibrium state is approached for the special
case of the step function input. If the process were a linear one, the solution of this
second problem (i.e. the determination of the step function response) would be sufficient
to characterise the response of the system and the outflow hydrograph, for any other
inflow pattern, could be calculated from it.
However, since the problem is inherently non-linear, the principle of superposition
cannot be used and each case of inflow must be treated on its merits. The third part
of the classical problem is that of determining the recession from the equilibrium condition
after the cessation of long continued inflow. Further problems that must be investigated are
the nature of the recession when the inflow ceases before equilibrium is reached, the case
where there is a sudden increase from one uniform rate of inflow to a second higher
uniform rate of inflow, and the case when a uniform rate of inflow is suddenly changed
to a second rate of uniform inflow, which is smaller than the first. The above
problems can be solved by numerical methods (Liggett and Woolhiser, 1967; Woolhiser,
1977) but such methods are outside the scope of the present discussion. Here we will be
concerned with simpler approaches to the problem and with attempts to find a simple
mathematical simulation or a simple conceptual model.
The first approach to the solution of overland flow in classical hydrology was based on
the replacement of the dynamic equation, given by equation (8.8) above, by an
- 149 -
- a ssumed relationship between the outflow at the downstream end of the plane and the
volume of storage on the surface of the plane. Because this method was first proposed by
Horton (1938) for overland flow on natural catchments and subsequently used by Izzard
(1946) for impermeable plane surfaces, it may be referred to as the Horton-Izzard
approach. It had been noted by hydrologists that for equilibrium conditions on experimental
plots, the relationship between the equilibrium runoff and the equilibrium storage could
be approximated by a power relationship. Such a relationship would indicate that the
outflow and the storage would be connected as follows
q( L, te ) qe a( S2 )c (8.9)
where qe, is the equilibrium discharge at the downstream end (x = L) after a lapse of time
te sufficient for equilibrium to occur, S e is the total surface storage at equilibrium
conditions, and a and c are parameters which could be determined from experimental data
by means of a log-log plot.
In the Horton-Izzard approach to the overland flow problem, the assumption is made
that such a power relationship holds, not only at equilibrium, but also at any time during
the period of unsteady flow, either during the rising hydrograph or during the recession.
Recession This assumption can be written as
q(L,t) = QL = qL = a(S)c (8.10)
where qL. is the discharge at the downstream end at any time t and S is the
corresponding storage on the surface of the plane of overland flow at the same time.
Izzard (1946) illustrates the nature of this approximate relationship on Figure 8.7 for two of
the experimental cases examined.
The equation of continuity in its lumped form for the whole plane can be written for
the case of constant input r as
dS
(8.11)
rL qL
dt
which is in reality an integrated form of equation (8.7). If the HortonIzzard assumption
given in equation (8.10) is made, then this equation of continuity can be written as
dS
qe aS c (8.12a)
dt
or in more convenient form as
dS
(8.12b)
dt
qe aS c
Equation (8. 12) can be integrated to give the time as a function of the storage
Se d (S / Se
(8.13)
qe 1 ( S / Se )c
t
Equation (8.13) can be solved analytically for values of c = 1 (linear case), c = 2, c = 3
or c = 4 and also for values of c which ratios of these integral values i.e. for c = 3/2 or
c = 4/3. It is interesting to note that the integral in equation (8.13) occurs also in the case
of non-uniform flow in an open channel (Bakmeteff, 1932) and in the case of the
relationship between actual and potential evapo-transpiration (Bagrov. 1953; Dooge, 1991).
Horton (1938) solved the equation of the rising hydrograph equation (8.13) for the
Rising hydrograph
case of c = 2, which he described as “mixed flow” since the value of c is intermediate
- 150 -
- between the value of 5/3 for turbulent flow and the value of 3 for laminar flow. If we define a time
parameter Ke by
Se 1 1
(8.14)
Ke c 1 1/ c ( c 1) / c
qe aSe a qe
the integration of equation (8.13) with c = 2 for zero initial condition gives
1 ( S / Se
2t
(8.15a)
loge
Ke 1 ( S / S e
which gives the time as a function of the storage. This can readily be rearranged to
give the storage as an explicit function of time:
t
S
(8.15b)
tanh
Se Ke
t
q
tanh 2 (8.15c)
qe Ke
Since the system is non-linear, the time parameter K e will depend on the
intensity of the inflow. Horton's equation as given above has been widely used in the
design of airport drainage systems in the period since he proposed it almost forty
year ago. The solution of equation (8.13) for the case of c = 3 i.e. for laminar flow was
presented by Izzard (1944) in the form of a dimensionless rising hydrograph. Izzard,
who appears to have followed the theoretical analysis of Keulegan (1944), uses as his
time parameter a time to virtual equilibrium, which is defined as twice the time parameter
given in equation (8.14) above.
For recession from equilibrium, the recharge in equation (8.11) becomes zero,
and the substitution for qL from equation (8.10) and a slight rearrangement gives us
the simple differential equation
dS
(8.16)
adt c
S
which can be solved for any value of c. For the linear case (c = 1) the storage and
hence the outflow shows an exponential decline. For all other values of e the recession of
equilibrium is given by
q 1
(8.17)
qe [1 (c 1)(t / K e ]c / ( c 1)
where q is the outflow at a time t after the start of recession i.e. after the cessation of
Partial recession
inflow. The special case of equation (8.17) for laminar flow (i.e. for c =3) was given by
Izzard (1944). Figure 8.6 shows the rising hydrograph and the recession for the Horton-
lzzard solution with c = 2 for a duration of inflow equal to twice the time parameter defined by
equation (8.14). The double curvature of the rising hydrograph is characteristic of the
shape of 1.e rising hydrograph for the Horton-Izzard solution for all values of c other
than c = 1. If the duration of inflow D is less than the time required to reach virtual
equilibrium, we get a partial recession from the value of the outflow qD which has been
reached at the end of inflow. It can be shown that the partial recession curve has the
same shape as the recession from equilibrium given by equation (8.17), except that
the recession curve from partial equilibrium starts at the point on the curve defined
- 151 -
- b y the appropriate value of qD/q i.e. the end of inflow does not correspond to the
time origin in equation (8.17).
If there is a change to a new rate of uniform inflow during the rising hydrograph
one of two cases can occur. If the new rate of inflow is higher than the rate of outflow
when the change occurs, the remainder of the rising hydrograph will follow the same
dimensionless curve as before; but since qe is equal to the inflow at equilibrium, the
value of q/qe will change as soon as the rate of inflow changes. Such a case is shown in
Figure 8.7. (run No. 138). If the new rate of inflow is less than the outflow at the time when
the change occurs, the hydrograph will correspond to the general falling hydrograph for the
case in question. An example of such a falling hydrograph taken from Izzard (1944) is
shown on Figure 8.7b (run No. 143). For the case of c = 2, the equation for the falling
hydro-graph can also be obtained. For the case of recession to equilibrium, the
integration given in equation (8.15.) above is not valid, because it results in the logarithm of
a negative number. However, it can be shown that the appropriate integration in this case is
S / Se 1
2t
(8.18a)
loge
Ke S / Se 1
which can be rearranged to give the discharge as an explicit function of time:
1
q
coth 2 (8.18b)
qc Kc
In practice, it would often be convenient to take advantage of the relationship
between tanh and coth and to write
t
qe
tanh 2 (8.18c)
q Ke
so that one could use either a single graph, or a single computer routine, to evaluate
the value of q/qe for the hydrograph, or the value of q/qe for the falling hydrograph.
The master recession curve defined by equation (8.18) applies to all cases
Master recession
where there is a uniform rate of lateral inflow and an initial storage on the plane
curve
which is higher than the equilibrium storage for that particular inflow. There will be a
similar master recession curve for any other value of the index of non-linearity c. The
only case to which this master recession curve will not apply is when the inflow drops to
zero. In the later case, the governing equation will be equation (8.16) and if allowance is
made for the discharge q0 at the cessation of inflow, we will have the relationship
t
[ qe / q]( c 1)/ c [qe / q0 ]( c 1)/ c (c 1) (8.19a)
Ke
which can be rearranged to give the discharge as a function of time
q 1
(8.19b)
qe [ q / q ]c /( c 1) (c 1) / (t / K ) c / ( c 1)
e 0 e
of which equation (8.17) is a special case.
The Horton-Izzard approach as described above clearly involves the use of a simple
conceptual model. In fact, the whole approach is based on treating the overland flow as a
lumped non-linear system, which can be represented by a single non-linear reservoir
whose operation is described by equation (8.10) above. The Horton-Izzard solution for the
- 152 -
- various cases arises from the application of the equation of continuity to such a single non-
linear reservoir. Even though this conceptual model is extremely simple in form, the fact that
it is non-linear makes it less easy to handle than some of the apparently more complex
conceptual models used to simulate linear or linearised systems. Thus the impulse
response for such a system no longer characterises the operation of the system because
the output in any individual case will depend on the form and intensity of the input.
The cumulants of the impulse response can no longer be added to the cumulants of the
input to obtain the cumulants of the output. The solution for a step function input (described
above as the rising hydrograph solution) cannot be used to predict the output for a varying
pattern of input. However, because the method involves only a single non-linear reservoir
the output for a complex pattern of input can be obtained by dividing the input into intervals of
uniform inflow and applying to each interval one or other of the three master solutions (the
rising hydrograph, the recession to equilibrium, the recession to zero) for the chosen
value of the index of non-linearity c.
The second simplified approach to the problem of overland flow, which appears in the
Kinematic wave
literature, is the kinematic wave solution (Henderson and Wooding, 1964). This also
solution
involves a power relationship between discharge and depth. In this case the relationship
is not a lumped one for the whole system, but a distributed relationship between the
discharge and the depth at each point. This basic relationship can be written as
q(x, t) = b[y(x, t)]c (8.20)
where q(x, t) is the discharge per unit width, y(x, t) is the depth of the flow, and a and
c are parameters of the system. The above basic assumption corresponds to neglecting all
of the terms in the equation for the conservation of linear momentum given by
equation (8.8) above in comparison with the terms for the bottom slope and the friction
slope, so that we can write
S0 – Sf = 0 (8.21)
value of c in equation (8.20) will be 3/2, whereas if it is taken according to the
Manning formula, the value of c will be 5/3.
As might be expected, the solution of the basic problems of overland flow for this
distributed model cannot be derived as easily as the solution for the lumped model
represented by the Horton-Izzard approach. The problem is to solve equations (8.7)
Method of and (8.20) subject to the appropriate boundary conditions. The method of
characteristics characteristics can be used (Eagleson, 1969; Woolhiser, 1977) to obtain the
solution, but the details of the derivation are outside the scope of the present chapter.
For the case of the rising hydrograph it can be shown that the solution is given by
c
t
q
(8.22)
qe tk
from t = 0 to t = tk where the latter time parameter is given by
1/ c
ye 1 qe
(8.23)
tk
r rb
where qe is the equilibrium discharge at the downstream end of the overland slope and ye
is the depth of flow at the downstream end for equilibrium conditions and small b is
the parameter in equation (8.20) above. For times greater than the kinematic time
- 153 -
- parameter, the discharge is equal to the equilibrium discharge. The recession from full
equilibrium for the kinematic wave solution can be shown to be
t
1 q / qe
(8.24)
c
( q / qe )( c 1)/ c tk
where t is the time since the cessation of inflow.
Care should be taken to distinguish the kinematic time parameter tk from the time
Horton-Izzard model parameter for the Horton-Izzard model Ke defined by equation (8.14) above. The relationship
between the two time parameters can be demonstrated as follows. The storage at
equilibrium is given by
L
Se ye ( x) dx (8.25)
0
where ye, is the equilibrium depth at a distance x from the upstream end. Since equation
(8.20) holds at equilibrium, this can be written
1/ c
L q ( x)
Se e (8.26)
dx
b
0
At equilibrium we have
(8.27)
qe ( x) rx
so that equation (8.26) can be written
1/ c
rx
L
(8.28)
Se dx
b
0
The latter equation can be integrated to give
1/ c
c r ( c 1) / c
(8.29a)
Se ( L)
c 1 b
which can be regrouped as
1/ c
c rL
(8.29b)
Se L
c 1 b
which by virtue of equations (8.23) and (8.27) is
c rL
(8.29c)
Se tk qe
c 1 b
Recalling the definition of the Horton-Izzard time parameter as given by
equation (8.14) we see that
Se c
(8.30)
Ke tk
qe c 1
which is the required relationship between the parameter for the Horton-Inard approach
and the kinematic time parameter.
The only sound basis for evaluating the above two models is by comparing the results
with solutions of the complete non-linear equation. The numerical solution of the
overland Sow problem has been tackled by Woolhiser and Liggett (1967). They reduce the
equations of continuity and momentum to dimensionless forms by expressing variables in
terms of the normal depth and the normal velocity at the downstream end of the plane for
the equilibrium discharge. When this is done, the continuity equation given above by
equation (8.7) becomes
- 154 -
- q y
(8.31)
1
x t
where q is the ratio of the discharge to the equilibrium discharge at the downstream
end and y is the ratio of the depth to the equilibrium depth at the downstream end of the
plane. The equation for linear momentum becomes
u u u S0 L u 2
y
F02 u (8.32)
1
x x t y y0 y
where u is the ratio of the velocity to the equilibrium velocity at the downstream end and
F0 is the Froude number for the equilibrium flow at normal depth. In both equations
(8.31) and (8.32) the independent variable x is the ratio of the distance from the upstream
end to the total length of overland flow. The independent variable t is the ratio of the
time to the characteristic time t0 obtained by dividing the length of the plane by the
velocity at the downstream end at equilibrium. This reference time to can be shown to be
equal to the kinematic time parameter defined by
equation (8.23). It appears from these dimensionless equations that there are only
two parameters governing the flow, namely the Froude number for normal flow at
equilibrium discharge F0 and the dimensionless length factor S0L/y0. W hen the ratio of
Dimensionless length
these two parameters defined by
S0 L
(8.33)
K
F02 y0
is appreciably greater than 1, all of the terms on the left-hand side equation (8.32), except the
first, can be neglected, thus reducing the momentum equation to convective-diffusion form. If
in addition the dimensionless length factor S0L/y0 is also appreciably greater than 1, the first
term on the left hand side of equation (8.32) can be neglected as well and the kinematic
wave approximation results.
Woolhiser and Liggett (1967) found that for values of K as defined by equation (8.33)
greater than 10, the kinematic wave solution was a good approximation to the rising
hydrograph, but that for values of K smaller than 10 the approximation was a poor one.
Figure 8.8 shows a solution obtained by Ligett and Woolhiser for a value of the Froude
number equal to 1.0 and a value of the parameter K equal to 3.0.
- 155 -
- 8.3 LINEARISATION OF NON-LINEAR SYSTEMS
If we know the input and the output for a non-linear system, there is nothing to
prevent us using the techniques discussed in Chapter 4 to obtain the apparent unit
hydrograph of the system. If we convolute the “unit hydrograph" thus derived with
the input, we should obtain an accurate reconstitution of the output. The difficulty is
that for the non-linear system, we cannot use such a derived "unit hydrograph" to
predict by direct convolution the output for any other input. Only in the case of a
linear time-invariant system is the unit hydrograph obtained by the solution of the
convolution equation independent of the particular event from which it is derived.
However, if we can derive such an apparent unit hydrograph from a number of
Apparent unit
storms with markedly different inputs of effective precipitation, we might be able to
hydrograph
make progress. We seek information on the manner in which key parameters (the
time to peak and the peak discharge, or the values of the lower moments) of the
apparent unit hydrograph vary with the characteristics of the input. And, we may be
able to predict the shape of the apparent unit hydrograph for a given input.
The approach outlined in the last paragraph can also be used where the apparent
unit hydrograph is assumed to be represented by a simple conceptual model. In this case
the parameters of the model will depend not only on the properties of the catchment, but
also on the characteristics of the particular input for a given event. The approach can be
illustrated readily for the case where it is assumed that the input shape remains
constant from event to event. The only variation is in the mean value of the input x. We
then seek the apparent unit hydrograph for the ith input-output event by solving the
equation
(8.34)
yi (t ) hi (t ) * x i f (t )
where yi (t) is the output for the ith event, f (t) is the constant input shape, xi is the
mean rate of input for the ith event, and hi(t) is the apparent unit hydrograph derived by
treating the data, as if the input were transformed to the out in this event, in a linear time-
invariant fashion.
If the analysis is done on the basis of a three-parameter conceptual model, we
can express the apparent unit hydrograph hi(t) as a function of time and of the three
parameter values which give the closest fit in equation (8.34) for the particular event
i.e.
(8.35)
hi (t ) (t , ai , bi , ci )
The next step in the linearisation approach is to seek the relationship between
each of the parameters a, b, c and the level of input x. If we can obtain such a
relationship, it is possible to predict the values of the parameters a, b, c for any
given level of x and hence to determine the apparent unit hydrograph h(t) for this
level of input. If this has been successfully accomplished, the apparent unit
hydrograph corresponding to the design level of input can be convoluted with that
input to give a prediction of the output based on the linearisation of the system and
the correlation of the response parameters with levels of input.
Linearisation of the Linearisation can also be applied to the basic non-linear equations of physical
equations for open hydrology. Solutions of these linearised equations can be used to study the general
channel flow behaviour of systems but have the disadvantage that certain phenomena, which
- 156 -
- o ccur in non-linear systems, do not appear in their linearised versions. The
linearisation of the Richards equation for unsaturated flow of soil moisture and the
linearisation of the Boussinesq equation for groundwater flow have already been
discussed in Section 7.6 above. Accordingly, attention will be concentrated here on
the linearization of the equations for open channel flow.
The general non-linear equation for unsteady free surface flow in a wide
rectangular channel was mentioned in Section 8.1 above. The equation of continuity for
the case of routing an upstream inflow (i.e. r = 0) is given by
q y
(8.3)
r ( x, t )
x t
and equation for the conservation of linear momentum for Chezy friction is given by
y q q g
( gy 3 q 2 ) 2qy y 2 S0 gy 3 2 q 2 (8.4b)
x x t C
This highly non-linear equation can be linearised (Dooge and Harley, 1967) by
considering a perturbation about a steady uniform flow q0 and the following linear
equation derived for the perturbed discharge
2q 2 q 2q q 2 gS0 q
(1 F02 ) gy0 (8.36)
2u0 2 3gS0
2
x xt t x u0 t
where the coefficients depend on the assumed reference discharge q0 about which
the perturbation is taken but do not depend on the variable discharge q. Strictly
speaking equation (8.36) is only valid for perturbations small enough that the
variations of the coefficients in the non-linear equation are not sufficient to produce
large errors in the result. In fact, we may consider (8.36) as a version of the non-
linear equation in which the coefficients are frozen at the value corresponding to the
reference discharge. The use of equation (8.36) for large perturbations is equivalent
to the acceptance remains linear at all levels of flow (see Appendix D).
Since equation (8.36) is linear, it is only necessary to determine the solution for
a delta function input, since the outflow of the downstream end for any other inflow
can be obtained by convolution. The impulse response of a channel obtained in this way
Linear channel
can be referred to as the linear c hannel response (LCR) of the channel reach in
response
question (Dodge and Harley, 1967). Any of the standard mathematical techniques for
the solution of linear partial differential equations can be used in order to find the solution to
equation (8.36); but the Laplace transform method is probably the most convenient.
When the latter method is used (Dooge, 1967), the system function, i.e. the Laplace
transform of the impulse response (or LCR), is found to be (see Appendix D)
(8.37a)
H ( x, s) exp(2 ( s ).x )
where
2 e.s f as 2 b.s c (8.37b)
and a, b, c, e and f are parameters which depend on the hydraulic characteristics of the channel.
Since the above system function is of exponential form, the cumulants of the response
function can be determined by repeated differentiation of the quantity inside the brackets in the
equation and evaluated at s = 0 (see Chapter 3.5, page 51). When this is done and
the values of the parameters are substituted in the result, we obtain the first three cumulants of
the linear channel response as follows
- 157 -
- x
k1 U1' (8.38a)
1.5u0
2
2 F 2 y0 x
(8.38b)
k2 U 2 1
4 S0 x 1.5u0
3
2 3
4 F 2 F 2 y0 x
(8.38c)
k3 U 3 1 1
2 S0 x 1.5u0
3 4
Dimensionless The higher cumulants can be made dimensionless by dividing by the appropriate power of the
cumulants lag and the resulting dimensionless cumulants or shape factors are seen to be functions of the
Froude number and the dimensionless length parameter D defined by
S0 x
(8.39)
D
y0
Consequently, even if we were unable to invert the transform function given by equation
(8.10), it would be possible to determine the cumulants of the solution and to plot the
solution of the linearised equation (8.36) on a shape factor diagram. It is clear from equation
(8.38) above, that there is a general relationship between the shape factor sR and the
dimensionless length D of the form
( F )
(8.40)
sR R 1
D
where ( F ) depends only on the Froude number. In particular, for a vanishingly small Froude
number we have
2
(8.41a)
s2
3D
4
(8.41b)
s3
3D 2
so the solution for a vanishingly small Froude number would plot on (s3, s2) shape-
factor diagram as the line
2
(8.42)
s3 3 s2
Further information about the behaviour of the solution can be obtained from an
examination of the results for the first moment about the origin given by equation
(8.38.) and for the second moment about the centre given by equation (8.38b). The
lag given by equation (8.38.) indicates that, for the linearised solution, the average rate of
propagation of the flood wave is 1.5 times the velocity corresponding to the reference
discharge. This corresponds to the value for the celerity c of the flood wave in a wide
Celerity c of the
rectangular channel with Chezy friction given by
flood wave
Q q 3 1/ 2 1/ 2
(8.43)
c CS y 1.5u
A y 2 0
which was first proposed by Kleitz (1877) and Seddon (1900). The lag given by
equation (8.38.) of course depends on the choice of reference discharge q0 and is seen from
that equation to vary inversely with the cube root of the reference discharge. If we
examine equation (8.38b) for the effect on the second moment about the centre of the
reference discharge we obtain an interesting result. Since the reference discharge refers to
uniform flow, we have the following relationships between the reference discharge q 0. the
reference velocity u 0, and the reference depth y0:
- 158 -
- (8.44a)
q0 u0 y0
u0 CS0 / 2 y0 2
1 1/
(8.44b)
q0 CS0/ 2 y0 2
1 3/
(8.44c)
When the reference discharge is altered, the variation of y0 in the numerator of
equation (8.38b) is compensated for by the variation of u20 in the denominator, so
that the second moment about the centre, which measures the dispersion of the
linear channel response, is independent of the reference discharge.
Actually, it is possible to invert the system function given by equation (8.37a) to
the time domain, but the result is a complicated one. The solution in the original
coordinates ( x, t) consists of two terms
(8.45)
q( x, t ) q1 ( x, t ) q2 ( x, t )
The term q 1 is of the following form
x
(8.46)
q1 ( x, t ) t exp( px )
c1
where c1 is the characteristic wave speed given by
(8.47)
c1 u0 gy0
and where p is a parameter, which depends on channel characteristics, being given by
(2 F0 ) ( S0 )
(8.48)
p
( F0 F02 ) (2 y0 )
We see from the above result, that, for the linear solution of the case of a delta function
input at the upstream end of a channel, we will have a delta function moving downstream at the
head of the wave at the characteristic speed given by equation (8.47). However, its volume
declines exponentially at a rate determined by the parameter defined by equation
(8.48). It is clear that for conditions resembling those in most natural rivers (i.e. for a
Froude number of 0.2 or less) the decline in the volume under the head of the wave
given by equation (8.46) will be very rapid.
The second term in equation (8.45) which represents the attenuated body of the
wave is given by
x x I [2ha]
q2 ( x, t ) h exp[rt sx] 1 (8.49)
c1 c2 a
where c2 is the upstream characteristic wave speed give by
Characteristic
wave speed (8.50)
c2 u0 gy0
and where I1[2ha] is a modified Bessel function and a is given by
(8.51)
a [(t x / c1 ][t x / c2 ]
and r, s and h are parameters depending on the hydraulic properties of the channel (see
Appendix D).
- 159 -
- The use of the linearised equation is illustrated in Figure 8.9 for a channe state
rating curve given by
3/2
(8.52)
q2 50 y 0
and for an upstream inflow between t = 0 and t = 96 given by
t
(8.53)
q(0, t ) 125 75cos
48
which corresponds to the inflow used by Thomas (1934) in his classical paper on unsteady
flow in open channels. Figure 8.9 shows a modification of the flood wave for a distance up
to 500 miles.
For any linearisation of the flood routing problem, it is necessary to choose a value of
the reference discharge qc, about which the discharge is perturbed. Since this value q0 is
used to evaluate the reference depth y0 from equation (8.44a) and hence the reference
velocity u0 from equation (8.44b) and finally determines the coefficients in the linearised
equation (8.36), the choice of the reference discharge will naturally affect the result.
- 160 -
- Figure 8.10 shows the effect of reference discharge on the outflow 50 miles
downstream for the channel characterised by equation (8.52) and the inflow given by
equation (8.53). In the latter equation the inflow varies from 50 cubic feet per second per foot
width to 200 cubic feet per second per foot width. Accordingly the reference
discharge was taken at values of 50, 100, 150 and 200 cubic feet per second per
foot width. It can be seen from Figure 8.10 that the resulting flood waves are
Flood waves
displaced in time, as might be expected, since the lag given by equation (8.38a)
depends on the reference velocity u 0 and hence on the reference discharge q 0.
However, it is interesting to note that the shape of the flood wave and the value of
the peak are very similar for the various reference discharges. This would appear to
be a reflection of the fact that the second moment about the centre, which measures
the dispersion, is given by equation (8.38b) and is independent of the reference
discharge. It would appear to indicate that the effect of the third and higher
cumulants, which are affected by the reference discharge, have little influence on the
shape of the outflow. It also suggests that the use of a linear model for flood routing
might be expected to show greater accuracy in prediction of the peak outflow than in
regard to the time to peak.
Subsequent work on this problem between 1967 and 1989 has been well reviewed
by Napiorkowski (1992). His account covers the work done on such topics as (a) the
generalisation of the approach from a wide rectangular channel with Chezy friction to the
general case of any channel shape and any friction law; (b) the extension from a semi-
infinite channel with an upstream condition to a finite channel with both upstream and
downstream control; (c) general expressions for the moments and cumulants of any
order of the linearised channel response; (d) the attenuation and phase shift for a
harmonic input to a linearised channel.
The general form of the linearised equations for the case where the friction law
Friction law
takes the form Q = kAm can be derived easily. Such an assumption is almost always
adopted in the study of both steady and unsteady flow in open channels. For Chezy friction
the value of m varies from m = 3/4 for a triangular cross section to m = 3/2 for a wide
rectangular channel. For Manning friction the corresponding values are m = 5/6 for a
triangular cross section and m = 5/3 for a wide rectangular channel. The parameter m
corresponds to the ratio of the kinematic wave velocity (Uk ) to the average velocity of flow
( u ). When this generalisation is made, the linearised equation becomes
2Q 2Q 2Q Q 2 gS0 Q
(1 F02 ) y 0 (8.54)
2u0 2 2mgS0
2
x xt t x u0 t
where y0 is the value at the reference condition of the hydraulic mean depth, i.e. the
ratio of the area of flow to the surface width and u0 is the velocity corresponding to
the reference conditions used as a basis for linearisation. The cumulants given by
equation (8.38) above are then generalized (Dooge, Napiorkowski and
Strupczewski, 1987b) to
x
(8.55a)
k1
mu0
- 161 -
- 2
y x
1
k2 (1 ( m 1) 2 F02 0 (8.55b)
m S0 x mu0
2 3
y x
3
(1 (m 1) 2 (1 (m 1) F02 0 8.55b)
k3
m2 S0 x mu0
In general, it can be shown that all higher cumulants can be expressed in the
form
R 1 R
y x
k R ( m, F0 ) 0 (8.56)
S0 x mu0
Closed form expressions can be developed for the function (m, F0) for any order of
cumulant (RomanoiAe., Dooge and Kundzewicz, 1988).
The above results relate to the classical problem of flood routing, i.e. downstream
Flood routing
movement of water with a given upstream control and negligible effect of the downstream
control, i.e. a semi-infinite channel. The relative effects of the upstream and controls for a
finite channel reach were first solved for the case of the diffusion analogy which corresponds
to the assumption that the Froude (F0) number is vanishingly small (Dooge and Napiorkowski,
1984). The solution for the case of any value of m and any value of Fo was derived by the
same authors (Dooge and Napiorkowski, 1987). The downstream and upstream waves
are each reflected at the other boundary and consequently the two single terms
representing the transformation of the upstream boundary condition and the upstream
transformation of the downstream boundary condition respectively are replaced by two
corresponding infinite series of similar form with strong convergence. The rate of
convergence is the same for these infinite series and is given by
2mS0 L
(8.57)
exp 2
y0 (1 F0 )
so that the convergence is very rapid except for small values of y0 and F0. In practical
usage, a close approximation is obtained by including only the first term in each of the
series corresponding to (a) the initial downstream propagation of the upstream
boundary condition, (b) the initial upstream propagation of the downstream boundary
condition, and (c) the first reflection of the downstream propagation at the downstream
boundary (Napiorkowski and Dooge, 1988).
An extension can be made from the analysis of the response to a delta function
input to the more general case, by seeking the response to an upstream inflow of an
infinite periodic train of sinusoids, and hence the response to any function which can
Sinusoids
be broken down into a series of such functions. (Kundzewicz and Dooge, 1989). On
the basis of the continuity equation alone, it is possible to establish that the
necessary and sufficient condition for the attenuation of a flood wave is identical to
the necessary and sufficient condition for the existence of a phase shift between the
flow rate and the area of flow and for the existence of a looped rating curve. The
common condition is that either the frequency w or the wave length β is complex. If
the upstream inflow is taken as a train of non-attenuating waves, then the frequency
must be real. If the wave number is also real then there can be no attenuation with
distance and we have the limiting case of kinematic wave propagation.
- 162 -
nguon tai.lieu . vn
