PyDebrisFlow2D: A 2D Conservative Python Multiphysics Solver for Variable-Density Debris Flows with Voellmy Resistance, Bed Exchange, and Grain Size Segregation

Reliable debris-flow modelling requires the coupled treatment of rapid free-surface propagation, variable mixture composition, wet–dry fronts, basal resistance, bed exchange, and grain size segregation. This work presents PyDebrisFlow2D, a depth-averaged Python (3.13.14) solver based on seven conservative variables describing mobile fluid, layered fine and coarse solids, and two horizontal mixture momentum components. Mixture density and composition-dependent Voellmy resistance are derived consistently from the transported constituent volumes, while conservative source operators account for erosion, entrainment, class-resolved deposition, and segregation opposed by diffusive remixing. The governing equations are solved using a Cartesian finite volume method with hydrostatic reconstruction, Harten–Lax–van Leer–Contact (HLLC) fluxes, optional Monotonic Upstream-Centered Scheme for Conservation Laws (MUSCL) reconstruction, and first- or second-order time integration. Local Harten–Lax–van Leer (HLL) fallback, stage repetition, positivity repair, and adaptive time-step rejection provide additional robustness near wet–dry interfaces and strong gradients. Controlled verification experiments confirmed flux consistency, well-balanced behaviour, constituent conservation, and stable convergence for different slope limiters. The verification programme also examines the reduced segregation formulation, an alternative concentration-aware resistance closure, and the applicability limits of the depth-averaged model through targeted real-topography analyses. Across the Ritter refinement experiment, the nominal second-order scheme required 1.91–3.13× the first-order (O1) wall time while reducing the depth error; at the finest grid, the relative L1 error was 13.35× smaller for a 2.309× runtime increase. A propagation-only application on the Marsicano digital elevation model reproduced the principal channelized path and terminal accumulation area while conserving the released volume. A complementary ablation analysis showed that erosion, deposition, and segregation produce physically coherent and conservative modifications of the simulated flow without introducing numerical mass imbalance, while CUDA speedup over a 24-thread CPU reached 67.39× for the largest Marsicano benchmark. These results provide process-level verification of the numerical implementation and its multiphysics coupling, although they do not constitute observational validation of the simulated event. The Marsicano comparison offers a qualitative cross-model consistency check against a published scenario produced with the established Rapid Mass Movement Simulation (RAMMS) framework, which has been extensively applied and evaluated in real-event studies. The solver, verification tests, and reproducible configurations are publicly available as free and open-source research software. Event-scale predictive applications nevertheless require site-specific calibration and independent validation against field observations.