Flow-driven branching in a frangible porous medium

Channel formation and branching is widely seen in physical systems where movement of fluid through a porous structure causes the spatiotemporal evolution of the medium in response to the flow, in turn causing flow pathways to evolve. We provide a simple theoretical framework that embodies this feedback mechanism in a multi-phase model for flow through a fragile porous medium with a dynamic permeability. Numerical simulations of the model show the emergence of branched networks whose topology is determined by the geometry of external flow forcing. This allows us to delineate the conditions under which splitting and/or coalescing branched network formation is favored, with potential implications for both understanding and controlling branching in soft frangible media.

Branching patterns, or arborization, in porous media are common in many natural settings that include both living and non-living matter [1]. The formation of arborized patterns in physical and chemical systems is driven by a variety of processes all of which involve a combination of erosion, transport and deposition. On the laboratory scale, these processes can involve the chemical dissolution of brittle matrices by a penetrating reactive fluid [2,3], the advective rearrangement of unconsolidated media, dielectric breakdown of conducting media [4,5], the formation of fingerlike protrusions in dense granular suspensions [6], formation of beach rills in natural drainage systems [7,8] etc. On planetary scales, melt transport in the mantle arises via branching morphologies that lead to localized channels of widths up to 100 m [9][10][11], and water-driven erosion and branching in glaciers arises on scales of the order of 10 m [12]. In biological systems, the best known arborized systems are vasculatures in plants and animals. These arise through morphogenetic mechanisms involving gradients and physical flows that arrange and rearrange matter through a variety of feedback mechanisms at the cellular, organismal, and societal level [13,14], and examples include slime molds [15], vascular networks [16], and nest architectures of social insects [17].
Models based on porous flow theory [18,19] are capable of describing flow through these branched networks. However, their formation requires nonlinear models with multiple evolving phase boundaries which are still only partially understood both theoretically and experimentally. Here, we propose a simple model via an effective continuum theory that links flow, permeability and pressure gradients by considering pore-scale grain dislodgement in a relatively brittle structure. Numerical solutions of the resulting governing equations show the emergence of branching morphologies through selective erosion and subsequent flow enhancement.  (4) and (5), eroded regions of low solid fraction φ(x, t) are blue (see Fig. 3 for colorbar limits). At a given point on the macroscale (black dot), φ is the fraction of a pore-scale integration volume (inset) occupied by rigid grain microstructure (red circles). Fluid-mediated forces on grains induce stresses over a macroscopic region (shaded green circle) characterized by the communication length ξ. The response fraction ϕ is the spatial average of φ throughout this region. (b) The erosion threshold function ψ(ϕ) = c0 tanh ω(ϕ − ϕ * ) + c1 ∈ [0, 1] represents resistance to grain dislodgement at response fraction ϕ.
Mathematical model: Our starting point is a fluid-filled porous domain Ω comprised of a rigid grain microstructure with characteristic pore size l, as in the Fig. 1(a) inset. The fluid is of viscosity η and density ρ. On length scales large compared to the pore size L l, we can define macroscopic continuum fields that include the solid fraction φ(x, t) and volumetric fluid flux q(x, t) as averages of microscopic quantities [20]. Pressure gradients over macroscopic lengths Γ ∼ |∇p| drive motion of the interstitial fluid V ∼ |q| relative to the pore structure, so that at a scaling level Γ ∼ ηV /l 2 , leading to individual grains feeling forces of magnitude Γl 3 . When these overcome the attractive forces providing microstructural integrity, grains are dislodged and the local permeability of the medium evolves. If we denote the magnitude of the local network breaking stress by B(x, t), which can be heterogeneous, the most general rate law consistent with this mechanism reads where e 0 is an erosion rate and f (α, β) is a nonnegative dimensionless function which vanishes for α < β. A previous model [21] accounts for the relative motion of the grains, fluid and the static porous medium via a threephase model of fluid, immobile solid and mobile grains.
Here, we focus on a simpler two-phase model assuming loose grains to be indistinguishable from fluid.
In terms of a characteristic breaking stress B 0 and time scale τ = 1/e 0 , we can define a characteristic length L = l(B 0 /ηe 0 ) and pressure gradient magnitude Γ = B 0 /l. Rescaling our variables and parameters accordingly, and assuming that that the solid is relatively stiff but brittle so that it does not deform, the volumetric fluid flux q is well described by Darcy's law, where the dimensionless permeability κ(φ) is the wellknown Carman-Kozeny relation [18,19]. Assuming the fluid is incompressible, conservation of mass implies where s(x, t) is the rate at which fluid is depleted due to processes such as bulk reaction or evaporation. By combining the previous two equations, q can be eliminated to obtain an elliptic equation for the pressure, Boundary conditions correspond to specified fluxes q in and q out on boundary regions of inflow ∂Ω in and outflow ∂Ω out (See Supplementary Information (SI) section SI.1 for details).
To close the system, we must relate the erosion rate f to the fields φ(x, t) and p(x, t). A minimal analytic form for f with a breaking threshold based on symmetry arguments [21] suggests f = max{0, ∇p · ∇p − B 2 }, where B is the breaking threshold. The breaking stress itself is a nonlocal function of the solid fraction, depending on the grain density within a region of size ξ, a stress communication length which may depend on the porosity. Here, we assume the following hierarchy of lengths l ξ L, consistent with frangible brittle solids. In this limit, we introduce a simple erosion threshold B 2 = ψ (ϕ), defining the response fraction ϕ(x, t) as the convolution of φ(x, t) with a Gaussian kernel of length scale ξ, representing a spatial average of the solid fraction as shown in Fig. 1(a) (see SI.2 for details). Thus, the erosion rate law (1) becomes For the functional form of the threshold, we consider a sigmoid ψ(ϕ) ∈ [0, 1] centered at a critical phase fraction ϕ * , where the behavior is roughly linear over a scale ∆ϕ ∼ 1/ω, where ω represents a sharpness parameter as shown in Fig. 1(b). See SI.3 for the exact form. We note that our functional choice satisfies ψ (ϕ) > 0, i.e. the medium becomes more resistant to erosion at larger ϕ.
Equations (4) and (5) together determine the evolution of the permeability of the porous medium, φ(x, t), and the pressure, p(x, t), once we specify an initial condition. Ignoring anisotropy in grain orientation and packing, we set φ(x, 0) = φ 0 + δφ(x), with φ 0 a constant and δφ a perturbed packing structure described as a random Gaussian thermal noise field with zero mean, variance σ 2 φ , and correlation length ζ l, such that Here, * r = Ω ( * )dxdy/vol(Ω) is a spatial average over all x, y ∈ Ω such that |x − y| = r. width w c . From (5), loss of solid material at a point reduces resistance to further erosion in a surrounding neighborhood of size ξ-qualitatively similar to descriptions of nonlocal damage accumulation in settings such as hydraulic fracturing [22]. Features in the φ-field, initially of size ζ, correspond to smoothed features in the ϕ-field. Thus, the channel width scaling satisfies ζ 2 < w 2 c < ξ 2 + ζ 2 , approaching the small limit for large values of the packing variance σ 2 φ and vice versa, consistent with results obtained using the three-phase model [21]. The width of a given channel scales with w c and varies with the amount of flux it conducts. See SI.4 for details.
Before considering the spatiotemporal evolution of the flow and permeability, we examine the local dependence of erosion on the threshold shape ψ(ϕ) and local flux q. Letting * = A ( * )dx/vol(A) denote a spatial average over a mesoscopic region A, we introduce the scalar fields Φ = φ , G 2 = |∇p| 2 , and Q = |q| . We see they satisfy Q = −κ(Φ)G, derived by averaging (2). Differentiating this relation and combining it with an averaged (5) yields a set of purely time-dependent equations describing trajectories through forcing-response phase space. For eroding states with G 2 > ψ(Φ), Sustained erosion does not occur ifQ = 0, in which case points on the threshold surface G 2 = ψ are stable equilibria of the system. Eroding states reach the threshold in finite time, as can be seen from (7a). ForQ/Q > 0, this is not the case. The quantity Φκ /κ < 0 is negative, so the squared-gradient decays or grows when the first or second term in (7b) respectively dominates the other.
The majority of the system's evolution takes place along a monotonically increasing slow manifold G s (Φ) where the two are balanced, corresponding to d(G 2 )/dt = 0, whence from (7b) a cubic with one real root. In Fig. 2(a) we show the trajectories and slow manifolds for varyingQ. In Fig. 2(b) we plot the rate of erosion on the manifold, f s = G 2 s − ψ, for thresholds of varying steepness at a particular flux rateQ. Theoretical bounds f 0 > f s > f 1 , corresponding to constant thresholds ψ = 0 and 1, are plotted as dotted black lines. Both are monotonically increasing, diverge as Φ → 1, and vanish as Φ → 0, so the rate of erosion slows over long times. This effect is mitigated by a transition from f ≈ f 1 to f 0 near Φ ≈ ϕ * . For sharp thresholds of large ω, this effect is dominant and erosion accelerates upon reaching the transition region. The relative difference between the bounding rates, (f 0 − f 1 )/f 1 , vanishes as Φ grows. It can be shown that for ϕ * > 1/2, sharper thresholds yield faster average erosion over the entirety of the system's evolution. See SI.5 for details.
Branching morphospaces: We now turn to the spatiotemporal evolution of the flow and permeability fields in two-dimensional simulations. We aim to understand when, how and what arborization motifs arise as a function of the boundary conditions, the dynamical rate of boundary fluxes, and the nature of the fragility/breaking threshold function. We integrate the coupled set of equations (4) and (5) on a square domain Ω = [−5, 5] 2 ∈ R 2 , employing a second-order forward Euler method with Richardson extrapolation for error estimation and adaptive time stepping [23]. See SI.6 for details. We adopt boundary conditions which ramp up the flux from zero over a duration q in (x, t) =q in r(t), q out (x, t) =q out r(t), where we have introduced a set of hatted constants corresponding to final magnitudes. This formulation yields a uniform bulk fluid sinkŝ evenly distributed throughout the domain. Similarly, the boundary fluxes are assumed to be uniform everywhere on the regions ∂Ω in and ∂Ω out , which we center on the bottom and top walls of the domain, respectively. We note the sign ofŝ may be reversed and the labels "in" and "out" swapped with no change to morphogenic pattern formation, because the erosion rule (1) is agnostic to the substitution ∇p → −∇p.
There are two feedback mechanisms through which erosion in the model promotes itself. The first, observed in the homogeneous system, is the threshold reduction due to previous erosion. The second is a direct effect of the coupling between flux and permeability. According to (2), flux is preferentially directed along paths of larger permeability, so that as it grows, flow from other parts of the domain is redistributed to eroded areas. In terms of the homogeneous phase space shown in Fig. 2(a), the resulting flux increase moves quickly eroding areas onto slow manifolds G 2 s of higherQ, speeding up erosion. Slowly eroding areas experience the opposite effect until so much flow is diverted thatQ ≤ 0, so erosion ceases. In this way, flow enhancement leads directly to selective erosion of high-κ channelized regions of width w c . For a given integrated fluid flux at the boundary F , the number of channels to form in the absence of geometric constraints will scale as N c ∼ F/w c ; in what follows, N c > 1.
In Fig. 3 we show the results of simulations with four different combinations of boundary conditions and the role of bulk evaporation. In the first four panels Fig. 3(a-d) we setŝ = 0 and consider the effect of variation in boundary flux width. Generically, if both the inlet width w in and outlet width w out are larger than the emergent channel size w c , boundary fluxes induce the formation of multiple channels, as seen in Fig. 3(a). If both are less than w c , a single channel is favored as in Fig. 3(b). (We note branching in these settings is possible- Fig. 1(a) shows a single channel split and consolidate-but only given conveniently located low-κ regions of the initial condition in the ξ ζ limit.) If the reverse is true, i.e. w in < w c < w out , then N c channels are created at the outlet and one at the inlet, as in Fig. 3(c) and (d). In Fig. 3(e), we show the effect of bulk-evaporation driven flow withŝ > 0, a single inlet and no outlet. Because the channel width w c L the system size, multiple channels form in the bulk, although their number and width is attenuated with distance from the inlet. These results may be summarized via a simple geometric argument suggesting a formula for reliable branch generation. If the number of channel heads distributed along the inlet and outlet are not the same, branching junctions arise in their linking, which is favored by flow continuity.
Finally, we consider the effects of varying the form of the erosion threshold function ψ(ϕ) via its steepness ω, and the rate of flux increase, via the ramp-up time T . Fig. 4 shows a grid of eroded patterns corresponding to combinations of these two parameters. Low rates of flux increase correspond to slow manifolds G 2 s close to the threshold ψ, so small drops in the pressure gradient can yield |∇p| 2 < ψ. Conversely, rapidly increasing fluxes induce large pressure gradients |∇p| 2 ψ before flow reorganization can occur, leading to large-scale washout, also seen in the three-phase model [21]. We conclude that T 1 is necessary for selective erosion. Increasing ω yields more erosion across the domain and appears to form sharper channel boundaries; this is consistent with the relationship between the erosion rate f = |∇p| 2 − ψ and ω as in Fig. 2(b). As discussed, sharper thresholds induce faster average erosion, so more solid is eroded overall. In particular, systems with sharper thresholds have relatively higher rates of erosion in the region ϕ < ϕ * . This speeds up flow enhancement and thus increases erosion selectivity.
Conclusions: Our minimal continuum model for the coupled dynamics of erosion, flow and permeability in a porous material shows how complex branching patterns can arise from simple causes. While the model and discussion are rooted in the language of fragile solids, our framework is broadly applicable beyond this setting, to branching patterns generated by local interactions subject to non-local flow constraints. Generalizing this to biological settings that feature non-linear couplings such as that between nutrient concentration and flow behavior, e.g. if portions of solid may be flow-seeking or flow-avoiding [14] is a natural next step.