Numerical Simulation of Debris Flow Behavior Over a Series of Groundsills

1. Introduction

Debris flows, which are highly concentrated mixtures of water, soil, and debris can cause devastating damage as they propagate downslope, making their mitigation and control a necessary concern (Hungr et al., 2005, Han et al., 2015). In natural environments, debris flows rarely propagate on a flat, uninterrupted trajectory. Instead, they encounter complex and uneven terrain characterized by various obstacles and irregular topography. These irregularities can take the form of bottom roughness, topographic features such as humps and hilly terrains, or even man-made structures like check dams, groundsills, or other mitigation measures (Iverson, 1997). One of the established techniques for mitigating debris flows involves the strategic placement of groundsills, which are structural barriers designed to dissipate energy, trap sediment, and regulate flow propagation (Choi et al., 2018; Ng et al., 2014). When debris flows interact with bottom irregularities or man-made structures, the changes in topography introduce additional drag forces that decelerate the flow velocity. The presence of obstacles and uneven surfaces disrupt the continuous flow, causing a deceleration in the progression of the debris flow (Jacob et al., 2005). Consequently, the flow dynamics become inherently unpredictable and complex, as the flow dynamics are influenced by the interactions between the debris material and the varying bottom irregularities (Sand-Jensen, 1998). Moreover, the presence of larger solid particles in the mixtures may interact differently with the obstacles, which can significantly influence the sediment trapping between them (Han et al., 2024).

Recently, several researchers investigated gravity currents over a smooth rigid boundary moving down the slope (Britter and Linden, 1980), flow over or through obstacles (Hatcher et al., 2000), and flow over permeable surfaces (Marino and Thomas, 2002). To better understand how surface roughness affects fluid flow, experimental and theoretical studies have been conducted on various artificial roughness elements, leading to the classification of different ’types’ of roughness with effects comparable to those observed in natural environments (Rouse, 1965). One such approach is to use beam-type groundsill arrays, which are transverse structures installed on river channels to guide the flow, reducing velocity, mitigating erosion, and adjusting bed slope (Shima et al., 2016; Lin et al., 2024; Ikhsan et al., 2009). Ground sills are engineered structures designed to reduce the sediment-laden flow velocities and increase sediment deposition upstream while allowing flow to propagate over them (Lin et al., 2008). However, the practical implementation of ground-sills has revealed challenges in both design and operation, where the intended functionality often differs from the observations in the real field (Lin et al., 2008). For instance, studies have demonstrated that ground sills enhance sediment deposition upstream and stabilize channel morphology, but real-world performance has varied depending on the design and context. Addressing the complexity of ground-sill projects requires a more integrated strategy rather than fragmented solutions (Bagli, 2019). In Taiwan, debris flows frequently occur due to the mountainous terrain, heavy rainfall, and frequent typhoons (Jan and Chen, 2005). This has necessitated the implementation of various mitigation measures, including check dams, slit dams, and groundsills (Jan and Chen, 2005). After the 1990 debris-flow event, a series of check dams were constructed along a debris-flow gully in Hualien County to manage debris transport and reduce damage (Figure 1). These structures function similarly to beam-type groundsill arrays, which are transverse structures installed on river channels to guide the flow, reduce velocity, mitigate erosion, and adjust bed slope.

Research on flow over groundsills has mainly focused on Newtonian fluids, especially regarding erosion control and riverbed stabilization (Hairani and Legono, 2016). Recently, however, interest has grown in non-Newtonian fluids due to their complex flow characteristics and practical applications (Bigham, 2020). Observing non-Newtonian flow behavior experimentally is challenging because of their intricate rheological properties and flow dynamics. Consequently, we employed numerical simulation in this study to investigate the flow behavior of debris flow approximated by the Bingham Fluid Model with two rheological parameters, namely the Bingham yield stress (${\tau }_B$) and viscosity (${\mu }_B$). Debris flows can be typically classified into three main categories: mudflows, which consist primarily of fine materials such as clay and silt (Blasio, 2011), and stone- or boulder-rich debris flows dominated by larger particles (Julien & O’Brien, 1997). Between these two extremes, intermediate debris flows exist with varying sediment compositions and concentrations (Iverson and Denlinger, 1987). This study simulates the propagation of sediment mixtures modeled as the debris flows that follow the Bingham Fluid Model, as they move over a beam-type array of groundsills positioned within an inclined channel set at a 15° slope. The groundsill array density refers to the spacing between consecutive groundsills, represented by the ratio between the spacing and their relative height (i.e., $d/h$). It is classified as high array density ($d/h\le 4$) and low array density ($d/h>4$). Our main goal is to investigate how varying this groundsill array density influences the flow behavior of sediment mixtures including flow depth, front position of the flow, and flow velocity.

2. Materials and Methods

2.1. Sediment Mixture Preparation

The sediment mixtures are highly concentrated clay-silt-water mixtures having a sediment fraction $\varphi$. A volume, ${\forall }_c$ of sediment was well mixed with tap water of volume, ${\forall }_w$ to form a sediment mixture with a sediment fraction, $\varphi$, as shown in Eq. 1.

$$\varphi ={\displaystyle \frac{{\forall }_c}{{\forall }_c+{\forall }_w}}$$

(1)

The clay-silt materials were used as the sediment in the present experiments, which were collected from a reservoir sediment deposition. After washing them multiple times to remove any organic debris, the materials were oven-dried and stockpiled for experiments. As shown in Figure 2, the particle size distribution of the clay-silt materials yielded a medium diameter ($d_{50}$) of 0.0036 mm and a particle density of 2.65 g/cm³. The sediment mixtures prepared for this study had a sediment fraction $\varphi =0.30$.

2.2. Rheometer Setup and Measurement of Rheological Parameters

To characterize the rheological properties of the sediment mixtures, a conventional Brookfield rheometer (model DV-III) was utilized. This rheometer operates by rotating a spindle (RV-6) within the sediment mixture and measuring the torque required to overcome the viscous resistance during the spindle’s rotational motion. The spindle is immersed into the mixture contained in a 600 ml glass beaker, and the rheometer is controlled by a computer running the Rheocal-32 software. A schematic representation of the rheometer test, depicted in Figure 3a, follows the same approach we used in our previous studies (Jan and Dey, 2022; Dey et al., 2021).

The measurements were conducted in a controlled shear rate mode, where the rotational speed of the spindle was regulated, and the corresponding torque required was recorded. To ensure the repeatability and reliability of the results, each test was performed at least three times, and the final reported values are the average of these repeated measurements. These measurements were carried out at room temperature, approximately 25°C. The recorded spindle speed and torque values were then converted to shear rate and shear stress, respectively, using Mitschka’s method (Mitschka, 1982), as described by Dey et al. (2021). Subsequently, the calculated shear stress and shear rate data were used to characterize the rheological properties of the sediment mixture, as illustrated in Figure 3b. The relationship between shear stress and shear rate follows the Bingham Fluid Model, which is expressed mathematically in Eq. 2.

$${\tau =\tau }_B+{\mu }_B\dot{\gamma }$$

(2)

in which ${\tau }_B=$ Bingham yield stress (Pa), ${\mu }_B=$ Bingham viscosity (Pa.s), and $\dot{\gamma }=$ shear rate (1/s). The measured values of ${\tau }_B$ and ${\mu }_B$ for sediment mixture used in this study with $\varphi =0.30$, were 12.83 Pa and 0.28 Pa.s, respectively. The rheological parameters (${\tau }_B$and${\mu }_B$) of these Bingham fluids are functions of sediment fraction ($\varphi$) and sediment size ($d_s$), expressed as ${\tau }_{B}or{\mu }_B=f(\varphi ,d_s)$. The parameter $d_s$ represents the influence of larger particles within the mixture by changing its rheological parameters (Major and Pierson, 1992; Dey et al., 2021; Jan and Dey, 2022).

At low shear rates ($0\le \dot{\gamma }\le 4$), shear stress of sediment mixtures used in this study first increases and then decreases, potentially causing data scatter (Dey et al., 2021; Ancey and Jorrot, 2001). In contrast, at higher shear rates ($4\le \dot{\gamma }\le 20$), the shear stress-shear rate relationship aligns well with the Bingham fluid model. The Bingham fluid model is widely used in numerical simulations of debris flow movement not only for its simplicity but most of the natural debris flows exhibit a shear rate within the range of 5 to 20 s⁻¹ (Major and Pierson 1992; Jan et al., 2019; Dey et al., 2021). Therefore, we focus on the rheological properties of sediment mixtures tested in this study within the shear rate range of $4\le \dot{\gamma }\le 20$ s⁻¹.

3. Numerical Simulation Methodology

In this study, we utilized the commercial CFD software Flow-3D to simulate the flow behavior and dynamics of sediment mixtures flowing over a beam-type array of groundsills placed in an inclined channel. The simulation employs the principles of mass conservation (the continuity equation) and momentum conservation (Navier-Stokes momentum equation) to comprehensively characterize the complex dynamics of three-dimensional fluid movement. The fluid phase interface is specified through the utilization of the Volume-of-Fluid (VOF) method, which accurately monitors the interface between different fluids or between a fluid and a solid. The Flow-3D User Manual provides detailed descriptions and mathematical formulations of the mass continuity and momentum equations governing the fluid flow dynamics (Flow Science, 2008).

In this study, the Flow-3D solver uses the Navier-Stokes equations combined with the Volume of Fluid (VOF) method to simulate sediment mixture flow on an inclined channel. The Navier-Stokes equations incorporating the Bingham Fluid Model can be expressed as:

$${\displaystyle \frac{d{u}_i}{dt}}+u_j{\displaystyle \frac{d{u}_i}{d{x}_j}}=-{\displaystyle \frac{1}{{\rho }_m}}{\displaystyle \frac{dp}{d{x}_i}}+{\displaystyle \frac{d{\tau }_{ij}}{d{x}_j}}+g_i$$

(3a)

$${\displaystyle \frac{d{u}_i}{d{x}_i}}=0$$

(3b)

where $u_i$ are the components of the velocity vector u, $t$ is the time, ${\rho }_m$ is the density of the mixture, $p$ is the pressure, ${\tau }_{ij}$is the stress tensor for the Bingham Fluid Model, $g_i$ is the i-th component of the gravitational acceleration vector.

The stress tensor ${\tau }_{ij}$ for the Bingham Fluid Model is defined as:

$${\tau }_{ij}=\left\{\begin{array}{c}0if‖\dot{\gamma }‖<{\displaystyle \frac{{\tau }_B}{{\mu }_B}}\\ \left({\displaystyle \frac{{\tau }_B}{‖\dot{\gamma }‖}}+{\mu }_B\right){\dot{\gamma }}_{ij}if‖\dot{\gamma }‖\ge {\displaystyle \frac{{\tau }_B}{{\mu }_B}}\end{array}\right.$$

(4)

where ${\tau }_B$ is the Bingham yield stress, ${\mu }_B$ is the Bingham viscosity, ${\dot{\gamma }}_{ij}$are the components of the rate of strain tensor, $‖\dot{\gamma }‖$is the magnitude of the shear rate tensor. The above Eqs. ($3-4$) modeled the sediment mixtures tested in this study as a single-phase Bingham fluid. When modeling debris flows as single-phase fluids, the effects of larger particles, such as pebbles and rocks, are typically captured through the density of the mixture (${\rho }_m$), and rheological parameters (i.e., yield stress, ${\tau }_B$ and viscosity, ${\mu }_B$). Although this approach simplifies the system, the single-phase assumption can still offer valuable insights into the large-scale behavior of debris flows (Uchida et al., 2019; Lin et al., 2022).

3.1. Inclined Channel and Groundsill Array Setup

A laboratory-scale inclined channel was designed using AutoCAD software, with dimensions of 2.00 m length, 0.10 m width, and 0.1 m height. The channel bottom and walls were constructed from 5-mm-thick plexiglass material. At the upstream end of the channel, spanning $L_0=0.2$ m in length and separated by a 5-mm-thick, 0.10-m-wide control gate to form a reservoir. The sediment mixtures used within the reservoir had a total volume of approximately 2 liters, characterized by dimensions $L_0=0.2$ m and $h_0=0.1$ m. The initial dimensions of the reservoir were characterized by length ($L_0$) and height ($H_0$), as shown in Figure 4. The beam-type array of groundsills was placed along the inclined channel at a distance of 0.8 m from the control gate (Figure 4b). This arrangement allowed the sediment mixtures to have sufficient flow length to develop before interacting with the groundsills.

.

The beam-type groundsills had initial dimensions of height ($h$), and length ($b$), as shown in Figure 4a. Notably, the length of these groundsills aligned with the width of the channel (W = 0.1 m). All the groundsills were positioned over a length spanning 0.5 m in the middle of the inclined channel, which is within the range of $1.0-1.5$m as shown in Figure 4b. This setup allowed us the investigate the flow behavior and dynamics of sediment mixtures as they interacted with the beam-type array of groundsills in the inclined channel. Field investigations have shown that the majority of debris flows occur in gullies with channel slopes ranging from 10° to 20° (Tsunetaka et al., 2019; Pierson, 2020). Therefore, in this study, a $15°$ channel slope was used to approximate the typical range observed in natural debris flows, ensuring the experimental results are relevant to real-world observations (Tsunetaka et al., 2019; Pierson, 2020). The design considerations, such as the channel dimensions, reservoir configuration, and groundsill placement, were carefully chosen to facilitate a comprehensive study of the fluid-sediment interactions and associated physical processes.

3.2. Groundsill Array Configuration

To simulate sediment flow over the array of groundsills, the beam-type groundsills were positioned at various array densities. The dimensions of the groundsill are illustrated in Figure 4a, with respective values of $b=0.02$ m for the base width and $h=0.02$ m for the height. The array density is characterized by the ratio ($d/h$), where $d$ is the spacing between two consecutive groundsills, and $h$ is the height of the groundsill. In this study, five different array densities were simulated, with the ($d/h$) ratio ranging from 2 to 10 (Table 1). Based on the density of the groundsills, the array density can be broadly categorized as low array density ($d/h>4$) and high array density ($d/h\le 4$), as illustrated in Figure 5. The categories of high and low array density can be defined based on the geometric relationships provided in Eqs. 5a and 5b:

$${\displaystyle \frac{d}{h}}\le {\displaystyle \frac{1}{\tan \theta }}=\cot \theta$$

(5a)

$${\displaystyle \frac{d}{h}}>{\displaystyle \frac{1}{\tan \theta }}=\cot \theta$$

(5b)

For a channel slope of $\theta =15°$, Eq. 5a defines a high array density with $d/h$ less than $\mathrm{c}\mathrm{o}\mathrm{t}(15°)=3.72\approx 4.0$, and a low array density with $d/h>4$. High array densities have smaller spacings, which lead to different flow behaviors compared to low array densities with larger $d/h$ values (Wu, 2022). Figure 5c shows a schematic diagram for the definition of groundsill array density illustrating the relationship between $d/h$ and channel slope, $\theta$.

3.3. Simulation Domain, Boundary, and Initial Conditions

After creating the models in AutoCAD 3D, the channel geometry, reservoir column, and control gate were seamlessly integrated into the Flow-3D software for numerical simulation. The non-slip boundary condition is considered within the Flow-3D framework for the inclined channel with a rigid bed and wall configuration. At the inlet and outlet boundaries of the computed domain, a “wall (W)” boundary condition was applied to emulate a virtual barrier (Figure 6a). This “wall” boundary condition was also extended to the bottom surface and lateral sides of the channel, creating a confined flow domain. Further, we applied a no-slip wall boundary condition on the bottom surface, which ensures that the flow velocity at the boundary is zero, thus simulating the frictional resistance imposed by the groundsill on the sediment mixtures. This boundary condition is integrated with the turbulence and sediment transport models within Flow-3D, allowing the frictional effects to be incorporated into the governing equations of the flow (Flow Science, 2008). To characterize the unobstructed free surface, the “symmetry (S)” condition was employed. This condition allows for the proper representation of the air-fluid interface and its dynamics. A geometry resolution of 0.002 mm was employed within the channel, gate, and fluid regions to better represent the boundaries of walls, bottom surfaces, and interfaces in the flow simulation (Figure 6b). The lifting speed of the sliding gate was set at 0.3 m/s, which corresponds to the observed lifting speed during experimental investigations by Jan and Dey (2024).

4. Results and Discussions

4.1. Sediment Mixture Propagation

After the removal of the control gate, the sediment mixture starts to move from the reservoir, propagate downward, and interact with the groundsills. For high array density (specifically, $d/h=2$) the propagation of sediment mixtures is significantly impeded. A portion of the sediment mixtures are trapped and accumulate behind the groundsills. In contrast, for the low array density case (specifically, $d/h=10$), the sediment mixtures propagate more freely over the widely spaced groundsills. With larger gaps between the groundsills, the sediment mixtures can flow with less obstruction and deposition compared to the high array density case.

4.2. Temporal Evolution of the Flow Front

Figure 7 shows the velocity magnitude of sediment mixtures propagating over a low array density of groundsills in the inclined channel at different time frames. At $t=0.2s$, after releasing the sediment mixture from the reservoir, the velocity of the mixture is relatively high near the flow front. As time progresses to $t=0.7s$, the sediment mixtures start interacting with the groundsills, obstructing the flow and slightly decreasing the velocity near the front. By $t=0.8s$, as the mixtures propagate further over the groundsills, the velocity near the front continues decreasing due to the increased resistance. At $t=0.9-1.2s$, the sediment mixtures have traveled even further, and the velocity near the front is lower. Finally, at $t=3.0s$, after the sediment mixtures have propagated over the entire array of groundsills, the velocity near the front is very low (Figure 7).

Figure 8 presents the evolution of the flow front position over time for sediment mixtures propagating with and without the presence of a groundsill array. The data points for involving groundsills demonstrate a slower propagation of the flow front compared to without groundsills. This observation indicates that the presence of groundsill obstructions impedes flow propagation and retards the propagation of the sediment mixture front. The influence of groundsill array density is evident in the flow front behavior. Higher groundsill array densities, characterized by smaller $d/h$ ratios, lead to more obstructed flow conditions with reduced flow depths and slower front propagation. However, these high-density arrays also result in higher localized velocities due to the constrained flow paths. On the other hand, lower array densities, represented by larger $d/h$ ratios, allow for smoother propagation with less obstruction but also lower overall velocities. Interestingly, it is observed that at times between $t=$ 1.2 and 1.5 s, for lower array densities ($d/h=$ 8 and 10), the flow front travels faster than the cases without groundsills (Figure 8). This unusual behavior might be attributed to the increased local velocities induced by the spacing between groundsill arrays.

4.3. Influence of Groundsill Array Density on the Flow Profiles

The contour plots in Figs. 9a-e presents the fraction of sediment mixture at three different simulation times: 0.7, 1.0, and 1.5 s for flow over groundsills at various array densities ranging from $d/h=2$ to $10$. Based on these figures, we can observe how the flow patterns and extent of sediment mixture propagation vary with the groundsill array density. For high array densities ($d/h=$2 and 4) at $t=0.7s$, the sediment mixture flows relatively smoothly over the closely spaced groundsills with minimal disturbance (Figure 9a and 9b). However, at later times ($t=1.0s$ and $1.5s$), the mixture exhibits a more wavy pattern as it propagates over the groundsills, indicating increased interaction and disturbance caused by the groundsills. The sediment mixture also appears to propagate slightly less far downstream for $d/h=$ 2 compared to $d/h=$ 4, suggesting a stronger influence of the more closely-spaced groundsills on the flow.

For low array densities ($d/h=$ 6, 8, 10), even at the early stage ($t=0.7s$), the sediment mixture flow exhibits a more irregular pattern, with pockets of higher and lower fluid fractions forming around the widely-spaced groundsills (Figure 9c-e). At later times ($t=1.1s$and 1.5 s), the flow patterns become increasingly complex, with noticeable wavy structures and fluctuations in the fluid fraction distribution. It can be observed that as the $d/h$ ratio increases, the flow patterns become more complex and disturbed, with more wavy structures and irregular fluid fraction distributions. Furthermore, the extent of sediment mixture propagation downstream reduces, indicating a stronger influence of the more closely spaced groundsills on the flow behavior.

Figure 10 presents the sediment mixture fraction contours at a time $t=15s$, displaying a snapshot of how the sediment deposits appear after overflowing the groundsills and settling. When $d/h=2$, the sediment mixture overflows the groundsills almost uniformly, with a significant portion of the mixture flowing over each sill. The fluid fraction contours reveal that the sediment is well-distributed across the array, with minimal accumulation between the sills, nearly filling the space between them (Figure 10a). As the $d/h$ ratio increases to 6, the overflow becomes less uniform, with visible gaps emerging between the areas of sediment overflow. The fluid fraction buildup between the sills roughly forms a triangular shape, with the highest deposit height near the base upstream of each groundsill, gradually tapering off further upstream (Figure 10b). At $d/h=10$, the sediment mixture propagates over the groundsills with even larger spacing. The overflow becomes concentrated near the sills, with minimal sediment accumulation in the spaces between them (Figure 10c). This indicates that at higher $d/h$ ratios, the sediment mixture tends to bypass the intermediate areas (Figure 9), resulting in sediment deposition primarily near the sills.

4.4. Streamwise Flow Depth and Depth-Average Velocity Profiles

Figure 11a-c presents quantitative profiles of the flow depth and depth-averaged velocity for sediment flow over different groundsill array densities. The plots compare scenarios with and without groundsills when the simulation time $t=1.5s$. For $d/h=$ 2, the flow depth profile exhibits oscillations corresponding to the locations of the closely spaced groundsills along the inclined channel. Compared to the case without groundsills, the presence of the groundsills introduces disturbances and fluctuations in the flow depth and velocity profiles. The velocity profile also shows variations, with localized reductions in velocity at the groundsill locations (Figure 11a).

The fluctuations in the flow depth and velocity profiles become more evident and regular due to the increased $d/h$ ratios (i.e., $d/h=$ 6 and 10). The magnitude of the variations in both depth and velocity is larger compared to the high array density ($d/h=$ 2). These quantitative profiles, as illustrated in Figure 11, support the observations made from the contour plots shown in Figure 10, demonstrating that decreasing the groundsill array density (increasing $d/h$ ratio) leads to more complex and disturbed flow patterns, with larger deviations from the unobstructed flow conditions without groundsills.

The results presented in Figure 10 and Figure 11 reveal an important insight: when the sediment mixture interacts with the groundsills, a dynamic interaction occurs, resulting in a hydraulic jump-like flow pattern. The formation of these hydraulic jump-like flows is a result of the interaction between the sediment mixture and the groundsills. This phenomenon demonstrates the significant impact of the groundsill obstructions on the total hydraulic behavior of the propagating sediment mixture flow.

One of the key observations from Figure 10 and Figure 11 is that the groundsills act as barriers, disrupting the continuous flow and causing a rapid change in the flow regime. This change is characterized by an increase in flow depth and a reduction in velocity. As the sediment mixture approaches each groundsill, the flow depth is relatively shallow, indicating a supercritical flow regime with higher velocity. Upon reaching the groundsill, the obstruction causes a sudden rise in the flow depth, producing a hydraulic jump-like flow pattern. This abrupt change in flow depth is accompanied by a dissipation of energy and a decrease in flow velocity, transitioning the flow to a subcritical regime.

4.5. Formation of Hydraulic Jump-Like Flows Near the Groundsill

The analysis of the flow depth (Figure 11a) and velocity profiles (Figure 11b) provides further insights into the formation of hydraulic jumps. Along the inclined channel at $x=0.7$ m, the average approach flow depth ($y_1$) for sediment flow approaching the first groundsill is observed to be $0.014$ m. Conversely, the average post-jump flow depth ($y_2$) found to be 0.030 m (Table 2). This increase in post-jump flow depth ($y_2$) is approximately $2.2$ times higher than the approach flow depth ($y_1$), revealing the presence of hydraulic jumps.

Alternatively, the depth-averaged velocity obtained from the simulation, allows us to determine the Froude number for the flow approaching the groundsill ($F_{r1}=V_{a1}/\sqrt{g{y}_1}$) and the flow after the hydraulic jump ($F_{r2}=V_{a2}/\sqrt{g{y}_2}$). The parameters obtained for approach flow and post-jump flow are termed the approach flow parameters ($y_1,V_{a1}$, and $F_{r1}$) and post-jump flow parameters ($y_2,V_{a2}$, and $F_{r2}$), as shown in Table 2. The approach flow parameters include the approach flow depth ($y_1$), approach velocity ($V_{a1}$), and Froude number ($F_{r1}$) of the flow approaching the groundsill. Similarly, the post-jump flow parameters include the post-jump flow depth ($y_2$), post-jump flow velocity ($V_{a2}$), and post-jump Froude number ($F_{r2}$). The terms “approach” and “post-jump” characterize the flow conditions before and after the hydraulic jump occurs.

The observed approach flow parameters and post-jump flow parameters are presented in Table 2, showing that the values of the average $F_{r1}$ and $F_{r2}$ are $1.67$and $0.79$, respectively. This significant decrease in the Froude numbers indicates a transition from a supercritical flow regime ($F_{r1}>1.0$) before the groundsill to a subcritical flow regime ($F_{r2}<1.0$) after the first groundsill. This transition is evidence of a hydraulic jump-like flow occurring at the location of the groundsill, as schematically illustrated in Figure 12.

The fact that sediment flow front flowing over low array density groundsills travels faster than the cases without groundsills, as can be seen in Figure 9, Figure 10 and Figure 11, can be attributed to the increased local velocities induced by the spacing between groundsill arrays. In the context of sediment mixture propagation over groundsill arrays, lower array densities correspond to higher $d/h$ ratios, indicating that the groundsills are spaced farther apart relative to their height. This configuration leads to less obstructed flow conditions, with larger spacing between the groundsills (Choi et al., 2018). While the overall flow is less impeded by these fewer obstructions, resulting in smoother propagation on average, the larger spacing between groundsills can create localized acceleration and higher velocities in certain regions (Ng et al., 2014). As the sediment flow encounters these larger openings, the lack of confinement can cause the flow to accelerate through these unobstructed sections, potentially resulting in faster front propagation over short time intervals. This phenomenon can be attributed to the principles of fluid dynamics and the conservation of mass and momentum. When the flow is less constricted, with larger cross-sectional areas available, the debris flow can accelerate to maintain continuity and conserve momentum, leading to localized increases in velocity (Jacob et al., 2005).

At a lower $d/h$ ratio, the sediment mixtures tend to accumulate uniformly between the two consecutive groundsills and this can be explained by the fact that lower array densities correspond to dissipate energy and promote more uniform sediment deposition (Takahashi, 2007). As the spacing between obstacles increases, sediment mixtures tend to bypass the intermediate areas and concentrate around the groundsills (Figure 9), resulting in non-uniform sediment deposition. When obstacles are spaced farther apart, the flow can develop enough momentum between them to maintain its speed. This allows the flow to jump from one obstacle to the next, resulting in sediment being deposited primarily near the obstacles rather than in the spaces between them, where the flow remains relatively undisturbed (Armanini and Gregoretti, 2000). The findings and analysis presented in this study are based primarily on numerical experiments of single-phase fluids, using a controlled setup with constant channel slope and groundsill dimensions. This setup allows us to gain valuable insights into the flow behavior of sediment mixtures flowing over varying groundsill array density, which can enhance our understanding of debris flow behavior in real-world scenarios. However, it did not include practical flume experiments or account for the interaction between solid phases like pebbles and the flow-bed dynamics. Future works are recommended to explore these aspects, including comprehensive simulations involving fluid-solid coupling and validation through real-world applications, to refine the design and operational strategies for groundsills.

5. Conclusions

In this study, we numerically investigate the behavior of sediment mixture with a sample volume of 2 liters flowing over arrays of groundsills installed in a laboratory-scale inclined channel with a fixed slope of $\theta =15°$. The sediment mixture is modeled as the debris flow that follows the Bingham Fluid Model, characterized by two rheological parameters: Bingham yield stress (${\tau }_B$) and Bingham viscosity (${\mu }_B$). Using the commercial CFD software Flow-3D, we simulate the sediment flow over groundsill arrays with varying densities, represented by the ratio of spacing ($d$) and a fixed groundsill height ($h$ = 0.02 m). Specifically, the array densities vary between $d/h=2$ and $10$. The key findings are:

• 1. Increasing groundsill array density significantly obstructed sediment mixture propagation, reducing the flow front velocity and causing a portion of the mixture to accumulate behind the closely spaced groundsills.
• 2. High array density ($d/h\le 4$) exhibited smoother flow over the groundsills initially but transitioned to more wavy flow profiles with more flow interaction at later times.
• 3. Low array density ($d/h>4$) allowed free propagation but induced irregular flow patterns and fluctuating depth and velocity profiles from the start due to larger spacing between groundsills.
• 4. The interaction between the propagating sediment mixture and the groundsill array resulted in a hydraulic jump-like flow pattern near the groundsills. The flow depths at the location of the first groundsill were observed to be approximately $2.2$ times higher than the approaching flow depths.

Overall, modifying the groundsill array density provides a potential mitigation strategy by obstructing and slowing debris flows while dissipating energy through induced turbulence and wave-like patterns. High array density best restricts initial propagation, while low array arrays promote dissipation over longer runouts. Further investigations are recommended to study the effects of varying sediment mixture volumes, and different channel slope conditions on the flow behavior when encountering obstructions like groundsills. The insights obtained from this study could guide the strategic design and placement of groundsill arrays for effectively reducing debris flow mobility and mitigating related hazards.