pith. sign in

arxiv: 2606.12800 · v1 · pith:YWJHID4Onew · submitted 2026-06-11 · 🧮 math.NA · cs.NA

Massively parallel flow routing and drainage area determination

Wolfgang Bangerth This is my paper

Pith reviewed 2026-06-27 06:35 UTC · model grok-4.3

classification 🧮 math.NA cs.NA
keywords flow routingdigital elevation modelsparallel algorithmsdrainage areamassively parallel computinghydrological modelingDEM processingdistributed memory
0
0 comments X

The pith

Parallel algorithm variations route flow on 1.88 billion point DEMs in 4 seconds using 12,288 processes.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper develops parallel algorithms for the flow routing problem on digital elevation models whose size exceeds what sequential methods can handle. It builds directly on the Richardson et al. (2014) approach and introduces several variations intended to reduce communication cost while preserving correctness when thousands of processes cooperate. A sympathetic reader would care because flow routing determines drainage areas and is a foundational step in flood prediction, erosion modeling, and other hydrological applications that now require continental-scale or high-resolution inputs. The reported timing shows that the best variation completes the task for the largest test case in seconds rather than hours.

Core claim

Variations on the Richardson et al. (2014) method can be implemented to solve the flow routing and drainage-area problems on massively parallel computers, with the strongest result being that a model containing 1.88 billion points can be processed in 4.0 seconds on 12,288 processes.

What carries the argument

A family of distributed flow-routing algorithms and their implementation variants that partition the DEM across processes and exchange only the minimal neighbor information required to resolve flow directions and accumulate drainage areas.

If this is right

  • Flow routing and drainage-area calculations become practical for DEMs containing billions of cells on current supercomputers.
  • Higher-resolution or larger-extent elevation data can be used directly in hydrological models without first down-sampling.
  • The same partitioning and neighbor-exchange strategy applies to any grid-based accumulation or path-following computation on distributed memory machines.
  • Runtime remains acceptable even when the number of processes grows into the low tens of thousands.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Similar parallelization tactics could be applied to other accumulation problems on unstructured meshes, such as contaminant transport or sediment routing.
  • The observed scaling behavior supplies a concrete benchmark for testing future load-balancing libraries aimed at irregular grid graphs.
  • If the communication pattern stays sparse, the same code base could support on-the-fly routing inside ensemble weather or climate simulations.

Load-bearing premise

The proposed algorithm variations can be coded so that communication overhead and load imbalance remain small enough at 12,288 processes to preserve both correctness and the reported wall-clock times.

What would settle it

An independent run on the same 1.88-billion-point model and 12,288-process configuration that either produces flow paths differing from the sequential reference or requires substantially more than 4 seconds because of communication stalls.

Figures

Figures reproduced from arXiv: 2606.12800 by Wolfgang Bangerth.

Figure 1
Figure 1. Figure 1: A simple 3 × 3 digital elevation model. Left: Nodes are numbered left-to-right and bottom-to-top. Right: Nodes are numbered high-to-low, as indicated by the ↓ symbol next to each index. about the area that lies upstream of a point. (Or perhaps more correctly: It asks for the area surrounding the nodes that drain through the area surrounding a given node.) The latter problem is a special case of the former:… view at source ↗
Figure 2
Figure 2. Figure 2: Left panel: A digital digital elevation model with 𝑁 = 897 × 513 = 460 161 points of the state of Colorado, USA, the home state of the author. The model is shown projected onto the surface of Earth in three dimensions and clearly shows the Rocky Mountains in the western half of the state, and the valleys of the South Platte river existing at top right, of the Arkansas river existing at bottom right, and of… view at source ↗
Figure 3
Figure 3. Figure 3: A visualization of the nonzero entries of the matrix that appears in the reformulation of flow routing as a system of linear equations (2)–(3). Diagonal matrix entries have a value of +1, all other nonzero entries are −1. Left: The matrix 𝐴 we obtain if we enumerate nodes as in the left panel of [PITH_FULL_IMAGE:figures/full_fig_p007_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: The matrices 𝐴 (left) and 𝐴𝑝𝑝↓ (right) that correspond to the Colorado elevation model of [PITH_FULL_IMAGE:figures/full_fig_p008_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Computed water flow in m3∕year based on the digital elevation model shown in [PITH_FULL_IMAGE:figures/full_fig_p014_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Number of Richardson iterations required to solve the flow routing problem. The three panels show results for digital elevation models of the Colorado topography with 7.3 × 106 , 117 × 106 , and 1.88 × 109 node points. The first three formulations result in nearly identical numbers of iterations. For the rightmost panel, computations are too large to be run with 𝑃 < 12 (matrix-based) or 𝑃 < 8 processes (fo… view at source ↗
Figure 7
Figure 7. Figure 7: Run times required to solve the flow routing problem. The three panels show results for digital elevation models of the Colorado topography with 7.3 × 106 , 117 × 106 , and 1.88 × 109 node points, respectively solved on a laptop with an Intel Core i9-13900H processor (6 performance plus 8 efficiency cores); a workstation with two AMD EPYC 7773X processors (a total of 128 physical cores); and a cluster with… view at source ↗
read the original abstract

Digital elevation models (DEMs) have reached resolutions and sizes that only parallel computaters can efficiently process. One important application of DEMs is predicting how much water flows where, the so-called ``flow routing problem'' (a variation of which is the problem of determining the drainage area upstream of a point in a DEM). The traditional algorithm for flow routing is sequential, and attempts to parallelize this method have so far only been moderately successful. Herein, we build on earlier work in Richardson et al. (2014) and propose an algorithm and several variations that can efficiently solve the flow routing problem on very large models with very large numbers of parallel processes. For the largest model we use, with 1.88 billion points, the best algorithm herein can route water in 4.0 seconds on 12,288 processes of a computer cluster.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

3 major / 1 minor

Summary. The manuscript proposes parallel algorithms and variations, extending Richardson et al. (2014), for solving the flow routing and drainage area problems on large digital elevation models (DEMs). The central empirical claim is that the best variant processes a DEM with 1.88 billion points in 4.0 seconds using 12,288 processes on a computer cluster.

Significance. If the reported timings are shown to reflect correct results with proper verification, the work would address a practical need for scalable processing of high-resolution DEMs in hydrology and geosciences. The absence of any algorithm description, pseudocode, validation data, or scaling analysis in the provided text leaves the significance of the performance claim unestablished.

major comments (3)
  1. [Abstract] Abstract: The headline performance claim (4.0 s on 1.88 B points, 12 288 processes) is presented without any reference to validation against analytic solutions, sequential reference implementations, error metrics, or baseline comparisons. This omission makes it impossible to determine whether the timings reflect correct flow routing or implementation artifacts.
  2. [Abstract] Abstract: No description is given of the proposed variations on the Richardson et al. (2014) method, nor of how boundary handling, halo exchanges, or load balancing are managed at 12 288 processes. Without these details the claim that linear scaling is preserved cannot be evaluated.
  3. [Abstract] Abstract: The manuscript supplies no quantitative data on communication volume, load-imbalance metrics, or global reduction costs, all of which are load-bearing for the parallel-efficiency claim at the reported scale.
minor comments (1)
  1. [Abstract] Abstract, line 2: 'computaters' is a typographical error and should read 'computers'.

Simulated Author's Rebuttal

3 responses · 0 unresolved

We thank the referee for the detailed and constructive report. We address each major comment below and will incorporate the suggested improvements in the revised manuscript.

read point-by-point responses

  1. Referee: [Abstract] Abstract: The headline performance claim (4.0 s on 1.88 B points, 12 288 processes) is presented without any reference to validation against analytic solutions, sequential reference implementations, error metrics, or baseline comparisons. This omission makes it impossible to determine whether the timings reflect correct flow routing or implementation artifacts.

    Authors: We agree that the abstract should reference the validation performed. We will revise the abstract to include a brief statement on validation against analytic solutions and sequential reference implementations, along with the error metrics achieved. The main text will be checked to ensure these validations are presented with sufficient detail and baseline comparisons. revision: yes

  2. Referee: [Abstract] Abstract: No description is given of the proposed variations on the Richardson et al. (2014) method, nor of how boundary handling, halo exchanges, or load balancing are managed at 12 288 processes. Without these details the claim that linear scaling is preserved cannot be evaluated.

    Authors: The abstract is concise by nature and omits these algorithmic details. We will add a short overview sentence to the abstract describing the variations on the Richardson et al. (2014) method and the approaches to boundary handling, halo exchanges, and load balancing. The main text already expands on these topics and will be reviewed for clarity. revision: yes

  3. Referee: [Abstract] Abstract: The manuscript supplies no quantitative data on communication volume, load-imbalance metrics, or global reduction costs, all of which are load-bearing for the parallel-efficiency claim at the reported scale.

    Authors: We acknowledge that the current manuscript does not provide these quantitative metrics. We will add data on communication volume, load-imbalance metrics, and global reduction costs, including new tables or figures in the results section to support the parallel-efficiency claims. revision: yes

Circularity Check

0 steps flagged

No circularity: empirical performance claims with no derivational reduction

full rationale

The paper presents an algorithmic extension of Richardson et al. (2014) for parallel flow routing on DEMs and reports measured wall-clock times (e.g., 4.0 s on 1.88 B points with 12,288 processes). These are direct experimental results, not predictions or derivations from equations. No self-definitional steps, fitted inputs renamed as predictions, or load-bearing self-citation chains appear; the cited prior work is external and the claims remain falsifiable against independent implementations or hardware runs. The derivation chain is therefore self-contained.

Axiom & Free-Parameter Ledger

0 free parameters · 0 axioms · 0 invented entities

Abstract supplies no information on free parameters, axioms, or invented entities; all technical content is deferred to the full manuscript and the 2014 reference.

pith-pipeline@v0.9.1-grok · 5664 in / 1062 out tokens · 19941 ms · 2026-06-27T06:35:20.354758+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

17 extracted references · 15 canonical work pages

  1. [1]

    GeoInformatica 7, 283–313

    Efficient flow computation on massive grid terrain datasets. GeoInformatica 7, 283–313. URL:http://doi.org/10.1023/A:1025526421410, doi:10.1023/A:1025526421410. Arndt, D., Bangerth, W., Bergbauer, M., Blais, B., Fehling, M., Gassmöller, R., Heister, T., Heltai, L., Kronbichler, M., Maier, M., Munch, P., Scheuerman, S., Turcksin, B., Uzunbajakau, S., Wells...

    work page doi:10.1023/a:1025526421410

  2. [2]

    The deal.ii library, version 9.7,

    The deal.II library, version 9.7. Journal of Numerical Mathematics 33, 403–415. doi:10.1515/jnma-2025-0115. Arndt, D., Bangerth, W., Davydov, D., Heister, T., Heltai, L., Kronbichler, M., Maier, M., Pelteret, J.P., Turcksin, B., Wells, D.,

    work page doi:10.1515/jnma-2025-0115 2025

  3. [3]

    The deal.ii finite element library: Design, features, and insights,

    The deal.II finite element library: design, features, and insights. Computers & Mathematics with Applications 81, 407–422. URL:https: //arxiv.org/abs/1910.13247, doi:10.1016/j.camwa.2020.02.022. Balay, S., Abhyankar, S., Adams, M.F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, ...

    work page doi:10.1016/j.camwa.2020.02.022 1910

  4. [4]

    Technical Report ANL-21/39 - Revision 3.24

    PETSc/TAO Users Manual. Technical Report ANL-21/39 - Revision 3.24. Argonne National Laboratory. doi:10.2172/2998643. Balay, S., Abhyankar, S., Adams, M.F., Benson, S., Brown, J., Brune, P., Buschelman, K., Constantinescu, E.M., Dalcin, L., Dener, A., Eijkhout, V., Faibussowitsch, J., Gropp, W.D., Hapla, V., Isaac, T., Jolivet, P., Karpeev, D., Kaushik, D...

    work page doi:10.2172/2998643

  5. [5]

    URL:https://petsc.org/

    PETSc Web page.https://petsc.org/. URL:https://petsc.org/. Bangerth,W.,Burstedde,C.,Heister,T.,Kronbichler,M.,2011. Algorithmsanddatastructuresformassivelyparallelgenericadaptivefiniteelement codes. ACM Transactions on Mathematical Software 38, 14/1–28. Barnes,R.,2016. Parallelpriority-flooddepressionfillingfortrillioncelldigitalelevationmodelsondesktopso...

    work page doi:10.1016/j.cageo.2013.04.024 2011

  6. [6]

    Geomorphology180-181,170–179

    A very efficient o(n), implicit and parallel method to solve the stream power equation governing fluvial incisionandlandscapeevolution. Geomorphology180-181,170–179. URL:https://www.sciencedirect.com/science/article/pii/ S0169555X12004618, doi:https://doi.org/10.1016/j.geomorph.2012.10.008. Burstedde, C., Wilcox, L.C., Ghattas, O., 2011.p4est: Scalable al...

    work page doi:10.1016/j.geomorph.2012.10.008 2012

  7. [7]

    European Space Agency,

    Upslope area—forming and solving the flow matrix, mathworks.http://blogs.mathworks.com/steve/2007/08/07/ upslope-area-flow-matrix/. European Space Agency,

    2007

  8. [8]

    doi:10.5069/G9028PQB

    Copernicus Global Digital Elevation Model, distributed by OpenTopography.https://doi.org/10.5069/ G9028PQB. doi:10.5069/G9028PQB. accessed: 2026-03-19. Freeman,T.,1991.Calculatingcatchmentareawithdivergentflowbasedonaregulargrid.Computers&Geosciences17,413–422.URL:https:// www.sciencedirect.com/science/article/pii/009830049190048I, doi:https://doi.org/10....

    work page doi:10.5069/g9028pqb 2026

  9. [9]

    Hysom, D., Pothen, A.,

    Efficient parallel computation of ilu(k) preconditioners, in: Proceedings of the 1999 ACM/IEEE Conference on Supercomputing,AssociationforComputingMachinery,NewYork,NY,USA.p.29–es.URL:https://doi.org/10.1145/331532.331561, doi:10.1145/331532.331561. Hysom, D., Pothen, A.,

    work page doi:10.1145/331532.331561 1999

  10. [10]

    SIAM Journal on Scientific Computing 22, 2194–2215

    A scalable parallel algorithm for incomplete factor preconditioning. SIAM Journal on Scientific Computing 22, 2194–2215. doi:10.1137/S1064827500376193. Jenson, S.K., Domingue, J.O.,

    work page doi:10.1137/s1064827500376193

  11. [11]

    Water Resources Research 58, e2021WR030948

    An improved D8-LTD for the extraction of total con- tributing area (TCA) by adopting the strategies of path independency and local dispersion. Water Resources Research 58, e2021WR030948. URL:https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021WR030948, doi:https://doi.org/ 10.1029/2021WR030948,arXiv:https://agupubs.onlinelibrary.wiley.com/doi/pdf/...

    work page doi:10.1029/2021wr030948

  12. [12]

    doi:10.5067/ASTER/ASTGTM.003

    ASTER Global Digital Elevation Model V003.https: //www.earthdata.nasa.gov/data/catalog/lpcloud-astgtm-003. doi:10.5067/ASTER/ASTGTM.003. accessed: 2026-03-19. O’Callaghan, J.F., Mark, D.M.,

    work page doi:10.5067/aster/astgtm.003 2026

  13. [13]

    Computer Vision, Graphics, and Image Processing 28, 323–344

    The extraction of drainage networks from digital elevation data. Computer Vision, Graphics, and Image Processing 28, 323–344. URL:https://www.sciencedirect.com/science/article/pii/S0734189X84800110, doi:https: //doi.org/10.1016/S0734-189X(84)80011-0. Quarteroni, A., Valli, A.,

    work page doi:10.1016/s0734-189x(84)80011-0

  14. [14]

    Clarendon Press, Oxford

    Domain Decomposition Methods for Partial Differential Equations. Clarendon Press, Oxford. Richardson,A.,Hill,C.N.,Perron,J.T.,2014. IDA:Animplicit,parallelizablemethodforcalculatingdrainagearea. WaterResourcesResearch50, 4110–4130. URL:https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2013WR014326, doi:10.1002/2013WR014326, arXiv:https://agupubs.onl...

    work page doi:10.1002/2013wr014326 2014

  15. [15]

    Philosophical Transactions of the Royal Society of London, Series A: Containing Papers of a Mathematical or Physical Character 210, 307–357

    The approximate arithmetical solution by finite differences of physical problems involving differential equations, with an application to the stresses in a masonry dam. Philosophical Transactions of the Royal Society of London, Series A: Containing Papers of a Mathematical or Physical Character 210, 307–357. URL:http://doi.org/10.1098/rsta.1911.0009, doi:...

    work page doi:10.1098/rsta.1911.0009 1911

  16. [16]

    Environmental Modelling & Software 25, 770–781

    Topotoolbox: A set of matlab functions for topographic analysis. Environmental Modelling & Software 25, 770–781. URL:https://www.sciencedirect.com/science/article/pii/S1364815209003053, doi:https://doi.org/10.1016/ j.envsoft.2009.12.002. Tarboton, D.G.,

    2009

  17. [17]

    Water Resources Research 33, 309–319

    A new method for the determination of flow directions and upslope areas in grid digital elevation models. Water Resources Research 33, 309–319. URL:http://doi.org/10.1029/96WR03137, doi:10.1029/96WR03137. Wallace,R.,Tarboton,D.,Watson,D.,Schreuders,K.,Tesfa,T.,2010. Parallelalgorithmsforprocessinghydrologicpropertiesfromdigitalterrain, in:Proceedingsofthe...

    work page doi:10.1029/96wr03137 2010