An evaluation of flow-routing algorithms for calculating contributing area on regular grids

Abstract. Calculating contributing area (often used as a proxy for surface water discharge) within a digital elevation model (DEM) or landscape evolution model (LEM) is a fundamental operation in geomorphology. Here we document the fact that a commonly used multiple-flow-direction algorithm for calculating contributing area, i.e., D∞ of Tarboton (1997), is sufficiently biased along the cardinal and ordinal directions that it is unsuitable for some standard applications of flow-routing algorithms. We revisit the purported excess dispersion of the multiple-flow-direction (MFD) algorithm of Freeman (1991) that motivated the development of D∞ and demonstrate that MFD is superior to D∞ when tested against analytic solutions for the contributing areas of idealized landforms and the predictions of the shallow-water equation solver FLO-2D for more complex landforms in which the water surface slope is closely approximated by the bed slope. We also introduce a new flow-routing algorithm entitled IDS (in reference to the iterative depth- and slope-dependent nature of the algorithm) that is more suitable than MFD for applications in which the bed and water surface slopes differ substantially. IDS solves for water flow depths under steady hydrologic conditions by distributing the discharge delivered to each grid point from upslope to its downslope neighbors in rank order of elevation (highest to lowest) and in proportion to a power-law function of the square root of the water surface slope and the five-thirds power of the water depth, mimicking the relationships among water discharge, depth, and surface slope in Manning's equation. IDS is iterative in two ways: (1) water depths are added in small increments so that the water surface slope can gradually differ from the bed slope, facilitating the spreading of water in areas of laterally unconfined flow, and (2) the partitioning of discharge from high to low elevations can be repeated, improving the accuracy of the solution as the water depths of downslope grid points become more well approximated with each successive iteration. We assess the performance of IDS by comparing its results to those of FLO-2D for a variety of real and idealized landforms and to an analytic solution of the shallow-water equations. We also demonstrate how IDS can be modified to solve other fluid-dynamical nonlinear partial differential equations arising in Earth surface processes, such as the Boussinesq equation for the height of the water table in an unconfined aquifer.