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 on 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 5 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, 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 10 impacts of these flows and the effects on the tributary-river junction connectivity.


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 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 60 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 65 (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 on 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 70 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 tributaryjunction alluvial fan connectivity changes during a debris flow event? To answer this question we used the FLO-2D commercial 75 software due to its reasonable computational cost and its broad applicability for hazard assessment projects.
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 80 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 this 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 85 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 flows 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 90 behaviours when artificially changing the fan river connectivity by means of a channel under different configurations of the debris flows surges.
Tributary-junction alluvial fans in the Huasco Valley, located between 1500 and 3000 m.a.s.l., supply sediments to the main river. They are so abundant that they can reach densities of 1.9 fans per km. During the 25M event, rainfall gauges in the valley 95 recorded precipitation ranging from 20 to 76 mm in three days, with a maximum intensity of 16 mm h −1 . 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 100 event, where a sequence of surges impacted the fans repeatedly (Cabré et al., 2020a). Mitigation and reconstruction works were not performed immediately in the Crucecita Alta fan, allowing a complete identification of the debris flow surges and geomorphological evidence. Therefore, the flood sequence of the 25M event in this fan was reported by Cabré et al. (2020a).

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. 1.a).

105
Its catchment has an area of ∼13 km 2 , 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 0.14 km 2 area; 8.6% mean slope; and hosts few houses in the northern area ( Fig. 1.b). The northern 110 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. Previous to the event, riverbank erosion trimmed the alluvial fan toe (Cabré et al., 2020a), disrupting the alluvial fan gradient ( Fig. 1.b). 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).

115
A post-event topography of a 51 km long segment of El Carmen River valley was acquired with 1×1 m 2 horizontal resolution between Feb-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 1.a 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.

120
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 previous to the junction were 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 3 s −1 during the fourth surge (Figure 2.b).

125
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 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 −1 during the last 6 − 10 M yr 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 15% up to 130 40% of the total sediment volume eroded, on average, every one hundred 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. 2.a). The different facies were named from F 1 to F 5, where F 1 was the first deposit during the storm, and F 5 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, 135 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), F 1 and F 2 correspond to debris flow deposits associated with matrix-supported flows with high fine sediment concentrations. F 1 is the thickest deposit (above 1 m) and reaching near 7000 m 3 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. F 2 is a wider but lower thickness deposit (10 -30 cm) that overlays F 1 and covers all the fan's length, from the 140 apex to the distal toe, with near 5000 m 3 of deposited sediment volume (Aguilar et al., 2020). F 1 and F 2 deposits overlain pre-event sediments, indicating that the associated flows had none or negligible erosion capacity (Cabré et al., 2020a). 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 where erosion and sedimentation processes become important (Church and Jakob, 2020). The debris flood surges generated the F 3, F 4, F 5 deposits, and were responsible for 145 the alluvial fan entrenchment and consequent generation of the feeder channel (Figure 2.b). F 3 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 F 3, 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 F 4 and F 5 facies further downstream. Feeder channel depth ranges from 100 to 350 cm in the mid zone and 150 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) 155 (F 1 to F 5) are fundamental evidence to determine the sequence of flows. Cabré et al. (2020a) interprets the significance of the mechanical classification of debris flows and classifies F 1 and F 2 as non-Newtonian flows. In contrast, F 4 and F 5 are classified as Newtonian flows, which is a similar differentiation, to some extent, to the viscous and inertial debris flows classification of Takahashi (2014). F 3 is a transitional flow between the observed non-Newtonian and Newtonian flows.

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, Eq. (1) and Eq. (2), known as Saint-Venant equations, for water and debris flows.
where h is the flow height, V x and V y are flow velocities in the x and y directions and S 0x and S 0y are the channel slopes in each transverse direction (x, y). Sf x,y denotes the frictional slope, accounting for flow resistance. In the case of water flows, Sf x,y is calculated using the Manning equation. In the case of mudflows, Sf x,y is calculated using the so-called quadratic rheology model, Eq. (3), which combines different terms accounting for yield/coulomb, viscous, and turbulent/dispersive 170 stresses (O'Brien and García, 2009).
In Eq. (3) where α 1,2 and β 1,2 are empirical coefficients that have to be calibrated. n is the conventional Manning coefficient, and b = 0.0538 and m = 6.0896.

180
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 185 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 timestep (O'Brien and García, 2009). However, no erosion/deposition processes can be included when using the mudflow module.
The sediment transport model accounts for the local deficit/excess of transported sediment, modifying the bed channel height 190 using the well-known Exner equation, Eq. (5).
In Eq. (5), h 0 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 eleven different expressions based on unique river conditions (O'Brien and García, 2009). For this study, we used the Parker-Klingeman-Mclean equation , Eq. (6). This equation is suitable for gravel and sandy bed material 195 and simulates sediment transport mechanisms from bedload to mixed suspended sediment loads.
where W * is the dimensionless sediment transport rate, Eq. (7), u * = τ 0 /ρ the frictional velocity, τ 0 the bottom shear stress, ρ s the density of the sediments, ρ the water density, s = ρ s /ρ the specific weight and g the gravitational acceleration. where ϕ 50 is the normalized Shields stress, also known as transport stage, Eq. (8): where τ * 50 is the Shield Stress, τ * r50 = 0.0876 is a reference Shields stress value, and D 50 the mean sediment size of the substrate.
Originally, Eq. (7) considered W * = 0 for ϕ 50 < 0.95, but it was modified by Parker (1990) to have a positive transport 205 rate in every flow condition.  did another modification for multiple sediment sizes, where the dimensionless fractional transport rate W * i for diameter Di is calculated as a function of Eq. (9).
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.   (Rickenmann, 1999).

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 225 the balance of the sediment transport capacity before reaching the fan apex, the LiDAR topography of Crucecita Alta was extended 550 m upstream using ALOS PALSAR DEM (12,5 x 12,50 m 2 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 represents 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 230 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 pre-event available imagery and 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 ten meter wide and ten meter depth rectangular channel was inserted in the pre-event topography (Fig. 3).
The resolution of the numerical grid was set at 5x5 m 2 , as a finer resolution results in numerical instabilities due to the high 245 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 −1/3 for the alluvial fan, which is within the range suggested by Rickenmann et al. (2006)   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 flows 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 where the results of the previous surge modify the topography for the next surge. The workflow steps are presented in Fig. 4.

260
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 F 1 and F 2 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 F 3, F 4, and F 5 show significant erosion and deposition and consist primarily of coarse sediment.

265
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 is presented in Zegers et al. (2020). Porosity p is set equal to 0.3 for facies F 1 and F 2 (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.
For the sediment transport model branch of the workflow (Fig. 4), on the other hand, the model is directly run since no 275 parameters need to be calibrated. The sediment transport model only needs the parameters D 50 and S CC (Eq. (10)) of the substrate. In Crucecita Alta alluvial fan, D 50 =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. The pixels flooded by the simulation were assigned as predicted positive (PP), whereas the not flooded pixels as predicted negatives (PN). Hence, the screenings filter the results based on thresholds θ 1 and θ 2 associated with true positives (TP= 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 θ 1 =0.7 and θ 2 =0.3 for surge 1 and θ 1 =0.9 and θ 2 =0.1 for surge 2. The deposited volume threshold was set as ±60% of the deposited volume estimated from field observations for both surges. 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 310 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.

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 Surge 1 does not reach the river due to its high C mean V and the absence of a channel (Fig. 5.b). In a fan without a feeder channel, viscous flows spread over the surface, increasing flow resistance, promoting flow detention, and consequent sediment deposition (Jakob and Hungr, 2005). Surge 2 reaches the river (Fig. 5.d), probably given to its more diluted nature (lower C mean V than surge 1). Thus, a small portion of the sediment transported by surge 2 reaches the river, but it was insufficient to 335 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 though 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. 5.b and 5.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.,  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 360 (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 C = D 75 /D 25 . Dang et al. (2009) proposed a critical index C = R M S C for dam formation with different partial and complete blockage thresholds. Stancanelli and Musumeci (2018) highlighted that the thresholds have to be calibrated 365 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. 6.b). Conversely, the river experiences 370 erosion in the downstream section for surge 4, i.e., a deficit of sediment load (Fig. 6.d). 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.

Morphological evolution of the tributary-junction alluvial fan 375
We characterized the tributary-junction alluvial fan evolution by analyzing the topographic change using four longitudinal sections ( Fig. 7.a). 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.
River avulsion is well defined in profile A-A' (Fig. 7.b), where the river channel shifts to the opposite valley side after  lobe section follows the river avulsion path. Therefore, profiles B-B', C-C', and D-D' are presented for T 2 , T 3 , and T 4 , which correspond to the status of the main channel after surge 2, 3, and 4, respectively (Fig. 7.c). Initial status and the status of the river after surge 1 remain the same as T 2 and are not presented here. Profile B-B' illustrates how the sediment yielded from the 395 tributary reaches the river forming the new lobe of the telescopic-like deposit with max deposited depths around 6 m for T 3 and 4 m for T 4 . Surge 3 is responsible for the first channel avulsion, represented in the C-C' profile for T 3 , and evidences the new lobe extent. Surge 4 increases the fan progradation, so profile C-C' exhibits deposition for T 4 in the lobe section. During surge 4, the river shifts again, and the D-D' profile follows the final river path where T 4 evidence the erosion and formation of the final river path.

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 (10m width and 10m 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 405 flow surges observed in 25M event, where two viscous debris flows are followed by 2 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 types combinations for the same mitigation work.

Scenario 1
In Scenario 1 (Fig. 8.a to 8.d), the proposed channel confines viscous debris flows, which allows the flow to reach the river for 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.

415
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 of 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 420 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 under risk must join these works.

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 435 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 where 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.

440
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 on the field that sediment sources influence 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 445 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 volumes 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 mean V = 0.52 and a volume of sediment around 9500 m 3 ). Conversely, for surge 2, 450 a still viscous but more diluted debris flow, the fan partially buffers the sediment discharge (C mean V = 0.27 and a volume of sediment around 4200 m 3 ). 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.

455
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. 8. 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 on the deposited volume of sediment (V ol 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 is available to estimate 470 this parameter, SD could be added to our novel Decision Support System. However, we calibrated a parameter set of 4 unknowns (β 1 , β 2 , C V , V ol sediment ), whereas considering 5 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 475 Zegers et al. (2017Zegers et al. ( , 2020. Despite these limitations, it is shown here that fan-river interaction studies can be performed with a commercially available software that simplifies the physics around the flow of solid-liquid mixtures.

Insights for a new modelling approach
High-intensity rainstorm events in the Atacama Desert have increased their occurrence in the last century, and it is expected that this trend continues due to climate change (Ortega et al., 2019). A higher frequency in high-intensity rain events in these 480 transport-limited (i.e., sediment-unlimited) catchments (Aguilar et al., 2020) result in a higher frequency in the debris flows events. Following this reasoning, the necessity of properly reproducing the debris flows mechanics should be a priority. Most of the models focus on run-out distance and affected area, but none cope with the notorious changes of 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°S -30°S). This 485 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 where 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 impacts considerably in the sediment budget analysis of a fluvial system. This work has presented new 490 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.

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); Cabré et al. (2020a) 495 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 500 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.

505
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/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 510 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/deposition processes and changing rheologies to better represent sediment transport processes in tributary-junction alluvial fans.