Add multiple flow direction (MFD) with adaptive exponent (Qin et al. 2007)

[brendancol]

Author of Proposal: brendan

Reason or Problem

xarray-spatial has D8 (single steepest neighbor) and D-infinity (continuous angle to steepest facet) for flow direction. Both send all flow in one direction, so they can't represent dispersive flow on hillslopes. Real water spreads out on convex slopes and converges in valleys; a single-direction model can't capture that.

Proposal

Add flow_direction_mfd, a Multiple Flow Direction algorithm based on Qin et al. (2007). MFD partitions flow from each cell to all downslope neighbors. The fraction going to each neighbor is controlled by an adaptive exponent that responds to local terrain. Steep convergent areas get a high exponent (concentrating flow, closer to D8 behavior), while gentle slopes get a low exponent (spreading flow out).

Design:

The function takes an elevation DataArray and returns a 3D DataArray of shape (8, H, W). Each band is the flow fraction to one of the 8 neighbors (E, SE, S, SW, W, NW, N, NE). Fractions sum to 1.0 at each cell. Pits and flats get all-zero fractions; edges and NaN neighborhoods get NaN.

Per cell:

1. Compute slope to each of the 8 neighbors, accounting for cardinal vs. diagonal distance.
2. Keep only downslope neighbors (positive slope).
3. Compute an adaptive exponent p from the ratio of maximum to mean downslope slope.
4. Weight each neighbor: w_i = (S_i * L_i)^p, where L_i is the contour length (cell size for cardinal, cell_size * sqrt(2) / 2 for diagonal).
5. Normalize: f_i = w_i / sum(w_j).

Supports NumPy, CuPy, Dask+NumPy, and Dask+CuPy backends.

Usage:

from xrspatial import flow_direction_mfd

fractions = flow_direction_mfd(elevation)
# fractions.shape == (8, H, W)
# fractions[0] = fraction of flow toward east neighbor
# fractions.sum(axis=0) == 1.0 for cells with downslope neighbors

Value:

MFD is common in distributed hydrological models and erosion/sediment work. D8, Dinf, and MFD are the main flow routing approaches used in practice, so this rounds out the set.

Stakeholders and Impacts

Primarily useful for distributed runoff or sediment transport modeling. No changes to existing functions. MFD-aware flow accumulation would be a separate follow-up.

Drawbacks

The 3D output (8 bands) uses 8x the memory of a single-band D8 result. Dask chunking helps, but it's still a lot more data per cell.

Alternatives

• Freeman (1991) MFD with a fixed exponent (simpler, less adaptive)
• Quinn et al. (1991) equal-weight MFD (no slope-dependent weighting)
• DEMON tube routing (Costa-Cabral & Burges, 1994)

Unresolved Questions

None for the core algorithm. MFD-aware flow accumulation is a separate issue.

Additional Notes or Context

Reference: Qin, C., Zhu, A.X., Pei, T., Li, B., Zhou, C., and Yang, L. (2007). An adaptive approach to selecting a flow-partition exponent for a multiple-flow-direction algorithm. International Journal of Geographical Information Science, 21(4), 443-458.