Yet another preprint is available now on ESSOAr. The original manuscript has been submitted to Water Resources Research and can be found here doi.org/10.1002/essoar.10504363.1.
Capturing complex infiltration dynamics in porous-fractured media remains a difficult challenge for any numerical method. While established Richards-based approaches are valid under certain conditions they often fail to capture the erratic and rapid flow dynamics in partially saturated fractures and due to their volume effective nature they may miss important spatial and temporal evolution of infiltration instabilities. Discrete modeling approaches for partially saturated fractures such as SPH on the other hand allow to model the intricate flow patterns based on the Navier-Stokes (NS) equation. Yet, in the presence of a porous matrix they are facing the high contrasts between average pore throat diameters and the dimension of the adjacent fracture void space (several order of magnitude). Resolving these two compartments with the same numerical resolution and underlying Navier-Stokes equation inevitably leads to immense computational demands that are currently out of scope even for modern HPC systems.
In this work we therefore coupled a NS-based SPH model (Kordilla et al. 2017) to a classical Richards equation representing a porous matrix. The code entirely relies on a particle-based meshless discretization of the fluid as well as the continuum pore space. Exchange of fluid between water particles and the pore space is realized via an efficient particle removal algorithm and a virtual saturation function. We validated the code via laboratory experiments and demonstrate the importance of the two flow domains for various examples including infiltration along simple fracture geometries as well as into complex rough fracture sets.