Kinematic wave

Surface routing

The main flow routing scheme available in wflow is the kinematic wave approach for river and overland flow, assuming that water flow is mostly controlled by topography. The kinematic wave equations are (Chow, 1988):

$$\begin{array}{c}{\displaystyle \frac{\partial Q}{\partial x}}+{\displaystyle \frac{\partial A}{\partial t}}=q_{\text{inflow}},\\ A=\alpha Q^{\beta }.\end{array}$$

These equations can then be combined as a function of flow only:

$${\displaystyle \frac{\partial Q}{\partial x}}+\alpha \beta Q^{\beta -1}{\displaystyle \frac{\partial Q}{\partial t}}=q_{\text{inflow}},$$

where $Q[{\mathrm{m}}^3{\mathrm{s}}^{-1}]$ is the surface flow in the kinematic wave, $x[\mathrm{m}]$ is the length of the flow pathway, $A[{\mathrm{m}}^2]$ is the cross-section area of the flow pathway, $q_{\text{inflow}}[{\mathrm{m}}^2{\mathrm{s}}^{-1}]$ is the lateral inflow per unit length into the kinematic wave, $t[\mathrm{s}]$ is the integration timestep and $\alpha$ and $\beta$ are coefficients. These coefficients can be determined by using Manning’s equation (Chow, 1988), resulting in:

$$\alpha ={\left(\frac{n{P}_{\mathrm{w}}^{\frac{2}{3}}}{\sqrt{c_{\mathrm{slope}}}}\right)}^{\beta }\text{and}\beta =0.6,$$

where $P_{\mathrm{w}}[\mathrm{m}]$ is the wetted perimeter, $c_{\mathrm{slope}}$ ($c_{\mathrm{land}\mathrm{slope}}$ for overland flow and $c_{\mathrm{river}\mathrm{slope}}$ for river flow) is the slope $[\mathrm{m}{\mathrm{m}}^{-1}]$, and $n$ ($n_{\mathrm{land}}$ for overland flow and $n_{\mathrm{river}}$ for river flow) is Manning’s coefficient $[\mathrm{s}{\mathrm{m}}^{-\frac{1}{3}}]$. The wetted perimeter $P_{\text{w}}$ for river flow is calculated by adding the river width ($w_{\text{river}}$) and 2 times half of the river bankfull depth ($h_{\text{bankfull}}$). A bankfull river depth map (default value is $1\mathrm{m}$) for the river can be provided as follows in the TOML configuration file:

[input.static]
river_bank_water__depth = "riverdepth"

For overland flow, $P_{\text{w}}$ is set equal to the effective flow width, determined by dividing the grid cell area by the flow length and subtracting $w_{\text{river}}$.

The kinematic wave equations are solved with a nonlinear scheme using Newton’s method. By default, the iterations are performed until a stable solution is reached ($ϵ<{10}^{-12}$). The internal time stepping is based on the Courant number or a fixed time step size $[\mathrm{s}]$. To enable a fixed internal time step for the kinematic wave the following lines can be inserted in the TOML file of the model:

[model]
# Adaptive (internal) time stepping for the kinematic wave
kinematic_wave__adaptive_time_step_flag = false    # optional, default is false
# Fixed internal timestep for river flow (river cells)
river_kinematic_wave__time_step = 900        # optional, default is model timestep
# Fixed internal timestep for overland flow (land cells)
land_kinematic_wave__time_step = 3600        # optional, default is model timestep

Kinematic wave adaptive (internal) timestepping can be set as follows in the TOML file of the model:

[model]
# Adaptive (internal) time stepping for the kinematic wave
kinematic_wave__adaptive_time_step_flag = true    # optional, default is false

Reservoir and lake models can be included as part of the river kinematic wave network.

External inflows

External inflows $[{\mathrm{m}}^3{\mathrm{s}}^{-1}]$ for example water supply or abstractions can be added to the kinematic wave for river flow routing, as a cyclic parameter or as part of forcing. For example, cyclic inflow can be provided in the TOML file as follows:

[input.cyclic]
"river_water_inflow~external__volume_flow_rate" = "river_inflow"

These inflows are added or subtracted from the upstream inflow before running the kinematic wave to solve the impact on resulting river flow $Q[{\mathrm{m}}^3{\mathrm{s}}^{-1}]$. In case of an abstraction (negative inflow), the abstraction is limited by a fixed scaling factor of $0.80$ applied to sum of upstream inflow, lateral inflow and river storage, all terms expresssed in unit ${\mathrm{m}}^3{\mathrm{s}}^{-1}$.

Abstractions

Internal abstractions $[{\mathrm{m}}^3{\mathrm{s}}^{-1}]$ from the river are possible when water demand and allocation is computed. The abstraction is set from the water demand and allocation module each time step. The abstraction is divided by the flow length of the river and subtracted from the lateral inflow of the kinematic wave routing scheme for river flow.

Subsurface flow routing

In the wflow_sbm model the kinematic wave approach is used to route subsurface flow laterally by default. Different vertical hydraulic conductivity depth profiles are possible as part of the SBM soil model concept, and these profiles (after unit conversion) are also used to compute lateral subsurface flow. The following profiles (see SBM soil model for a detailed description) are available:

• exponential (default)
• exponential_constant
• layered
• layered_exponential

For the profiles exponential and exponential_constant, the saturated store $S_{\text{sat}}$ is drained laterally by saturated downslope subsurface flow for a slope with width $w[\mathrm{m}]$ according to:

$$Q_{\text{subsurface}}=K_{\mathrm{h}0}{c}_{\mathrm{land}\mathrm{slope}}w\{\begin{array}{ll}\frac{1}{f}(e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,watertable}}}-e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,exp}}})+e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf, exp}}}(z_{\text{ssf,soil}}-z_{\text{ssf,exp}})& \text{if }{z}_{\text{ssf,watertable}}<z_{\text{ssf,exp}}\\ \\ e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,exp}}}(z_{\text{ssf,soil}}-z_{\text{ssf,watertable}})& \text{if }{z}_{\text{ssf,watertable}}\ge z_{\text{ssf,exp}},\end{array}$$

where $c_{\mathrm{land}\mathrm{slope}}[\mathrm{m}{\mathrm{m}}^{-1}]$ is the land slope, $Q_{\text{subsurface}}[{\mathrm{m}}^3{\mathrm{d}}^{-1}]$ is subsurface flow, $K_{\mathrm{h}0}[\mathrm{m}{\mathrm{d}}^{-1}]$ is the saturated hydraulic conductivity at the soil surface, $z_{\text{ssf, watertable}}[\mathrm{m}]$ is the water table depth, $z_{\text{ssf,soil}}[\mathrm{m}]$ is the soil depth, $f_{\mathrm{ssf},\mathrm{Kv}}[{\mathrm{m}}^{-1}]$ is a scaling parameter that controls the decrease of $K_{\mathrm{h}0}$ with depth and $z_{\text{ssf,exp}}[\mathrm{m}]$ is the depth from soil surface for which the exponential decline of $K_{\mathrm{h}0}$ is valid. For the exponential profile, $z_{\text{ssf,exp}}$ is equal to $z_{\text{ssf,soil}}$.

Combining with the following continuity equation: $$({\theta }_s-{\theta }_r)w\frac{\partial h}{\partial t}=-\frac{\partial Q_{\text{subsurface}}}{\partial x}+wR$$

where $h[\mathrm{m}]$ is the water table height, $x[\mathrm{m}]$ is the distance downslope, $R[\mathrm{m}{\mathrm{d}}^{-1}]$ is the net input rate to the saturated store, ${\theta }_s[-]$ and ${\theta }_r[-]$ are the saturated and residual soil water contents, respectively. Substituting for $h(\frac{\partial Q}{\partial h})$, gives:

$$\frac{\partial Q_{\text{subsurface}}}{\partial t}=-c\frac{\partial Q_{\text{subsurface}}}{\partial x}+cwR$$

where celerity $c$ is calculated as follows: $$c=\frac{K_{\mathrm{h}0}{c}_{\mathrm{land}\mathrm{slope}}}{{\theta }_s-{\theta }_r}\{\begin{array}{ll}{e}^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,watertable}}}+e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,exp}}}& \text{if }{z}_{\text{ssf,watertable}}<z_{\text{ssf,exp}}\\ \\ e^{-f_{\mathrm{ssf},\mathrm{Kv}}{z}_{\text{ssf,exp}}}& \text{if }{z}_{\text{ssf,watertable}}\ge z_{\text{ssf,exp}}.\end{array}$$

For the layered and layered_exponential profiles the equivalent horizontal hydraulic conductivity $K_{\text{h}}[\mathrm{m}{\mathrm{d}}^{-1}]$ is calculated for water table height $h=z_{\text{ssf,soil}}-z_{\text{ssf,watertable}}[\mathrm{m}]$ and lateral subsurface flow is calculated as follows:

$$Q_{\text{subsurface}}=K_{\text{h}}h{c}_{\mathrm{land}\mathrm{slope}}w,$$

and celerity $c$ is given by:

$$c=\frac{K_{\text{h}}{c}_{\mathrm{land}\mathrm{slope}}}{{\theta }_s-{\theta }_r}.$$

The kinematic wave equation for lateral subsurface flow is solved iteratively using Newton’s method.

Note

For the lateral subsurface flow kinematic wave the model timestep is not adjusted. For certain model timestep and model grid size combinations this may result in loss of accuracy.

Multi-Threading

The kinematic wave calculations for surface - and subsurface flow routing can be executed in parallel using multiple threads. In the model section of the TOML file, a minimum stream order can be provided to define subbasins for the river (default is 6) and land domain (default is 5). Subbasins are created at all confluences where each branch has a minimal stream order. Based on the subbasins a directed acyclic graph is created that controls the order of execution and which subbasins can run in parallel.

[model]
river_streamorder__min_count = 5 # minimum stream order to delineate subbasins for river domain, default is 6
land_streamorder__min_count = 4  # minimum stream order to delineate subbasins for land domain, default is 5

Subbasin flow

Normally the the kinematic wave is continuous throughout the model. By using the pit__flag entry and basin_pit_location__mask in the model and input sections of the TOML file all flow is at the subbasin only (upstream of the pit locations, defined by the netCDF variable wflow_pits in the example below) and no flow is transferred from one subbasin to another. This can be convenient when connecting the result of the model to a water allocation model such as Ribasim.

[input]
# these are not directly part of the model
basin_pit_location__mask = "wflow_pits"

[model]
pit__flag = true

Limitations

The kinematic wave approach for river, overland and lateral subsurface flow, assumes that the topography controls water flow mostly. This assumption holds for steep terrain, but in less steep terrain the hydraulic gradient is likely not equal to the surface slope (subsurface flow), or pressure differences and inertial momentum cannot be neglected (channel and overland flow). In addition, while the kinematic wave equations are solved with a nonlinear scheme using Newton’s method (Chow, 1988), other model equations are solved through a simple explicit scheme. In summary the following limitations apply:

• River flow, and to a lesser degree overland flow, may be unrealistic in terrain that is not steep, and where pressure forces and inertial momentum are important.
• The lateral movement of subsurface flow may be very wrong in terrain that is not steep.

References

• Chow, V., Maidment, D. and Mays, L., 1988, Applied Hydrology. McGraw-Hill Book Company, New York.