Flow-induced channelization in a porous medium

We propose a theory for erosional channelization induced by fluid flow in a saturated granular porous medium. When the local fluid flow-induced stress is larger than a critical threshold, grains are dislodged and carried away so that the porosity of the medium is altered by erosion. This in turn affects the local hydraulic conductivity and pressure in the medium and results in the growth and development of channels that preferentially conduct the flow. Our multiphase model involves a dynamical porosity field that evolves along with the volume fraction of the mobile and immobile grains in response to fluid flow that couples the spatiotemporal dynamics of the three phases. Numerical solutions of the resulting initial boundary value problem show how channels form in porous media and highlights how heterogeneity in the erosion threshold dictates the form of the patterns and thus the ability to control them.

💡 Research Summary

This paper presents a comprehensive theoretical and numerical framework for flow‑induced channelization in saturated granular porous media. The authors start by highlighting the prevalence of erosion‑driven channel formation in natural systems such as river networks, karst caves, and engineered contexts like groundwater remediation. Existing models typically treat the porous matrix as static, thereby neglecting the feedback between fluid flow, grain mobilization, and evolving permeability. To overcome this limitation, the authors formulate a three‑phase continuum model that simultaneously tracks the fluid phase, a mobile grain phase, and an immobile grain phase. The key state variable is the porosity field φ(x,t), which evolves dynamically as grains are dislodged by shear stresses exceeding a spatially heterogeneous critical threshold τc(x).

The fluid flow is described by Darcy‑Darcy–Forchheimer type equations, with the hydraulic conductivity k expressed as a Kozeny‑Carman function of φ and the mobile grain fraction θm. The erosion rate R is taken proportional to the excess of the local shear stress τ over τc, i.e., R = α·max(τ−τc,0), where α is a material‑specific erosion coefficient. Mobile grains are advected by the seepage velocity and may re‑deposit at a rate S(θm), completing a mass‑conserving exchange among the three phases. The governing equations therefore consist of: (1) fluid mass conservation, (2) mobile grain mass balance, and (3) immobile grain mass balance, all coupled through k(φ,θm) and the stress‑dependent source terms.

Numerical implementation is carried out on a two‑dimensional rectangular domain with an imposed inlet flux and a fixed outlet pressure. The critical stress field τc(x,y) is generated as a Gaussian random field characterized by a mean τ̄c and a standard deviation σc, allowing the authors to explore the impact of material heterogeneity. Spatial discretization uses a second‑order finite‑difference scheme, while temporal integration employs an implicit backward‑Euler method to ensure stability in the presence of stiff source terms.

Simulation results reveal three distinct regimes. When τc is low and spatially uniform, a single dominant channel rapidly forms, with porosity inside the channel rising from the initial 0.35 to values exceeding 0.6, and hydraulic conductivity increasing by an order of magnitude. If τc exhibits strong spatial variability, multiple channels nucleate simultaneously, compete, and sometimes coalesce, producing a complex network whose geometry mirrors the underlying τc landscape. Finally, for inlet fluxes below a critical value, erosion remains localized and the system retains a near‑uniform porosity, whereas fluxes above this threshold trigger a nonlinear acceleration of channel growth, leading to abrupt pressure drops and large‑scale permeability enhancement.

Parametric sweeps of the erosion coefficient α, the mean and variance of τc, and the inlet flux Qin quantify the conditions for channel initiation and propagation. Larger α accelerates channel widening, higher τ̄c suppresses erosion, and greater σc promotes multi‑channel patterns. The authors also demonstrate that modest changes in Qin can shift the system from a stable, low‑erosion state to an unstable, high‑erosion regime, highlighting the sensitivity of channelization to hydraulic forcing.

In the discussion, the authors connect their findings to field observations of river incision, subsurface conduit development, and engineered erosion‑controlled remediation strategies. They acknowledge that the present model neglects grain size distributions, chemical dissolution/precipitation, and compressibility effects, suggesting these as avenues for future refinement. The conclusion emphasizes that a dynamically evolving porosity field, coupled with a stress‑controlled erosion law, captures the essential positive feedback that drives channelization. Moreover, by manipulating the spatial pattern of τc (e.g., through material layering) or the imposed flow rate, one can potentially steer the morphology of the emergent channels. The paper proposes extensions to three‑dimensional geometries, non‑Newtonian fluids, and reactive transport to broaden applicability and to validate the model against laboratory and field data.