  Algorithms for extracting hydrologic features and properties from digital
elevation models (DEMs) are challenged by large datasets, which often cannot
fit within a computer's RAM. Depression filling is an important preconditioning
step to many of these algorithms. Here, I present a new, linearly-scaling
algorithm which parallelizes the Priority-Flood depression-filling algorithm by
subdividing a DEM into tiles. Using a single-producer, multi-consumer design,
the new algorithm works equally well on one core, multiple cores, or multiple
machines and can take advantage of large memories or cope with small ones.
Unlike previous algorithms, the new algorithm guarantees a fixed number of
memory access and communication events per subdivision of the DEM. In
comparison testing, this results in the new algorithm running generally faster
while using fewer resources than previous algorithms. For moderately sized
tiles, the algorithm exhibits ~60% strong and weak scaling efficiencies up to
48 cores, and linear time scaling across datasets ranging over three orders of
magnitude. The largest dataset on which I run the algorithm has 2 trillion
(2*10^12) cells. With 48 cores, processing required 4.8 hours wall-time (9.3
compute-days). This test is three orders of magnitude larger than any
previously performed in the literature. Complete, well-commented source code
and correctness tests are available for download from a repository.
