The determination of flow directions is an essential step for drainage network extraction, and flat surfaces are common features in flow direction determination. With the challenge of a massive volume of digital elevation models (DEMs), to reduce the running time and memory usage, there is a growing need to develop parallel algorithms to calculate flow directions over flat surfaces. We propose an efficient parallel algorithm for flow directions over flat surfaces based on the existing serial algorithm and three-step parallel framework. The proposed algorithm assigns pre-divided tiles to consumer processes to build local graphs. Then the producer process builds global graphs based on all the local graphs. Finally, consumer processes update the local graphs based on the global graphs and determine flow directions over flat surfaces. For all tested DEMs, the speed-up ratios are greater than 5 with 11 consumer processes. The strong scaling efficiencies are greater than 40% with 11 consumer processes. The proposed algorithm can run generally faster, use less memory, and process massive DEMs that cannot be successfully processed using the existing serial algorithm. This study shows that the proposed algorithm is an ideal parallel algorithm for determining flow direction over flat surfaces in massive DEMs.