Geoscience Reference
In-Depth Information
increase in roughness) or where there is flow convergence. A number of models that purport to use
two-dimensional kinematic wave equations (in plan) have been proposed, going all the way back to the
original grid-based model of Bernard in 1937. More recent examples are the catchment surface flow
model of Willgoose and Kuczera (1995) and the subsurface stormflow model of Wigmosta
et al.
(1994).
Mathematically, this is not really a good idea. Wherever there is flow convergence, there is the possibility
of two different kinematic waves meeting and a shock front developing. This was, in fact, explicitly
recognised by Bernard in the 1930s (see Hjelmfelt and Amerman, 1980) but there appears to be nothing
in most of the more recent models to explicitly deal with shocks. It seems that they generally rely on
numerical dispersion to smear out any effects of shock fronts. The models work and have been shown to
give good results in some test cases, but if shocks are dispersed in this way they are not true solutions of
the original kinematic wave equations. I am aware of this problem because I once spent several months
trying to develop a finite difference two-dimensional kinematic wave model. The solution kept going
unstable for some test runs on a real catchment topography. It took a long while to realise that this was not
just a program bug, or a matter of getting the numerics of the approximate solution right but a problem
inherent in the mathematics of a kinematic wave description. With hindsight, the reason is fairly obvious!
There are, however, a number of models in the literature that claim to solve the kinematic wave equations
in two dimensions, without any explicit shock-handling routines. It is not always clear how this has
been achieved.
The problems of kinematic shocks can, in fact, largely be avoided in many situations by solving each
node separately in the manner shown in Box 5.7, as if it was a node on a one-dimensional plane, by taking
advantage of the fact that kinematic waves only move downslope. Thus, in theory, there is no dependence
of the solution on downslope conditions. If the kinematic wave equation is solved for flow,
q
, at a node,
then all the inputs from upslope can be lumped together as one input, even if they converge from more than
one node. The solution is made and the resulting inputs can then be dispersed to form the inputs for one
or more downslope nodes as necessary. This is the approach used in DVSHM, the subsurface kinematic
model of Wigmosta
et al.
(1994), which is based on a two-dimensional raster elevation grid, solving
for a depth of saturation for every grid element. Palacios-Velez
et al.
(1998) provide an algorithm for
constructing such a cascade of kinematic solution elements for both TIN and raster spatial discretisations
of a catchment. The only cautionary note that should be added here is that this approach still adds
numerical dispersion in a way that is not well controlled and may still be subject to stability problems
under rapidly changing conditions. The user should also always remember that a kinematic wave model
cannot model backwater effects that arise in channels because of weirs or restrictions to the flow, or in
subsurface flow as a result of groundwater ridging or the drawdown effect of an incised channel in low
slope riparian areas. At least a diffusion analogy model is required to simulate such effects.
Having raised the issue of kinematic shocks and numerical dispersion, it should be added that it is
not necessarily a problem that should worry us unduly. Problems of parameter calibration might well
dominate any physical and theoretical approximations inherent in the numerical solution of kinematic
wave equations. We can in fact use the numerical dispersion to our advantage, if the resulting model
is actually more realistic. An example here is the still widely used Muskingum-Cunge channel flow
routing model. The Muskingum method was originally developed as a conceptual flow routing model
in the 1930s. Cunge (1969) showed that the Muskingum method was equivalent to a four-point, explicit
finite difference solution of the kinematic wave approximation for surface flow. It follows that, for this
interpretation, any dispersive and peak attenuation effects of the Muskingum-Cunge routing model come
from numerical dispersion associated with the finite difference approximation. In applying the model,
it is normal to fit its two parameters so as to match the observed peak attenuation, which allows some
control over the numerical dispersion by parameter calibration. This is an interesting historical example,
but the details of the Muskingum-Cunge model will not be given here as, in applications to long reaches,
it suffers another serious defect of not properly allowing for the advective time delays in the channel (i.e.
any change in the upstream inflow has an immediate effect on the predicted reach outflow, regardless
