Headward growth and branching in subterranean channels

We investigate the erosive growth of channels in a thin subsurface sedimentary layer driven by hydrodynamic drag toward understanding subterranean networks and their relation to river networks charged by ground water. Building on a model based on experimental observations of fluid-driven evolution of bed porosity, we focus on the characteristics of the channel growth and their bifurcations in a horizontal rectangular domain subject to various fluid source and sink distributions. We find that the erosion front between low- and high-porosity regions becomes unstable, giving rise to branched channel networks, depending on the spatial fluctuations of the fluid flow near the front and the degree to which the flow is above the erodibility threshold of the medium. Focusing on the growth of a network starting from a single channel, and by identifying the channel heads and their branch points, we find that the number of branches increases sublinearly and is affected by the source distribution. The mean angles between branches are found to be systematically lower than river networks in humid climates and depend on the domain geometry.

Figure 1

(a) Schematic of the experimental system. Fluid enters the porous granular bed through the top boundary from a reservoir either all across or near the two corners. An image of the initial medium corresponding to the area that can be visualized in the experiments is also shown. The fluid is injected into the inlet reservoir with a prescribed flow rate and is allowed to exit from an outlet at the center into a reservoir which can also collect the eroded material. (b, c) The calculated magnitude of the fluid velocity for the prescribed boundary conditions along with streamlines (black and blue lines). The gray scale bar and the numbers corresponds to the velocity $V$ scaled by the injection flux $J$.

Figure 2

(a) The distribution of the local granular packing fraction ${\varphi }_g$ measured using the Voronoi cell size. The mean volume fraction ${\varphi }_m=0.41$. (b) The ratio of the mean grain speed $v_g$ and the prescribed mean fluid flux $J$ as a function of $J/J_c$. The vertical dashed lines indicate range of critical flow required to erode a grain at rest on the substrate. The error bars indicate the percentage range of measured $v_g$ from grain to grain.

Figure 3

(a) The fluid velocity $V$ obtained using numerical simulations normalized by the maximum velocity $V_{\max }$ at $y=0.67L$ (uniform injection, $\mathrm{\Delta }\varphi =\pm 0.04$). (b) $V/V_{\max }$ at distance $r=0.1L$ from outlet (uniform injection, $\mathrm{\Delta }\varphi =\pm 0.04$). (c) $V/V_{\max }$ at $y=0.67L$ (nonuniform injection, $\mathrm{\Delta }\varphi =\pm 0.04$). (d) $V/V_{\max }$ at $r=0.1L$ from outlet (nonuniform, $\mathrm{\Delta }\varphi =\pm 0.04$). (e) $V/V_{\max }$ at $y=0.67L$ (nonuniform, $\mathrm{\Delta }\varphi =0.0$). (f) $V/V_{\max }$ at $r=0.1L$ from outlet (nonuniform, $\mathrm{\Delta }\varphi =0.0$). Velocity $V/V_{\max }$ is observed to be large along multiple directions when ${\varphi }_g$ fluctuates, but the velocity is peaked along the symmetry axis when fluctuations are absent.

Figure 4

A contour map of the front between volume fraction region corresponding to ${\varphi }_g=0.35$ observed in the experiments as the flux $J/J_c$ is increased in (a) the uniform injection case and (b) the nonuniform injection case. The channels develop at a lower value of $J$ compared with $J_c$ because of the convergence of the flow near the outlet. Several channels develop near the outlet, but only a few advance upstream toward the inlet. The bounding box corresponds to the viewing area shown in Fig. 1.

Figure 5

(a–d) Evolution of channels observed in the simulations for the uniform injection case. (a) $t_s=250$, $J/J_c=0.13$, (b) $t_s=350$, $J/J_c=0.18$, (c) $t_s=600$, $J/J_c=0.29$, and (d) $t_s=800$, $J/J_c=0.37$. (e–h) Evolution of channels observed in the simulations for the nonuniform injection case. The ramp rate ${\alpha }_J=0.1\mathrm{mm}{\mathrm{s}}^{-2}$. (e) $t_s=250$, $J/J_c=0.13$, (f) $t_s=500$, $J/J_c=0.24$, (g) $t_s=800$, $J/J_c=0.38$, and (h) $t_s=900$, $J/J_c=0.42$. (i–l) Evolution of channels observed in the simulations for the nonuniform injection case with $\mathrm{\Delta }\varphi =0$. The ramp rate ${\alpha }_J=0.1\mathrm{mm}{\mathrm{s}}^{-2}$. (i) $t_s=250$, $J/J_c=0.13$, (j) $t_s=500$, $J/J_c=0.24$, (k) $t_s=800$, $J/J_c=0.38$, and (l) $t_s=900$, $J/J_c=0.42$. The channel is not observed to bifurcate under the same driving conditions when the volume fraction is uniform.

Figure 6

The length $l_e$ and the width $w_e$ of the main channel in the experiments shown in Fig. 4 as a function of flux $J/J_c$ in the case of uniform and nonuniform injection. (b) The length $l_s$ of the main channel in the five different simulations each shows stick-slip motion as in the experiments, but increases faster than in the experiments (see text). (c) The length and width of the channels averaged over the five simulations for uniform and nonuniform injection cases grow similarly, although somewhat systematically faster in the uniform injection case as the length grows.

Figure 7

Evolution of erosion patterns starting from various fluid flux $J_o$. (a) $J_o/J_c=0$; $t=640\mathrm{s}$, 1600 s, 3200 s, and 4000 s. (b) $J_o/J_c=0.22$; $t=320\mathrm{s}$, 500 s, 1440 s, and 1920 s. (c) $J_o/J_c=0.44$; $t=320\mathrm{s}$, 480 s, 800 s, and 1120 s. (d) $J_o/J_c=0.66$; $t=320\mathrm{s}$, 380 s, 800 s, and 1120 s. Higher initial flux gives rise to more branched patterns.

Figure 8

(a) The perimeter of the region where erosion has occurred and ${\varphi }_g<0.25$ compared to its area ($L=28.5\mathrm{cm}$). The ratio is observed to decrease corresponding to growth of linear structures, i.e., channels. (b) The total number of channel heads as a function of length of the longest channel is observed to grow with increasing $J_o/J_c$. The rate of growth of number of channels is observed to decrease as the channel length grows.

Figure 9

(a) The positions of the main channel heads over time corresponding to $J_o/J_c=0.44$ shown in Fig. 7 plotted in time increments of 12 seconds ($L=28.5\mathrm{cm}$). The tracks show that the channels grow radially outward initially. The ones closer to the water source are observed to grow further and bend toward the water source. (b) The evolution of the corresponding angle $\alpha$ between the direction that the channels grow and the $y$ axis as a function of the length $l_s$ of the longest channel.

Figure 10

Evolution of a branched channel when grains are only allowed to leave near the center of the bottom side boundary but fluid is allowed to flow out everywhere along the bottom boundary in the case of uniform injection (a–d) and nonuniform injection (e–h). The imposed flow corresponds to (a) $J/J_c=0.45$, (b) $J/J_c=0.55$, (c) $J/J_c=0.61$, and (d) $J/J_c=0.69$. (h) (a) $J/J_c=0.45$, (b) $J/J_c=0.57$, (c) $J/J_c=0.69$, and $J/J_c=0.72$.

Figure 11

(a) The length of the channels $l_s/L$ as a function of imposed flow rate $J/J_c$ for the five example simulations. The lengths increase steadily above a critical $J/J_c$ with small variations from run to run with the different medium initializations. (b) The averaged $\langle l_s/L\rangle$ over the five runs plotted after a time when the channel is initiated is observed to nearly overlap for uniform and nonuniform injection. (c) The channel length growth rate is observed to increase initially and then remain more or less constant as the channel length increases before decreasing somewhat.

Figure 12

(a) The number of channel heads $N_h$ as a function of length of the longest channel $l_s$ in the case of uniform and nonuniform injection. (b) Probability distribution function (PDF) of the angle between branches $\beta$ in the case of uniform and nonuniform injection.

Figure 13

(a) An example of an eroded channel observed in a simulation. (b) Corresponding filtered image where isolated eroded regions have been removed. (c) The skeleton of the channel network with bifurcation and channel heads identified. (d) Segments are used to join the branch points to the next branch points in the network or channel head. The angle between adjacent segments are use to identify branch angle $\beta$ as shown.