...
Storage Routing links offer a choice of several different hydrologic routing methods:
- linear Muskingum routing,
- non-linear Muskingum routing (using a power function),
- variable parameter Muskingum routing
Scientific provenance
The storage (hydrologic) routing and lagged flow techniques are all well established; an excellent review of storage routing techniques is provided by Koussis (2009) and discussion on this is provided by Perumal (2010). The Muskingum technique dates back to a 1938 study on the Muskingum River in the US (McCarthy, 1938); it has been used in literally thousands of practical applications and discussed in thousands of journal papers since. Non-linear routing also has a long history (eg. Linsley et al, 1949; Laurenson, 1959; Mein et al, 1974). The variable parameter Muskingum routing is from the work of Close (1996). The calculations of reach fluxes are all based on well-established hydrologic principles that are discussed in hydrology text books and handbooks.
Dependencies
Each routing link requires a node at its upstream end to define inflows, and a node at the downstream end to process outflows as appropriate.
...
The Storage Routing link can be partitioned into divisions of equal length. Other assumptions and constraints are listed in the Table below. Further details on assumptions and constraints are given in the Theory section below.
No | Assumption/Constraint |
|---|---|
1 | The routing methods are based on lumped conceptual storages (ie. they are not spatially distributed). This is appropriate in river reaches where the friction slope is approximately equal to the bed slope. |
2 | Hence, it is assumed the cross-sectional geometry for a routing link can be represented via a single cross-section with a single representative bed/friction slope between the upper and lower ends of the link. Results from the routing methods may be very sensitive to the assumptions made about this representative cross-section. |
3 | The velocity is uniform, the water surface is horizontal across any section perpendicular to the longitudinal flow axis, and flows are gradually varying. |
4 | The relationship between storage and discharge is a monotonically increasing function so that the storage theory is applicable. Also, the same relationship is assumed to apply to both the rising and falling limbs of hydrographs. |
Storage Routing Theory
This section covers aspects of the theory relevant to routing
...
The law of conservation of mass then lets us say that for a reach:
| Equation 1 |  |
|---|
where:
S is the reach storage in (m3)
...
For streamflow routing Source uses an implicit Euler scheme which can be expressed as (Clark and Kavetski, 2010):
| Equation 2 |  |
|---|
Which is to say the flux over a time period n to n+1 is a function of the storage at the end of the time-step (Sn+1). Other hydrological models use an implicit Heun scheme (also known as the Crank-Nicholson method):
| Equation 3 |  |
|---|
In Equation 3, the flux over a time period n to n+1 is an average of the functions of storage at the start of the period and the end. Clark and Kavetski (2010) describe the implicit Euler scheme as being ‘ubiquitous in "industry-standard" engineering software’. The truncation error for the implicit Heun scheme is more sensitive to the length of the time-step than the Euler scheme. The implicit Euler scheme produces hydrologically sensible results for a wider range of parameter values than the implicit Heun scheme.
...
To make use of the continuity equation (Equation 1) it is necessary to relate the storage in the reach to the flow entering and leaving that reach. Source uses the Muskingum storage function (Koussis, 2009):
| Equation 4 |  |
|---|
where:
S is the reach storage (m3)
...
To simplify the following description an index flow (q) can be defined as:
| Equation 5 |  |
|---|
From which the Muskingum storage function (Equation 4) can be rewritten as:
| Equation 6 |  |
|---|
Source offers a form of variable parameter Muskingum routing in which the value of x is a constant but the value of K is allowed to vary with flow (Koussis, 1978). From the Kleitz-Seddon rule (Brutsaert, 2005: p.190) the speed that a small monoclinal wave will move through a reach is:
| Equation 7 |  |
|---|
where:
cw is the wave speed (m/s)
...
Equation 7 can be rearranged to give a relationship for the time this wave would take to travel the length of a reach by noting that the volume of water in the reach is equal to the reach length times the cross-sectional area:
| Equation 8 |  |
|---|
where:
Tw is the time the wave takes to pass through (s)
...
For the convenience of the model users (it being easier to think in terms of the wave travel times) Source defines a variable K which uses the wave travel time interpretation. Combined with the assumption that the flow rate in the reach is represented by the index flow q:
| Equation 9 |  |
|---|
To convert this equation into a storage function it is necessary to integrate it:
| Equation 10 |  |
|---|
Where S(x) is the link storage at flow rate x.
...
The storage function has the following form:
| Equation 11 |  |
|---|
where:
a is an arbitrary constant
...
Integrating Equation 11 gives a storage function (Equation 10) of:
| Equation 12 |  |
|---|
Defining:
| Equation 13 |  |
|---|
and
| Equation 14 |  |
|---|
allows Equation 12 to be rewritten in the more familiar form:
| Equation 15 |  |
|---|
where:
k is a storage-delay constant
...
In terms of k and m, Equation 11 is:
| Equation 16 |  |
|---|
By choosing the appropriate functional form and appropriate values of m and x the following routing methods can be replicated:
...
If x=0 then for this formulation the routing is unconditionally stable. If any other value of x is used then for hydrological stability there are two criteria that must be met, and these are:
| Equation 17 |  |
|---|
| Equation 18 |  |
|---|
where:
is q
Rearranging Equation 18, the second criterion is that:
...
- For the power function with m=1, the function K(
) is in fact a fixed value of K and the original version of the Muskingum method is being used. A check is made to enforce the requirement that K ≤ dt/x; or - For the power function with m≠1 (and x>0), the value of K() can range between zero and infinity, presenting a potential stability problem. The flow/storage relationship must be modified to avoid this problem and the point at which this should be done is at the largest allowed slope dt/x. This gives:
| Equation 19 |  |
|---|
The corresponding storage value is S(limit).
If m < 1, a linear section is inserted into the start of the flow/storage curve starting from the point (0,0). For this case, the required equation of the linear section is:
| Equation 20 |  |
|---|
where:
Kmax is the routing constant required for the linearised section, and:
| Equation 21 |  |
|---|
However, if the linear section is inserted from the point (0,0) to the point (limit, S(limit), then from Equation 15 and Equation 19, the resultant slope is:
| Equation 22 |  |
|---|
This means the slope of the linear segment from the point (0,0) to the point (limit, S(limit)), K>Kmax, which is too steep, and further adjustment is required. The approach adopted is to offset the power curve so that it meets a linear segment starting from the point (0,0) with a slope of Kmax. The offset required is then:
| Equation 23 |  |
|---|
Therefore:
| Equation 24 |  |
|---|
The storage function then becomes:
Equation 25 |  | when |  |
|---|---|---|---|
 | when |  |
As the routing characteristics of the storage methods are a function of the slope of the storage curve, shifting the curve like this will not affect results.
...
Table 1. Summary of parameters requiring user input
Parameters | Routing- Power Function (Linear and Non-linear) | Routing- Piecewise |
|---|---|---|
Length of routing reach | Yes | Yes |
Initialisation Type (governs whether an initial value of storage or an initial value of flow is entered) | Yes | Yes |
Initial value of storage | Yes | Yes |
Initial value of flow (note this is mandatory for Lagged Flow) | Yes | Yes |
Lag time | No | No |
Number of routing divisions in the link | Yes | Yes |
Representative flow rate (to calculate link delivery time) | Yes | Yes |
x , the weighting factor in routing | Yes | Yes |
The storage delay constant, k, if m≠1 or the Muskingum K, if m=1, in the routing equation | Yes | No |
m, the exponent in the routing equation | Yes | No |
Flow versus travel time relationship | No | Yes |
Governing data for fluxes (reach processes) - details below | Yes | Yes |
Data needs for lateral fluxes (reach processes)
...
Clark, M.P. and Kavetski, D. (2010) Ancient numerical daemons of conceptual hydrological modeling: 1. Fidelity and efficiency of time stepping schemes, Water Resources Research, 46, W10510, doi:10.1029/2009WR008894.Close, A.F. (1996) A new daily model of flow and solute transport in the River Murray. Proc. 23rd Hydrology and Water Resource Symposium. Hobart, 21-24 May: 173-178. Institution of Engineers, Australia.
Koussis, A.D. (1978) Theoretical estimation of flood routing parameters. J. Hydraul. Div. Am. Soc. Civ. Eng., 104(HY1): 109-115.
...