A modeling methodology to study the tributary-junction alluvial fan connectivity during a debris flow event

Traditionally, interactions between tributary alluvial fans and the main river have been studied in the field and in the laboratory, giving rise to different conceptual models that explain their role in the sediment cascade. On the other hand, numerical modeling of these complex interactions is still limited because the broad debris flow transport regimes are associated with different sediment transport models. Even though sophisticated models capable of simulating many transport mechanisms simultaneously exist, they are restricted to research purposes due to their high computational cost. In this article, we propose a workflow to model the response of the Crucecita Alta alluvial fan in the Huasco Valley, located in the Atacama Desert, Chile, during an extreme storm event. Five different deposits were identified and associated with four debris flow surges for this alluvial fan. Using a commercial software application, our workflow concatenates these surges into one model. This study depicts the significance of the mechanical classification of debris flows to reproduce how an alluvial fan controls the tributary–river junction connectivity. Once our model is calibrated, we use our workflow to test if a channel is large enough to mitigate the impacts of these flows and the effects on the tributary–river junction connectivity.

1 Introduction

In arid and semi-arid regions, extreme storm events commonly trigger debris flows with high potential to modify the landscape (Mather and Hartley, 2005; Michaelides and Singer, 2014). The sediment stored in the catchments is transported towards tributary-junction alluvial fans and main rivers during these events. Tributary-junction alluvial fans play a crucial role in the transference of sediment from source areas to the main river (Fryirs, 2013; Aguilar et al., 2020). The efficiency of this sediment transference depends on the degree of connectivity within the fluvial system (Brunsden and Thornes, 1979), in which tributary-junction alluvial fans can modulate the sediment transference towards the main river by buffering or bypassing the sediment discharge, creating different coupling conditions that depend on the combination of the fan's geomorphological configuration and the flow properties (Heckmann et al., 2018; Savi et al., 2020). The degree of connectivity between tributaries and the main river determines the down system transmission of sediments and water and, consequently, its sensitivity to any environmental change (Mather et al., 2017; Cabré et al., 2020a). The concept of “sediment connectivity” or “(dis-)connectivity” implies a balance between top-down and bottom-up processes. Top-down processes are controlled by climatic and geological variabilities that directly impact the sediment supply. On the other hand, bottom-up processes are controlled by morphological characteristics such as base-level-driven processes (Mather et al., 2017; Heckmann et al., 2018). Understanding the fan dynamics and fan–river interactions is necessary to comprehend the sediment cascade during debris flow events. However, these processes are still not properly modeled in hazard assessment projects (Savi et al., 2020).

The alluvial fan's geomorphic configuration may help us foresee their coupling conditions. For example, mild slope fans or the absence of a feeder channel connecting the tributary with the main river can promote sediment deposition, thus buffering the sediment discharge (Mather et al., 2017; Zegers et al., 2017). Hence, alluvial fans can act as a sediment decoupler feature in the sediment cascade, preventing lateral sediment discharges from reaching the main river (Fryirs et al., 2007). In contrast, alluvial fans with a trimmed fan toe, known as “toe cutting” (Leeder and Mack, 2001), enhance alluvial fan trenching and consequent generation of lobes at the fan toe. Telescopic-like deposit morphologies result from a local increase in the sediment transport in the alluvial fan, followed by an alluvial fan progradation (Colombo, 2005). In this case, the lateral sediment discharge becomes completely coupled with the river. However, when this progradation causes a river blockage, it acts as a barrier by decoupling the longitudinal sediment discharge along the river (Fryirs et al., 2007). The alluvial fan's geomorphic configuration and its connectivity can also be affected by anthropogenic activities and infrastructure. For example, in debris-flow-prone areas, hydraulics works are used to retain sediments (check dams) and/or to safely conduct the flows (channels). These works artificially change the coupling condition between the tributary fan and the river.

Debris flows are influenced by forces arising from particle–particle, fluid–particle, and fluid–fluid interactions (Iverson, 1997). When the interstitial fluid is a slurry composed of a high proportion of fine sediments, the density difference between the fluid and particles is small, allowing coarse sediment to flow at the same velocity as the slurry (Takahashi, 2014). These flowing mixtures are known as viscous debris flows and are well represented by a single-phase viscoplastic rheological model such as Bingham or Herschel–Bulkley rheologies (Takahashi, 2014; Naef et al., 2006; Montserrat et al., 2012). The turbulent stress dominates the flow behavior for higher velocities, and Manning- or Chézy-type relations provide good results (Naef et al., 2006; Takahashi, 2014). Conversely, when the interstitial fluid is mainly water, coarse sediment and water move at different relative velocities (von Boetticher et al., 2016). In this case, grain collisions, or dispersive stress, dominate the bulk flow energy dissipation (Takahashi, 2014; Naef et al., 2006). Similar to turbulent stress, collisional stresses are proportional to the square velocity (inertial stress); thus, pseudo-Manning or Chézy-type approaches have been used to estimate bulk stress (Naef et al., 2006; Rickenmann et al., 2006). In these approaches, pseudo-Manning or Chézy coefficients could be functions of grain types and size and volume concentration, among others (Naef et al., 2006; Rickenmann et al., 2006). For slow and/or highly concentrated flows, instead of collisions, particles can experience long-lasting contacts. In this case, the Coulomb plasticity model has been shown to be a good approach for modeling interparticle frictional stress (Ancey, 2007; Montserrat et al., 2012). Coulomb stress operates similarly to yield stress in Bingham or Herschel–Bulkley rheological models; thus, it can be used as a surrogate. In numerical models, Coulomb and yield stresses allow debris flows to stop, replicating the generation of alluvial deposits. Existing flow resistance relationships for debris flow modeling basically combine viscous, yield and coulomb, and turbulent and dispersive stress in a single equation accounting for bulk frictional losses (Naef et al., 2006).

One-phase models are preferred when representing bed material entrainment and deposition in numerical models due to computational costs. The incorporation of sediment in the flow produces changes in the rheological properties and an increase in the volume of the mixture (Iverson, 1997; Naef et al., 2006; Zegers et al., 2020). Entrainment and deposition continue to be pivotal components for modeling and reproducing the topographic adjustments during debris flow events (Cao et al., 2004). However, there is still no consensus on the parameters that govern the sediment entrainment. For example, Takahashi et al. (1992) proposed a method that incorporates sediment into the flow until reaching an equilibrium concentration C_∞. Cao et al. (2004) proposed entrainment and deposition relationships as a function of the shear stress, volume concentration, size of the particles, and the flow height and velocity, together with three coefficients that need to be calibrated. McDougall and Hungr (2005) represent sediment entrainment with an exponential growth parameter E, which is independent of the flow velocity.

The interaction between tributary alluvial fans with the main river has been studied in the field (Mather et al., 2017; Cabré et al., 2020a) and in the laboratory (Savi et al., 2020), giving rise to different conceptual models that explain its role in the sediment cascade (Mather et al., 2017; Cabré et al., 2020a; Savi et al., 2020). On the other hand, numerical simulation is limited because available models cannot tackle all the possible flows types and sediment transport processes occurring in an alluvial fan and its interaction with a river. Some approaches, like r.avaflow, have integrated many essential physical phenomena into one generalized debris flow model (Pudasaini, 2012; Mergili et al., 2018). However, the lack of sound field evidence has limited the validation of these models.

From an engineering perspective, this work addresses the following question. Are we able to numerically model tributary-junction alluvial fan connectivity changes during a debris flow event? To answer this question we used the FLO-2D commercial software application due to its reasonable computational cost and its broad applicability for hazard assessment projects.

Figure 1 Study area. (a) Crucecita Alta catchment (28.895569^∘ S, 70.449925^∘ W). The gray polygon is the river segment that includes the studied fan topography section used in the FLO-2D numerical model. The unfilled polygon depicts the Crucecita Alta catchment (13 km²). (b) Crucecita Alta alluvial fan main geometric features where the feeder-channel, generated during the 25M event, is marked with dotted lines. Map data: Google Earth Pro (CNES/Airbus) taken in October 2016 (post event).

In March 2015, an extreme rainstorm event occurred over a large area of the Chilean Atacama Desert (Bozkurt et al., 2016; Wilcox et al., 2016; Jordan et al., 2019; Cabré et al., 2020b), hereinafter called the 25M event. This study benefits from a detailed sedimentological field characterization of the debris flow deposits on the Crucecita Alta fan for the 25M event (Cabré et al., 2020a). The Crucecita Alta catchment is a tributary of the El Carmen River in the Huasco Valley (∼29^∘ S, 70^∘ W). The conceptual connectivity model proposed by Cabré et al. (2020a) is based on the characteristic storm signature of the 25M event, which resulted in characteristic debris flow surges. We used these data to differentiate two main sediment transport modes and assigned them to different sediment transport models for each surge. We developed a workflow that includes the calibration of the flow's rheological parameters and a Python routine to concatenate different numerical models for the different debris flow surges. We used the geomorphic changes that occurred in Crucecita Alta to test the suitability of our workflow in the reproduction of the 25M event on this alluvial fan.

Inhabitants of the Huasco Valley tend to dwell in the alluvial fans because of their gentle slopes and because alluvial fans are safe places during river floods. However, high-magnitude debris flow events directly impact these populated areas. To reduce the impacts of debris flows, mitigation works in such environments are divided into sediment retention (pools, barriers, etc.) and/or bypass strategies (channel, levees, etc.). After reproducing the 25M event, we used our workflow to simulate different behaviors when artificially changing the fan river connectivity by means of a channel under different configurations of the debris flow surges.

2 Study area and data

Tributary-junction alluvial fans in the Huasco Valley, located between 1500 and 3000 m a.s.l. (above sea level), supply sediments to the main river. They are so abundant that they can reach densities of 1.9 fans per kilometer. During the 25M event, rainfall gauges in the valley recorded precipitation ranging from 20 to 76 mm in 3 d, with a maximum intensity of 16 mm h⁻¹. The high elevation of the snow line in March 2015 (3200 m a.s.l.), 400 m higher than its average altitude (Lagos and Jara, 2017), exposed areas usually covered by snow, thus increasing run-off (Wilcox et al., 2016; Jordan et al., 2019).

Many of the tributary catchments in the Huasco Valley were activated during the 25M event causing casualties and great economic losses (Izquierdo et al., 2021). The alluvial fans present a characteristic sedimentological response to the 25M event, in which a sequence of surges impacted the fans repeatedly (Cabré et al., 2020a). The fieldwork was performed before any mitigation or reconstruction works in this fan were done, so the data were unaffected by anthropogenic activity. Therefore, the flood sequence of the 25M event in this fan was reported by Cabré et al. (2020a).

2.1 The Crucecita Alta fan

The studied fan is situated at the river junction of an ephemeral catchment with the main river (∼28.9^∘ S, 70.4^∘ W) (Fig. 1a). Its catchment has an area of ∼13 km², a length of 6 km, and a mean slope of 30 %. The river junction is at 1044 m a.s.l., and the maximum height of the catchment is 3129 m a.s.l.; consequently, the whole catchment was under the snow line elevation during the 25M event. Catchment lithology is dominated by volcanic rocks (andesites), conglomerates, sandstones, mudstones, and by large accumulations of unconsolidated sediments in hillslopes and alluviated channels (Cabré et al., 2020a). The Crucecita Alta fan has an area of 0.14 km² and a mean slope of 8.6 %, and it hosts a few houses in the northern area (Fig. 1b). The northern area of the fan was not affected during the 25M event presumably due to the presence of a 5 m high deflection barrier made of unconsolidated debris. Before the event, riverbank erosion trimmed the alluvial fan toe (Cabré et al., 2020a), disrupting the alluvial fan gradient (Fig. 1b). The difference in elevation between the main river and the alluvial fan surface was 13.5 m. This geomorphic configuration promotes the fan entrenchment and the formation of new fan lobes at its toe during debris flow events. These lobes may act as barriers because the valley is narrow (200–300 m wide).

A post-event topography of a 51 km long segment of the El Carmen River valley was acquired with 1 m × 1 m horizontal resolution between February and March 2017 by the Chilean Ministry of Public Works as part of a debris flow mitigation project in the area (IDIEM, 2019). The vegetation and buildings were removed from this available topography. Figure 1a shows the river segment where the Crucecita Alta fan is located. In the pre-event satellite imagery, retrieved from © Google Earth Pro v.7.3.3.7786 (CNES/Airbus images), there is no evidence of a feeder channel connecting the fan apex with the main river.

A calibrated HEC-HMS hydrological model (USACE, 2015) of the entire El Carmen River basin is also available from the same mitigation project (IDIEM, 2019; Zegers et al., 2020). Water flow discharge at the catchment outlet and at the river before the junction was obtained from this hydrological model. The obtained hydrograph shows that the flood event consisted of four main surges, with a maximum peak flow of ∼7 m³ s⁻¹ during the fourth surge (Fig. 2b).

2.2 Characteristics of the 25M event

There has been a consensus that storms like the 25M event in the Atacama Desert occur once a century (Ortega et al., 2019). About the magnitude of the event, Aguilar et al. (2020) estimated the mean erosion rate of the Huasco basin during the 25M event to be equal to 1.3 mm for the entire catchment. On the other hand, Aguilar et al. (2014) estimated erosion rates of 0.03–0.08 mm yr⁻¹ during the last 6–10 Myr for the same catchment. Assuming that an event such as the 25M event occurs once a century, we can say that the 25M event has a “high magnitude” because this single event contributed from 15 % up to 40 % of the total sediment volume eroded, on average, every 100 years.

Cabré et al. (2020a) characterized the flood event sedimentology and spatial distribution of the five detected deposits in the Crucecita Alta alluvial fan. These deposits were mapped and classified into facies (Fig. 2a). The different facies were named from F1 to F5, in which F1 was the first deposit during the storm, and F5 was the last one. “Sedimentary facies” or “facies” is a concept widely used in sedimentological studies because it assists in summarizing the grain shape, textural parameters, the relief, and the stratification type onto a single area or zone that can be mapped (Wells and Harvey, 1987). As described by Cabré et al. (2020a), F1 and F2 correspond to debris flow deposits associated with matrix-supported flows with high fine sediment concentrations. F1 is the thickest deposit (above 1 m) reaching near 7000 m³ of deposited sediment volume (Aguilar et al., 2020). The sediment deposition occurred in a narrow area in the upper portion of the fan without reaching the main river. F2 is a wider but less thick deposit (10–30 cm) that overlays F1 and covers all the fan's length, from the apex to the distal toe, with nearly 5000 m³ of deposited sediment volume (Aguilar et al., 2020). F1 and F2 deposits overlay pre-event sediments, indicating that the associated flows had no or negligible erosion capacity (Cabré et al., 2020a).

Figure 2 Available data for the 25M event in Crucecita Alta alluvial fan. (a) The facies F1 (magenta), (b) F2 (blue), (c) F3 (green), and (d) F4 and F5 (orange) retrieved from Cabré et al. (2020a) and lidar topography surveyed by IDIEM (2019). (e) Flow hydrograph obtained from the hydrologic model performed by IDIEM (2019). The colors used to identify the facies in (a–d) depict their correlation with the surges in (e).

Following facies were formed by inertial debris flows. These flows can be subcategorized into debris floods, which are water-driven floods with high bedload transport of gravel to boulder-size material in which erosion and sedimentation processes become important (Church and Jakob, 2020). The debris flood surges generated the F3, F4, and F5 deposits and were responsible for the alluvial fan entrenchment and consequent generation of the feeder channel (Fig. 2b). F3 deposits are associated with cohesionless transitional flows (Cabré et al., 2020a). These flows were channelized below the middle portion of the fan and formed the feeder channel. During the formation of F3, sediment deposition mainly occurs in the fan's upper part and after the fan toe in the El Carmen river. After this phase, more dilute flows erode the previously deposited materials and enable the deposition of F4 and F5 facies further downstream. Feeder channel depth ranges from 100 to 350 cm in the mid zone and reaches up to ∼500 cm in the fan toe (distal zone). The local base level, controlled by the main river, was reached during the incision. Therefore, the subsequent flows were deposited in a new lobe at the fan toe, exhibiting a telescopic-like morphology with open framework boulders and lobate gravel lenses. These deposits correspond to inertial debris flows with clast-supported fabrics and a matrix-free top (Cabré et al., 2020a).

In our analysis of the 25M event, the cross-cutting relationships of the facies mapped in Crucecita Alta by Cabré et al. (2020a) (F1 to F5) are fundamental evidence to determine the sequence of flows. Cabré et al. (2020a) interpret the significance of the mechanical classification of debris flows and classify F1 and F2 as non-Newtonian flows. In contrast, F4 and F5 are classified as Newtonian flows, which is a similar differentiation, to some extent, to the viscous and inertial debris flow classification of Takahashi (2014). F3 is a transitional flow between the observed non-Newtonian and Newtonian flows.

3 Methodology

3.1 Governing equations

We use the two-dimensional FLO-2D model to solve the flood wave progression for water and debris flows in complex topographical terrains (O'Brien and García, 2009). FLO-2D solves two-dimensional depth-averaged continuity and momentum equations, Eqs. (1) and (2), known as Saint-Venant equations, for water and debris flows.

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial h{V}_x}{\partial x}}+{\displaystyle \frac{\partial h{V}_y}{\partial y}}=\mathrm{0},\text{(2)}& {\displaystyle }\begin{array}{rl}& S_{{\mathrm{f}}_x}=S_{{\mathrm{0}}_x}-{\displaystyle \frac{\partial h}{\partial x}}-{\displaystyle \frac{V_x}{g}}{\displaystyle \frac{\partial V_x}{\partial x}}-{\displaystyle \frac{V_y}{g}}{\displaystyle \frac{\partial V_x}{\partial y}}-{\displaystyle \frac{\mathrm{1}}{g}}{\displaystyle \frac{\partial V_x}{\partial t}}\\ & S_{{\mathrm{f}}_y}=S_{{\mathrm{0}}_y}-{\displaystyle \frac{\partial h}{\partial y}}-{\displaystyle \frac{V_x}{g}}{\displaystyle \frac{\partial V_y}{\partial x}}-{\displaystyle \frac{V_y}{g}}{\displaystyle \frac{\partial V_y}{\partial y}}-{\displaystyle \frac{\mathrm{1}}{g}}{\displaystyle \frac{\partial V_y}{\partial t}},\end{array}\end{array}$$

where h is the flow height, V_x and V_y are flow velocities in the x and y directions, and $S_{{\mathrm{0}}_x}$ and $S_{{\mathrm{0}}_y}$ are the channel slopes in each transverse direction (x,y). $S_{{\mathrm{f}}_{x,y}}$ denotes the frictional slope, accounting for flow resistance. In the case of water flows, $S_{{\mathrm{f}}_{x,y}}$ is calculated using the Manning equation. In the case of mudflows, $S_{{\mathrm{f}}_{x,y}}$ is calculated using the so-called quadratic rheology model, Eq. (3), which combines different terms accounting for yield and coulomb, viscous, and turbulent and dispersive stresses (O'Brien and García, 2009).

$$\begin{array}{}\text{(3)}& S_{\mathrm{f}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{y}}}{{\mathit{\gamma }}_{\mathrm{m}}h}}+{\displaystyle \frac{K\mathit{\eta }V}{\mathrm{8}{\mathit{\gamma }}_{\mathrm{m}}{h}^{\mathrm{2}}}}+{\displaystyle \frac{n_{\mathrm{td}}^{\mathrm{2}}{V}^{\mathrm{2}}}{h^{\mathrm{4}/\mathrm{3}}}}\end{array}$$

In Eq. (3), τ_y is the yield stress, γ_m the specific weight of the solid–liquid mixture, h the element (cell) flow depth, K a resistance parameter for laminar flow, η the dynamic viscosity of the fluid phase, V the element (cell) flow velocity, and n_td a pseudo-Manning coefficient corrected by the sediment concentration. O'Brien and García (2009) suggest the following empiric relationships, Eq. (4), for estimating the parameters η, τ_y, and n_td as functions of sediment volume concentration, C_V (O'Brien et al., 1993; O'Brien and García, 2009).

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathit{\eta }={\mathit{\alpha }}_{\mathrm{1}}{e}^{{\mathit{\beta }}_{\mathrm{1}}{C}_{\mathrm{V}}}\\ & {\mathit{\tau }}_{\mathrm{y}}={\mathit{\alpha }}_{\mathrm{2}}{e}^{{\mathit{\beta }}_{\mathrm{2}}{C}_{\mathrm{V}}}\\ & n_{\mathrm{td}}=nb{e}^{m{C}_{\mathrm{V}}},\end{array}\end{array}$$

where α_1,2 and β_1,2 are empirical coefficients that have to be calibrated, n is the conventional Manning coefficient, b=0.0538, and m=6.0896.

The model also used a detention coefficient (SD) that controls flow detention. D'Agostino and Tecca (2006) suggest that SD works as a minimum physically plausible flow depth. In addition, Zegers et al. (2020) reported that SD is one of the two more sensitive parameters in the model, and, even though it is not part of the rheological model, it can be a surrogate for the flow rheology.

FLO-2D is also capable of solving sediment transport equations to simulate mobile bed and topographic adjustments. In the FLO-2D model, the coupling of the sediment transport computation and the flow hydraulics are solved in a two step procedure; i.e., water flow characteristics are computed first, and then sediment transport and bed geometry changes, due to sediment erosion or deposition, are computed for that time step (O'Brien and García, 2009). However, no erosion/deposition processes can be included when using the mudflow module.

Figure 3 Topographic data modifications. Post-event topography corresponds to the available lidar topography while the pre-event topography is a restitution based on satellite images and the available topography. The synthetic channel (dashed brown line) attached to the lidar topography is 550 m long, while feeder channel (dashed orange line) is 450 m long. In the post-event topography, the feeder channel is the result of the inertial debris flow incisions on the alluvial fan. In the mitigation works' topography, the feeder channel was replaced by a straight rectangular channel.

The sediment transport model accounts for the local deficit/excess of transported sediment, modifying the bed channel height using the well-known Exner equation, Eq. (5).

$$\begin{array}{}\text{(5)}& {\displaystyle \frac{\partial h_{\mathrm{0}}}{\partial t}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-p}}\left({\displaystyle \frac{\partial q_{{\mathrm{s}}_x}}{\partial x}}+{\displaystyle \frac{\partial q_{{\mathrm{s}}_y}}{\partial y}}\right)=\mathrm{0}\end{array}$$

In Eq. (5), h₀ is the height of the channel bed, q_s is the sediment load, and p is the porosity. For the sediment load (q_s), FLO-2D has 11 different expressions based on unique river conditions (O'Brien and García, 2009). For this study, we used the Parker–Klingeman–Mclean equation (Parker et al., 1982), Eq. (6). This equation is suitable for gravel and sandy bed material and simulates sediment transport mechanisms from bedload to mixed suspended sediment loads.

$$\begin{array}{}\text{(6)}& q_{\mathrm{s}}={\displaystyle \frac{W^{*}{u}_{*}^{\mathrm{3}}{\mathit{\rho }}_{\mathrm{s}}}{(s-\mathrm{1})g}},\end{array}$$

where W^* is the dimensionless sediment transport rate, Eq. (7), $u_{*}=\sqrt{{\mathit{\tau }}_{\mathrm{0}}/\mathit{\rho }}$ the frictional velocity, τ₀ the bottom shear stress, ρ_s the density of the sediments, ρ the water density, $s={\mathit{\rho }}_{\mathrm{s}}/\mathit{\rho }$ the specific weight, and g the gravitational acceleration.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& W^{*}=\\ & \left\{\begin{array}{ll}\mathrm{11.2}{\left(\mathrm{1}-\frac{\mathrm{0.822}}{{\mathit{\varphi }}_{\mathrm{50}}}\right)}^{\mathrm{4.5}}& {\mathit{\varphi }}_{\mathrm{50}}>\mathrm{1.65}\\ \mathrm{0.0025}\exp \left[\mathrm{14.2}\left({\mathit{\varphi }}_{\mathrm{50}}-\mathrm{1}\right)-\mathrm{9.28}{\left({\mathit{\varphi }}_{\mathrm{50}}-\mathrm{1}\right)}^{\mathrm{2}}\right]& \mathrm{0.95}\le {\mathit{\varphi }}_{\mathrm{50}}\le \mathrm{1.65}\\ \mathrm{0.0025}{\mathit{\varphi }}_{\mathrm{50}}^{\mathrm{14.2}}& {\mathit{\varphi }}_{\mathrm{50}}<\mathrm{0.95}\end{array}\right.\end{array},\end{array}$$

where ϕ₅₀ is the normalized Shields stress, also known as transport stage, Eq. (8).

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{\varphi }}_{\mathrm{50}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{50}}^{*}}{{\mathit{\tau }}_{{\mathrm{r}}_{\mathrm{50}}}^{*}}}\\ & {\mathit{\tau }}_{\mathrm{50}}^{*}={\displaystyle \frac{u_{\mathrm{2}}^{*}}{(s-\mathrm{1})g{D}_{\mathrm{50}}}},\end{array}\end{array}$$

where ${\mathit{\tau }}_{\mathrm{50}}^{*}$ is the Shields stress, ${\mathit{\tau }}_{{\mathrm{r}}_{\mathrm{50}}}^{*}=\mathrm{0.0876}$ is a reference Shields stress value, and D₅₀ is the mean sediment size of the substrate.

Originally, Eq. (7) considered $W^{*}=\mathrm{0}$ for ϕ₅₀<0.95, but it was modified by Parker (1990) to have a positive transport rate in every flow condition. Parker and Klingeman (1982) did another modification for multiple sediment sizes, in which the dimensionless fractional transport rate $W_i^{*}$ for diameter D_i is calculated as a function of Eq. (9).

$$\begin{array}{}\text{(9)}& {\mathit{\varphi }}_i={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{50}}^{*}}{{\mathit{\tau }}_{{\mathrm{r}}_{\mathrm{50}}}^{*}}}{\left({\displaystyle \frac{D_{\mathrm{50}}}{D_i}}\right)}^{\mathrm{0.018}}\end{array}$$

The FLO-2D model allows the user to specify the sediment gradation coefficient, S_CC, Eq. (10), which is a metric of the sediment size distribution.

$$\begin{array}{}\text{(10)}& S_{\mathrm{CC}}={\displaystyle \frac{{\left(D_{\mathrm{30}}\right)}^{\mathrm{2}}}{D_{\mathrm{10}}{D}_{\mathrm{60}}}}\end{array}$$

Cabré et al. (2020a) reported that viscous debris flows observed in the Crucecita Alta fan show negligible erosion. Therefore, the mudflow model of FLO-2D is enough to represent the flow characteristics of viscous debris flows. On the other hand, Cabré et al. (2020a) indicate that inertial debris flows observed in the Crucecita Alta fan significantly modified its morphology, so erosion and deposition processes can not be neglected. Since FLO-2D can not use the mudflow and sediment transport models simultaneously, we used the following reasoning. For viscous debris flows, the first two terms of Eq. (3), i.e., yield and viscous stresses, dominate flow resistance. Conversely, for inertial debris flows, the third term of Eq. (3), i.e., turbulent and dispersive stresses, becomes more important for estimating S_f. Because of this behavior, inertial debris flows can be modeled using the same Manning approach for water flows but considering an appropriate Manning coefficient. This approach has been validated in previous studies, in which the Manning coefficient was estimated by back-calculation (Rickenmann et al., 2006; Naef et al., 2006). Typically used values for the Manning coefficient for modeling turbulent Newtonian debris flows are around 0.1 (Rickenmann, 1999).

3.2 Model configuration

The FLO-2D model built for this tributary-junction alluvial fan considers a 1400 m long river section and a 450 m long section of the alluvial fan, from the fan apex to the river junction (Fig. 3). To have an adequate flow development and achieve the balance of the sediment transport capacity before reaching the fan apex, the lidar topography of Crucecita Alta was extended 550 m upstream using the ALOS PALSAR DEM (digital elevation model) (12.5 m × 12.50 m resolution: © JAXA/METI ALOS PALSAR L1.0 2007). Secondly, because the ALOS PALSAR topography does not have the resolution to represent the creek topography, a synthetic channel was created in the extended zone to better represent the channel morphology. The synthetic channel dimensions were mapped based on the Google Earth Pro scenes. Thus, the model is 1 km long from the Crucecita Alta catchment (from the inflow element to the river junction), where 450 m corresponds to the lidar topography and 550 m to the synthetic channel (Fig. 3). No pre-event high-resolution topography is available for the Crucecita Alta fan. This is a common problem when studying the geomorphic consequences of low-frequency, high-magnitude debris flow events in arid areas (Zegers et al., 2017). To overcome this limitation, similarly to the approaches used by McMillan and Schoenbohm (2020), the post-event raster was modified to reproduce the pre-event topography based on available pre-event imagery and field evidence. The pre-event geomorphological configuration of the Crucecita Alta fan (Fig. 3) consists of an undisturbed fan surface without evidence of erosion (caused by the incision of a feeder channel or headward erosion). The modifications performed over the post-event lidar topography to obtain the pre-event geometry can be summarized into the following: (i) the incisions caused by the 25M flood event in the alluvial fan were removed from the DEM, and then the surface was smoothed; and (ii) facies F1 and F2 were subtracted from the DEM based on field thickness measurements of the deposits and assuming that cohesive debris flows have little to negligible bed erosion (Cabré et al., 2020a).

We also study the effects that a bypass channel, a typical debris flow mitigation work in alluvial fans, has on the fan–river connectivity. To this purpose, a 10 m wide and 10 m deep rectangular channel was inserted in the pre-event topography (Fig. 3).

The resolution of the numerical grid was set at 5 m × 5 m, as a finer resolution results in numerical instabilities due to the high flow velocities. In contrast, a coarser grid resolution results in a loss of terrain information. Manning number, n, was set equal to 0.1 s m${}^{-\mathrm{1}/\mathrm{3}}$ for the alluvial fan, which is within the range suggested by Rickenmann et al. (2006) (0.07–0.16 s m${}^{-\mathrm{1}/\mathrm{3}}$) for debris flows. Although the incisions measured in the field reach up to 5 m deep, erosion depth was limited to 4 m as greater depths affect the model stability. No sediment rating curve was set for the inlet because the 550 m synthetic channel allows the model to find its own sediment load equilibrium before reaching the fan apex.

3.3 Modeling approach and workflow

The flood sequence of the 25M event consisted of four different debris flow surges with different sediment loads. The differences in the rheology of the flows resulted in multiple geomorphological adjustments within the fan. Viscous debris flows mainly deposited sediments in the fan surface. In contrast, inertial debris flows strongly incised and eroded the fan. Consequently, our modeling approach divides the debris flow event between viscous and inertial debris flow surges and correlates them to the field evidence (Fig. 2). Therefore, if the facies (or any other sedimentological description) indicates that a viscous debris flow formed the deposit, the associated surge is modeled with the mudflow model. Conversely, if the facies suggests an inertial debris flow, the surge is computed using the water flow model with n=0.1 and including sediment transport. For this purpose, we developed a Python routine to concatenate the four surges into one model. Each surge is modeled in a sequence in which the results of the previous surge modify the topography for the next surge. The workflow steps are presented in Fig. 4.

Figure 4 Surge modeling workflow. The main challenge is to find a correlation between surges and field data that allows the user to discriminate between viscous and inertial debris flow surges. The developed routine is able to concatenate the surges by updating the model inputs for the next surge, based on the results of the previous one.

In the first step of this workflow, “subdivision into surges” (Fig. 4), we set the number of surges i=1 … N and the rheology model of each one, i.e., non-Newtonian (viscous debris flow) or Newtonian (inertial debris flow). For example, facies F1 and F2 have high fine sediment concentration, present evidence of a laminar flow, and show negligible erosion. Therefore, surges 1 and 2 are assigned as viscous debris flows. In contrast, the following surges 3 and 4 are assigned as inertial debris flows since F3, F4, and F5 show significant erosion and deposition and consist primarily of coarse sediment.

For the mudflow branch of the workflow (Fig. 4) the next step is the “calibration process and result screening”. This step finds the best set of rheological parameters, as explained in the next section. Multiple model runs may pass the screening. Therefore, in the step “model selection”, the user must visually inspect all the filtered runs and choose the best fit for the characterized facies. This step could not be automated because a critical analysis is needed. From the selected run, the deposit depth h_d is estimated for each grid element according to Eq. (11), where h_f is the resulting final flow depth reported by the numerical model, and $C_{\mathrm{V}}^{\mathrm{mean}}$ is the mean volumetric concentration of the surge. The sediment rating curve function C_V(Q,t) is presented in Zegers et al. (2020). Porosity p is set equal to 0.3 for facies F1 and F2 (Nicolleti and Sorriso-Valvo, 1991). Finally, h_d is added to the original topography at each cell by the Python routine, which updates the model for the next surge.

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& h_{\mathrm{d}}=h_{\mathrm{f}}{\displaystyle \frac{C_{\mathrm{V}}^{\mathrm{mean}}}{\mathrm{1}-p}}\\ & C_{\mathrm{V}}^{\mathrm{mean}}={\displaystyle \frac{\underset{t=\mathrm{0}}{\overset{T}{\int }}Q\left(t\right)C_{\mathrm{V}}(Q,t)\mathrm{d}t}{\underset{t=\mathrm{0}}{\overset{T}{\int }}Q\left(t\right)\mathrm{d}t}}\end{array}\end{array}$$

For the sediment transport model branch of the workflow (Fig. 4), on the other hand, the model is directly run since no parameters need to be calibrated. The sediment transport model only needs the parameters D₅₀ and S_CC (Eq. 10) of the substrate. In Crucecita Alta alluvial fan, D₅₀=7.8 mm and S_CC=1.79 (Cabré et al., 2020a). We assumed that the sediment transport equation and the grain size distribution of the substrate do not change between surges. After running the sediment transport model, our routine updates the model topography according to the computed erosion and deposition for each pixel; this updated topography is used for simulating the next flow surge.

3.4 Decision support system (rheology calibration)

Zegers et al. (2020) performed a sensitivity analysis for the FLO-2D numerical model and established the most sensitive parameters in the mudflow model. In their study, two alluvial fans located in the same valley were tested. They concluded that equifinality was present in the model, mainly related to the rheological parameters. Thus, to reduce over parameterization and restrict the parameter search, three parameters (α₁, α₂, and SD) have been left out of the calibration process. Consequently, the calibrated parameter set corresponds to the parameters β₁, β₂, $C_{\mathrm{V}}^{\max }$, and Vol_sediment (i.e., a reduction from seven to four parameters in the calibration process). Parameters α₁ and α₂ were fixed since they are the least sensitive parameters of the model (Zegers et al., 2020) and were set equal to 0.0075 poises and 0.152 dynes per cubic centimeter, respectively, using the values obtained by Zegers et al. (2017). Conversely, the SD parameter has a considerable influence on the results (Zegers et al., 2020). This parameter sets a minimum physically plausible flow depth, allowing the mixture to stop on the alluvial fan (i.e., the mudflow stops if the local flow depth is lower than SD). Based on estimated deposit height measurements by Cabré et al. (2020a), we tested different values for SD to reproduce the event by trial and error. SD was finally set equal to 1 m for the alluvial fan and 0.03 m for the river valley, in which the value in the valley is the default value in FLO-2D for water flows.

We used the algorithm of Zegers et al. (2020) as the basis to develop a novel decision support system (DSS) that allows us to run the model N times and automatically apply screenings to determine the model's best fit. Zegers et al. (2020) found that the calibration by flood area and deposited volume also constrained the flow velocity and flow height. Therefore, our DSS has incorporated filters for the affected area and the deposited volume. Similar to Mergili et al. (2018), we assigned the pixels within the observed deposit as observed positives (OPs) and pixels outside the deposit area as observed negatives (ONs). The pixels flooded by the simulation were assigned as predicted positives (PPs), whereas the non-flooded pixels were assigned as predicted negatives (PNs). Hence, the screenings filter the results based on thresholds θ₁ and θ₂ associated with true positives ($\mathrm{TP}=\mathrm{OP}\wedge \mathrm{PP}$) and false positives ($\mathrm{FP}=\mathrm{ON}\wedge \mathrm{PP}$) (Eq. 12). The aim of this thresholds is to find a manageable amount of filtered runs for visual inspection rather than finding specific threshold values. Moreover, these thresholds vary due to the distribution of the results, even for surges 1 and 2 in the same alluvial fan. In particular, the affected area thresholds were set equal to θ₁=0.7 and θ₂=0.3 for surge 1 and θ₁=0.9 and θ₂=0.1 for surge 2. The deposited volume threshold was set as ±60 % of the deposited volume estimated from field observations for both surges.

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\sum p{x}_i^{\mathrm{TP}}}{\sum p{x}_i^{\mathrm{OP}}}}>{\mathit{\theta }}_{\mathrm{1}}{\displaystyle \frac{\sum p{x}_i^{\mathrm{FP}}}{\sum p{x}_i^{\mathrm{OP}}}}<{\mathit{\theta }}_{\mathrm{2}}\end{array}$$

Figure 5 Viscous debris flow surges. Panels (a, b) correspond to the results for surge 1, whereas panels (c, d) correspond to the results for surge 2. Dashed polygons show the extent of mapped facies F1 (blue) and F2 (red). Panels (a, c) present the maximum flow depth, whereas panels (b, d) present the final flow depth at the end of each surge. Topographic changes for subsequent surges are updated based on the final flow depth of the previous surge.

For the selected parameters, we choose the range of each parameter based on the ranges defined by Zegers et al. (2020). To generate the combination of parameters for each simulation, we used the Latin hypercube sampling (LHS) method (Olsson and Sandberg, 2002). Since we simulated surges 1 and 2 with the mudflow model, both surges needed a calibration process to find their rheological parameters. We tested the model with sets of 50, 100, and 200 runs. Due to the constraint of calibrated parameters from seven to four, 100 model runs were sufficient to find reasonable results that fit the mapped deposits. For 50 runs, no result fits the screenings, whereas 200 runs show no significant improvements. The screenings applied, Eq. (12), returned four cases for each surge which were visually inspected to select the best fit. Considering the modified topography after the selected case for F1, we ran the same algorithm for the second surge. The minimum and maximum values for the calibrated parameters and the parameters adopted for surges 1 and 2 after the calibration process are presented in Table 1.

Table 1 Calibrated parameters and their range for surges 1 and 2. Parameters for both surges were expected to be different because of the equifinality and their different sedimentologic characteristics.

4 Results

In this section, we analyze if the proposed workflow can reproduce the debris flow event and the geomorphological adjustments that this alluvial fan experienced. Once the workflow capabilities are proven, we include a rectangular channel that connects the alluvial fan apex with the river. This channel is tested for the same chain of debris flow surges as in the 25M event and also for two inertial debris flow surges. Then, we study the effect of this channel on the lateral and longitudinal coupling status of the river junction.

Figure 6 Inertial debris flow surges. Panels (a, b) correspond to the results for surge 3, whereas panels (c, d) correspond to the results for surge 4. Dashed polygons delimit the telescopic-like deposit, i.e., the incision in the alluvial fan and the new lobe at the fan toe. Panels (a, c) present the maximum flow depth, while panels (b, d) present the topographic change for erosion (negative values) and deposition (positive values).

4.1 The 25M debris flow event reconstruction

The maximum and final flow depths of the viscous debris flow surges 1 and 2 are presented in the flood inundation maps of Fig. 5. In a first visual inspection of Fig. 5, the simulated flooded areas are in good agreement with the facies mapped. Integrating the inflow hydrograph and considering the variable C_V, for surge 1, the inflow volume is 15 450 m³ of water and 16 991 m³ of sediment resulting in $C_{\mathrm{V}}^{\mathrm{mean}}=\mathrm{0.52}$. According to the simulation, the volume of sediment deposited on the alluvial fan is 9482 m³ (Fig. 5a and b). In contrast, the volume estimated from field measurements is ca. 7000 m³ (F1). The integration of the inflow hydrograph for surge 2, on the other hand, results in an inflow volume of 28 218 m³ of water and 10 453 m³ of sediment, i.e., $C_{\mathrm{V}}^{\mathrm{mean}}=\mathrm{0.27}$. The lower value for $C_{\mathrm{V}}^{\mathrm{mean}}$ for surge 2 indicates a lower flow resistance than surge 1. According to the simulation, the deposited sediment volume is 4167 m³ (Fig. 5c and d), whereas the volume estimated from field measurements is 5000 m³ (F2).

Surge 1 does not reach the river due to its high $C_{\mathrm{V}}^{\mathrm{mean}}$ and the absence of a channel (Fig. 5b). In a fan without a feeder channel, viscous flows spread over the surface, increasing flow resistance and promoting flow detention and consequent sediment deposition (Jakob and Hungr, 2005). Surge 2 reaches the river (Fig. 5d) probably given to its more diluted nature (lower $C_{\mathrm{V}}^{\mathrm{mean}}$ than surge 1). Thus, a small portion of the sediment transported by surge 2 reaches the river, but it was insufficient to change the river geometry. Therefore, the alluvial fan acts as a sediment buffer, preventing the interaction with the main river. Consequently, the river's sediment discharge remains almost undisturbed even with the occurrence of surges 1 and 2.

The SD parameter controls the flow depth at the end of the simulation for the viscous debris flow surges and, consequently, their resulting deposit thicknesses (Fig. 5b and d). We expected this condition because the simulated deposited volume is very sensitive to SD when the alluvial fan operates as a buffer preventing sediments from reaching the river (Zegers et al., 2020). However, SD=1 m overestimates the sediment deposit thicknesses for surges 1 and 2. We also tested lower SD values of 0.5 and 0.7 m, associated with deposit thicknesses of 0.2 and 0.3 m according to Eq. (11), obtaining unsatisfactory results regarding the observed affected area. The problem with lower SD values is that the flow spreads over the surface, losing any correlation with the mapped facies. Therefore, we prioritize matching the flooded area instead of the deposit thickness for the model calibration as measuring the flooded area is more reliable than characterizing the deposit thickness distribution.

Our simulation results for the following inertial debris flow surges 3 and 4 are presented in Fig. 6. These surges have a higher amount of water than the previous viscous surges. Surge 3 has a total volume of water of 61 388 m³ and surge 4 of 143 568 m³. The maximum flow depth maps of Fig. 6a and c are evidence of the river avulsion triggered by surge 3, where the original river path and the newly formed river path are observed. The river is pushed to the opposite valley side due to surge 4 with flow depths over 2 m deep. The resulting topographic changes for surges 3 and 4 (Fig. 6b and d) are coherent with the primary erosion and deposition zones identified during fieldwork. Moreover, it is possible to locate the new lobe formed at the fan toe showing a telescopic-like pattern responsible for the river avulsion and the partial river blockage. On the other hand, the new scoured channel on the alluvial fan is wider than the channel observed in the field possibly due to the numerical model resolution. These results show that the fan–river interactions become important for surges 3 and 4 because the tributary sediment discharge experiences a complete coupling with the main river through alluvial fan trenching and progradation. Savi et al. (2020) reported similar behavior in their laboratory experiments.

We also analyzed the potential river blockage or dam formation. The obstruction ratio represents the space occupied by sediment deposits compared to the available space for the river to flow (Stancanelli and Musumeci, 2018). This ratio is useful to discriminate between three states (no blockage, partial blockage, full blockage). The obstruction ratio r_b is defined as the ratio between the obstruction lobe width and the flooded valley width, both along the orthogonal direction to the river flow (Stancanelli and Musumeci, 2018). The value r_b=1 means a complete river blockage. River blockage or dam formation is dominated by the momentum ratio R_M and the unevenness of the grain sizes S_C (Dang et al., 2009). The momentum ratio is calculated as the product of the flow rate ratio R_Q, the velocity ratio R_V, the bulk density ratio R_γ, and the confluent angle θ, while $S_{\mathrm{C}}=\sqrt{D_{\mathrm{75}}/D_{\mathrm{25}}}$. Dang et al. (2009) proposed a critical index C=R_MS_C for dam formation with different partial and complete blockage thresholds. Stancanelli and Musumeci (2018) highlighted that the thresholds have to be calibrated depending on the material adopted. For gravel deposits, they proposed a threshold value of C=9.

For surges 3 and 4, we obtained C indexes of 2.15 and 0.85, respectively. These values indicate that both surges do not have the potential to block the river. The C indexes are consistent with the obstruction ratios r_b of 0.75 (surge 3) and 0.86 (surge 4). These r_b values indicate a river's partial blockage due to the lateral input of sediment. Interestingly, the river experiences deposition in its downstream section for surge 3, i.e., an excess of sediment load (Fig. 6b). Conversely, the river experiences erosion in the downstream section for surge 4, i.e., a deficit of sediment load (Fig. 6d). This change indicates that a river obstruction r_b=0.75 was not able to act as a barrier. In contrast, an obstruction ratio r_b=0.86 was able to act as a barrier decoupling the river's sediment discharge in the river junction. For both inertial surges, the upper section of the river is a deposition zone because of the river's partial obstruction and, consequently, pounded water.

4.2 Morphological evolution of the tributary-junction alluvial fan

We characterized the tributary-junction alluvial fan evolution by analyzing the topographic change using four longitudinal sections (Fig. 7a). All longitudinal profiles are presented in the downstream direction. We established the A-A' profile of the alluvial fan according to the main incision observed in the field after the 25M event. In contrast, sections B-B', C-C', and D-D' correspond to the main river, which shifted during the event.

Figure 7 Morphological evolution of the tributary-junction alluvial fan. (a) Location of the longitudinal sections presented in (b) and (c). (b) Longitudinal profile A-A' of the alluvial fan topography after each surge. (c) Longitudinal profiles B-B', C-C', and D-D' of the river's topographic evolution.

River avulsion is well defined in profile A-A' (Fig. 7b), where the river channel shifts to the opposite valley side after surges 3 and 4 due to the formation of a new lobe. Profile A-A' shows a convex shape after surge 1, typical for viscous debris flows with sharp frontal boundaries. The following surge 2 stacks over surge 1, maintaining this convex shape and raising the topography. On the other hand, inertial debris flow surges 3 and 4 result in a characteristic concave shape for profile A-A' due to an alluvial fan incision. Analog laboratory experiments of Savi et al. (2020) reproduced similar sediment mobilization and geomorphic changes. The results of our simulation show less scour than the observed incision in the field possibly because of the fixed 4 m of maximum erosion. Another explanation for the lack of erosion is the difference between eroding the early event deposits (F1 and F2) and the pre-event old debris flow deposits. In our routine, the new deposits are accounted part of the topography instantly when our routine updates the topography between surges. However, the early event sediments should be eroded easier than old pre-event deposits.

The B-B' profile follows the river channel observed in the lidar topography. Profile C-C' corresponds to the main flow path of the river after the new lobe associated with surge 3. Profile D-D' corresponds to the main final flow path at the end of surge 4 and, therefore, at the model's end. The upper and lower sections are the same for these three profiles, whereas the lobe section follows the river avulsion path. Therefore, profiles B-B', C-C', and D-D' are presented for T₂, T₃, and T₄, which correspond to the status of the main channel after surge 2, 3, and 4, respectively (Fig. 7c). Initial status and the status of the river after surge 1 remain the same as T₂ and are not presented here. Profile B-B' illustrates how the sediment yielded from the tributary reaches the river forming the new lobe of the telescopic-like deposit with max deposited depths around 6 m for T₃ and 4 m for T₄. Surge 3 is responsible for the first channel avulsion, represented in the C-C' profile for T₃, and is evidence of the new lobe extent. Surge 4 increases the fan progradation, so profile C-C' exhibits deposition for T₄ in the lobe section. During surge 4, the river shifts again, and the D-D' profile follows the final river path where T₄ is evidence of the erosion and formation of the final river path.

4.3 Simulation of different scenarios considering mitigation works

In hazard assessment projects, numerical models are also used to design mitigation works. However, these models do not always consider the broad debris flow types that a creek can experience. Moreover, the same hydraulic works should be studied under the broad flow typology. As an example, we tested a rectangular channel (10 m width and 10 m depth) that connects the alluvial fan apex with the river junction under two different scenarios: (1) Scenario 1 consists of the same sequence of debris flow surges observed in the 25M event, in which two viscous debris flows are followed by two inertial debris flows. (2) Scenario 2 consists only of the two inertial debris flows. With both scenarios, we prove the importance of studying different debris flow type combinations for the same mitigation work.

4.4 Scenario 1

In Scenario 1 (Fig. 8a–d), the proposed channel confines viscous debris flows, which allows the flow to reach the river for surges 1 and 2. For surge 1 (Fig. 8a), deposit depth in the channel is up to 7 m and over the river junction up to 4 m, indicating a partial channel obstruction. Surge 2, on the other hand, spreads over the river junction but does not form a new lobe. Surges 3 and 4, shown in Fig. 8c and d, are deviated and forced to deposit due to the previous deposits. With surge 4, avulsion is present, inundating the southern portion of the fan. We first tested channels with smaller cross-sections where overflow also occurred, whereas a larger channel cross-section would be cost-inefficient and, therefore, not realistic.

Figure 8 Two future scenarios are modeled to understand the possible effects of a channel (mitigation works): first scenario (a–d) and second scenario (e, f).

This case shows that a channel is not sufficient to convey the debris flows safely. However, deposits observed in the river junction are smaller than deposits observed in the original case. Surprisingly, the presence of a channel does not directly mean an increase in lateral connectivity. In the 25M event, alluvial fan trenching creates a local increase in the sediment load, which is immediately deposited at the fan toe due to the change in the slope. From this test, we learned two important features when assessing this debris flow hazard. First, trimmed fan toes create a base level drop that increases the sediment load that reaches the river, enhancing river avulsion. Second, transport-limited catchments such as Crucecita Alta produce sediment discharges during extreme storm events, resulting in cost-inefficient works which are not realistic. To sum up, channels may reduce the impact of debris flows, but they do not entirely solve the problem. Therefore, targeted population education to avoid people settling in areas at risk must join these works.

4.5 Scenario 2

Scenario 2 (Fig. 8e and f), consisting of two inertial debris flow surges, also demonstrates that the presence of a channel avoids the telescopic-like deposit. The channel avoids an abrupt base level drop due to the trimmed fan toe and, therefore, avoids headward erosion and the local increase in sediment load. Consequently, no local sediment load increase takes place in the alluvial fan for surges 3 and 4. For this scenario, only a mild river avulsion without river obstruction occurs for surge 3 and surge 4.

5 Discussion

5.1 Dynamic response of the fan–river connectivity

The field “forensic” analysis of the debris flow event in the Crucecita Alta fan resulted in a clear differentiation of the transport mechanisms for each surge (Cabré et al., 2020a). Our workflow benefits from this differentiation to select a specific model to represent the flow of the water–sediment mixture in each surge. Our results show that the combination of these numerical models reproduces the main processes described by Cabré et al. (2020a). For the 25M event in Crucecita Alta alluvial fan, the coupling conditions changed within the event, but our workflow can reproduce this dynamic response. For the 25M event in the Huasco basin, Cabré et al. (2020a) reported 49 catchments with a characteristic response in which seven catchments are described in detail. This specific response of the catchments consists of viscous debris flow surges, followed by inertial debris flow surges. This catchment response in the Atacama Desert also occurred in 2017 and 2020 but has not been reported yet.

Characteristic stratigraphic and geomorphic features observed in the Huasco River basin can also be used to understand changes in the water–sediment ratio and the coupling degree between fans and the main river. We observed in the field that sediment sources influence the fine sediment vs. coarse sediment ratio present in the debris flows. On the one hand, viscous debris flows travel short distances because of their high flow resistance. Therefore, these surges come from sources close to the fan apex and arrive first. Conversely, inertial debris flows can travel longer distances and come from more distant sources (i.e., they arrive later at the alluvial fan). We expect this behavior to keep occurring in these arid zones; therefore, our workflow could be used to study future scenarios. Moreover, our workflow could assist in the volume estimations of water and sediment needed to modify the sediment cascade dynamics under viscous and inertial debris flows, which are different, as seen in this study.

Using our workflow on a past event in Crucecita Alta, we determined the following characteristics. The fan can completely buffer the viscous debris flow surge 1 ($C_{\mathrm{V}}^{\mathrm{mean}}=\mathrm{0.52}$ and a volume of sediment around 9500 m³). Conversely, for surge 2, a still viscous but more diluted debris flow, the fan partially buffers the sediment discharge ($C_{\mathrm{V}}^{\mathrm{mean}}=\mathrm{0.27}$ and a volume of sediment around 4200 m³). Moreover, for both later inertial debris flows, the sediment discharge of the tributary is totally coupled. Alluvial fan trenching for surge 3 increases the river's sediment load downstream, leading to total connectivity of the sediment cascade. Conversely, for surge 4, the greater volume of sediment acts as a barrier that decouples the sediment discharge of the river.

Our methodology is helpful to study hypothetical scenarios, such as a channel designed to mitigate affected areas but increasing the coupling status of the alluvial fan. The results of scenario 1 (Fig. 8a–d) show that viscous debris flows can reach the river since they remain confined to the channel. However, sediment deposition occurring during the flow reduces the channel transport capacity for the following surges. Thus, when the inertial debris flow surges occur, we observe that the channel can convey surge 3 but not surge 4. Compared to the original 25M event, smaller erosion areas result in lower sediment loads, reducing the size of the new lobe at the fan's toe. The results of scenario 2 (Fig. 8e and f), which consider only the inertial debris flows, show that the channel avoids the headward erosion observed in the 25M event. Consequently, the amount of sediment deposited at the new lobe at the fan toe is reduced considerably compared to the original case, even though we expected that the presence of the channel would increase the structural connectivity.

5.2 Model limitations

Even though our workflow successfully reproduces the main processes of interest, we observed two main limitations. First, we recognize that neglecting bank erosion may affect the sediment cascade dynamics. However, this process is rarely considered because most numerical models do not account for bank erosion. Second, the surface detention parameter (SD) directly impacts the deposited volume of sediment (Vol_sediment) in the calibration process because it controls the final flow depth. Therefore, estimating the SD parameter based on field data is advisable. If no data of the deposit depth are available to estimate this parameter, SD could be added to our novel decision support system. However, we calibrated a parameter set of four unknowns (β₁, β₂, C_V, Vol_sediment), whereas considering five unknowns increases the computational costs exponentially. In the absence of a proper estimation of SD, we recommend that the user chooses a value between 0.5 and 1 m for viscous debris flows in steep alluvial fans in the Atacama Desert instead of adding SD to the DSS. This will maintain computational costs reasonably. The range [0.5, 1] m is based on this study's findings and other previous studies in the Atacama Desert, such as Zegers et al. (2017, 2020). Despite these limitations, it is shown here that fan–river interaction studies can be performed with a commercially available software application that simplifies the physics around the flow of solid–liquid mixtures.

5.3 Insights for a new modeling approach

High-intensity rainstorm events in the Atacama Desert have increased their occurrence in the last century, and it is expected that this trend will continue due to climate change (Ortega et al., 2019). A higher frequency in high-intensity rain events in these transport-limited (i.e., sediment-unlimited) catchments (Aguilar et al., 2020) results in a higher frequency in the debris flow events. Following this reasoning, the necessity of properly reproducing the debris flow mechanics should be a priority. Most of the models focus on run-out distance and affected area, but none cope with the notorious changes in rheological macroscopic behaviors between surges. Based on the data for a specific event, our workflow studies the tributary-junction alluvial fan's response to a characteristic chain of flows with different rheological behaviors in the arid region of Chile (27–30^∘ S). This study highlights the effects of an alluvial fan on the lateral and longitudinal connectivity of the sediment cascade. Our workflow can be used to study future scenarios in which the catchment has a specific sedimentological response to rainfall events. Moreover, we show the decisive and dynamic role of an alluvial fan in the tributary–river junction (dis-)connectivity and, consequently, in the transference of the sediment down-system. To our knowledge, this important interaction is not yet considered in hazard assessment projects but has a considerably impact in the sediment budget analysis of a fluvial system. This work has presented new insights on how debris flow hazards should be simulated in a tributary-junction alluvial fan. In the future, we expect to propose a modeling protocol to guide modelers on how to study these important but neglected effects on catchments in northern Chile.

6 Conclusions

This study addresses how a tributary-junction alluvial fan results in local perturbations over the sediment cascade dynamics. We created a novel methodology based on the conceptual models of Mather et al. (2017), Savi et al. (2020), and Cabré et al. (2020a) to numerically model the connectivity development in an alluvial fan. The proposed methodology reproduces sequential mass flow surges that reach the Crucecita Alta fan and its morphological adjustments. Our approach copes with the challenge of coupling different simulation models by integrating cascading mass flow processes in one integrated modeling workflow.

In Crucecita Alta alluvial fan, we identified an extreme rainfall event that could be separated into four different surges with different flow characteristics. Our model reproduced the observed fan aggradation and the tributary's sediment discharge buffering for viscous debris flow surges 1 and 2 (Fig. 5). On the contrary, inertial debris flow surges 3 and 4 carve and prograde the alluvial fan (Fig. 6). In this study we showed that coupling conditions for the tributary-junction alluvial fan are dominated by the viscous/inertial debris flow nature due to their capacity to travel shorter/longer distances, respectively. Thus, the significance of the mechanical classification of debris flows is critical to reproduce how an alluvial fan controls the tributary–river junction connectivity and the effects of building a channel as mitigation work.

Our methodology provides new insights into how to improve hazard assessment projects from an engineering perspective in arid environments such as the Atacama Desert. Studies considering only the mudflow model, i.e., not taking into account erosion and deposition processes, underestimate affected areas due to sediment load changes and the avulsion of the river. On the other hand, studies considering only the sediment transport model in water flows do not consider the effects of non-Newtonian flows such as channel plugs due to highly viscous debris flows. The complex combination of different debris flow types makes it impossible to know beforehand how the fan and the mitigation works will behave, so the best approach is to study the solution under different scenarios by combining viscous and inertial debris flow surges. Our methodology is a good solution to represent debris flows and their effects on the tributary-junction alluvial fan connectivity, taking advantage of the capabilities of existing one-phase models. However, numerical models should improve how they solve erosion and deposition processes and changing rheologies to better represent sediment transport processes in tributary-junction alluvial fans.