Hydrological information analyses based on Digital Eleva-tion Models (DEM) provide hydrological properties derived from high-resolution topographic data represented as an el-evation grid. Flow direction is one of the most computa-tionally intensive functions in the current implementation of TauDEM, a broadly used high-performance hydrological analysis software in hydrology community. Hydrologic flow direction defines a flow field on the DEM that directs flow from each grid cell to one or more of its neighbors. This is a local computation for the majority of grid cells, but becomes a global calculation for the geomorphologically motivated procedure in TauDEM to route flow across flat regions. As the resolution of DEM becomes higher, the computational bottleneck of this function hinders the use of these DEM data in large-scale studies. This paper presents an efficient parallel flow direction algorithm that identifies spatial fea-tures (e.g., flats) and reduces the number of sequential and parallel iterations needed to compute their geomorphologi-cally motivated flow direction. Numerical experiments show that our algorithm outperformed the existing parallel D8 algorithm in TauDEM by two orders of magnitude. The new parallel algorithm exhibited desirable scalability on Stam-pede and ROGER supercomputers.