Add groundwater flow primitives (Darcy flow, particle tracking) #1147
Description
Problem
The existing hydrology module handles surface water routing (flow direction, flow accumulation, watershed delineation) but has nothing for subsurface flow. Groundwater modeling is a gap in most open-source GIS tools, not just xarray-spatial, so there's an opportunity to provide basic building blocks that hydrogeologists and environmental consultants currently have to code from scratch or pull from specialized (and often expensive) software.
Proposed scope
Darcy flow field
Compute flow direction and magnitude from a water table or piezometric surface raster, given a hydraulic conductivity raster (or scalar). This is analogous to surface slope/aspect but applied to the potentiometric surface, with flow proportional to the hydraulic gradient times conductivity.
Particle tracking
Trace flow paths through the Darcy velocity field, both forward (where does water from this cell go?) and backward (where did the water at this cell come from?). This is the subsurface equivalent of surface flow path tracing.
Capture zone delineation
Identify which cells contribute flow to a given well location by backward-tracking from the well through the velocity field. This is critical for wellhead protection area mapping.
Drawdown estimation
Rasterize analytical solutions (Thiem for steady-state, Theis for transient) around well locations. Given pumping rate, transmissivity, and storativity, produce a drawdown surface. These are textbook equations but having them produce georeferenced rasters directly is useful for overlay with other spatial data.
Implementation notes
- The audience is more specialized than for surface hydrology, but these are operations that people actually need and can't easily get elsewhere in the Python geospatial stack.
- The flow routing infrastructure from the surface hydrology module (particularly the velocity field traversal and path accumulation logic) could be adapted for particle tracking through the Darcy field.
- Darcy flow computation is cell-local (just gradient calculation), so it parallelizes trivially with dask. Particle tracking is sequential along each path but embarrassingly parallel across starting points.