Routing

After calculating the different runoff components, the cell-specific total runoff $(Q_{Tot})$ is calculated by adding these different runoff components. Depending on the modules being switched on, the different runoff components are i) rainfall runoff (RRo), (ii) snow runoff (SRo), (iii) glacier runoff (GRo), and iv) baseflow (BF). Rainfall runoff is the sum of surface runoff (RO, Section 2.7.4) and lateral flow from the first soil layer ($LF_1$, Section 2.7.5). If the groundwater module is not used, then baseflow is calculated as being the lateral flow from the second soil layer $Q_{Tot}$ is eventually calculated according to:

Equation 69

$$QTot=RRo+SRo+GRo+BF$$

with $Q_{Tot} (mm)$ the cell-specific total runoff, RRo (mm) rainfall runoff, SRo (mm) snow runoff, GRo (mm) glacier runoff, and BF (mm) baseflow from the third soil layer or lateral flow from the second soil layer. In order to obtain river discharge, $Q_{Tot}$ needs to be routed through a flow direction network. SPHY allows the user to opt between the use of a simple routing scheme (Section 2.9.1) or a more complex routing scheme (Section 2.9.2) that involves the calculation of lake outflow through $Q(h)$ relations. Both methods require a flow direction network map, which can be obtained by delineating a river network using PCRaster or GIS software in combination with a digital elevation model (DEM).

Runoff routing

In hydrology, streamflow routing is referred to as the transport of water through an open-channel network. Since open-channel flow is unsteady, streamflow routing often involves solving complex partial differential equations. The St. Venant equations (Brutsaert 1971; Morris and Woolhiser 1980) are often used for this, but these have high data requirements related to the river geometry and morphology, which are unavailable for the spatial scale SPHY is generally applied on. Additionally, solving these equations requires the use of very small time steps, which result in large model calculation times. The use of very small time steps in the St. Venant equations is required to provide numerical stability. Other models, such as, e.g., SWAT (Neitsch et al. 2009), use the Manning equation (Manning 1989) to define the rate and velocity of river flow in combination with the variable storage (Williams 1975) or Muskingum (Gill 1978) routing methods to obtain river streamflow. But, the Manning equation also requires river bed dimensions, which are generally unknown on the spatial scale that SPHY generally is applied on.

Therefore, SPHY calculates for each cell the accumulated amount of water that flows out of the cell into its neighboring downstream cell. This can easily be obtained by using the accuflux PCRaster built-in function, which calculates for each cell the accumulated specific runoff from its upstream cells, including the specific runoff generated within the cell itself. If only the accuflux function is used, then it is assumed that all the specific runoff generated within the catchment on one day will end up at the most downstream location within one day, which is not plausible. Therefore, SPHY implements a flow recession coefficient (kx (–)) that accounts for flow delay, which can be a result of channel friction. Using this coefficient, river flow in SPHY is calculated using the three equations shown below:

Equation 70

$$QTOT_t^*=\frac{QTot_t \cdot 0.001 \cdot A}{24 \cdot3600}$$

Equation 71

$$Q_{accu,t}=accuflux(F_{dir},QTot^*_t)$$

Equation 72

$$Q_{rout,t}=(1-kx)\cdot Q_{accu,t}+kx \cdot Q_{rout,t-1}$$

with $QTot_t^*(m^3s^{-1})$ the specific runoff on day $t$, $QTot_t$ the specific runoff in mm on day $t$, $A(m^2)$ the grid-cell area,$Q_{accu,t}(m^3s^{-1})$ the accumulated streamflow on day $t$ without flow delay taken into account, $Q_{rout,t}(m^3s^{-1})$ the routed streamflow on day $t$,$Q_{rout,t-1}(m^3s^{-1})$ the routed streamflow on day $t-1$, $F_{dir}$ the flow direction network, and kx (–) the flow recession coefficient. The coefficient has values ranging between 0 and 1, where values close to 0 correspond to a fast responding catchment, and values approaching 1 correspond to a slow responding catchment. This coefficient is typically used for model calibration.

The user can opt to route each of the four streamflow contributors separately, which may be useful if one wants to evaluate, for example, the contribution of glacier melt or snowmelt to the total routed runoff. However, this increases model run time substantially, because the accuflux function, which is a time-consuming function, needs to be called multiple times, depending on the number of flow contributors to be routed.

Lake/reservoir routing

Lakes or reservoirs act as a natural buffer, resulting in a delayed release of water from these water bodies. SPHY allows the user to choose a more complex routing scheme if lakes/reservoirs are located in their basin of interest. The use of this more advanced routing scheme requires a known relation between lake outflow and lake level height ($Q(h)$ relation) or lake storage.

To use this routing scheme, SPHY requires a nominal map with the lake cells having a unique ID, and the non-lake cells having a value of 0. The user can supply a Boolean map with “True” for cells that have measured lake levels, and “False” for lake cells that do not have measured lake levels. This specific application of SPHY is discussed in detail in Section 3.3.

Four different relations can be chosen to calculate the lake outflow from the lake level height or lake storage, being (i) an exponential relation, (ii) a first-order polynomial function, (iii) a second-order polynomial function, and (iv) a third-order polynomial function. The user needs to supply maps containing the coefficients used in the different functions.

The lake/reservoir routing scheme simply keeps track of the actual lake storage, meaning that an initial lake storage should be supplied. Instead of the simple accuflux function described in the previous section, the lake/reservoir routing scheme uses the PCRaster functions accufractionstate and accufractionflux. The accufractionflux calculates for each cell the amount of water that is transported out of the cell, while the accufractionstate calculates the amount of water that remains stored in the cell. For non-lake cells, the fraction that is transported to the next cell is always equal to 1, while the fraction that is transported out of a lake/reservoir cell depends on the actual lake storage. Each model time step, the lake storage is updated by inflow from upstream. Using this updated storage, the lake level and corresponding lake outflow can be calculated using one of the four relations mentioned before. The lake outflow can then be calculated as a fraction ($Q_{frac}$(-)) of the actual lake storage. Instead of using Equation 71, $Q_{frac}$ is then used in Equation 73 and Equation 74 to calculate the accumulated streamflow and updated storage, respectively:

Equation 73

$$Q_{accu,t}=accufractionflux(F_{dir},S_{act,t},Q_{frac,t})$$

Equation 74

$$S_{act,t+1}=accufractionstate(F_{dir},S_{act,t},Q_{frac,t})$$

with $S_{act,t}(m^3)$and $S_{act,t+1}(m^3)$ the actual storage and updated storage to be used in the next time step, respectively, and $Q_{accu,t}(m^3d^{-1})$ the accumulated streamflow on day $t$, without flow delay taken into account. Since $Q_{frac}$ is always equal to 1 for the non-lake cells, the accufractionflux function becomes equal to the accuflux function used in the previous section. This actually means that for the river network, the same routing function from Section 2.9.1 is used, and that Equation 73 and Equation 74 only apply to lake/reservoir cells.

In order to account for non-linearity and slower responding catchments, the same kx coefficient is used again. This involves applying Equation 72 as a last step after Equation 73 and converting the units from $m^3d^{-1}$ to $m^3s^{-1}$. Since the accufractionflux and accufraction state functions are more complex to compute, the use of these functions increases model run time.

Sediment transport

Sediment is transported using a routing scheme that takes into account both the transport capacity $TC (ton \space ha^{-1})$ of the accumulated runoff and the trapping efficiency of the reservoirs TE (-). The latter only applies when the reservoir module is used. The transport capacity TC (Figure 5a) of the accumulated runoff is based on Prosser & Rustomji (2000):

Equation 75

$$TC=flow_{factor}q_{surf}^{\beta}S^{\gamma}$$

Where $flow_{factor}$ is a spatially distributed roughness factor (-), $q_{surf}$ accumulated runoff per unit width $m^2day^{-1}, S$ the local energy gradient (°), approximated by the slope, and $\beta$ and $\gamma$ are model parameters (-). As suggested by Prosser & Rustomji (2000) $\gamma = 1.4$ and $\beta$ is used for model calibration.

The roughness factor $flow_{factor}$ is determined as follows:

Equation 76

$$flow_{factor}=\frac{\upsilon_{actual}}{\upsilon_{bare}}$$

Where $\upsilon_{actual}$ is the actual flow velocity $(ms^{-1})$ and $\upsilon_{actual}$ is the flow velocity for bare soil conditions $(ms^{-1})$. The actual flow velocity $\upsilon_{actual}$ is obtained from Equation 64-67, applying a water depth $d_{actual}$ of 0.25 m, which coincides with deeper rills from Morgan & Duzant (2008). The flow velocity for bare soil conditions $\upsilon_{actual}$ is obtained from Equation 64, applying values for $n'=0.015 \space s \space m^{-1/3}$ and $d_{bare} = 0.005 m$ (Morgan and Duzant, 2008).

Reservoir sediment trapping efficiency TE (Figure 5b), the percentage of sediment trapped by the reservoir, is calculated according to Brown (1943):

Equation 77

$$TE=100 \bigg [1-\frac{1}{1+0.0021D \frac{C}{A_{basin}} }\bigg]$$

Where D is a constant (-) within the range 0.046-1, depending on the reservoir operation characteristics that we set at 0.1, C the reservoir capacity $(m^3)$, and $A_{basin}$ the drainage area of the subcatchment $(km^2)$

Figure 5. Sediment routing: (a) transport of sediment through the catchment (Equation 75), and (b) trapping efficiency at the reservoirs (Equation 77).