/
Link storage routing - SRG

Note: This is documentation for version 5.0 of Source. For a different version of Source, select the relevant space by using the Spaces menu in the toolbar above

Link storage routing - SRG

Background

The aspects of a river reach that influence the movement of water are cross-section shape, length, slope and roughness (Figure 1). In addition, surface area and depth influence net evaporation and groundwater fluxes within the reach, respectively. The cross-section can be divided into the river channel and floodplain as shown in Figure 2.

Figure 1. Schematic of a reach river


Figure 2. Cross section of a river reach

The rate at which flow enters the upstream end of a river reach is not necessarily the same as the flow rate that exits the downstream end. As water passes through a reach there may be additional inflows from rainfall on the water surface and from groundwater. Alternatively, water may be lost due to evaporation or seepage to groundwater. The bed friction and slope within a river reach influence the amount of water stored in the reach and the shape of hydrographs. Inflowing water may be trapped in dead storage anywhere along the reach (eg. in pools or depressions). Consequently the flow rate into the reach is not always the same as the flow rate out of the reach even after adjusting for travel time (Figure 3).

Figure 3. Longitudinal section of a river reach showing dead storage

Description and rationale

Storage Routing links model the storage and movement of water through a length of river using a hydrologic routing method. They can represent the travel time of water through a reach, the attenuation of flow rates due to channel shape and roughness and reach processes such as lateral fluxes. Examples of lateral fluxes include net evaporation from the water surface and exchanges between groundwater and surface water.

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

Storage Routing links can also be configured to represent lag flow, equivalent to a Lagged Flow link, but with the capacity to include lateral fluxes. See the Link Routing user guide for how to configure the Storage Routing link to achieve this, and the section Modelling lagged flow routing with storage routing for information on the solution procedure.

Scale

Routing links represent river reaches and these can have lengths ranging from zero to perhaps hundreds of kilometres and widths ranging from zero to many kilometres. Links are used at every model time-step.

Principal developer

eWater CRC is the principal developer of the versions of software in Source.

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.

Version

Source v2.10

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.

Availability/conditions

Automatically included with Source.

Assumptions and Constraints

The Storage Routing link can be partitioned into divisions of equal length. Other assumptions and constraints are listed in Table 1. Further details on assumptions and constraints are given in the Theory section below.

Table 1. Assumptions and constraints applying to storage flow routing methods

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.

Storage Routing Theory

This section covers aspects of the theory relevant to routing and also:

  • Loss and gain fluxes;
  • Dead storage;
  • Implementation methodology with loss and gain fluxes and dead storage;
  • Stability criteria; and
  • Lagging flow when loss and gain fluxes are significant.
Note: Further details on all computations described in this section are given in the section on Solution Methodology, including requirements for the initialisation phase and flow phase. Requirements for the order phase are discussed in Rules-Based Ordering - SRG.

Background

Streamflow routing in Source is carried out using the Muskingum method, which is a lumped kinematic approach to flow routing (Brutsaert, 2005: p.224). It is based on mass conservation and the assumption of a monotonically increasing storage (stage) and discharge relationship in a link. The latter is a simplification of the full momentum equation for river reaches and assumes that diffusion and dynamic effects are negligible. In this approach, the reach is considered to be compressed to a single, representative point or, in the case where the reach has been subdivided, a number of representative points.

The law of conservation of mass then lets us say that for a reach:

Equation 1

where:

S is the reach storage in (m3)

t is time (s)

Put another way: the volume of water within a reach changes at a rate equal to the difference between the inflow and outflow rates to and from the reach.

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.

Form of the Routing Equation

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)

I is the inflow (m3/s)

O is the outflow (m3/s)

K is a constant (s)

x is a weighting factor denoting the importance of inflow relative to outflow

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)

Q is the flow rate (m3/s)

Ac is the channel's cross-sectional area (m2)

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)

From Equations 6 and 8, it can be seen for a fixed K that Tw=K. That is, the time a wave takes to pass through a reach is K. This means that the K in the linear Muskingum routing scheme has two interpretations: as a constant in the storage relationship or as the wave travel time.

Extending the method for variable K

 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.

There are two choices in Source for defining the relationship in Equation 9:

  • Using a functional relationship (a power function); or
  • Using a lookup table (ie. variable parameter Muskingum routing).
Using a Functional Relationship

The storage function has the following form:

Equation 11

where:

a is an arbitrary constant

b is an arbitrary exponent

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

As S has units of m3 and q has units of m3/s, k has the units m3(1-m)sm.

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:

  • Linear Muskingum routing (power function with m=1)
  • Non-linear Muskingum routing (power function with m≠1).
Using a look-up table

The lookup table option provides for the relationship between K and q to be specified as a series of pairs in a lookup table and to use linear interpolation for intermediate values of q. This results in a series of linear relationships:

Equation 16A

where:

ai is half the slope of the segment appropriate for q

bi is the intercept of the segment appropriate for q.

Integrating these straight line segments will result in a storage function in the form of a number of quadratic segments:

Equation 16B

where:

ci is a integration constant for the segment appropriate for q.

Stability Criteria

The routing formulation implemented in Source uses an implicit Euler scheme with the assumptions that:

  • the inflow and outflow rates are constant over the time-step; and
  • the storage function relates the storage at the end of the time-step to the average flow rates during the time-step.

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:

 

Where K() is the slope of the storage curve at .

These criteria are illustrated in Figure 4, and they apply to each reach division. It should be noted that these constraints are much more liberal than the equivalent for the implicit Heun scheme.

Figure 4. Stability criteria for routing methods using average flow formulation

The criterion expressed in Equation 18 is enforced is enforced as follows:

  • For the lookup table option, stability is a function of the largest value of travel time in the lookup table (Kmax). A check is made to ensure that the number of divisions is at least xKmax/dt. To reproduce the behaviour of BigMod (Close, 1996), x must be set to 1 and the number of divisions to twice the minimum required for stability;
  • 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.

If m>1, a linear section replaces the section of the flow/storage curve after the point (limit,S(