Debris Flows, Boulders and Constrictions: A Simple Framework for Modeling Jamming, and Its Consequences on Outflow

Large boulders can jam in constrictions, both in canyons and man‐made structures (e.g., slit‐dams). Boulder jamming occurs stochastically, and partially governs the outflow rate of debris flow material. Explicit hydro‐mechanical models of boulder‐laden flows are too computationally demanding to study this stochasticity. This study presents a new framework implemented into a numerical program. This program encapsulates basic hydraulic equations and criteria predicting boulder blockage at narrow sections, along with a simple statistical way to compute boulder sizes and numbers. The program applies empirical estimates by experts of the average number of boulders in a deposit to evaluate how equivalent multi‐phase flows would interact with a constriction. The model stochastically generates boulders that can obstruct the constriction(s). The program outputs simple descriptions of the upstream flow level over time, as well as the downstream outlet discharge rate and volume. The model is fast (5–30 s per run), allowing uncertainty propagation analyses of interactions between flows, boulders and constrictions. Interaction between both low‐ and high‐risk flows can be quickly evaluated for different scenarios. For uncertainty propagation, we use possibility theory, which accommodates uncertainties relevant to debris flows. The framework gives practitioners a much‐needed generalized approach for understanding the long‐term behavior of boulder‐laden flows through canyons, or for designing engineered structures. Finally, we apply the model to a field case – the design of a proposed North American slit‐dam – thence elucidating the design requirements for effective flow control. We found that horizontal bars are necessary for dependably controlling outflow from this structure.

1. At the larger scale, there is typically a main constriction that is narrower than the channel. The main constriction does not necessarily extend to the base of the channel. 2. At the smaller scale, the main constriction may optionally be "subdivided". These "sub-openings" are created by partially obstructing the main constriction using concrete, steel beams and/or grillage bars. Subdivision causes smaller boulders or wood chunks to become trapped.
It should be noted that certain details vary between structures such as slit-dams and sabō dams, as well as the properties of the flows with which they interact. By extension, the intended outcomes of installing such structures varies too. Nonetheless, a common point is that these structures are all designed with both hydraulic and mechanical control in mind. We also emphasise that the phenomenon of boulders becoming jammed is of strong interest to scientists looking at long-term geomorphological changes in natural topographies.
For simple constrictions that extend to the base of the channel, and where only hydraulic control is important, we can use the slit equation. The slit equation gives the rate of unsteady outflow from a constriction: where Q is the discharge rate [m 3 /s]; t is time; μ is a dimensionless opening coefficient taken as 0.65 (assuming that the main constriction has not been subdivided; Zollinger, 1983); g is gravitational acceleration [m/s 2 ]; w is the slit width [m]; h is the flow depth over the bottom of the channel [m]; and y(t) is the base of the constriction [m]. For purely hydraulic flows, one can trivially compute the range of possible outflow rates for given assumed hydrographs. Note that it is likely that different flow rheologies can be captured by modifying the coefficient μ, potentially describing the distinct behaviors of (say) water, mudflows and debris flows. However, the relationship between μ and flow rheology remains an open scientific question for future research to address.
In contrast to purely hydraulic flows, if flows contain large boulders, jamming can occur at constrictions. This jamming can have a large influence on the outflow rate. Trying to account quantitatively for jamming (i.e., mechanical control) greatly complicates matters. The number, size and spatial distribution of large boulders within a flow can vary enormously, affecting whether and when jamming occurs. Furthermore, "wood jams" can occur. Large wood chunks, which are prone to fluvial transport, can jam at constrictions before large flow events occur (Horiguchi et al., 2021;Piton & Recking, 2016c;Piton, Horiguchi, et al., 2020;Tateishi et al., 2020). These types of jamming affect the terms w and y(t) in Equation 1, that is, the slit width and the position of the base of the constriction. Furthermore, the same constriction might behave differently for two events with a similar overall volume of material, but with a substantially different number and/or size of boulders. Indeed, for constrictions 8-10 m in height, one or two very large boulders can cause the entire basin to become filled; without the boulders, only a small upstream deposit would occur (Piton, Charvet, et al., 2020). Furthermore, two flows with identical volumes and compositions may still behave differently given non-identical spatio-temporal distributions of boulders, because jamming depends on when boulders of particular sizes arrive! Currently, the expected degree of jamming comes from empirical estimates by expert engineers or scientists (hereafter termed "expert estimates"). Unfortunately, expert estimates have difficulties dealing with the high degree of stochasticity inherent to jamming problems. Past studies on geophysical flows and constrictions generally address either purely hydraulic processes, or mechanical trapping, but seldom both simultaneously. For instance, Armanini and Larcher (2001) investigated hydraulic sediment control of sand in water flows, whilst Armanini et al. (2006) performed experiments with mud and constrictions. Examples of studies on mechanical control cover sabō dams (e.g., Ishikawa & Mizuyama, 1988;Osanai et al., 2010;Shima et al., 2016;Watanabe et al., 1980), slit dams (e.g., Choi et al., 2016;Goodwin & Choi, 2020;Leonardi et al., 2019;Shrestha et al., 2012) and baffle arrays (e.g., Goodwin et al., 2021). Small-scale physical studies involving both fluids and large solids exist (e.g., Chen et al., 2020;Sun et al., 2018;Sun et al., 2021), but have issues with scaling (see Iverson, 2015).
In contrast, numerical solvers for interaction between flows and constrictions can be run at any scale. Furthermore, some numerical works explicitly consider both an interstitial fluid and boulders. One method is to use fully 2D or 3D fluid mechanics solvers coupled with the Discrete Element Method (e.g., Canelas et al., 2017;Horiguchi & Komatsu, 2018;Leonardi et al., 2014;Shan & Zhao, 2014). Other methods such as HyperKANAKO (e.g., Nakatani et al., 2016;Yanagisaki et al., 2016) which use shallow-water formulations together with relationships for the deposition and trapping of sediment of two specific sizes, also exist. These depth-averaged methods are relatively numerically efficient, although there are open questions relating to whether the assumptions underpinning models such as HyperKANAKO still hold in debris basins with very thick layers of debris, which are pertinent to slit-dams or canyons that can be more than 20 m high. It is these very high constrictions which are of primary interest in our study.
In sum, both depth-averaged and non-depth-averaged coupled methods have a strong theoretical basis for encapsulating both hydraulic and mechanical control of geophysical flows. However, both remain unsuitable for analyzing uncertainty propagation given their high computational demands, which stem from the explicit modeling of the flow and/or boulders. Furthermore, even relatively efficient coupled numerical methods such as HyperKANAKO are not yet able to account for the mechanical trapping of a mixture of arbitrary sizes of boulder.
In short, we lack numerical models that can explore the broad possibilities for hydrograph characteristics, boulder characteristics, and initial wood-jams over a meaningful number of trial runs. These explicit models help understand the underlying mechanisms of interaction between flows and constrictions, and are best suited as complementary tools for statistical investigations of best-and worst-case scenarios for highly uncertain multivariate problems.

Summary of Scientific Challenges & Scope of Present Work
The overarching problem to be solved is: "Given that boulders within debris flows can become jammed stochastically in constrictions, how can we most efficiently model the physical processes involved for a statistically significant number of trials?" An overarching condition for the scientific challenges is avoiding using explicit, computationally expensive numerical methods (such as the Discrete Element Method for the solid phase). Under this main condition, the scientific challenges include: 1. How can we model boulders within flows, as well as their jamming in constrictions? 2. Given limited information on the specifics of boulders within debris flows, what is the basis for the number, size distribution and spatio-temporal distribution of boulders that we model? 3. How can we combine the results from n simulations to give a meaningful overview of the possibilities for interactions between flows, boulders and structures?
To address these challenges, this manuscript describes a new computationally-light possibilistic framework for an entirely novel application -boulder jamming. The framework amalgamates well-established principles for (a) fluid mechanics and (b) statistics. We implement the framework into a simple numerical program. The program allows engineers and scientists to quickly understand how constrictions may moderate boulder-laden debris flows, for a wide set of plausible input parameters, and over a large number of simulations.
In this manuscript, we first discuss the theoretical basis for the framework and program, specifically (a) the governing equations for hydraulic control using constrictions and (b) the conditions required for mechanical jamming of boulders at constrictions. Thereafter, we present the possibilistic statistical method implemented for modeling boulders arriving at and interacting with constrictions. We then discuss possibilistic analyses for uncertainty propagation using the model. Thereafter, we use the new model to analyze the preliminary design of the Cheekye slit-dam. We use newly acquired information about the boulder size distribution at the site of the proposed dam.

Theoretical Basis of Novel Framework
In this section, we first detail the equations used for modeling hydraulic control of flows, and then the novel routine for considering mechanical control (i.e., random generation of boulders and determining when they jam). Figure 2 defines the geometry of the constrictions and sub-openings considered in this study. Figure 3 shows a conceptual schematic of how the individual components of the model developed in this paper fit together, with inputs in blue boxes and outputs in yellow boxes; each of the components are discussed in the text hereafter.

Hydraulic Control of Flows Caused by Constrictions
Analytic hydraulic equations can be applied to flows interacting with constrictions. We define the area upstream of the constriction as the "basin". We present here the slit equation (Equation 2), the Grand Orifice equation (Equation 3) and the weir equation (Equation 4) (Piton & Recking, 2016b). These equations are adapted for the possibility of the geometries of the constrictions changing with time, given the presence of boulders, by writing the slit width as w i (t) and the base of the flow as y i (t):  (c) show three possible geometric configurations of man-made structures constricting the channel and its sub-openings, including both vertical and horizontal openings, as well as spillways. Parameters are also displayed: flow depth over the barrier base level is denoted h, discharge, basal clogging, and flow width at opening number i are denoted Q i , y i and w i , respectively. If the opening is an orifice, its height is denoted a i . The spillway wing angle with the horizontal is denoted as Φ.
variables are as defined as for Equation 1. The index i denotes particular sub-openings. The variable y i is also the height of the bottom of the sub-opening above the base of the channel; a i is the slot height [m], whilst ϕ is the angle between the weir wing and the horizontal [°]. The variables are defined graphically in Figure 2.
For each flow event, we define an input hydrograph representing the supply of water and sediment into the basin Q sup (t). The input Q sup (t) and the mass outflow Q(h(t)) from the constriction must be calculated. A mass conservation equation links the input and output: where Q is the discharge from the constriction [m 3 /s], Δt is the time step [s] and ΔV b is the variation in the basin storage [m 3 ]. The term ΔV b can be used to find the local depth of debris within the basin. This requires a relationship between the volume stored V b and the depth h (see the Supporting Information), which depends on the basin geometry.
Equations 2-4 have been extensively evaluated based on experimental work dealing with water. We adopt them to model non-Newtonian fluids in this study, noting that no values are available for μ for non-Newtonian fluids. Nonetheless, debris flow materials are often shear-thinning, and the shear rate as a flow passes through the constriction is likely to be very high, suggesting that applying the value of 0.65 directly to non-Newtonian fluids may not be unreasonable. Previously, μ has been reduced to account for diminished outflow caused by partial obstructions (e.g., Piton, Horiguchi, et al., 2020). However, in reality, the obstruction develops with time during the events, and there is no basis for determining a single representative μ value that approximates the outflow reduction.
Fully accounting for boulder jamming in Equations 2-4 requires: (a) defining the spatio-temporal distribution of boulders reaching the constriction; and (b) using relevant criteria to define whether the boulders are jammed and obstructing the constriction, based on their spatio-temporal distribution within the flow. These topics are addressed in the following sections.

Granular Dynamics: Conditions for Mechanical Jamming
Conditions conducive to jamming of constrictions by boulders -boulders having a non-zero friction angle and being able to rearrange relative to one another -are known from experimental, numerical and field survey works (e.g., Choi et al., 2016;Goodwin & Choi, 2020;Osanai et al., 2010;Shima et al., 2016;Shrestha et al., 2012;Watanabe et al., 1980). There are two main types of jamming: "lateral jamming" and "vertical jamming". In this manuscript, the smallest class of boulder considered are boulders with a diameter above 0.5 m. This diameter was chosen based on the size of the constriction to be studied.  Figure 2c). The condition is:

Lateral Jamming
where D A and D B are the diameters of the two largest boulders passing at the same time in a particular constriction; and w i (t) is the free section width of the constriction. Note that the parameter w i (t) can also be modified to account for trapezoidal shapes and progressive widening with increasing height, as in the canyon in Figures 1a and 1b or the spillway of Figure 2c. We emphasize that lateral jamming causes a vertical increase of the base level 7 of 27 of the flow section. This is considered numerically as though the bottom part of the opening has been elevated by a height equal to the diameter of the biggest boulder blocked (e.g., Figures 2a and 2c). Furthermore, if more boulders than are required for a lateral jam try passing at the same time, they are assumed to be either ahead or behind the two biggest boulders, and thus do not pile up (see also Figure 4).

Vertical Jamming
"Vertical jamming" is where boulders do not pass through a constriction because they are higher than the height of the opening (Figures 1c and 1d). It is only relevant at sections where an element restrains the flow height (e.g., an overhanging cliff; a natural canyon bridge; a bridge deck; steel bars in sabō dams; or concrete bars in slot dams). A good example of this type of jamming is shown by the large boulders jamming the bottom slot on Figures 1c and 1d and 2b. The condition is: Jamming occurs for all boulders where − ( ) where D is the boulder diameter; a i is the height of opening number i; and y i (t) is the height of a lateral boulder jam at the base. When a group of boulders tries to pass at the same time through a given constriction, the condition of Equation 7 is tested for each of them. The constriction width is then decreased by the sum of all boulder diameters satisfying Equation 7. Thus, vertical jamming causes a lateral constraint on the flow section. Combined jamming states may emerge in constrictions after successive jamming events throughout a flow . Wide . Part (b) shows a simplified version of the same deposit, but where only the class-1 boulders are shown as being highlighted in yellow. White polygons correspond to mud or boulders of other class. Part (c) shows the same concept for class-2 boulders. Note that since class-2 boulders are smaller than class-1 boulders, their maximum virtual number N j is much bigger in the same reference volume. It should be emphasized that the flow itself is not explicitly modeled in the coupled method presented in this paper -the flow characteristics are based on various possible hydrographs. By extension, neither the upstream deposit nor the distribution of particular boulder classes within a deposit are explicitly modeled. However, a deposit profile similar to that shown within this Figure -including the potential location of boulders of particular classes -can be reconstructed based on information output by the program. and low slots will typically start being jammed vertically, thus restraining the effective section width until it becomes sufficiently narrow to become laterally jammed (e.g., the central orifice on Figure 2b). Narrow and tall slots will also typically start being jammed laterally until the flow section becomes so low that vertical jamming may appear (e.g., the bottom orifices on Figure 2c).
These jamming conditions are simple to implement numerically Shima et al., 2011), and probabilistic approaches exist too (Katsuki et al., 2015). However, an open challenge is to define a way to estimate how many boulders -and of which sizes -try passing at the same time for each of the sub-openings.

Implementation of Novel Numerical Method
We use a binominal statistical distribution for simulating the stochastic arrival of boulders at a constriction. This distribution provides a randomized sample of how many boulders exist at any given location near the constriction and its sub-openings at a given time. Previous works studying boulder jamming in sabō dams have used probabilistic approaches, but did not continue the analysis for modeling the downstream discharge of material. This was probably because a couple of boulders were sufficient to jam the structure being studied Shima et al., 2011), which is not necessarily always the case. We now introduce the concepts underpinning our novel method. An example is given in Section 5 for a large structure where many boulders are required to fully jam the opening.

Sub-Routine for Boulder Jamming
As discussed previously, the size of the boulders relative to the constriction size governs whether they might become jammed. Any meaningful predictive method must thus account for the natural variability of boulder sizes. For simplicity, we adopt classes of boulder sizes, rather than a continuous distribution. We denote Class-1 as the biggest (coarsest) class, so boulder size decreases as the class index increases ( Figure 4). Practically, the incidence rates of particular classes of boulder was estimated from a field study near the site of the proposed Cheekye slit-dam; the results of this study are given in Table 1.
The variable J gives the total number of classes of boulder size used. The subscript for a given class is j. For a given class-j, boulder diameters lie between a lower bound "min" and an upper bound "max". The range of boulder diameters for a given class-j are D j, min − D j, max . For each class-j, the average boulder diameter is the mean value of the minimum and maximum diameters in each class: Log jam height (case # 2) m 5 10 10 15 No. Of boulders per 0.1 Mm 3 0.5-1.0 m -1,500 3,000 3,000 6,000 Note. Only the volume and Q peak were different for the 0.8 Mm 3 event event. (9) our model performs successive jamming analyses for each class of boulders. Jamming analyses always start from class 1 (the largest) and finish with the smallest. One binomial distribution is defined for each class.
A binomial distribution includes two parameters -N, the number of independent trials, and p, the probability of success (so the probability of failure is 1 − p). By definition, the binomial distribution gives the discrete probability distribution of the number of successes in a sequence of N independent trials.
The binomial distribution is used when randomly generating the number of class-j boulders (n j, boulder ) in the debris flow volume (V control ) passing through an opening. The volume V control is divided into a set of smaller volumes, each of which have the same size as the average boulder volume for a particular class, that is, V j, boulder . These smaller volumes are hereafter termed "packets". The number of packets is N j, packet , corresponding to the maximum number of boulders of class-j that could appear within the "main" volume of V control . The equation giving the number of packets N j is: Note that the quantity N j, packet also represents the number of virtual trials required for a particular class of boulders at each time-step. In other words, at each time-step we inspect the probability of each packet including a boulder belonging to class-j. The "successes" -where the random number from the binomial distribution for a packet corresponds to a boulder belonging to class-j exists -give us the value of n j, boulder (i.e., the number of boulders belonging to class-j in V j, control ). The remaining packets of material, given by N j, packet − n j, boulder , represent a volume of mud or boulders of other classes. This volume is smaller than the original control volume V control . Whether the remaining volume contains mud or boulders of other classes is determined from further loops considering other classes of boulders. These additional loops involve discretizing the remaining control volume into ever-smaller packets. We assume that packets that are "unoccupied" after having considered all boulder classes comprise mud and/or boulders that are sufficiently small to be neglected in this simple modeling of the overall process.
The objective of the above routine is to systematically compute n j, boulder (i.e., the number of boulders of class-j in a control volume by a constriction at a given time-step), achieved by performing N j, packet randomized samples on a binomial distribution of known probability p j .
To obtain suitable binomial distributions, we suggest using field data where boulder counts are available. The mean value of the binomial distribution ⟨ boulder ⟩ is: Estimating p j requires identifying a reference deposit of a debris flow with a large enough volume V ref to assume that there are a representative number of boulders of each class-j (i.e., where boulder ≈ ⟨ boulder ⟩ ). Figure 4 depicts the concept. Counting the numbers of real class-j boulders in a large volume V ref is equivalent to physically measuring ⟨ ⟩ . Then Equation 11 can be inverted to define the probability of finding a boulder of class-j within a packet: This operation must be done class-by-class because N j, packet changes for each representative boulder diameter ⟨ boulder ⟩ , as shown in Equation 10. Figures 4b and 4c can be helpfully compared here. Physically counting the number of boulders belonging to each class in a large volume of a real debris flow V ref enables estimation of the parameter p j for each class. Then, parameter p j can be reused for any given volume of debris flow V control smaller than V ref . (Reminder: in our numerical implementation, V control always refers to a control volume of material interacting with a constriction.) The number of boulders of each class n j, boulder in the volume V at a particular instant can be simulated by performing N j, packet trials on a binomial distribution of parameter p j . The probability mass function P, that is, the probability that n j, boulder is exactly equal to a given value k (which is a dummy variable), is defined by: Random sampling on this distribution is easy to code. For our implementation, we used the function rbinom of the package "stats" in R (R Core Team, 2020).
As noted previously, boulder classes cover a range of diameters [D min , D max ]. Each time a class-j boulder is successively generated, its diameter must be defined. Rather than using the mean diameter value computed previously for the packet size ⟨ boulder ⟩ , we randomly select a value in the range [D j,min , D j,max ] assuming a uniform distribution. This allows consideration of both the inherent variability in boulder shape, and the impreciseness associated with using classes of sizes.

Overall Routine for Each Run
The overall routine for the program starts by defining all variables, including the constriction, hydrograph and boulder probability distributions. This is described further in the next section. At each timestep -which is set at a constant 1 s -the following sub-routines are carried out ( Figure 5): 1. Update the volume and depth of the inflow in the basin using the hydrograph. 2. Find the outflow from the constriction based on the hydraulics equations (including conservation of mass) and the jamming status at the previous time step, and update the volume and depth of the material in the basin again. 3. Run the sub-routine for boulder generation and jamming, and update the constricted opening geometries if necessary.
This process is looped until the outflow rate approaches zero (either due to jamming or due to all material having flowed out).

Statistical Basis for Dealing With Uncertainty Over Many Runs
The approach presented above is much more advanced than basic engineering practices such as those based on direct estimation of clogging rates by experts. An advantage of our approach is that we can statistically consider several uncertain parameters in outputting possible outcomes of events: 1. The flow volume, which is the main parameter for describing the event magnitude for a given event return period. Nonetheless, the volume for a specific event is often uncertain, even on sites where comprehensive case studies have been undertaken (Jakob & Friele, 2010). 2. The peak discharge and overall hydrograph, which are subject to many uncertainties, even for events that have already occurred. Furthermore, there is no well-defined relationship between the peak discharge and the flow volume (Nagl et al., 2020). 3. The deposit slope in the field, which varies spatially across the deposit, hampering characterization (Hungr et al., 1984;Piton & Recking, 2016a). This uncertainty greatly affects the "stage-volume" capacity curve of the basin: the total capacity for retention of material behind the constriction without overflow occurring. (As flow material drains through the constriction, this storage capacity can be partially recovered). 4. The extent of filling before an event, which depends on past events and maintenance policy (Carladous et al., 2022;Chahrour, 2021). It is often impractical to regularly assess sites, so the basin condition before a large event is generally uncertain. 5. Initial wood jams. Within debris flows, logs and trunks are often crushed and destroyed by impacting boulders, justifying the exclusion of wood from the model for mechanical jamming. However, surges of large chunks of wood gathering at, or ahead of, the front of some field debris flow events have been reported (Friele et al., 2020; Roberti et al., 2017). Clogging of constrictions by frontal woody surges are also regularly reported (Shima et al., , 2016, although the processes involved remain poorly understood. 6. The number of boulders in each class, which is tricky to assess and varies between events (Shima et al., 2016).
All of these sources of uncertainty are dealt with in our model using the possibilistic distributions described hereafter.
The parameters listed above are subject to both epistemic and stochastic uncertainty. Epistemic uncertainty is represented by a knowledge gap for the process, and difficulty in measuring the parameter values. Any set of these input parameters for any model is thus somewhat uncertain. Stochastic uncertainty is represented by natural variation of parameter values between events. Our model accounts for stochasticity in the spatio-temporal distribution of boulders through random generation. For any model accounting for both epistemic and stochastic uncertainties, any two model runs with identical input parameters will generate different results. The next question relates to how to model these uncertainties. We review here the Monte-Carlo method, imprecise probabilities and possibility distributions.

Monte-Carlo Method
Stochasticity is often considered using the Monte-Carlo approach, which requires a Probability Density Function (PDF) defined for each input variable. A Monte-Carlo model is called many times (typically ∼10,000 times). Examining the relationship between the input PDFs and the output from the model enables study of the output variability for inputs of known variability. The Monte-Carlo approach involves defining each parameter using an appropriate probability law (e.g., normal, log-normal, or gamma distributions). Fitting parameters (e.g., the mean or standard deviation) are used to further constrain PDFs.
For debris flow hazards, there are multiple parameters that could be selected for investigation using Monte-Carlo analyses. These parameters can have different degrees of generality. Broad studies that consider the most general attributes of debris flows, such as the debris flow volume and runout length (e.g., Paola et al., 2017) may be good candidates for Monte-Carlo treatments if sufficient data from sufficiently homogeneous data sources are available. However, Monte-Carlo analyses are certainly less well-suited to data that is more specific, such as the number, size distribution and spatio-temporal distribution of boulders within a flow, or the height of the jam that wood chunks will create in a constriction -there is almost no field data that could be used to support such an analysis. Accurately defining suitable probability laws for such specific parameters is impossible. Indeed, making such definitions would require vast amounts of instrumentation and data -not to mention time -which is not feasible.

Imprecise Probabilities
Alternatively, one could adopt imprecise probabilities. This involves specifying probability distributions governing specific parameters, whilst also explicitly accounting for uncertainty in the fitting parameters (i.e., the mean and standard deviation). This is infeasible because the type of PDF is even less well-known than the values of the fitting parameters.

Possibility Distributions
In practice, after field visits and site investigation experts can only give provide educated guess on their best estimates such as "The number of class-1 boulders is likely 1-2 in a reference volume of 200,000 m 3 of debris flow", along with bounds, for example, "In the same volume, the number of class-1 boulders cannot be less than 0 and it is virtually impossible that it is more than 5". The expert thus provides information on the most probable values (a single value or an interval) and on the minimum and maximum possible values. Importantly, the type of probability law is neglected, making classical uncertainty methods that use only probabilities (normal or imprecise) irrelevant.
Fortunately, another statistical framework exists for such limited input data: possibility distributions. Several studies have used possibility distributions for geosciences and natural hazards given their suitability for the accuracy of input data used in modeling (Baudrit et al., 2006;Guyonnet et al., 2003;Rohmer & Verdel, 2014). Presenting the full fundamental mathematical basis of possibility distributions and hybrid-uncertainty propagation, that is, merging possibilities and probabilities, is beyond the scope of this paper. We refer the reader to Guyonnet et al. (2003) and Baudrit et al. (2006) for the full story. Notwithstanding, we will present a simple synthesis of the method (see also Figure 6). The statistical framework of possibility distributions was implemented in this study using the software library HYRISK, developed by Rohmer et al. (2018), which is available for the programming language R. Our implementation was used to perform uncertainty propagation using possibility distributions in our model.
The objective of uncertainty propagation using hybrid-or possibility distributions is to define probability boxes ("p-boxes"), that is, upper and lower cumulative distribution functions (CDF). These bounds are defined mathematically such that the "true" CDF, that is, the one that would be obtained using a perfect Monte-Carlo analysis, is bounded by the box (the gray area in Figure 6, bottom panel), but remains unknown due to insufficient information.
To build p-boxes, one of three types of possibility distribution is developed based on an expert's assessment ( Figure 6, second panel): 1. Given maximum and minimum values, and a single-valued best estimate, a triangular distribution is adopted ( Figure 6, first column). 2. Given maximum and minimum values, and a range of best estimates, a trapezoidal distribution is adopted ( Figure 6, second column); 3. Given only minimum and maximum values, a rectangular distribution is adopted ( Figure 6, third column); Possibility distributions are defined between zero and unity. The interval between the minimum and the maximum parameter values, called the "support" in the specialist literature, are the outer bounds (written hereafter as [ ; ] ). A possibility of zero is assigned to the bounds. Mathematically, a possibility distribution means that a parameter's value cannot be outside of the bounds. Meanwhile, the best estimate, either a single value or a range, is called the "core" by statisticians and is written hereafter as [ ; ] . The best estimate is assigned a possibility of unity, and mathematically represents the most probable range of values of the parameter.
In the following explanation we use four generic variables. The generic input variables are called X 1 , X 2 , X 3 . The generic output from a model is called Y. We assume that Y is positively correlated with X 1 , X 2 and X 3 . If the output Y is negatively correlated with X, the upper bound of Y is computed with the lower bounding values of X, and inversely for the lower bound. If Y does not monotonically vary with X, additional sampling should be performed to get the local minimum or maximum values, although is computationally costly.
To compute the p-box using monotonic models, the upper and lower bounds are computed independently (the orange and red curves in Figure 6). Here, we explain only the process for the upper bound, which is nearly identical to that for the lower bound. The objective is to identify the worst-case distribution, which is given by the upper bound of the distributions for each variable (Figure 6, second panel).
For illustrative purposes, a few plausible PDFs are drawn in Figure 6b -second panel. They all have a mode, that is, a most probable value at the upper limit of the best estimate. The most conservative distribution, maximizing the parameter value, is the uniform distribution between − (and between − for the lower bound, respectively). The associated CDF is thus the upper bound of all possible CDFs consistent with the expert opinion, that is, with all possible information available to the designers (the fourth panel of Figure 6).
We stress that this approach of "combined worst cases" is a (sometimes very) conservative approach: although the expert may be able to give an appraisal of the most probable value of parameter ∈ , extreme cases could still occur, for example, where ⪆ or ⪅ . We will later explain our suggestions for interpreting the results and use several assessments to produce a type of sensitivity analysis.
Mathematically, to ensure that all plausible cases consistent with the expert opinion are covered, one must consider an upper bound. Any value between the upper best guess and the maximum (i.e., within − ) is In this conceptual figure, Y is assumed to be positively correlated with X 1 , X 2 and X 3 , each of them being a generic variable. equally probable. When one single best estimate value is provided, the upper and lower bounds will be sampled in two distributions close to each other ( Figure 6 second panel, left column). If a range of best estimates is provided, the sampling is done in two separate distributions (Figure 6 second panel, central column). If no best estimate values were provided, the conservative hypothesis is to systematically use the maximum value for the upper bound and the minimum value for the lower bound ( Figure 6 second panel, right column).
The CDFs are now defined for all parameters for both the lower and upper bounds, where pairs of Monte Carlo analyses can be performed. Parameter values are randomly sampled defining a set of input values used for the lower bound (the orange dots in Figure S1b in Supporting Information S1), and a set of input values are used for the upper bound (maroon dots in Figure S1b in Supporting Information S1). A sufficiently high number of paired calls enables identification of the upper and lower bounds of Y, and thus also the p-box (Figure 6 fourth panel). Using these conservative mathematical hypotheses enables the definition of a box which definitely contains the "true" CDF, given a reasonable expert opinion. (Side-note: see Figure S1 in Supporting Information S1 for more details on the basics of constructing probability distributions.) We now detail the computations with respect to building p-boxes to facilitate subsequent interpretation. At least two sets of simulations are required to construct these p-boxes, adopting (a) maximizing assumptions, namely namely uniform over the range − , and (b) minimizing assumptions namely uniform over the range − , (see Figure S1 in Supporting Information S1). The expert assessment provides all assumptions.
As mentioned, the entire upper bound corresponds to the pessimistic assumptions. We suggest that data interpretation can safely avoid focusing on the upper tail of the upper bound (e.g., the 75th percentile or above), because it represents a confluence of all pessimistic assumptions, on top of the least-ideal randomness for boulder generation. These cases are extremely unlikely to occur. Instead, we strongly tend toward the best estimates (i.e., the bottom of the upper bound and top of the lower bound), whilst also considering the progressive shift toward the central part of the distribution. Furthermore, the shape of the distribution gives additional information: the output parameter may undergo a smooth, continuous increase between the best estimates and the tails, or conversely may have "steps", that is, sudden changes in the output variable. Steps can indicate that parameters considered by the model have threshold effects.
Finally, using possibilistic distributions to simultaneously propagate all the uncertainties given in Section 4 is useful as a global reference, but often results in p-boxes that are so wide that there is limited insight into the specifics of interaction between flows and constrictions. Consequently, we adopted "what if?" scenarios where we independently changed a subset of the hypotheses, for example, "What if the initial woody jam is rather small?" Comparing these "what if?" p-boxes with the reference p-box highlights the dependency of output variables on the optimism/pessimism of the original hypothesis. Individual p-boxes can also aid in checking the robustness of the constriction under worst-case scenarios, for example, "what if there is no initial woody jam at all?" Combining simple outflow equations with possibility distributions requires <30 s per run on a laptop -fast enough for uncertainty propagation analyses using multiple parameter sets.

Modeling the Cheekye Slit-Dam
We now apply the newly developed numerical model to the proposed Cheekye slit-dam, a dam proposed to protect the Squamish community, British Columbia, Canada. The approach proposed in this paper was actually developed specifically for this case study. Boulder size distributions were estimated by BGC Engineering to feed into this analysis. The model and the approach is nonetheless equally applicable to other debris basins or to natural constrictions such as canyons (see Figure 1).
The Cheekye River is part of a small catchment with high debris activity, which has been thoroughly described by (Friele & Clague, 2005) and by (Jakob & Friele, 2010). It erratically experiences debris flows of magnitudes with volumes over one million m 3 which pose serious hazards for the local population, especially in the area of Brackendale (see Figure 7a and . The concentration of these flows may be 40% or more, which falls into the category of debris flows. When interacting with large volumes of dense flows, the Cheekye slit-dam is designed to trap such events with a basin capacity of 2.4 Mm 3 . It would be, to the best of our knowledge, the biggest debris flow basin in North America. The goal in designing the Cheekye slit-dam is for it to become blocked by large boulders during large, unsteady flow events, thus minimizing outflow downstream. It should be noted that in Canada, urban areas tend to be much less dense than places such as Japan or the Hong Kong SAR. This means that the some creeks may experience debris flow events of small or moderate magnitude without threatening assets. The erosion processes in the Cheekye River are very active (Friele & Clague, 2005;Jakob & Friele, 2010), and so its channel often transports a lot of sediment. The channel has the capacity to accommodate even moderate debris flow events without threatening assets (up to a volume of 0.4 Mm 3 , i.e., the 50 years return period event, as demonstrated by . Slit-dams such as the Cheekye slit-dam are installed to allow the passage of sediments during routine (safe) events, in addition to the passage of normal fluvial flows. This minimizes expensive and environmentally unfriendly activities for clean-up after routine events. Designing such structures is not straightforward because allowing even moderate debris flows to pass while trapping large ones is necessarily uncertain, hence the development of the approach presented in this paper.

Input Data
Here we discuss the hydrographs, the basic barrier design, the basin capacity, and binomial probabilities for boulder classes.

Hydrographs: Frequency-Magnitude Curve of Supplied Volumes and Peak Discharges
Consistent with the Cheekye Expert panel recommendation (Clague et al., 2014), the frequency-magnitude data assessed by Jakob and Friele (2010) are used with events classified according to seven volumes and associated peak discharges (see also Table S1 in Supporting Information S1). A comprehensive quantitative risk assessment showed the smallest dangerous event has a volume 0.8 Mm 3 and a peak inlet discharge of 3,400 m 3 /s. The largest plausible event that poses an unacceptable risk has a volume of 2.8 Mm 3 and a peak inlet discharge of 15,000 m 3 /s. Giving further explanation on rational ways to approach quantitative risk assessment on debris flow creeks is beyond the scope of this paper (but see Jakob & Holm, 2012;Strouth & McDougall, 2021 for more information).

Barrier Definition
A digital mockup of the Cheekye slit-dam is shown in Figures 7b and 7c, for perspectives of the rear and front of the dam, respectively. The full details of the constrictions of the slit-dam design are given in the Supporting Information (Table S2 in Supporting Information S1). The width of the slit is 6 m, whilst the crest is 100 m wide. The structure is between 24.6 and 28.6 m tall, since there is a weir-like feature above the slit (Figures 7b and 7c). It should be noted that the beams for the Cheekye dam will be manually movable, enabling "adaptive management" of flow hazards, but only between events. This is because the duration of the debris flows that the slit-dams are intended to arrest are rather short, relative to fluvial transport problems.
The slit dam has a concrete and bedrock base below an elevation of 275.4 m.a.s.l., and is thus unerodible. However, the Cheekye River longitudinal profile has a mean level about 1 m above this elevation. The Cheekye River has an alluvial bed in this section, so it is expected that an alluvial blanket will cover the bottom of the constriction. Depending on fluctuations in sediment supply and the local hydraulic conditions, the level of the alluvial blanket will undergo spatio-temporal variations. Expert assessment of the alluvial blanket thickness in the constriction estimated a blanket depth between 0 and 2 m, with a best estimate value of 1 m.
During major debris flow events, we assume that this alluvium blanket remains intact: the area slopes relatively gently, and the constriction promotes jamming of boulders. We reduce the height of the constriction according to the assumed alluvial bed. For lesser events, the alluvium blanket may vary substantially, but within the range of 0-2 m if bedload transport occurs. Given large wood jams and routine events, scouring of the alluvium blanket is expected; this is a type of beneficial self-cleaning behavior. In addition to the slit structure, the barrier is equipped with a spillway for events larger than the projected event, but modeling the spillway and associated overflow is beyond the scope of this paper. Here, we only describe the structure behavior for the 2.8 Mm 3 event, not for larger and less frequent events.

Basin Capacity
BGC analyzed LiDAR data for the basin area and appraised the storable volume for a given flow level at the constriction and the deposition slope. (See Table S3 in Supporting Information S1.)

Assessing Binomial Probabilities of Each Class of Boulders
Although it is rarely done, it is possible to estimate the quantity and size of boulders conveyed by debris flows at a given site, of course with some uncertainties. BGC provided such an estimate for the Cheekye barrier site to support modeling of debris flows with the barrier's slit outlet. The boulder count estimate was organized around the probable number of boulders of a given size class in a 100,000 m 3 reference volume. The estimate was developed by recording the size and quantity of boulders upstream and downstream of the barrier site, using aerial imagery, field surveys, and test pits in debris flow deposits. The test pit data included the estimated percentage of boulders within the overall deposit (2%-5%); these percentages were used to scale the boulder counts to the reference volume. The shape of the boulder grain size distribution was also cross-checked against a generalized debris flow boulder grain size distribution that was developed in Japan (Horiguchi & Richefeu, 2020). Although imprecise, we believe the boulder count estimate, with the wide uncertainty bounds provided, is an adequate representation of the range of possibilities. Table 1 in Section 5.2.2 summarizes information on these boulders. Figure 8 shows an example of model results for the 2.8 Mm 3 event. Each call to the code for a given set of input parameters yields this kind of results (see also other examples in Figure S2-S4 in Supporting Information S1). Note that Opening #1 always refers to the bottom opening; the index increases with height. For these examples, the height of the initial wood jam was randomly sampled in the distribution [2|4|6] m ([min | best estimate | max]).

Example of Design Event Modeling (With Horizontal Grillage Bars)
The first panel shows the discharge time series. The black line shows material entering the basin, and the gray line shows discharge through the constriction. The second panel shows the rate of vertical clogging of the sub-openings. The lowest initially available slit (Opening #2) progressively becomes jammed, with increasingly higher openings becoming jammed thereafter. The third panel shows the rate of clogging in slots. The fourth panel shows a time series of the flow level and the level of jamming in the slit; the slit does not become jammed all the way up to the crest. The last panel shows the storage volume time series for the basin. Peak storage occurs at the end of the inlet hydrograph. Discharge continues after the input hydrograph has finished, leading to continued progressive stochastic jamming. This continued discharge constitutes self-cleaning, which can continue beyond the data shown on the graph.

Uncertainty Propagation Analyses (Both With and Without Bars)
The run that is presented in Figure 8 is interesting, but the main value of the model is to analyze uncertainty propagation using many simulations. Table 1 gives the input parameter distributions for the possibilistic modeling of 2.8 Mm 3 events (see a graphical form in Figure S5 in Supporting Information S1).
The normalized peak lag, that is, the delay of the peak discharge compared to the event duration, is unknown.
These values are adopted: The deposition slope adopted is 4.7% following a site study. Site investigations showed that the deposits in the proposed debris basin are probably associated with the last major debris flow event of the Cheekye River (a so-called "garbage dump" event in Jakob & Friele, 2010). The deposit surface slope is about 7.1%, which is important because it demonstrates that this location is prone to naturally storing geological material. Furthermore, although an extended deposition was detected on the alluvial fan by Jakob and Friele (2010), that is, on a slope of about 4%-5%, a significant part of the debris flow material is likely to be able to deposit at a steeper angle. As such, the 4.7% slope is probably conservative.
Constrictions are prone to becoming clogged by woody debris, either before large debris flow events, or due to logs clustering at the surge front. The height of this initial log jam is a key unknown, and may vary between a few meters to more than 20 m in height depending on stochastic movements of woody debris based on observations in the field. This range is an expert estimate by the authors based on (a) the high height of the trees standing in the debris basin (see Figures 7b and 7c), as well as along the upstream confined debris flow channel (Jakob & Friele, 2010), and (b) evidence from field observations that large wood chunks (with sizes of several to more than ten m) are systematically transported at the front of large scale collapses transforming into debris flows (see pictures in Guthrie et al., 2012;Roberti et al., 2017;Friele et al., 2020). For jammed woody debris, evidence is lacking for a best estimate range, and even bounds are difficult to define. We thus proposed two scenarios to highlight the sensitivity of the outflow to this parameter: a "low" scenario regarding the clogging of the barrier where the large wood jam would be rather small (a range of 2-6 m, with a most probable value of 4 m: ] ); and a "high" scenario with a rather high jam: [5|10|15] m. (The hypothesis will be further discussed in §5.2.2.) To extend the analysis, worst-case scenarios were computed without any large wood bottom jam, although such events are considered highly unlikely.
The number of boulders of each class is uncertain, and is defined according to the data provided in Table 1. Each class is considered as an independent parameter, so, for instance, a relatively high number of class-1 boulders can be selected for one run along with a relatively lower number of class-2 boulders (but always separating an optimistic range and a pessimistic range from the best estimate to the upper-or lower bound respectively, as explained in Figure 6).
To highlight the uncertainty on the output parameters associated with particular hypotheses regarding the number of boulders of certain classes -for instance, the effects on the released volume of material -two p-boxes were computed for each run in Figure 9. One p-box accounts for the full uncertainty on the number of boulders (shown using continuous lines in all p-boxes) and another one considers only the best-estimate numbers of boulders (shown using dotted lines in all p-boxes). Several sets of parameters were also considered independently for this analysis: the constriction with and without horizontal grillage bars (gray vs. blue p-boxes in Figure 9); non-existent, small and large initial wood jams (rows of Figure 9); and 0.8 and 2.8 Mm 3 events (columns in Figure 9). Figure 9 show the p-boxes of the released volume for these cases.
We chose to focus on these event classes for two reasons: 1. An event with a volume of 0.8 Mm 3 is the most likely event that is assessed as posing a threat to lives. Events of volume lower than 0.4 Mm 3 are non-lethal "routine events" sensu Camiré et al. (2019), intermediate values represent a transition. 2. An event with a volume of 2.8 Mm 3 is the biggest plausible design event. Event class-7 is the biggest plausible event, but risk associated with it is tolerable because it is so unlikely. An event with a volume of 2.8 Mm 3 is the design event because it is the largest event posing an unacceptable risk. The intended function of a slit dam is to reduce the outflow volume of design events to be lower than 0.8 Mm 3 and if possible than 0.4 Mm 3 .
The p-boxes with continuous lines in Figure 9 are much wider than the p-boxes with dotted-lines, testifying to the importance of the number of boulders for the functioning of the structure. We first focus on the efficacy of the dam without bars (blue p-boxes). The released volume computed for event 0.8 Mm 3 event is likely relatively high, with the middle part of the p-box with the continuous line ranging from about 0.2 to 0.7 Mm 3 . This suggests that it should be partially or fully self-cleaning. The middle range of the p-box with the continuous line ranges from 0.5 to 0.6 Mm 3 for a low large wood jam, implying that it releases more material than 0.4 Mm 3 even if the boulders are present according to the best estimates. For the hypothesis of a high wood jam (bottom row), the lower bound of the release volume is slightly lower, but remains around a 0.4 Mm 3 event. Essentially, too much flow material passes through the structure if the slit is not equipped with bars, even for relatively frequent events.
The conclusion is similar for the event with a volume of 2.8 Mm 3 releasing 0.4-0.9 Mm 3 of debris flow for boulder numbers according to the best estimates (p-box with the dotted line) and up to a range of 0.2-1.5 Mm 3 when accounting for the full uncertainty on the boulder numbers. Here again, the objective of the barrier is to moderate the debris flow volume of 2.8 Mm 3 events such that the released volume is not life-threatening, that is, such that the outflow remains significantly below 0.4 Mm 3 . These simulations shows that this objective is not achieved using this design option, that is, a 6-m wide slit without bars.
As for the case with horizontal grillage bars (gray p-boxes), the behavior of the barrier was drastically changed. Each interval between bars was (and should be) treated as a single sub-opening that can become progressively obstructed (see second and third panels of Figure 8). The self cleaning is consequently dramatically lower and also slower.
The bottom panels show that acceptably low release volumes are probable for both events of volume 0.8 and 2.8 Mm 3 , if horizontal grillage bars are installed. The hypothesis regarding the number of boulders in the case of high deposits of initial wood jams (>5 m) plays a role only for the 2.8 Mm 3 event, decreasing the value of the worst-case results (the upper tail of the upper bound). Almost no differences appear for the 0.8 Mm 3 event event; indeed, the continuous-and the dotted-lines p-boxes are virtually the same.
The central panels show that these conditions should satisfy the design requirements. For the events of both volume of 0.8 and 2.8 Mm 3 , the released volume of flow material is very probably lower than or equal to 0.4 Mm 3 , providing that the boulder numbers are at least as numerous as the best estimate (the dotted-line p-box). In case their number is closer to the lower bound, the released volume could approach the upper bound of the gray p-box (the red line). Its median would still satisfy the design conditions since it is also closest to 0.4 Mm 3 . However, its upper tail approaches a volume of material released through the constrictions equal to the supply for the 0.8 Mm 3 event event, and is possibly higher for the 2.8 Mm 3 event.
For the worst-case scenario without large chunks of wood becoming trapped at the bottom of the constriction, one can only rely on boulders to jam the constriction. For this case, the released volume is still probably acceptable providing that the number of boulders is around the best estimate(s) or higher (see the p-box with the dotted line). However, given fewer boulders, the volume of material released through the constriction remains substantially higher than the acceptable value, that is, much more than 0.4 Mm 3 . We assume that this last scenario (with no initial wood jam and very few boulders) is very improbable, and suggest neglecting it for design work. Notwithstanding, this computation demonstrates that for this unlikely case, the release volume approaches or possibly becomes higher than the 0.8 Mm 3 event event, but that a part of the supplied material is still trapped.

Value of Our Approach
Our new model requires quantifying the number and diameters of boulders within a given volume of deposited material, to be used as an input parameter for the generation of binomial distributions describing boulders. The advantages of this type of boulder counting include that the data are factual, observable and relatively simple to measure, and can be collected during typical post-event dredging operations. Nevertheless, the average number of boulders per volume of debris flows can vary between events. Thus, even in convenient cases where a precise count of boulders is available, our new model enables the sensitivity of the model results to the variability of these input data to be assessed. Indeed, it is the stochasticity of the process and the epistemic uncertainty on the input parameters that are the underlying motivation for the implementation of the possibilistic approach presented in this manuscript.
It is also worth highlighting that despite their clear value for problems involving many unknowns that are difficult to characterize, possibilistic and hybrid probabilistic-possibilistic approaches have yet to become widely used in the geosciences or in the field of natural hazards (Baudrit et al., 2006;Guyonnet et al., 2003;Rohmer & Verdel, 2014). Although they are well-suited to such fields of application, where some model input data are systematically imperfectly known, they remain seldom used. We believe this is partially because they are not easily implemented from scratch, and so we emphasize the immense utility of the HYRISK R-package of Rohmer et al. (2018). Another reason for the lack of widespread adoption of possibilistic approaches is likely that existing explanations of their functioning tend to be insufficiently intuitive. We believe this to be particularly true for their output results, for example, the p-boxes, as well as how they are to be used and interpreted. We hope that this paper will highlight other potential applications of possibilistic approaches to other areas of geosciences and natural hazard mitigation.
The design of large debris basins, as in our case-study of the Cheekye barrier, is challenging our knowledge and understanding of interactions between debris flows and structures. Systematically taking conservative decisions would have resulted in severe over-sizing of the structure, which would then have been prohibitively expensive (Strouth & McDougall, 2021). We developed the approach presented in the present paper specifically for the Cheekye barrier because we needed an approach that was (i) capable of modeling the main processes that affect the functioning of the barrier, whilst (ii) being sufficiently fast to run to perform error propagation analyses. As hitherto stressed, it can be used to study other cases where debris flows pass through and interact with constrictions. The approach enabled the study of the barrier efficacy in trapping a sufficient volume of debris for the project design events. Considering the complexity of the processes and the many remaining uncertainties, precise results are not obtainable. As shown in our interpretation of the p-boxes of Figure 9, the approach is welladapted for providing ranges of estimates of the release volumes, enabling engineers and scientists to take design decisions.
In short, our new program is adept at quantifying the range of possibilities for flow-structure interaction. We now supplement a short discussion on how to apply this range of possibilities for the design of engineered constrictions. Our experience in France, confirmed by many informal discussions with colleagues from other countries, is that given all these uncertainties, for engineered channel constrictions, it is better to err on the side of having constrictions and/or sub-openings that are too large. Excess discharge can later be controlled with extra bars and/ or flexible barrier nets. The option for these additions can be "baked into" a modular initial design of a slit dam. This is preferable to having a slit which is too narrow and persistently clogs up even for small events, causing large ongoing maintenance costs (Ballesteros-Cánovas et al., 2016), or requiring a complicated cutting operation (Chiu et al., 2021). For instance, around a third of debris basins designed in France lead to either excessive maintenance costs or even abandonment of the structure (Carladous et al., 2022).
In addition to such analyses, one can record other outputs from the model: statistical distributions of the maximum flow level helping to design the barrier crest, sizes of the boulders jammed to refine and adapt the bar spacing (e.g., Figure S6 in Supporting Information S1), or the hydrograph released downstream of the constriction for later reuse in a residual hazard study (e.g., Figure S7 in Supporting Information S1).

Possible Extensions to the Method
A key novelty proposed in the present paper is the simple approach adopted for randomly computing the presence of boulders in debris flows. It is used in our case to update the geometry of constrictions using simple jamming laws. Nevertheless, extended applications could also be envisaged. For instance, our new model could be coupled with a depth-averaged model to compute propagation of debris flows. Our new model could also be extended to enable the computation of impact forces against barriers or buildings, including the effect of boulders (Hubl et al., 2017;Nagl et al., 2020). A further step could also be to include large wood chunks in the approach to compute compound jams made of boulders and wood chunks. This would require the development of slightly more sophisticated jamming conditions to account for (a) their elongated shapes and (b) their tendency to float at the free surface. Existing numerical models explicitly modeling the transport of large wood coupled with flows can be a source of inspiration in this direction (Ruiz-Villanueva et al., 2014. Since the interactions between constrictions and flows laden with boulders and/or woody chunks are fundamentally stochastic, our novel approach of performing several simulations to capture the stochastic variance of the results will remain necessary.

Assumptions and Limitations
The newly developed program involves certain assumptions and limitations of which potential users ought to be aware.
First, we assume that stable arching only occurs for one or two boulders with a cumulative diameter greater than the smaller opening size. In fact, arches can also form when three or more boulders interact with an opening simultaneously (Choi et al., 2016;Marchelli et al., 2019). These cases are neglected, since the probability of stable arch formation is low; the vast majority of longer arches are, at best, quasi-static.
Second, the statistical random generation of boulders is performed a number of times equivalent to the number of packets. It is thus computed based on the total volume passing through an opening at a given time step, which was set to one second. This was assumed to be a relevant approximation of the time it takes for boulders to pass the constriction. Boulders are thus only "eligible" for clogging at a single timestep, regardless of their size or the rate of outflowing material. This means that jams involving large boulders that may take several seconds to pass through the openings could be underestimated.
Third, the numerical analysis considers the sub-openings of the constriction independently. This assumption is reasonable for debris flow breakers which include many parallel openings, which can behave independently depending on the spacing of the sub-openings (see Goodwin and Choi (2020)). However, if openings are stacked vertically, as with the Cheekye barrier, the model can produce physically implausible results: upper slits can be clogged whilst lower slits remain open. This phenomenon may generate questionable outcomes during the long-lasting self-emptying phase. We thus advise against extrapolating too much from the output for self-cleaning, especially since other basic features of the self-cleaning of basins are neglected in our reservoir model. These processes are intrinsically 2D and thus cannot be well captured in fewer dimensions, as stressed by Armanini and Larcher (2001).
Fourthly, the input parameters for the numerical model, including the estimated bounds and most likely values of any given parameter, are subject to expert assessment, rather than directly considering field data. The newly developed numerical model only gives an aggregate output of these estimates for a wide range of values considered simultaneously. In particular, an important next step will be to improve the characterization and modeling of the wood jam, perhaps in a way analogous to the way that boulders are currently handled. Relevant literature on wood-jams includes NILIM (2007) Fifthly, the probability p j that each packet is a boulder is assumed constant for the entire flow. Migration and segregation effects (Zhou & Ng, 2010) are thus neglected, as is runup of frontal surges over the barrier. Direct-or post-event observations of debris flow deposits suggest that debris flow fronts, laden with large boulders when reaching debris basins, tend to become deposited when losing the confinement of the upstream channel (Piton et al., 2018). We thus hypothesize that boulders approach man-made constrictions as random groups rather than as a well-defined bouldery front. This approximation may not apply for canyons or barriers which are not located downstream from wide deposition basins.
Finally, our new program does not model the flow explicitly. Other numerical methods such as Hyper KANAKO Yanagisaki et al., 2016), which implements depth-averaged equations for flows, could be used in a complementary fashion. This type of complementary modeling would be especially helpful for checking the edge cases (i.e., "best" and "worst"-case scenarios). Of course, it should be noted that the definition of "sediment" is not necessarily the same across different regions; in particular, Hyper KANAKO only models two sizes of sediment within a particular flow, whereas our approach allows the consideration of an arbitrary number of sizes (via the boulder classes).

Conclusions
Boulders becoming jammed in structures such as slit-dams is a concern to scientists and engineers in mountainous regions worldwide. Ideally, boulders in major flows should jam the constriction, limiting outflow, whereas smaller flows should pass harmlessly, enabling self-cleaning. Nevertheless, small variations in the characteristics of the flow, boulders or the constriction can lead to very different outflow behavior. Existing methods are unable to dynamically account for potential changes in the geometry of constrictions due to a range of boulder sizes, despite their governing role. This paper has detailed an entirely new method for determining the possibility of boulders becoming jammed in constrictions. The method uses analyses of flow and basin properties from experts with a new framework coupling statistics and analytic outflow equations for flows interacting with slit-dams. We implemented this framework into a simple numerical model. Boulder jamming was modeled using a binomial law which randomly generates boulders. Input for the binomial law was obtained from new field data in which the distribution of boulder sizes was characterized. We dealt with uncertain parameters, such as the initial clogging of the constriction due to woody debris, using an uncertainty propagation approach.
The new model was used to analyze a proposed slit-dam, where field data for boulder size distributions was available. The analysis considered cases both with and without horizontal grillage bars. The new model showed that, absent significant jamming of large chunks of wood, a 6-m wide slit would probably release an excessive amount of debris flow material during projected events. The analysis suggests that adding bars across the constriction would drastically decrease the released volume during large events. The analysis also quantified the probability of a particular volume of material being over the risk threshold. We found the probability to be reasonably low for designs incorporating bars. A designer could use our model to further refine the design of a slit-dam, either to minimize the possibility of exceeding the risk threshold for all cases, or to accept higher risk as a trade-off for improved long-term self-cleaning.

Data Availability Statement
Data for the debris basin geometry, the proposed slit-dam geometry and the boulder size distributions in the field were provided by BGC Engineering. They can be found is the Supporting Information of this manuscript. The source code and data files can also be found at https://github.com/GuillaumePitonInrae/CheekyeDebrisFlowBarrier.git.