Microbial mats promote surface water retention in proglacial streams

The retreat of glaciers opens up large proglacial areas which become available for colonization and primary succession. Yet, factors that contribute to habitability during early succession in proglacial areas remain poorly understood. In proglacial streams, biofilms, which are matrix-enclosed microbial communities, colonize the streambed and grow into millimeter thick mats. Particularly in proglacial streams draining relatively flat and stable lateral terraces, these biofilms may augment habitability by reducing water scarcity through clogging of the streambed. To quantitatively address this phenomenon, we performed streamside flume experiments and conceived the idealized terrace model, which models stream length elongation as a function of microbially induced clogging, sediment hydraulic properties, stream roughness, slope, width and inflow. Significant stream elongation, and hence habitabilization, occurs when clogging suffices to induce unsaturated conditions below the streambed. Considering multiple terrace configurations with educated parameter bounds, we found a wide range of possible elongation, ranging from none to 100-fold. Sensitivity analysis suggests that sediment hydraulic properties mostly contribute to variability in stream elongation due to biofilm induced clogging. Taken together, we here show that microbial communities can significantly extend the habitability of proglacial stream ecosystems by inducing streambed clogging and retaining water. This is relevant in light of the rapid glacier retreat.

1 Introduction

Globally, glaciers are projected to lose 20 % to 60 % of their actual surface by the end of the century exposing large areas to colonization and primary succession (Bosson et al., 2023). Understanding the formation of newly deglaciated ecosystems is critical for anticipating consequences for their ecohydrology and ontogeny (e.g., Gurnell et al., 2000; Eichel, 2019). It has been postulated that geomorphological and ecological processes in proglacial terrain were mostly independent and that primary succession progressed linearly with time since deglaciation (Cooper, 1939; Matthews and Whittaker, 1987). This notion has recently been redressed, acknowledging the coupling of geomorphological and ecological processes (Ficetola et al., 2021). It is well established that disturbance regimes, soil moisture, sediment size and microclimate, for instance, affect the distribution of species within glacier forefields (e.g., Siegfried et al., 2023; Burga et al., 2010; Roncoroni et al., 2023b). However, there is also building evidence that certain organisms can augment ecosystem habitability by altering the physical environment (Eichel, 2019; Miller and Lane, 2019; Roncoroni et al., 2019). Miller and Lane (2019) synthesized previous models (from Matthews (1992) and Corenblit et al. (2007)) into a theoretical framework of primary succession accounting for environmental stress and the balance between biotic and abiotic controls in glacier forefields (Fig. 1a). This model identifies four ecosystem phases: (1) in the geomorphic phase, abiotic factors preclude life settlement; (2) in the pioneer phase, abiotic factors dominate and limit life to pioneer species; (3) in the biogeomorphic phase, biotic and abiotic factors contribute equally; and finally (4) in the ecological phase biotic factors dominate. Depending on severity and frequency of environmental stress, primary succession can remain limited to the geomorphic phase, or reset to the pioneer phase. Under limited stress, primary succession in proglacial landscapes can transit towards biogeomorphic and ecological controls.

Following incision–aggradation mechanisms (Maizels, 2002), terraces often form between the active floodplain and lateral moraines (both typical elevated stress environments) and offer relatively stable terrain compatible with primary succession. Surface water retention is very limited on these terraces making water scarcity the main stress. Nascent soils are characterized by high permeability as they are composed of poorly sorted sediments of heterogeneous sizes (Burga et al., 2010). Furthermore, the elevated location of terraces prevents capillary water rise (Siegfried et al., 2023). High permeability coupled with disconnection results in large infiltration rates (Brunner et al., 2009), up to of a few meters per hour (Müller et al., 2022), strongly limiting surface water extent. Consequently, water accessibility exerts a strong control on spatial distribution of vegetation. For instance, Siegfried et al. (2023) found that drought tolerant pioneer species could colonize the whole terrace albeit sparsely, while early successional plant species were confined to the vicinity of surface water or to shallow water table locations.

Besides ponding water following rain events, surface water is mainly found within both groundwater-fed and glacier-fed channels (Robinson et al., 2015) (Fig. 1b), which can cover up to 10 % of proglacial terrace areas (e.g., Roncoroni et al., 2023b). Intermittency (that is, zero-flow days) is inherent to these terrace streams (i.e., draining proglacial terraces) with expansion–contraction cycles arising from the superposition of daily and seasonal melting cycles and stochastic rain events (e.g., Paillex et al., 2020; Malard et al., 1999; Llanos-Paez et al., 2025). In addition to the diversity of water sources, the heterogeneity of water storage and transmission times within different proglacial landforms contribute to the complexity of the spatiotemporal dynamics of these terrace streams (Langston et al., 2011).

Figure 1 (a) Miller and Lane (2019) succession model within proglacial forefields. Ecosystem phases (colored regions) are positioned relative to the importance of (bio-)geomorphic and ecological controls. Different succession pathways (solid black arrows) correspond to different environmental stress. Dashed black arrows indicate disturbances. (b) Conceptual representation of the effect of microbial mats on terrace streams. Expansion–contraction cycles are illustrated by the solid and dashed stream lines, which correspond to minimal and maximal inflows for bio-clogged and unclogged streams, respectively. (c) Terrace stream succession pathway.

It was recently suggested that microbial mats (or biofilms) coating the streambed may promote the habitability of these streams via fertilization, stabilization and impermeabilization (Miller and Lane, 2019; Roncoroni et al., 2019). After snow melt, microbial mats rapidly grow up to a few millimeters in thickness and accumulate organic matter and nutrients on an otherwise mineral streambed. This fertilization effect benefits biofilm development (Vincent and Howard-Williams, 1986), but also pioneering plants adjacent to the streams (Ciccazzo et al., 2016; Frey et al., 2013). Furthermore, microbial mats are rich in extracellular substances which increase the cohesiveness of fine sediments (e.g., Gerbersdorf and Wieprecht, 2015). Therefore, benthic microbial mats are thought to contribute to the stability of these streams, particularly in the face of flow fluctuations. Finally, microorganisms can contribute to streambed clogging through active growth but also trapping of fine sediments (Dubuis and De Cesare, 2023). Because of water scarcity caused by highly permeable sediments, clogging by microbial mats may extend the wetted perimeter of streambeds and even flow length, therefore expanding their spatial cover. For instance, experimental work showed that the growth of microbial mats significantly reduced infiltration in flumes that mimicked a small proglacial tributary stream (Roncoroni et al., 2023a). Because the impermeabilization effect strongly depends on the substrate, resolving the permeability of clogging layers is needed to predict surface water extension following biofilm development. However the permeability of microbial mats in streams remains poorly studied.

The aim of this study was to quantify potential effect of microbial mats on the linear extension of streams draining proglacial terraces. We asked the following questions: (1) What is the microbial mat permeability and does it vary spatially and temporally? In particular, how is it affected by intermittent desiccation? (2) Is the permeability small enough to extend the flow length of terrace streams ? If so, what would be the effect of the underlying sediments? (3) Can we quantify surface water extension at the scale of the terrace despite the complexity of intermittent inflows? To address these questions, we first carried out a series of flume experiments in two different Alpine catchments to monitor microbial mat growth and permeability. We also studied how permeability responds to mat desiccation under controlled experiments. Finally, we conceptualized a simplified model (referred to as terrace model, hereafter) to estimate clogging effect on surface water retention. Experimental results are used as input parameters in the terrace model.

2 Methods

2.1 Experiments

A streamside flume (1 m long) was designed to grow microbial mats under near natural conditions and monitor mat permeability, k_c, non-destructively. A middle section of the flume is drained by a box of side length L=20 cm filled with sediments. As streamwater was fed into the flume, microbial mats started colonizing the sediment interface. The mat thickness, w_c, was measured with a micrometer screw gauge (0.1 mm accuracy) by averaging over 15 points regularly covering the flume. Monitoring the infiltration rate q and the pressure profile throughout the box (with two columns of nine piezometers each), the mat hydraulic conductivity, K_c, and the sediment hydraulic conductivity, K_s, are obtained by solving the double Darcy equation

$$\begin{array}{}\text{(1)}& q=K_{\text{c}}\left(\mathrm{1}+{\displaystyle \frac{h^{\text{up}}-h_{\mathrm{I}}}{w_{\text{c}}}}\right)=K_{\text{s}}\left(\mathrm{1}+{\displaystyle \frac{h_{\mathrm{I}}-h^{\text{dn}}}{L}}\right),\end{array}$$

where h^up is the water depth in the flume (above the box) and h^dn the water column imposed at the box outlet (Fig. 2b). Assuming uniformity of both sediments and mat, an algorithm was developed to find the best value of the interface pressure head h_I (between mat and sediments) with conservative uncertainty range based on the piezometer data (see Sect. S2; Paccolat, 2025a). Finally, to cancel out the effect of temperature on hydraulic conductivity, K – water viscosity almost doubles in the range of temperature considered –, it is converted to permeability, k, according to the monitored temperature. Five experimental runs were carried out: v1 (17 d) and v2 (34 d) in an alpine catchment (Valsorey, Switzerland) and r1 (9 d), r2 (19 d), and r3 (20 d) in a subalpine catchment (Réchy, Switzerland). Sediments ranging from 0.1 to 1 mm (fine to coarse sand) were used for v1, v2, r1, and r2, while sediments ranging from 1 to 2 mm (very coarse sand) were used for r3. The two categories are referred to as fine and coarse sediments, respectively. The flow velocity above the sediments was maintained close to 1.5 cm s⁻¹ during the different runs.

How mat permeability responded to desiccation was investigated on nine sediment cores (10–15 cm, 3 cm ∅) sampled at the end of r2. In each core, a millimeter thick mat overlaid the sediment column. After being re-saturated, the cores were dried (30 °C) and triplicates were withdrawn during phase I of evaporation (after 4 d), three more at the transition to phase II (after 15 d), and the last three during phase II (after 25 d) (see Sect. S3 or Lehmann et al. (2008) for a definition of evaporation phases). Immediately after withdrawal, the cores were re-saturated and the column (sediment + mat) permeability, k(t), was measured during eight hours using the falling head method (Klute and Dirksen, 1986). Because the saturation procedure was applied on wet sediments, saturation was incomplete and air entrapment altered permeability. A linear drift, $k_{\text{sed}}\left(t\right)=\mathrm{47}\left(\mathrm{8}\right)\times {\mathrm{10}}^{-\mathrm{6}}+\mathrm{0.6}t$ mm², was observed on a series of barren sediment cores (see Sect. S3). The evolution of the mat resistivity, $r_{\mathrm{c}}=w_{\mathrm{c}}/k_{\mathrm{c}}$, after desiccation is thus given by

$$\begin{array}{}\text{(2)}& r_{\mathrm{c}}\left(t\right)={\displaystyle \frac{l}{k\left(t\right)}}-{\displaystyle \frac{l}{k_{\text{sed}}\left(t\right)}},\end{array}$$

where l is the core height.

2.2 Terrace model

The effect of clogging on surface water retention in terrace streams was investigated with the following idealized terrace model. The terrace is assumed infinitely long such that streams eventually fade out because of infiltration. Each channel is assumed uniform in width, W, slope, s, and roughness, n_M (Manning coefficient). To account for the variability of topographic and pedologic conditions, the following (educated) bounds are used: 0.5 m < W < 2 m, ${\mathrm{10}}^{-\mathrm{4}}<s<{\mathrm{10}}^{-\mathrm{1}}$ and 0.1 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$ < n_M < 0.1 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$. They are chosen to reflect the variability of proglacial streams. Because no specific prior knowledge is available, we consider that channel width and Manning roughness follow uniform distributions; whereas to recognize the important fraction of flat conditions, the slope is assumed log-uniformly distributed. The hydraulic properties of the sediments underlying a channel are characterized by their hydraulic conductivity, K_a, and the van Genuchten–Mualem (vG–M) shape and scale parameters, n and h_g (van Genuchten, 1980). Two distinct joint probability distribution are considered: the sand distribution of Carsel and Parrish (1988) and the floodplain distribution. The latter is introduced because the reported hydraulic conductivity values in glacial floodplains (Müller et al., 2022) are mainly above the upper range of the sand distribution. For the floodplain distribution, K_a is taken log-uniformly between $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ and $\mathrm{8}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹ (i.e., range reported in Müller et al. (2022)) and the associated vG–M parameters are computed according to the model of Peche et al. (2024).

2.2.1 Stream length

For a constant discharge W Q_in (in m³ s), the stream length is given by

$$\begin{array}{}\text{(3)}& L=\underset{\mathrm{0}}{\overset{Q_{\text{in}}}{\int }}{\displaystyle \frac{\mathrm{d}Q}{q\left(h^{\text{up}}\left(Q\right)\right)}},\end{array}$$

where q is the infiltration rate (in m s⁻¹) which depends on the ponding depth h^up itself found from the Manning rating curve $Q\left(h^{\text{up}}\right)=h^{\text{up}}{R}^{\mathrm{2}/\mathrm{3}}\sqrt{s}/n_{\mathrm{M}}$, with $R=h^{\text{up}}W/(\mathrm{2}{h}^{\text{up}}+W)$ the hydraulic radius of a rectangular channel. The stream inflow condition is parametrized by the inlet stage, h_in, taken randomly between 0 and 0.3 m. The above integral solution relies on the assumption of steady-state conditions: the inflow must be constant during a period longer than the time required for the surface front to propagate to the stream end (where it fades out) and the transient time of infiltration. It is also assumed that the tributary streams are disconnected from the groundwater, yielding an infiltration rate independent of the water table position (Brunner et al., 2009). Note that transient groundwater movement, such as mounding, is admissible as long as it does not alter stream disconnection.

The infiltration rate through a channel uniformly clogged by a layer of thickness, w_c, and hydraulic conductivity, K_c, is given by the one-dimensional equation (Fox and Durnford, 2003)

$$\begin{array}{}\text{(4)}& q_{\mathrm{c}}=K_{\mathrm{c}}\left(\mathrm{1}+{\displaystyle \frac{h^{\text{up}}-h_{\mathrm{I}}}{w_{\mathrm{c}}}}\right)=K_a(-h_{\mathrm{I}}),\end{array}$$

in case of unsaturated streambed and by

$$\begin{array}{}\text{(5)}& q_{\mathrm{c}}=K_{\mathrm{c}}\left(\mathrm{1}+{\displaystyle \frac{h^{\text{up}}-h_{\mathrm{I}}}{w_{\mathrm{c}}}}\right)=q_{\mathrm{u}}\left(h_{\mathrm{I}}\right),\end{array}$$

in case of saturated streambed, where h_I is the pressure head at the interface between both layers (negative for Eq. 4 and positive for Eq. 5) and K_a(⋅) is the unsaturated hydraulic conductivity of the terrace sediments. The unsaturated one-dimensional infiltration (Eq. 4) is a good approximation as long as $h^{\text{up}}<h_{\mathrm{c}}=w_{\mathrm{c}}(K_{\mathrm{c}}/K_a-\mathrm{1})$. For larger ponding depths (or equivalently weaker clogging), an inverted water table is formed beneath the stream (Brunner et al., 2009; Xian et al., 2017) and lateral flow must be considered. The equation $q_{\mathrm{u}}\left(h^{\text{up}}\right)=K_a(\mathrm{1}+(\mathrm{2.7}{h}^{\text{up}}/W{)}^{\mathrm{0.77}}{)}^{\mathrm{1.3}}$ is (a good fit of) the disconnected infiltration rate without clogging (Swamee et al., 2000). The numerical integration of Eq. (3) coupled with Eqs. (4) and (5) is referred to as the exact stream length L_c. In the absence of clogging the notation L_u is used.

It is shown in Paccolat and Battin (2025) that the solution to Eq. (4) is well approximated by the linear expression

$$\begin{array}{}\text{(6)}& q_{\mathrm{c}}\left(h^{\text{up}}\right)\approx q_{\mathrm{0}}+h^{\text{up}}/R_{\mathrm{c}},\end{array}$$

where $R_{\mathrm{c}}=w_{\mathrm{c}}/K_{\mathrm{c}}$ is the clogging resistance (related to the resistivity by $R_{\mathrm{c}}=r_{\mathrm{c}}\mathit{\nu }/g$, with g the gravitational acceleration and ν the water kinematic viscosity) and q₀=q_c(0) is the channel infiltrability. This work also reveals asymptotic solutions associated to three clogging categories, namely

$$\begin{array}{}\text{(7)}& q_{\mathrm{0}}\to \left\{\begin{array}{ll}{K}_a,& \mathrm{\text{for }}{R}_{\mathrm{c}}\ll R_{-}\\ & \text{(negligible clogging)}\\ K_a{\left(R_{\mathrm{c}}/R_{-}\right)}^{-\mathit{\xi }},& \mathrm{\text{for }}{R}_{-}\ll R_{\mathrm{c}}\ll R_{+}\\ & \text{(soft clogging)}\\ K_{\mathrm{c}},& \mathrm{\text{for }}{R}_{\mathrm{c}}\gg R_{+}\\ & \text{(hard clogging)}\end{array}\right.,\end{array}$$

with $\mathit{\xi }=\frac{\mathrm{5}n-\mathrm{1}}{\mathrm{5}n+\mathrm{1}}$. The transition resistances $R_{-}=h_g/K_a$ and $R_{+}=R_{-}(K_a/K_{\mathrm{c}}{)}^{(\mathrm{5}n+\mathrm{1})/(\mathrm{5}n-\mathrm{1})}$ delimit the three regimes: increasing the clogging resistance, the infiltration is first controlled by the underlying sediments (negligible clogging), then depends on both layers (soft clogging) and is finally controlled by the clogging layer only (hard clogging). Using the approximate Eq. (6) and assuming wide channel conditions (W≫h^up), Eq. (3) integrates to the approximate clogged stream length (see Appendix A)

$$\begin{array}{}\text{(8)}& L_{\mathrm{c}}^{\text{app}}={\displaystyle \frac{Q_{\text{in}}}{q_{\mathrm{0}}}}{g}_{\mathrm{c}}\left({\displaystyle \frac{h_{\text{in}}}{q_{\mathrm{0}}{R}_{\mathrm{c}}}}\right).\end{array}$$

The hypergeometric function g_c decreases from 1 to 0.5 for arguments growing from 0 to 2. By definition, this solution only applies if h_in<h_c¹. Similarly, assuming wide channel conditions (W≫h^up) for the unclogged configuration, Eq. (3) integrates to the approximate unclogged stream length (see Appendix A)

$$\begin{array}{}\text{(9)}& L_{\mathrm{u}}^{\text{app}}={\displaystyle \frac{Q_{\text{in}}}{K_a}}{g}_{\mathrm{u}}\left({\displaystyle \frac{h_{\text{in}}}{W}}\right),\end{array}$$

where the hypergeometric function g_u decreases from 1 to 0.5 for arguments growing from 0 to 0.5.

2.2.2 Elongation factor

Clogging may increase stream length by reducing infiltration. This effect is quantified by the elongation factor, the ratio of clogged to unclogged stream length, i.e., $\text{EF}=L_{\mathrm{c}}/L_{\mathrm{u}}$. At the scale of the terrace, we introduce the multi-stream elongation factor

$$\begin{array}{}\text{(10)}& {\text{EF}}_m={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^m{L}_{\mathrm{c}}^{\left(i\right)}}{{\sum }_{i=\mathrm{1}}^m{L}_{\mathrm{u}}^{\left(i\right)}}},\end{array}$$

where m is the number of streams draining the terrace. $L_{\mathrm{c}}^{\left(i\right)}$ and $L_{\mathrm{u}}^{\left(i\right)}$ are the clogged and unclogged stream lengths associated to the ith inlet of ponding depth $h_{\text{in}}^{\left(i\right)}$. Confluent and divergent stream branching will reduce, respectively, increase the overall streams length (see Appendix A). However, because clogging only slightly alters ponding depth, the latter corrections almost cancel out from the elongation factor and are thus discarded.

2.2.3 Input parameter sets

To account for the variability of spatiotemporal factors such as inflows, pedology and topography, stream properties (i.e., W, s, n_M and h_in) are drawn from independent distributions introduced in the previous subsubsections. Clogging properties (i.e., w_c and K_c) are drawn from independent distributions characterized from experiments discussed in Sect. 3. The sediment hydraulic properties (i.e., K_a, h_g and n) are drawn from joint distributions, either the sand or the floodplain distribution introduced previously. Altogether, the distribution of the nine input parameters is referred to as the sand parameter distribution (SPD) or the floodplain parameter distribution (FPD), depending on the considered distribution of sediment hydraulic properties, as summarized in Fig. 4a. Note that despite the randomness of input parameters each stream is described by a fixed set of parameters (uniform conditions).

3 Experimental results

Similar microbial mat evolution was observed during the five flume experiments (Fig. 2a): sediment interstices were first filled within a few days before the mat grew homogeneously over the sediment bed to reach a maximal thickness of roughly 5 mm within 15 d. Because of technical issues, r1 was stopped during the growth phase. Mat permeability remained within a relatively small range, well captured by a 10-fold window ($\mathrm{6}\times {\mathrm{10}}^{-\mathrm{6}}$ mm² < k_c < $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ mm²), with an inter-run variability comparable to intra-run variability. The initial reduction of permeability suggests that the mat grew denser while colonizing the sediments, to stabilize later to a more constant structure. Intense phototrophic activity resulted in the formation of oxygen bubbles within and on top of the microbial mat forming vertical, tube-like connections through the mat. Particularly large bubbles (up to a few millimeters) were observed during v1, v2 and r3. The buoyancy force resulted in large bubbles escaping and perforating the mat on day 5 of v2, which corresponds to the observed sudden reduction of permeability. Note that only measurements taken at dawn were considered to discard diel permeability changes. Roncoroni et al. (2023a) conducted a similar experiment with faster flow and coarser sediments. The authors report a similar microbial mat growth reaching a maximal thickness of 8(3) mm within 10 d. With a simplified flow path calculation, we estimated, from the reported infiltration rates, that the clogging layer permeability was $\mathrm{3.2}\left(\mathrm{1.3}\right)\times {\mathrm{10}}^{-\mathrm{7}}$ mm² (see Appendix B).

Figure 2 (a) Evolution of the mat thickness on coarse (top) and fine (middle) sediments, and of the mat permeability (bottom), during all flume experiments. The solid line and the shaded area respectively correspond to the best fit of the logistic function and the uncertainty arising from varying the fit parameters within two standard deviations. The dashed lines delimit the representative 10-fold window. The yellow band corresponds to the range of values obtained from Roncoroni et al. (2023a). (b) Sediment box and pressure head profile. (c) Evolution of the mat permeability after desiccation. Colors indicate mat dryness before resaturation.

Permeability measurements rely on sediment uniformity, which is corroborated by the algorithm: despite a reduction in sediment permeability (by 10 %–30 % for fine and 80 % for coarse sediments), the goodness-of-fit improved over time for all experimental runs (see Sect. S2). Permeability reduction is explained by the formation of trapped bubbles in the sediment matrix due to outgassing from air saturated water. Also, we were unable to quantify organic matter (as ash free dry mass) within the cores sampled at the end of r2 and r3.

Dried mats were characterized by the evaporation depth (integrated evaporation rate), E, at withdrawal normalized by the evaporation depth, E_c, at the transition to the second phase of evaporation (when no capillary water is left in the mat). The permeability of dry mats (E>E_c) was significantly (Welch's t-test: p=0.00, t=4.29, df=18) larger than the one of wet mats (E<E_c), but all values remained within the 10-fold window obtained from the flume experiments (Fig. 2c). Crack formation in dry mats probably explain the permeability rise. Permeability reduction of wet mats following re-saturation may be attributed to remobilized particulates.

4 Terrace model results

4.1 Stream length

Stream length, L, is defined in Sect. 2.2 as the distance over which water flows before it infiltrates underground. Figure 3 illustrates how L increases as clogging reduces infiltration for different inlet depths and clogging hydraulic conductivities. The hydraulic properties of the terrace sediments matches the mean of the floodplain distribution (i.e., $K_a=\mathrm{1.9}\times {\mathrm{10}}^{-\mathrm{3}}$ m s⁻¹, h_g=0.14 m and n=2.74), while the channel slope, width and roughness are respectively set to $s={\mathrm{10}}^{-\mathrm{3}}$, W=1 m and n_M=0.022 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}/\mathrm{3}}$. For all combinations, three phases can be identified. When clogging is too small to cause unsaturated condition in the underlying sediments (h_in>h_c), the stream length remains close to its unclogged reference value. For intermediate resistances, it scales as $L\sim R_{\mathrm{c}}^{\mathit{\xi }}$, with $\mathit{\xi }=(\mathrm{5}n-\mathrm{1})/(\mathrm{5}n+\mathrm{1})$, as predicted by the approximate solution for soft clogging (Eq. 8). Finally, for resistances larger than R₊ (hard clogging), the stream length converges to a maximal value which depends on K_c.

Figure 3 Stream length as function of clogging resistance for different inlet depth (colors) and clogging hydraulic conductivities (line styles). The exact and approximate clogged channel length L_c and $L_{\mathrm{c}}^{\text{app}}$ are respectively represented by the thick colored lines and the thin black lines. The saturation transition, h_in=h_c, is marked by white dots and curves are shaded in the saturated region. The color triangles correspond to the values associated to the unclogged limit. The color dots indicate the limit of application of the approximate formula. The crosses mark the transition point R₊. The predicted scaling is outlined by the triangle slope (ξ≈0.86).

4.2 Single-stream elongation factor

To evaluate the extension of terrace streams due to mat induced clogging, the single stream (i.e., m=1) elongation factor (Eq. 10) was computed for N=10 000 terrace configurations (0.03 (SPD) and 0.9 (FPD) margin of error on the mean elongation factor at 95 % confidence interval) drawn from each input parameter distribution introduced in Sect. 2.2, i.e., SPD and FPD. The results from the previous section suggest the following bounds on the mat thickness and its hydraulic conductivity: 4 mm < w_c < 10 mm and $\mathrm{1.3}\times {\mathrm{10}}^{-\mathrm{6}}$ m s⁻¹ < K_c < $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{5}}$ m s⁻¹. Values were respectively drawn from the uniform and log-uniform distributions. The input distributions are summarized in Fig. 4a and the output EF distributions are shown on Fig. 4b. For the SPD, the median (IQR) is 1.3 (1.1–2.2) and only 35 % of the configurations are unsaturated. For the FPD, the median (IQR) is 7.2 (2.4–24.2) and 85 % of the configurations are unsaturated. The fraction of streams shorter than a given (cutoff) terrace size is represented in Fig. 4c for both input distributions. Considering finite terrace sizes, the EF distributions are shifted toward smaller values, as exemplified on Fig. 4b.

Figure 4 Single stream elongation factor (EF) for both SPD (violet) and FPD (yellow). (a) Input distributions. Bounds on the (log)-uniform distributions of the stream parameters and joint distributions of the sediment hydraulic properties. (b) Distribution of the elongation factor. Solid lines correspond to infinite terraces, while dashed and dash-dotted lines correspond to finite terraces (resp. of length 1000 m and 100 m). The pie charts indicate the prevalence of unsaturated configurations, as well as of the different clogging categories. (c) Cumulative distribution function of stream lengths. (d) Sensitivity analysis (2¹² samples). S₁ and S_T are respectively the first and total order Sobol indices. Error bars display the 95 % confidence intervals.

The exact (numerical integration) and approximate (analytical) stream length solutions are compared in Fig. 3. The approximation only applies to unsaturated conditions as it relies on the approximate linear formula for unsaturated infiltration. Using the analytical solutions, the elongation factor takes the simple expression

$$\begin{array}{}\text{(11)}& \text{EF}\approx {\displaystyle \frac{K_a}{q_{\mathrm{0}}}}{\displaystyle \frac{g_{\mathrm{c}}(h_{\text{in}}/(q_{\mathrm{0}}{R}_{\mathrm{c}}\left)\right)}{g_{\mathrm{u}}(h_{\text{in}}/W)}},\end{array}$$

which is compared to the exact solution in Fig. 5. Consistently, it is only reliable for unsaturated conditions. In particular, the relative error vanishes with increasing EF values, since the linear infiltration approximation improves for larger clogging resistances (Paccolat and Battin, 2025).

Figure 5 Validation of the approximate (analytical) EF expression for both SPD (violet) and FPD (yellow). Dashed lines indicate the identity function. Saturated configurations (red dots) are discarded in the relative error computation.

For both input parameter distributions, the two most influential factors are the sediment hydraulic properties (grouped as a single factor due to the variable dependence) and the clogging hydraulic conductivity: for the FPD (SPD), the sediments explain 43.5 % (16.2 %), K_c 17.6 % (50.4 %) and the interaction between both factors 30.7 % (16.4 %) of the EF variance. First- and total-order Sobol indices of the different factors are shown on Fig. 4d. The inlet depth strongly influences stream length (see Fig. 3), but its effect on the elongation factor is highly reduced as it almost factorizes out. Equation (11) reveals that h_in only impacts EF through the 𝒪(1) term $g_{\mathrm{c}}/g_{\mathrm{u}}$. In particular, the sign and magnitude of the effect depend on the stream width and clogging. Also, the elongation factor is independent of the slope and the roughness coefficient since they factorize out from the stream length integral (i.e., $L\sim \sqrt{s}/n_{\mathrm{M}}$).

4.3 Multi-stream elongation factor

The multi-stream system also exhibits the three phase behavior as illustrated for the 10-stream elongation factor in Fig. 6. However, the soft clogging regime no longer converges to a power law. Indeed, in general ${\text{EF}}_m\sim {\sum }_{i=\mathrm{1}}^m{A}_i{R}_{\mathrm{c}}^{{\mathit{\xi }}_i}$, with A_i random variables depending on the stream, clogging and sediment properties. However, due to the small support of ξ, the apparent scaling $\text{EF}\sim R_{\mathrm{c}}^{{\mathit{\xi }}_m}$ is observed. The distribution of ξ_m is different from the one of ξ because of two reasons. First, only a fraction of the considered configurations falls into the soft clogging category (see Fig. 4b). And second, even if considering larger clogging resistances corresponding to the maximal logarithmic slope (gray distributions in Fig. 6), the asymptotic behavior is not fully established because $K_a/K_{\mathrm{c}}$ is not diverging, yielding smaller effective exponents. Nonetheless, with 81 % of soft clogging configurations and a $K_a/K_{\mathrm{c}}$ median (IQR) of 126 (39–394), the FPD results in a fairly good scaling and exponents mostly between $\mathrm{2}/\mathrm{3}$ and 1. Besides altering the scaling behavior, the presence of multiple streams averages out the elongation factor. It follows from Jensen's inequality that $\mathbb{E}\left[{\text{EF}}_m\right]\ge \mathbb{E}\left[L_{\mathrm{c}}\right]/\mathbb{E}\left[L_{\mathrm{u}}\right]=\mathbb{E}\left[{\text{EF}}_{\mathrm{\infty }}\right]$. Numerical computations (data not shown) indicate further that $\mathbb{E}\left[{\text{EF}}_m\right]>\mathbb{E}\left[{\text{EF}}_{m+\mathrm{1}}\right]$. In other words, the clogging effect decreases with the number of inlets in the stream network.

Figure 6 Elongation factor for 1000 multi-streams systems (m=10) as function of the clogging resistance. All parameters are drawn from the terrace model input distributions (Fig. 4a), but w_c which scans a wider range. The thick intervals indicate the mat thickness values from the terrace model. The slope distributions are highlighted in the inset plots (max slope in gray and mean thick interval slope in color). The red band corresponds to the range of ξ values from the considered sediment distribution.

4.4 Model assumptions

The above results rely on the assumption of stream disconnection and steady-state. Rough estimates of the sum of the mound and capillary zone heights, discussed in Appendix C, indicate a median (IQR) of 15 (10–26) cm for the FDP and 18 (26–37) cm for the SDP. Similarly, the transient time (surface front propagation plus transient infiltration time) is quantified in Appendix C, yielding a median (IQR) of 1.1 (0.7–1.8) h for the FPD and 0.4 (0.2–0.8) h for the SPD.

5 Discussion and conclusion

Combining experimental observations with idealized modeling, we here show how microbially induced clogging contributes to the habitability of proglacial landscapes. While stream length elongation may be considerable (i.e., up to 100 fold), sedimentary properties typically dictate the magnitude of this effect (i.e., geomorphic control). However, our work clearly indicates that microbially induced clogging can be ecologically relevant in proglacial landscapes – by retaining water available for pioneer and early successional plants growing adjacent to proglacial streams, and increasing habitat size and connectivity available for aquatic organisms. While the importance of sediment biostabilization in intertidal zones and fluvial systems has already been recognized (Paterson et al., 2018), the relevance of this impermeabilization process during early successional stages of proglacial ecosystems remained largely overlooked.

5.1 Microbial mat permeability

Our experimental results show that the microbial mat permeability varies between $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{6}}$ and $\mathrm{6}\times {\mathrm{10}}^{-\mathrm{5}}$ mm² independent of sediment type and experiment context (e.g., water source, meteorological conditions), and compares well to values estimated from Roncoroni et al. (2023a), i.e., $\mathrm{3.2}\left(\mathrm{1.3}\right)\times {\mathrm{10}}^{-\mathrm{7}}$ mm². These experiments provide a first sensible range of the permeability of microbial mats in streams, which improves on the eight orders of magnitude variability reported from model biofilm systems (Kurz et al., 2023). These values are however limited to artificial setups and field measurements are critically lacking. Our results on the resilience of the permeability of microbial mats to desiccation corroborate infiltration observations from Roncoroni et al. (2023a). Also the relatively short growth timescales (up to 15 d) and widely reported cases of fast mat recovery following droughts (e.g., Hawes and Howard-Williams, 1998; Romaní et al., 2013; Touchette et al., 2025) support the first order approximation of constant mat properties used in the terrace model.

5.2 Stream elongation

It was postulated that microbial mat clogging can increase surface water availability (Miller and Lane, 2019; Roncoroni et al., 2019). The terrace model provides a simplified description of this effect. Clogging is shown to significantly increase stream length if large enough to unsaturate the underlying sediments, i.e., if $h_{\text{in}}<w_{\mathrm{c}}(K_a/K_{\mathrm{c}}-\mathrm{1})$ (Brunner et al., 2009), which generally applies to the FPD but not to the SPD. For such configurations, the elongation factor is well approximated by the analytical solution (Eq. 11), which highlights the effect of the different stream and clogging properties. In case of soft clogging (prevalent for both parameter distributions), the elongation factor (but also the stream length) scales with the clogging resistance, i.e., $\text{EF}\sim (w_{\mathrm{c}}/K_{\mathrm{c}}{)}^{\mathit{\xi }}$, where $\mathit{\xi }=(\mathrm{5}n-\mathrm{1})/(\mathrm{5}n+\mathrm{1})\in [\mathrm{2}/\mathrm{3},\mathrm{1}]$ and n increases with the uniformity of the sediment sizes. For a network of m streams of different properties, an apparent scaling subsists (due to the small exponent range) with an effective exponent ξ_m<ξ. This observation offers a rule of thumb to predict surface water cover response to clogging: for a x-fold increase of the clogging resistance, the overall length of the streams increases by $x^{{\mathit{\xi }}_m}$, which may help anticipate expansion–contraction dynamics.

5.3 Uncertainty and sensitivity

Because of the inherent heterogeneity of the physical parameters involved in the terrace model (i.e., h_in, W, s, n_M, K_a, h_g, n, w_c and K_c) and their values being poorly constrained from field observations, the uncertainty analysis indicates elongation factors ranging from one (no effect) to more than one hundred. The wider range and larger values found with the FPD than with the SPD are explained from the hydraulic conductivities being larger for the former distribution. More precisely, they are sufficiently large to trigger unsaturation for most configurations, resulting in large clogging effect, while this is not the case for the SPD. Fixing the input parameter distribution, the sensitivity analysis showed that the sediment hydraulic properties are the first (FPD) or second (SPD) most influential parameters. Gathering both distributions in the same sensitivity analysis (e.g. using a random trigger to select the distribution (Puy et al., 2020)) would undoubtedly result in a strong dominance of the latter parameters on the EF variability. Reducing uncertainty on the sediment hydraulic properties is consequently of prime importance to improve clogging effect estimates.

Because the stream length is proportional to $\sqrt{s}/n_{\mathrm{M}}$, the elongation factor is independent of these two parameters. In reality, streambed roughness is altered by microbial mat growth such that EF should be corrected by a factor $n_{\mathrm{M},\mathrm{u}}/n_{\mathrm{M},\mathrm{c}}$, where the indices refer to unclogged and clogged conditions, respectively. For reduced roughness, as observed by Roncoroni et al. (2023a), the clogging effect would hence be accentuated.

Intermittent network dynamics are observed in glacial forefields (Paillex et al., 2020; Malard et al., 1999) in response to the multiple water inputs (i.e., from kryal, krenal and stormflow sources) and variable landform transmissivities. Resolving this complex behavior is, however, not necessary to study the effect of clogging since the elongation factor is only slightly influenced by the inflow. Nonetheless, the effect is slightly larger during low flow periods as quantified by the ratio $\text{EF}{|}_{h_{\text{in}}=\mathrm{20}\mathrm{cm}}/\text{EF}{|}_{h_{\text{in}}=\mathrm{2}\mathrm{cm}}$, whose median (IQR) is 0.7 (0.8–1.0) for the SPD and 0.7 (0.8–0.9) for the FPD (data not shown).

5.4 Validity of the terrace model assumptions

The terrace model relies on several idealizations (Sect. 2.2). First, groundwater disconnection is motivated by the relative elevation of the terrace with respect to the main glacier-fed streams, which control the position of the water table, as well as diverse field observations (e.g. Miller and Lane, 2019; Malard et al., 1999; Siegfried et al., 2023). Estimates of vadose zone and water mounding, for predicted infiltration rates, in Sect. 4.4, corroborate further this intuition – terraces elevation of a few decimeters would be sufficient for disconnection. However, these values strongly depend on the arbitrarily chosen terrace geometry and must be taken cautiously. Second, the steady-state assumption requires stream response to be fast compared to inflow variations. Time to complete stream inundation is estimated around one hour in Sect. 4.4 for the parameters used in the study, such that response to inflow variations – even faster due to the smaller amount of water involved – may safely be assumed instantaneous in front of diurnal and seasonal timescales. Finally, parameters describing flow and infiltration (i.e., w_c, K_c, K_a, h_g, n, n_M, s and W) are assumed homogeneous along each stream. This is an arguably necessary first order approximation in front of our ignorance regarding spatial distributions. Also it permits an explicit formulation of stream length (i.e., Eq. 3) from the conservation equation, $\mathrm{d}Q/\mathrm{d}x=-q\left(h^{\text{up}}\right)$, by separation of variables. While specific studies could address the effect of in-stream heterogeneity, it is generally expected that the variability of stream lengths will be reduced if part or all of the between stream variability is present at the single stream scale.

5.5 Conclusion

Although highly idealized, the terrace model provides a mechanistic description of stream elongation, which refines our view on the bio-engineering potential of microbial mats and shows their transitional role linking the pioneer and the biogeomorphic ecosystem phases. The extent to which primary succession benefits from microbial mats remains, however, poorly quantified, in view of the variability of parameters and processes involved. Nevertheless, the terrace model indicates that for given stream conditions, impermeabilization is a relevant mechanism if (1) clogging is sufficient to unsaturate the underlying sediments, i.e., $h_{\text{in}}<w_{\mathrm{c}}(K_a/K_{\mathrm{c}}-\mathrm{1})$, and (2) the terrace is longer than the unclogged stream length, $L_{\mathrm{u}}\approx Q_{\text{in}}/K_a{g}_{\mathrm{u}}(h_{\text{in}}/W)$. Under these conditions, the effect can be quantified by the clogged stream length formula, $L_{\mathrm{c}}\approx Q_{\text{in}}/q_{\mathrm{0}}{g}_{\mathrm{c}}(h_{\text{in}}/q_{\mathrm{0}}{R}_{\mathrm{c}})$. Resolving the dynamics of streams draining proglacial terraces would require more realistic models accounting for additional processes (e.g., mat dynamics, coupling with water flow and nearby vegetation). However, model calibration is delusive with the current state of knowledge. Site specific, data-driven studies could provide field evidence for our predictions and set premises for further model development. In particular, drone-based imagery has shown promising results to simultaneously map flow and biofilms (Roncoroni et al., 2023b).

Appendix A: Approximate stream length and branching

A1 Stream length integral

The stream length integral (Eq. 3) can be solved analytically for the clogged unsaturated (h_in<h_c) and the unclogged configurations under the wide channel assumption (W≪h_in). In both cases, we make use of the equality

$$\begin{array}{}\text{(A1)}& {}_{\mathrm{2}}{F}_{\mathrm{1}}(\mathrm{1},b,b+\mathrm{1},z)=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}u}{\mathrm{1}-z{u}^{\mathrm{1}/b}}},\end{array}$$

which is a particular case of the Euler integral expression of the hypergeometric function. For clogged unsaturated profiles, injecting the approximate linear infiltration equation $q_{\mathrm{c}}\left(h^{\text{up}}\right)\approx q_{\mathrm{0}}+h^{\text{up}}/R_{\mathrm{c}}$ into the stream length integral yields

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{c}}^{\text{app}}& {\displaystyle }=\underset{\mathrm{0}}{\overset{Q_{\text{in}}}{\int }}{\displaystyle \frac{\mathrm{d}Q}{q_{\mathrm{0}}+h^{\text{up}}\left(Q\right)/R_{\mathrm{c}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{q_{\mathrm{0}}}}\underset{\mathrm{0}}{\overset{Q_{\text{in}}}{\int }}{\displaystyle \frac{\mathrm{d}Q}{\mathrm{1}+\frac{{\left(n_{\mathrm{M}}Q/\sqrt{s}\right)}^{\mathrm{3}/\mathrm{5}}}{h^{\star }}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_{\text{in}}}{q_{\mathrm{0}}}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}u}{\mathrm{1}+\frac{h_{\text{in}}}{h^{\star }}{u}^{\mathrm{3}/\mathrm{5}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_{\text{in}}}{q_{\mathrm{0}}}}{}_{\mathrm{2}}{F}_{\mathrm{1}}(\mathrm{1},\mathrm{5}/\mathrm{3},\mathrm{8}/\mathrm{3},-h_{\text{in}}/h^{\star })={\displaystyle \frac{Q_{\text{in}}}{q_{\mathrm{0}}}}{g}_{\mathrm{c}}(h_{\text{in}}/h^{\star }),\end{array}$$

where $h^{\star }=q_{\mathrm{0}}{R}_{\mathrm{c}}$ and $h_{\text{in}}={\left(n_{\mathrm{M}}{Q}_{\text{in}}/\sqrt{s}\right)}^{\mathrm{3}/\mathrm{5}}$. The wide channel assumption was used to simplify the Manning formula in the second line and the third line follows from the change of variable $u=Q/Q_{\text{in}}$. It is shown in Swamee et al. (2000) that seepage from rectangular channels into unconfined aquifer is well fitted by the function $q_{\mathrm{u}}\left(h^{\text{up}}\right)=K_a{\left[\mathrm{1}+{\left(a^{\prime }{h}^{\text{up}}/W\right)}^{\mathrm{0.77}}\right]}^{\mathrm{1.3}}\approx K_a\left[\mathrm{1}+a{\left(h^{\text{up}}/W\right)}^{\mathrm{0.77}}\right]$, where the approximation follows from the wide channel assumption and $a=\mathrm{1.3}{a}^{\prime \mathrm{0.77}}\approx \mathrm{2.79}$, with $a^{\prime }=\mathit{\pi }(\mathrm{4}-\mathit{\pi })$. This expression can be used to solve the unclogged channel length integral, by following the same steps as in the previous integral, namely

$$\begin{array}{ll}{\displaystyle }{L}_{\mathrm{u}}^{\text{app}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{K_a}}\underset{\mathrm{0}}{\overset{Q_{\text{in}}}{\int }}{\displaystyle \frac{\mathrm{d}Q}{\mathrm{1}+a\left(h^{\text{up}}\right(Q)/W{)}^{\mathrm{0.77}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_{\text{in}}}{K_a}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{d}u}{\mathrm{1}+a{\left(\frac{h_{\text{in}}}{W}\right)}^{\mathrm{0.77}}{u}^{\mathrm{0.462}}}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_{\text{in}}}{K_a}}{}_{\mathrm{2}}{F}_{\mathrm{1}}\left(\mathrm{1},\mathrm{2.164},\mathrm{3.164},-a(h_{\text{in}}/W{)}^{\mathrm{0.77}}\right)\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{Q_{\text{in}}}{K_a}}{g}_{\mathrm{u}}(h_{\text{in}}/W).\end{array}$$

The functions g_c and g_u are plotted in Fig. A1a.

A2 Branching

For a given set of inflows generating a stream network, the total network length depends on the presence and positions of stream branchings. If two streams of respective surface flow Q₁ and Q₂ merge into a single stream of surface flow $Q=Q_{\mathrm{1}}+Q_{\mathrm{2}}$, the length of the resulting stream L is smaller than the sum of the two independent streams (L₁ and L₂). To estimate this effect on the elongation factor, we define the (two streams) branching ratio as the ratio between the elongation factor with and without branching, namely

$$\begin{array}{}\text{(A2)}& B={\displaystyle \frac{L_{\mathrm{c}}^{\left(\mathrm{1}\right)}+L_{\mathrm{c}}^{\left(\mathrm{2}\right)}}{L_{\mathrm{u}}^{\left(\mathrm{1}\right)}+L_{\mathrm{u}}^{\left(\mathrm{2}\right)}}}\times {\displaystyle \frac{L_{\mathrm{u}}}{L_{\mathrm{c}}}}\approx {\displaystyle \frac{Q^{\left(\mathrm{1}\right)}{g}_{\mathrm{c}}^{\left(\mathrm{1}\right)}+Q^{\left(\mathrm{2}\right)}{g}_{\mathrm{c}}^{\left(\mathrm{2}\right)}}{Q^{\left(\mathrm{1}\right)}{g}_{\mathrm{u}}^{\left(\mathrm{1}\right)}+Q^{\left(\mathrm{2}\right)}{g}_{\mathrm{u}}^{\left(\mathrm{2}\right)}}}\times {\displaystyle \frac{g_{\mathrm{u}}}{g_{\mathrm{c}}}},\end{array}$$

where the approximation only applies to unsaturated conditions and $g_{\mathrm{c}}=g_{\mathrm{c}}(h^{\text{up}}/h^{\star })$, $g_{\mathrm{u}}=g_{\mathrm{u}}(h^{\text{up}}/W)$, $g_{\mathrm{c}}^{\left(i\right)}=g_{\mathrm{c}}(h_i^{\text{up}}/h^{\star })$ and $g_{\mathrm{u}}^{\left(i\right)}=g_{\mathrm{u}}(h_i^{\text{up}}/W)$ for $i=\mathrm{1},\mathrm{2}$. The ponding depth h^up and $h_i^{\text{up}}$ are computed with the Manning equation and the respective inflows Q and Q_i. Considering both parameter distributions, the branching ratio essentially remains between 1 and 1.1 (Fig. A1b). A stream may also split in two, resulting in a increased total length. In that case the branching ratio is inverted (i.e., B⁻¹) such that the branching effect is negligible at the multi-stream scale.

Figure A1 (a) Hypergeometric functions associated to the clogged (blue) and unclogged (red) streambeds. (b) Construction of the branching ratio and its distribution for both sand and floodplain parameter sets.

Appendix B: Estimation of microbial mat permeability from Roncoroni et al. (2023a)

During the flume experiment carried out by Roncoroni et al. (2023a), the infiltration flux dropped from Q_i=0.18(0.01) L s⁻¹ to Q_f=0.025(0.005) L s⁻¹ as a consequence of microbial mat growth. The sediment chamber length, height and width were respectively L_x=2 m, L_z=7 cm and W=30 cm, while the ponding depth was h^up=11.5 cm. The infiltrated water was drained through three pipes of diameter D_pipe=7 mm and unknown length L_pipe=10–20 cm (range estimated by the experimental designer upon request). In order to estimate the sediment permeability, k_s, and the clogging resistivity, r_c, we first compute the pressure at the pipe entrance according to the Darcy-Weisenbach equation

$$\begin{array}{}\text{(B1)}& h_{\text{out}}={\displaystyle \frac{L_{\text{pipe}}}{D_{\text{pipe}}}}{\displaystyle \frac{f_{\mathrm{D}}\left(Re\right)u^{\mathrm{2}}}{\mathrm{2}g}},\end{array}$$

where u is the mean flow through each pipe obtained from $Q_{\text{pipe}}=Q/\mathrm{3}=\mathit{\pi }(D/\mathrm{2}{)}^{\mathrm{2}}u$, ρ is the water density and f_D is the Darcy friction factor. The Reynolds numbers of the pipe flows are respectively Re_i=7200(400) and Re_f=1000(200) for the initial and final rates. The friction factors f_D(Re_i) and f_D(Re_f) are thus respectively computed from the turbulent smooth-pipe and the laminar expressions. We then compute the total flux Q according to the chamber dimensions and boundary conditions. It is given by the integral

$$\begin{array}{}\text{(B2)}& Q=W\underset{\mathrm{0}}{\overset{L_x}{\int }}q\left(x\right)dx,\end{array}$$

where q(x) is the flow rate along the streamline starting a distance x from the chamber downstream end. We solve this integral by assuming that the streamlines follow a right-angle path of length $L=x+L_z+w_{\mathrm{c}}$, where w_c is the clogging layer thickness (Fig. B1a). Consequently the flow rates are given by

$$\begin{array}{}\text{(B3)}& q\left(x\right)=K{\displaystyle \frac{h^{\text{up}}-h_{\text{out}}+L_z}{x+L_z+w_{\mathrm{c}}}}=K_{\mathrm{s}}{\displaystyle \frac{h^{\text{up}}-h_{\text{out}}+L_z}{x+L_z+r_{\mathrm{c}}{K}_{\mathrm{s}}}},\end{array}$$

where $K_{\mathrm{s}}=g/\mathit{\nu }{k}_{\mathrm{s}}$ is the sediment hydraulic conductivity and K is the average hydraulic conductivity, which are related by $K/(x+L_z+w_{\mathrm{c}})=K_{\mathrm{s}}/(x+L_z)+r_{\mathrm{c}}$.

The above integral hence yields

$$\begin{array}{}\text{(B4)}& Q(k_{\mathrm{s}},r_{\mathrm{c}})\approx {\displaystyle \frac{gb{k}_{\mathrm{s}}}{\mathit{\nu }}}\left(h^{\text{up}}-h_{\text{out}}+L_z\right)\log \left(\mathrm{1}+{\displaystyle \frac{L_x}{L_z+k_{\mathrm{s}}{r}_{\mathrm{c}}}}\right).\end{array}$$

The sediment permeability is found by solving this equation for the initial condition, i.e., $Q(k_{\mathrm{s}},\mathrm{0})=Q_i$. The clogging resistivity is then obtained by solving it for the final condition, i.e., $Q(k_{\mathrm{s}},r_{\mathrm{c}})=Q_f$. The results depend on the unknown pipe length L_pipe. Varying it from 10 to 20 cm results in a variability similar to the one arising from the uncertainty on the initial and final fluxes as shown in Fig. B1b. Altogether we find $k_{\mathrm{s}}=\mathrm{3.0}\left(\mathrm{0.6}\right)\times {\mathrm{10}}^{-\mathrm{4}}$ mm² and $r_{\mathrm{c}}=\mathrm{2.5}\left(\mathrm{0.3}\right)\times {\mathrm{10}}^{\mathrm{7}}$ mm⁻¹. With respect to the reported mat thickness of w_c=8(3) mm, we hence find $k_{\mathrm{c}}=\mathrm{3.2}\left(\mathrm{1.3}\right)\times {\mathrm{10}}^{-\mathrm{7}}$ mm².

Figure B1 (a) Sketch of the flume experiment from Roncoroni et al. (2023a) illustrating the considered simplified flow trajectory. (b) Dependence of the sediment permeability and the clogging resistance on the unknown pipe length. The shaded areas correspond to the uncertainty propagated from the flow rate uncertainties.

Appendix C: Assessment of the terrace model assumptions

C1 Groundwater mound

The water table rises below the stream to evacuate the flux laterally. For a two-dimensional symmetrical configuration, the linearized Boussinesq equation predicts the following water table level (Brunner et al., 2009)

$$\begin{array}{}\text{(C1)}& \mathrm{\Delta }{z}_{\text{WT}}\left(x\right)={\displaystyle \frac{q}{\mathrm{2}T}}(W{L}_a-W^{\mathrm{2}}-x^{\mathrm{2}}),\mathrm{\text{for }}\mathrm{0}<x<{\displaystyle \frac{W}{\mathrm{2}}},\end{array}$$

where q is the infiltration rate, x is the lateral distance from the centre of the stream of width W and L_a is the distance to the reference water table level. Considering an aquifer of depth B_a and transmissivity T=K_aB_a, the height of the saturated mound below the stream is $M_{\text{sat}}=qW(L_a-W/\mathrm{4})/\mathrm{2}{K}_a{B}_a$. On top of the mound, a capillary zone extends the groundwater zone of influence. It is estimated from the following linear approximation. At steady-state, the infiltration rate within the unsaturated zone is equal to the one at the water table, yielding

$$\begin{array}{}\text{(C2)}& q=K_a(-h_{\mathrm{I}})=K_a\left(\mathrm{1}-{\displaystyle \frac{\partial h}{\partial z}}{|}_{\text{WT}}\right),\end{array}$$

where h_I<0 is the constant pressure head below the clogged streambed. Assuming the pressure gradient to remain constant within the capillary zone, the size of the latter is approximated as

$$\begin{array}{}\text{(C3)}& w_{\text{cz}}\approx {\displaystyle \frac{-h_{\mathrm{I}}}{\partial h/\partial z{|}_{\text{WT}}}}={\displaystyle \frac{K_a^{-\mathrm{1}}\left(q\right)}{\mathrm{1}-q/K_a}},\end{array}$$

where $K_a^{-\mathrm{1}}(\cdot )$ is the inverse of K_a(⋅), such that $-h_{\mathrm{I}}=K_a^{-\mathrm{1}}\left(q\right)$. Considering B_aq=8 m and L_aq=15 m, the distribution of M_sat and w_cz are shown in Fig. C1a for both parameter distributions (SPD and FPD). Note that only the unsaturated configurations (i.e., 85 % for the FPD and 35 % for the SPD) are considered. For the saturated configurations, an inverted water table is formed below the stream which may significantly increase the required water table depth for disconnection (Xie et al., 2014; Xian et al., 2017).

Figure C1 Distribution of (a) the saturated mound height and the capillary zone height, and (b) the surface and infiltration transient times obtained from the input parameter distributions of Fig. 4a (FPD in yellow and SPD in violet).

C2 Transience

The problem of surface water propagation downward from a constant inflow Q_in was studied by Turbak and Morel-Seytoux (1988). A power law rating curve is used which corresponds to the Manning formula in the large width limit, namely $Q=\sqrt{s}(h^{\text{up}}{)}^{\mathrm{5}/\mathrm{3}}/n_{\mathrm{M}}$. The position of the front at time t and the time required for the front to reach the distance x are respectively given by

$$\begin{array}{ll}{\displaystyle }{x}_F\left(t\right)& {\displaystyle }={\displaystyle \frac{\sqrt{s}}{n_{\mathrm{M}}I\left(t\right)}}\left[(h_{\text{in}}{)}^{\mathrm{5}/\mathrm{3}}-{\left(h_{\text{in}}-{\displaystyle \frac{\mathrm{3}}{\mathrm{5}}}I\left(t\right)t\right)}^{\mathrm{5}/\mathrm{3}}\right],\\ {\displaystyle }{t}_F\left(x\right)& {\displaystyle }={\displaystyle \frac{\mathrm{3}}{\mathrm{5}I\left(x\right)}}\left[h_{\text{in}}-{\left((h_{\text{in}}{)}^{\mathrm{5}/\mathrm{3}}-{\displaystyle \frac{n_{\mathrm{M}}}{\sqrt{s}}}xI(x)\right)}^{\mathrm{3}/\mathrm{5}}\right],\end{array}$$

where I(t) and I(x) correspond to the averaged infiltration rate up to time t, respectively up to distance x. Here we merely observe that when the front has propagated to its maximal distance L after a time T_sur, one can set x_F(T_sur)=L, t_F(L)=T_sur and $I\left(T_{\text{sur}}\right)=I\left(L\right)=I_{\text{end}}$ in the above two equations. Solving them yields

$$T_{\text{sur}}={\displaystyle \frac{\mathrm{5}{h}_{\text{in}}}{\mathrm{3}{I}_{\text{end}}}}\mathrm{\text{and}}L={\displaystyle \frac{\sqrt{s}(h_{\text{in}}{)}^{\mathrm{5}/\mathrm{3}}}{n_{\mathrm{M}}{I}_{\text{end}}}},$$

which results in the relation $T_{\text{sur}}=\mathrm{5}{n}_{\mathrm{M}}L/\mathrm{3}\sqrt{s}(h_{\text{in}}{)}^{\mathrm{2}/\mathrm{3}}=\mathrm{5}L/\mathrm{3}{u}_{\text{in}}$. The transience time is increased by the infiltration relaxation time T_inf required for the infiltration rate to reach its maximal value. It is computed by solving the Richards equation for h_in=0 (i.e., stream end) until the infiltration rate q(t) stabilizes. The simulations is stopped when $\partial q/\partial t<\mathit{ϵ}q$, with $\mathit{ϵ}={\mathrm{10}}^{-\mathrm{3}}{\mathrm{s}}^{-\mathrm{1}}$. The distribution of T_sur and T_inf are displayed in Fig. C1b for both parameter distributions.