## 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. |

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 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.

#### 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 (m^{3})

*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 (*S ^{n+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 (m^{3})

*I* is the inflow (m^{3}/s)

*O* is the outflow (m^{3}/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:

*c _{w} *is the wave speed (m/s)

*Q* is the flow rate (m^{3}/s)

*A _{c}* is the channel's cross-sectional area (m

^{2})

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:

*T _{w}* is the time the wave takes to pass through (s)

From Equations 6 and 8, it can be seen for a fixed *K* that *T _{w}=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 m^{3} and *q* has units of m^{3}/s, *k* has the units m* ^{3(1-m)}*s

*.*

^{m}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:

*a _{i}* is half the slope of the segment appropriate for

*q*

*b _{i} *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:

*c _{i}* 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 (
*K*). A check is made to ensure that the number of divisions is at least_{max}*xK*. To reproduce the behaviour of BigMod (Close, 1996),_{max}/dt*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:

*K _{max}* 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>K_{max}*, 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

*K*. The offset required is then:

_{max}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(_{ limit})) which continues the slope (*K _{max}*) at this point.

#### Dead storage

Source has the ability to model the storage remaining in the reach after the stream ceases to flow (called dead storage here - see Figure 3). If the reach dries out this storage has to be refilled before the stream can start to flow again. The total volume in a reach is defined as:

Equation 26 |
---|

where:

*S* is the active storage (m^{3}) and related to through the storage function

*D* is the dead storage (m^{3})

*V* is the total storage (m^{3})

There are three possible outcomes for modelling a reach with dead storage:

- The reach is below dead storage and the fluxes would be insufficient to fill it above dead storage;
- The fluxes in the reach over the time-step would put the reach into dead storage; or
- There is flow in the reach above dead storage.

Evaluation of the state of the reach is discussed further in the section on "Routing process for solution" (in the Process section) below.

#### Reach subdivision

As noted in the Assumptions and Constraints section, above, Source can divide the storage representing the flow routing in a routing link into a number of equal divisions. This means that the flow in a routing link passes through a cascade of storages. The form of the *K() *relationship used for each of the divisions depends upon which of the two options is chosen.

- For the power function relationship, the user inputs a value of the storage-delay constant,
*k*, which applies to each division on the link for use in the K() relationship (Equation 16) for each division; or - When using the lookup table option, in the
*K()*relationship the values of*K*from the input lookup table are divided by the number of divisions used; that is, wave travel times for the entire link are entered into the table.

#### Lateral fluxes (reach processes) modelled

The following lateral fluxes are modelled to determine the volume of water gained or lost in each division during a model time-step:

- Net evaporation
*Flux*;_{NE} - Groundwater
*Flux*;_{GW} - General purpose (function)
*Flux*; or_{flow} - General purpose (time series)
*Flux*._{TS}

By convention, positive numbers indicate losses, negative numbers indicate gains.

Net evaporation and flow based fluxes require iterative solutions as they are functions of index flow. For example, net evaporation changes water surface area, impacting the flow that determines surface area. The general form of each type of flux equation is described below.

#### Net evaporation

The net evaporation flux* Flux _{NE}* is a combination of the evaporation and precipitation that occurs at the water’s surface, and is a function of surface area. Surface area can be determined by multiplying stream width by reach length. The equation below describes these relationships:

Equation 27 |
---|

where:

*E* is rate of evaporation typically specified in mm/day (converted to m/s)

*P* is a rate of precipitation typically specified in mm/day (converted to m/s)

*Width ()* is the relationship between the flow rate and an average top width of the stream (m) (which when multiplied by the reach length will give the surface area of the link).

*Length* is the reach length (m)

*ndiv* is the number of divisions in the routing link

#### Groundwater

Groundwater flux *Flux _{GW}* is a function of head/level which, in turn, is a function of flow. A head or depth from the previous time-step is provided to a linked groundwater model and this returns a corresponding flux:

Equation 28 |
---|

where:

*GroundwaterModel* is a separate model that returns a flux for a given river level (m^{3}/s).

*fsLevel* is a piecewise monotonically increasing function between head and index flow (m)

This approach relies on the assumption that the change in head between time-steps is small. Refer to the Groundwater - SRG for a description of the available groundwater models and their configuration.

#### General purpose function

A general purpose flux *Flux _{flow}* is provided to allow users to specify a flux as a function of flow. This calculation can be represented as:

Equation 29 |
---|

where:

*fFlow *is a piecewise monotonically increasing relationship between flux and index flow (m^{3}/s).

#### General purpose time series

A general purpose time series flux *Flux _{TS}* is provided to represent known losses or gains of particular owners from the river reach. When the time series flux is a loss (positive values) it is limited to reach storage.

Equation 30 |
---|

where:

*TimeSeries* returns a flow rate for a particular owner, *o*, from a data set containing date/time and associated flux (m^{3}/s) for the current time-step *t*.

*o* is the owner

*no* is the number of owners in the system

#### Governing equations

There are two forms of the finite time-step version of the continuity equation (Equation 1) for a reach division depending on what state the division is currently in:

- Flowing: There was enough water in the reach during the last time-step that when the current time-step’s inflow and lateral flows are considered a value of outflow will be able to be found that balances the continuity equation in conjunction with the routing storage function (Equation 10); or
- Ceased to flow: A solution for division outflow will be unable to be found that also satisfies the routing storage function. In this case the reach division has to be treated as if it were a reservoir with no outflow.

The mass balance for a reach division that is flowing can be expressed as:

Equation 31 |
---|

where:

*V _{t-1}* was the volume of water in the reach division during the last time-step (m

^{3}).

*I* is the inflow rate to the reach division this time-step (m^{3}/s).

*O* is the outflow rate from the reach division this time-step (m^{3}/s).

*Flux _{NE}* is net evaporation flux (m

^{3}/s) for storage associated with index flow (from Equation 27)

*Flux*_{GW }is groundwater flux (m^{3}/s) for _{t-1} (from Equation 28)

*Flux _{flow}* is a general purpose flow based flux (m

^{3}/s) for index flow (from Equation 29)

*Flux _{TS }*is a general purpose time series flux (m

^{3}/s) - from Equation 30

*S() *is the storage function that relates the active storage in the reach division to the index flow (m^{3}).

*D _{max}* is the capacity of the reach division dead storage (m

^{3}) and is the dead storage for the entire routing link divided by the number of divisions,

*ndiv*.

*dt *is the model time-step (s).

Note that when the power function is being used the storage function, *S()*, used is Equation 15:

Whereas when the lookup table option (ie. variable parameter Muskingum routing) is being used the equation is:

where:

*a, b, c *are coefficients of the quadratic equation used to interpolate between two points in the piecewise storage versus flow rate relationship.

For the case where a division has ceased to flow the mass balance is:

Equation 32 |
---|

In this case* Flux _{NE}* and

*Flux*are defined as being related to the reach storage,

_{flow}*V*, and

_{t}*Flux*is defined as being related to the head associated with

_{GW}*V*.

_{t-1}## Solution Methodology

The solution methodology for Storage Routing links is discussed below

#### Ordering phase

In the ordering phase, the wave travel time must be estimated before the flows for the current time step are known. The estimated travel time is called “Order Travel Time” and its calculation depends on the type of routing used:

- For Storage Routing links with linear Muskingum routing (), the wave travel time is constant and known. The Order Travel Time is hence equal to the wave travel time and is calculated automatically by Source.
- For Storage Routing links with non-linear and variable parameter Muskingum routing, the travel time depends on the flow. The user must supply an “Average Regulated Flow”, which Source uses to calculate the corresponding
*order*travel time.

Further information on how ordering, resource assessment and ownership are handled in links is provided in the chapters on Rules-Based Ordering - SRG, Resource Assessment - SRG and Ownership - SRG.

#### Flow phase

It should be noted that if *x=1* then a direct solution for O can be made (as =) when the division is flowing. Otherwise an iterative solution is required, for which a modified Newton-Raphson technique is used with two convergence criteria. If either one is satisfied the procedure is considered to have converged. The criteria are a change in index flow of less than 1.0e^{-8} m^{3}/s (equivalent to 0.864 litres per day) and mass balance error less than 0.001 m^{3} (equivalent to 1 litre).

- Assume = 0 then
. Determine the storage associated with the minimum index flow_{min}= x•*S(*and the mass balance error based on Equation 31_{min})*MB*. If_{min}=f_{MB}(_{min},0)*MB*>-0.001 then the stream has ceased to flow and the storage is based on Equation 32;_{min } - Determine the maximum possible index flow:

- If max ≤ min then the stream has ceased to flow as fluxes exceed inflow and storage.
- Determine the storage associated with the maximum index flow
*S(*and the mass balance error based on Equation 31_{max})*MB*. If_{max}=f_{MB}(_{max},O)*MB*< 0.001 then the outflow is the maximum possible outflow;_{max} - Initial estimate .

Determine the storage associated with the index flow*S()*,*O*and the mass balance error based on Equation 31*MB = f*(, ). If |MB| < 0.001 a solution has been found;_{MB} - Use a modified Newton-Raphson solver limited to 20 iterations
- If
*MB*> 0.0 then

max = and MB_{max} = MB

else

min = and MB_{min} = MB

- If |-
_{halve}| < 1.0e^{-8}a solution has been found =_{halve}.

- If |-
_{bisect}| < 1.0e^{-8}a solution has been found =_{bisect}. - Determine the storage associated with the halve index flow S(
_{halve}), and the mass balance error*MB*(_{halve}= f_{MB}_{halve},). If |*MB*| < 0.001 a solution has been found =_{halve}_{halve}. - Determine the storage associated with the bisect index flow S(
_{bisect}), and the mass balance error*MB*(_{bisect}= f_{MB}_{bisect},). If*|MB*|<0.001 a solution has been found =_{bisect}_{bisect}. - Determine the slope of the mass balance index flow function

- Determine the Newton-Raphson correction =
- NR = -
*correction* - If
_{min }< NR <_{max}then - If correction < 1.0e
^{-8}then a solution has been found =_{NR}. - Determine the storage associated with the index flow S(
_{NT}) and the mass balance error MB_{NR}=fMB(_{NR}, ). If |MB_{NR}| < 0.001 a solution has been found =_{NR}. - If MB
_{NR }> 0.0 and_{min}<_{NR }<_{max }then_{max}=_{NR}and*MB*else_{max }= MB_{NR}_{min }=_{NR}and*MB*_{min }= MB_{NR} - If MB
_{bisect}> 0.0 and_{min}<_{bisect}<_{max}then_{max }=_{bisect}and*MB*else_{max }= MB_{bisect}_{min }=_{bisect }and*MB*_{min }= MB_{bisect}. - If MB
_{halve}> 0.0 and_{min}<_{halve }<_{max}then_{max }=_{halve}and*MB*else_{max }= MB_{halve}_{min }=_{bisect }and*MB*_{min }= MB_{halve} - If
_{min}has not been changed

= _{max}

else if _{max} has not been changed

= _{min}

else if |MB|_{min }≤ MB_{max}

= _{min}

_{max} and MB_{max} set to previous max

else

= _{max}

_{min} and MB_{min} set to previous min.

- Determine the storage associated with the index flow
*S()*, and the mass balance error*MB = f*(, )._{MB} - If
*S() ≤ D*a solution has been found with zero outflow_{max} - If |MB| < 0.001 a solution has been found
- Calculate the outflow:

In reaches with multiple divisions the outflow from the reach becomes the inflow into the downstream division. The outflow from the final division is the outflow for the routing reach.

## Data

Refer to the Source User Guide for detailed data requirements and formats.

#### Input data

The information that users may provide is summarised in Table 1. Note that there is also a global requirement to specify the model time-step, *dt*. In addition, of course, data on the inflow to the routing link is required for each model time-step.

###### 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 |

| Yes | Yes |

The storage delay constant, | Yes | No |

| 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)

For the net evaporation flux, *Flux _{NE}*, the additional data requirements are:

- A monotonically increasing rating curve relating a representative average stream width to flow. The average width should be such that when multiplied by the length gives the surface area of the reach;
- A time series of evaporation; and
- A time series of precipitation.

If groundwater exchanges, *Flux _{GW}*, are to be modelled then a monotonically increasing rating curve relating head (or water level) to flow is required. Data for the groundwater model being used is also required (see the chapter on Groundwater modelling for details).

For the general purpose (function) flux, *Flux _{flow}*, a piecewise monotonically increasing relationship between flux and flow is required.

For the general purpose (time series) flux, FluxTS, a time series of losses and/or gains is required for each Owner in the system being modelled. Data required for Owners is discussed in the chapter on Ownership.

## Parameters or settings

Information on the meaning and function of each parameter, whether it is a "physical" parameter or otherwise, and its units can be found in the Theory section above. Where applicable, information on default values and the range of valid values can also be found in the Theory section.

#### Sensitivity of parameters

Clearly, which parameters are the most sensitive depends on the routing option being used. For Storage Routing, the parameters that are relevant from the point of view of sensitivity depend on whether the stream is flowing or it has ceased to flow and, if it is flowing, whether the power function or piecewise routing is being used.

When the stream has ceased to flow, results can be sensitive to the total volume of dead storage specified and to values of parameters defining the fluxes affecting the draw down of this storage. Depending on the fluxes being modelled, sensitive parameters include the adopted routing link length and the adopted shape of the representative cross section of the routing link (ie. as expressed in relationships for width versus flow and/or depth versus flow). Results will also be influenced by other input data, including evaporation and precipitation data if net evaporation is modelled, and the configuration of the groundwater model if groundwater fluxes are modelled.

When the stream is flowing, irrespective of whether the power function or piecewise routing is being used, results can also be sensitive to the adopted routing link length and the adopted size and shape of the representative cross section of the routing link (noting that deciding on this can be problematic when the geometry of the cross section along the link is highly variable). When the power function is being used, results can be sensitive to the choice of values of *k* (or *K*), *m*, and *x*, which in turn are affected by the routing link length, the characteristics of the representative cross section and the friction characteristics of the routing link, and assumptions made about these. When the lookup table option (variable parameter Muskingum routing) is being used, results can be sensitive to the adopted values of travel time, especially the maximum travel time, and travel time is also affected by routing link length, cross section size, shape and variability, and friction characteristics.

As for the case where the stream has ceased to flow, when the stream is flowing results may also be sensitive to values of the parameters defining the fluxes being modelled. The degree of sensitivity will depend on the magnitudes of these fluxes relative to the magnitudes of the flows and reach storage. A circumstance where the results could be expected to be most sensitive to these fluxes would be near cease to flow.

#### Valid ranges of values of parameters

Reiterating the section on "**Form of the Routing Equation**", by choosing the appropriate functional form and appropriate values of m and x the following hydrologic routing methods can be replicated:

- Linear Muskingum routing (power function with
*0 ≤ x < 1*and*m=1*) - Non-linear reservoir routing (power function with
*x=0*and*m>0*) - Non-linear Muskingum routing (power function with
*0 < x ≤ 1*and*m>0*) - Variable parameter Muskingum (BigMod) routing (lookup table)

Values of *m* such that *m>1* are possible for natural channels (see for example, Bates and Pilgrim, 1982; however, such values would typically only occur for overbank flows). Other studies (eg. Wong and Laurenson, 1983, 1984) show that values of *m* such that *m<1* also occur for natural channels, and highlight that the validity of assuming constant values of m and k apply over the full flow range is fairly tenuous for natural channels. The assumption is reasonable for regular geometric shapes though (eg. trapezoidal and triangular cross sections), as found by Mein et al (1974).

From the above list it may be seen that the valid range of values of *x* is *0 ≤ x ≤ 1*. More information on the valid range of values of *K* is given in the section on "**Stability Criteria**".

It is also worth noting that as the flow rate approaches zero, travel time is going to get longer (potentially approach the infinite). However, when the power function is used the maximum travel time (encapsulated in *K*) allowable in the routing method is governed by the stability criteria, where *K _{max} = dt/x* for

*0 < x ≤ 1*, as discussed in the section on "

**Stability Criteria**" (also recalling particularly that when

*m=1*, travel time (

*K*) is constant over the full flow range).

For the lookup table option (ie. variable parameter Muskingum routing), where the maximum travel time in the lookup table input by the user is used to determine the number of routing divisions needed to ensure stability, it may often be appropriate to truncate the travel time in the low flow region by specifying a constant value of travel time for all flows below a minimum threshold (ie. effectively use linear Muskingum routing for low flows).

Information on evaluating *k* (or travel time or *K*), *x* and *m* is given in a number of sources such as Pilgrim (1987). It is worth noting that when discussing storage-discharge relationships, Pilgrim (1987: eg equation 7.4) refers to the term being a "representative discharge for a reach" and this is the same as the index flow used here. Pilgrim (1987: Figure 7.21) also recasts the wave speed-discharge relationship developed by Wong and Laurenson (1983, 1984) into a relationship between time of travel of flood peak and discharge, which is more relevant to the way routing in Source would be typically applied. However, note that the discussion of aspects such as stability criteria by Pilgrim (1987) applies to the implicit Heun scheme and not the implicit Euler scheme used in Source.

## Output data

Outputs include time series of link outflows and modelled fluxes.

## References

Bates, B.C. and Pilgrim, D.H. (1982) Investigation of storage-discharge relations for river reaches and runoff routing models. Proc. Hydrology and Water Resources Symposium. Melbourne, 11-13 May: 120-126. Institution of Engineers, Australia.

Brutsaert, W. (2005) Hydrology - an introduction. Cambridge University Press, Cambridge.

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.

Koussis, A.D. (2009) Assessment and review of the hydraulics of storage flood routing 70 years after the presentation of the Muskingum method. Hydrological Sciences Journal, 54(1): 43-61. February.

Laurenson, E.M. (1959) Storage analysis and flood routing in long rivers. Journal of Geophysical Research, 64(12): 2423-2431, doi:10.1029/JZ064i012p02423.

Linsley, R.K., Kohler, M.A. and Paulhus, J.L.H. (1949) Applied Hydrology. McGraw Hill, New York.

McCarthy, G.T. (1938) The unit hydrograph and flood routing. Manuscript presented at a conference of the North Atlantic Division, US Army Corps of Engineers, 24 June 1938 (unpublished)

Mein, R.G., Laurenson, E.M. and McMahon, T.A. (1974) Simple non-linear model for flood estimation. J. Hydraul. Div. Am. Soc. Civ. Eng., 100: 1507-1518.

Perumal, M. (2010) Discussion of "Assessment and review of the hydraulics of storage flood routing 70 years after the presentation of the Muskingum method". Hydrological Sciences Journal, 55(8): 1427-1430, doi:10.1080/02626667.2010.491260.

Pilgrim, D.H. (1987) Flood routing. Chapter 7 in Pilgrim, D.H. (ed): Australian Rainfall and Runoff - a guide to flood estimation. Vol 1. Institution of Engineers, Australia, Barton, ACT. ISBN: 085825 434 4.

Wong, T.H.F. and Laurenson, E.M. (1983) Wave speed-discharge relations in natural channels. Water Resources Research, 19(3): 701-706, doi:10.1029/WR019i003p00701.

Wong, T.H.F. and Laurenson, E.M. (1984) A model of flood wave speed-discharge characteristics of rivers. Water Resources Research, 20(12): 1883-1890, doi:10.1029/WR020i012p01883.

## Bibliography

Gill, M.A. (1978) Flood routing by the Muskingum method. Journal of Hydrology, 36: 353-363.

Ladson, A.R. (2008) Hydrology: an Australian introduction. Oxford University Press, South Melbourne, Vic., Australia. 304 p. ISBN: 0195553586.

Linsley, R.K., Kohler, M.A. and Paulhus, J.L.H. (1982) Hydrology for Engineers. 3rd Ed., McGraw Hill, Auckland.