The mudflow disaster at Villa Santa Lucía in Chilean Patagonia: understandings and insights derived from numerical simulation and post-event field surveys

Abstract. The evaluation of potential landslides in mountain areas is a very complex process. Currently, event understanding is scarce due to information limitations. Identifying the whole chain of events is not a straightforward task, and the impacts of mass-wasting processes depend on the conditions downstream of the origin. In this paper, we present an example that illustrates the complexities in the evaluation of the chain of events that may lead to a natural disaster. On 16 December 2017, a landslide occurred in the Yelcho mountain range (southern Chile). In that event, 7 million m3 of rocks and soil fell on the Yelcho glacier, depositing 2 million m3 on the glacier terminal, and the rest continued downstream, triggering a mudflow that hit Villa Santa Lucía in Chilean
Patagonia and killing 22 people. The complex event was anticipated in the
region by the National Geological and Mining Survey (Sernageomin in Spanish). However, the effects of the terrain characteristics along the
run-out area were more significant than anticipated. In this work, we
evaluate the conditions that enabled the mudflow that hit Villa Santa Lucía.
We used the information generated by Sernageomin's professionals after the
mudflow. We carried out geotechnical tests to characterize the soil. We
simulated the mudflow using two hydrodynamic programs (r.avaflow and Flo-2D)
that can handle the rheology of the water–soil mixture. Our results indicate that the soil is classified as volcanic pumices. This
type of soil can be susceptible to the collapse of the structure when
subjected to shearing (molding), flowing as a viscous liquid. From the
numerical modeling, we concluded that r.avaflow performs better than Flo-2D.
The mudflow was satisfactorily simulated using a water content in the
mixture ranging from 30 % to 40 %, which would have required a source of about 3 million m3 of water. Coupling the simulations and the soil
tests that we performed, we estimated that in the area scoured by the
mudflow, there were probably around 2 800 000 m3 of water within the soil. Therefore, the conditions of the valley were crucial to enhancing the impacts of the landslide. This result is relevant because it highlights the importance of evaluating the complete chain of events to map hazards. We suggest that in future hazard mapping, geotechnical studies in combination with hydrodynamic simulation should be included, in particular when human lives are at risk.



1
Introduction 50 Landslides processes are particularly dangerous in areas close to human settlements. They can affect nearby villages, directly destroying houses and taking human lives (Gariano and Guzzetti, 2016) or indirectly affecting the connectivity of remote areas (Winter et al., 2016). The impacts of landslides are a function of the size of the event but also of the conditions downstream. For example, glaciar lakes susceptible to overflow, as well as unstable valleys that, given the right soil matrix and water content, can mobilize and produce mudflows (Carey et al., 2011;Haeberli et al., 2013). 55 Areas where glaciers are receding worsen this situation because they expose unstable hillslopes that can collapse as well as potentially create glacier lakes. Currently, baseline information availability still critical in austral zones of South America, especially in Northern Patagonia, with a low population density that has not encouraged rigorous evaluation. Moreover, in recent years landslides events have increased due to anthropic and climatic effects (Aldunce and González, 2009). Parallelly, Northern Patagonia shows an increase in the population (INE, 2018) increasing the 60 risk. Therefore, a better understanding of the mudflows chain of events triggered by landslides is urgent.
In our contribution, we will evaluate the generation of a cascade of events associated with the Villa Santa Lucia mudflow in Northern Patagonia. In this study, we will evaluate the mechanisms that enable a landslide of 7x10 6 m 3 to

Geological setting
The study area consists mainly of 9 geological units ( Figure 2). The dominant geological unit corresponds to an intrusive rock (Cretaceous and Miocene age) composed of tonalites, granodiorites, granites, diorites, and others.
Moreover, the basement material is an old metamorphic rock. This unit is composed of micaceous shales, amphibolite.

Event Background
The event in Villa Santa Lucía was triggered during a rainfall event front that passed over the area in a year where the entire Los Lagos Region. The event generated a significant record of rainfall surplus. Therefore, the hydrometeorological condition could constitute a condition of soil saturation that is favorable to natural hazards (Nguyen et al., 2018). 105 According to the information provided by the "General Water Directorate" (DGA in Spanish), the annual average rainfall in this area is 3,420 mm. At the time of the event, the total annual rainfall was 3,650 mm, and in the 30 hours prior to the event the rainfall reached 124.8 mm, with a maximum intensity of 10.6 mm/h at 16:00 hours on December 15, 2017. This hydrometeorological event exceeds 99% of the historical precipitation events in the area (CECS, 2017).
The precipitation events in Villa Santa Lucía occurred after two weeks with maximum daily temperatures that exceeded 22°C on at least nine days between December 1 and 15, with December 2 and 5 being the hottest days with temperatures that exceeded 27°C. During December, before the event, the air temperature reached an average of 14.9°C (CECs, 2017). The Center for Scientific Studies (CECs in Spanish) determined through the data of a radiosonde located in Puerto 73.09ºW) that the isothermal level on December 15 was 2,771 meters above sea level (m a.s.l.) for the Villa Santa Lucía coordinates. This implies that in the study area, there were only 115 liquid precipitations in the days before the event, even in the glaciers of the Yelcho range, considering that the maximum heights of the peaks do not exceed 1,800 m a.s.l. approximately.
Previous studies shows that the surroundings of the Villa Santa Lucia presents a high and moderate risk of being affected by mass wasting processes, highlighting possible debris flows at the bottom of the valleys and active channels (Sernageomin, 2008). Areas that are not classified in high or moderate risk have low or no risk, which should be 120 assessed at greater detail. High flood hazard would affect only the flood plain adjacent to Villa Santa Lucia, although some degree of risk in the town itself must be evaluated on a more detailed scale. This shows that the area of Villa Santa Lucía and its surroundings were and are prone to mass wasting phenomena. Although the town itself was not in risk, it was mentioned that more detailed studies are needed in the area.

Description of the event 125
On December 16, 2017, a 7 million cubic meter volcanic rock slide detached ( Figure 3, point 1), falling next to the Yelcho glacier toe that sits on an intrusive formation that has a drop of about 80 meters. The glacier shows a retreat of at least 250 meters during the last decade. Sernageomin (2018c) indicated that there might have been a small lake and water available on top of the intrusive formation, but there is no conclusive information. Satellite images before the events do not indicate a lake above the intrusive unit. Also, there were indications that ice contributed to the 130 mudflow. Remote sensing data do not show that the landslide felt on top of the glacier unless detritus had covered the glacier, although we could see ice on top of the flank that collapsed. However, we do not have data to prove or disregard the covered ice on the intrusive and how much water it contributed to the mix. Two million cubic meters of the landslide material stayed above the intrusive, and five million cubic meters of material continued downstream, sliding on the intrusive unit at a slope of more than 70 degrees. At this point, the mudflow reaches the Burrito River 135 with high topographic differences. The flow continued along with the drainage network at high speed ( Figure 1 and  Additionally, the soil was almost saturated, which added water to the mix, transforming the landslide flow into a mudflow. Moreover, dense forest was present, which added a significant amount of biomass to the mix. Then the mudflow reached an area with low slopes at a distance of 8.6 km from its origin through the Burritos River to the east. In its trajectory, the mudflow crossed Route 7 (Carretera Austral) in a stretch estimated to be 2 km and filling an old wetland. In that sector, the flow was channeled in a canyon oriented north-south toward Villa Santa Lucia in a section 145 of 1.5 km (Figures 1 and 3  3. Non-channeled flow: This section has a length of 2.3 km in a west to east direction until the Burritos River 160 changes its orientation to a north-south direction. The elevation goes from 600 to 380 m a.s.l. The slopes decrease, slowing down the flow which deposited sediments. The non-channeled flow width reaches 1.4 km. In this sector, the flow goes over the road (Carretera Austral) to the east in 1.3 km, affecting an old wetland (Figure 3 and 4c-d).

4.
North-south channeled flow: The flow descends through the enclosed channel at the Burritos River canyon from 380 to 250 m a.s.l. on a 2 km long section ( Figure 3 and 4e-f). 165

5.
The alluvial fan on Villa Santa Lucía: Due to the loss in confinement, the flow expands with a radial distribution in a southeast-southwest direction of 600 to 1,000 m (

Fieldwork, Topography and Soil Characteristics
A fieldwork campaign was carried out in January 2019 to map in high resolution the area using an unmanned aerial vehicle (UAV). We used aerial photogrammetry to produce a high-resolution DEM for some parts of the study area.
Aerial photogrammetry creates 3D models from 2D images, obtaining geometric characteristics of the objects they 185 represent. We used a UAV Inspire II with a Zenmuse X4 camera to capture aerial photos ( Figure 5). We did not have differential GPS at the moment of this survey to take control points to correct the DEM generated. Therefore, we just used this DEM to corroborate that the DEMs with lower resolution and freely available were able to capture specific features such as canyons or small changes in slopes that could affect the hydrodynamic simulation and the path of the flood. We also generated a high-resolution mosaic to observe details of the mudflow deposition to improve the 190 delimitation of affected area. We used the software Agisoft PhotoScan Professional 1.4 desktop to process the images.
The software performs the restitution of images by spatial coincidence among the elements represented in each image. It was not possible to carry out the extraction of unaltered geotechnical samples from the upper part, mainly due to the 5 kilometers of steep terrain, including a section that needed to be climbed to transport the material by foot.
However, we extracted soil samples in the middle-low part of the first section of the flow ( Figure 6). The geotechnical soil sample represents the material added to the mudflow during the event. Therefore, the soil properties become 205 relevant to reproduce the mudflow flow of Villa Santa Lucía.
The tests performed to determine the geotechnical properties of the primary materials were: direct shear test, unconfined compression, density, and soil classification following the American Society of Testing Materials ASTM D3080 (ASTM, 2017a). The parameters obtained in each laboratory test will serve as constraints for the numerical modeling of the mudflow of Villa Santa Lucía. For the soil moisture, a wet portion of the soil was extracted and allowed to dry in an oven at 60°C for two days.
Following, we estimated mass soil water content (by weight) using the relation 1.96 gwater/gsolid 215 Then, a fraction of the soil sample was covered with paraffin to seal all pores and thus determine its density through the submerged weight difference. The wet soil density was 1.24 grams per cm 3, and the dry soil's apparent density was 0.419 grams per cm 3 -the specific weight Gs=2.65 gr. Therefore, the void ratio was 5.324. Using these values, it turned out that the soil porosity was 84.18%, and the water saturation of the soil sample 97.4%; therefore, 81.95% of the total unaltered soil volume corresponded to water. 220 Then a direct shear test was performed in three probes. Each probe was loaded with 2, 4, and 8 kg, respectively, allowing them to consolidate for 24 hours. We recorded the deformations versus time in the first sample to determine the speed for the direct shear test using the square root scale method (Table 1). In order to carry out the classification of the soil, we used two tests. First, we determined the soil granulometry using the American Society of Testing Materials ASTM D2487 (ASTM, 2017a), and then we determined soil consistency limits (liquid and plastic) using the American Society of Testing Materials ASTM D4318-17e1 (ASTM, 2017b). For the granulometry, Table 2 shows what percentages pass through the different sieve sizes. We omitted the larger sieves 230 since 100% of the material passed through them. Finally, the liquid and plastic limits of the soil are 50% and 27%, respectively. According to the USCS classification, the soil corresponds to a CH, an inorganic clay of high plasticity, whereas for the AASHTO classification, it is 235 considered to be clay.

Numerical Modeling
The modeling of the mudflow was carried out using two pieces of software: r.avaflow and Flo-2D, which we describe below. According to the antecedents of the event, the landslide fell on the intrusive formation. A fraction of the material stayed deposited between the glacier and the intrusive formation. The difference (about five million cubic 240 meters of sediment) continued downstream, initiating the scouring process and mudflow that led to the Santa Lucia event. Therefore, the domain of the models starts at the base of the intrusive formation with an initial volume of 5,000,000 cubic meters.

R.Avaflow
R.avaflow is an open-source model that offers a practical, innovative and unified solution to simulate granular and 245 debris flows produced in high mountains around the world. The model can handle rapid mass flows, including avalanches and two-phase flows (Mergili, M., Pudasaini, 2019). R.avaflow calculates the propagation of mass flows from one or more given release areas on a defined baseline topography until (I) all material has been stopped and deposited; (II) all the material has left the area of interest; or (III) the maximum simulation time has been reached.
R.avaflow is developed in two formats for its environment and operation, r.avaflow expert and r.avaflow professional. 250 The latter is an autonomous version with reduced functionalities. It operates through a graphical user interface (GUI).
In this study, we used the professional version (Mergili et al., 2017). This model represents a complete open-source computational framework based on a geographic information system (GIS) that offers a two-phase flow model. The  (Julien and Lan, 2007;O'Brien et al., 1993;O'Brien and Julien, 1988), which describes the dynamic viscosity and the shear stress of the mixture as an exponential function of the sediment content 270 (FLO-2D, 2018).

Model calibration
For the calibration of the models we used a trial and error aproach seeking too match the following three pieces of data availablse from Sernagomin: (1) Inundated area, which was mapped after the event using aerial imagery (SAF,

2017); (2) Flow heights estimated by Sernageomin in Villa Santa Lucía and at the beginning of the canyon; and (3) 275
Flow velocities in the canyon curve and at the beginning of Villa Santa Lucia (Figure 7). See Table 3 for the values.

R.Avaflow
To simplify the calibration, we divided the process into two. First, we set the sediment concentration by volume of the mudflow in 50% and change the entrainment coefficient, basal friction angle, ambient drag coefficient and fluid friction coefficient. For the description of the parameters in r.avaflow see Mergili et al. (2017). Table 4 shows the first 290 set of parameters used.

295
The best set of parameters is 10 -5.75, 2, 0.022 and. 0.0005 for the entrainment coefficient, basal friction angle, ambient drag coefficient and fluid friction coefficient, respectively, and the corresponding map result is in Figure 8.  (Sernageomin, 2018). However, the mixture is still largely fluidized, and it does not follow the edges of the inundation. For this reason, we performed simulations with the same input parameters but a varying percentage of initial water. We varied the percentage of water between 20% and 70%.
The error for the heights, speeds calculated in each model are in Figure 9. Therefore, we propose that a mudflow with a 30% water volume could reproduce best the VSL event ( Figure 10). 305

Flo2d
The main results from Flo-2D for the mudflow flow in Villa Santa Lucía are presented below. For the laminar 315 resistance k, we used the default number 3000. The influence of the value of does not affect simulations significantly compared to other parameters related to flow resistance (Hsu et al., 2010). The specific gravity of sediments Gs is equal to 2.65 g cm -3 . For the sediment concentration by volume, we used values between 40% to 50%, which correspond to mudflow flows. The error for the heights, speeds calculated in each model are in Figure 11. Figure 12 shows the extension of the best model using Flo2D.

Geotechnical sample
Our geotechnical data shows that soil correspond to clay with high plasticity (CH). The soil has an internal friction angle of 23.8 ° and a cohesion value of 22.5 kPa, characteristic of a consolidated clay. However, the soil extracted had various types of soils and minerals in its composition, which made this a rare soil not described or classified by USCS.
For the granulometry, 72% of the soil passed through sieve No. 200. Mainly sand and volcanic soils did not pass sieve 330 No 200. Similar results were founded by Gonzalez-Pulgar (2012) in volcanic soils that have an internal soil friction angle of 25° and a cohesion value of 2.9 kPa. Moreover, our results were consistent with their high soil water content, greater than 150%, and a very low dry density of solids with values between 0.4 and 0.7 kg cm -3 . The values obtained by (Gonzalez-Pulgar, 2012) are comparable to the results of this study.
The results of the dry density and natural moisture determination tests help to ascertain the natural state of the unaltered 335 soil sample obtained in the mudflow zone of Villa Santa Lucía. Additionally, these tests help to calculate intrinsic soil parameters, such as the relative density of solid particles (2.65 g cm -3 ), the void ratio, the degree of saturation and porosity. The humidity estimated under ASTM D2216-19 (ASTM, 2019) showed a value of 195.85% and a dry soil density of 48.1 kPa. The soil saturation was 97.4%. Therefore, 81.58% of the soil volume was water.
The previous result is relevant because the liquid and plastic limits are 50 and 27%, respectively. These parameters exaggerating the runout. The sensitivity analysis showed that r.avaflow is also sensitive to changes in the basal friction 355 angle. These strongly condition the rheology of the flow, determining the height and speeds. The results obtained in the modeling in the FLO-2D software were diverse. We could not achieve the same level of r.avaflow's results. The inundation area was the most sensitive result in the concentration of sediments by volume and the rheology model ( Figure 11). The best combination of parameters was the Glenwood 1 rheology model and a 40 % concentration sediment, which has an error of 32%. We concluded that we could reproduce the mudflow satisfactorily using 360 r.avaflow with water content in the mixture ranging from 30 to 40 % (Figure 11).
For the calibration, we calculated the standard deviation from the modeling results and the parameters from Table 3.
The parameter combination that results in less standard deviation is considered the best parameterization for the particular software used. Our procedure, as described in Saltelli and Annoni (2010), corresponds to a one-factor-at-atime type of calibration where we changed one parameter at a time manually, trying to match Table 3. The limitations 365 of this procedure, according to Saltelli and Annoni (2010), are: Its efficiency is poor. The method does not capture the interaction of the factors because that would require the movement of more than one element at a time. Therefore, the approach assumed that all the variables are independent, and the processes are linear, which is not usually the case.
Another limitation of our calibration is the potential presence of equifinality. We did not test if different combinations of parameters would provide similar performance. We selected one best combination of parameters without checking 370 the degrees of freedom that these variables may have had to replicate the observations (Beven, 1996).

Source of water
From the observation and information provided by Sernageomin, our soil samples analyses and numerical simulation, we found that the primary source of water was in the soil. Our analyses determined that in the area scoured, the soil has a porosity of 84.18%. Similar values were reported by Cuevas et al. (2013) in volcanic soils in the south of Chile. 375 From our observations, we also notice the ground is saturated all the time. We did fieldwork after a couple of weeks of no rain, and the whole area around the river was saturated; indeed, our soil sample had a saturation of 97.4%. Therefore, using the results from r.avaflow, we estimated that a volume of 3.402.100 m 3 was scoured. Consequently, there was the potential of adding roughly 2.8 million m 3 of water and 600,000 m 3 of soils, which evolved to a mudflow due to the water added to the event. Additionally, extra water was added, not quantified in this study, from the ice Our results highlight the importance of considering the potential chain of events in risk analysis and hazard mapping.
In the practice of predicting a mudflow hazard, we seek areas where evident water sources are available such as glacial lakes. However, in this work, we showed that it is possible to have the catastrophic event that hit Villa Santa Lucia only with the water within the soil in the valleys downstream of the landslide. Additionally, the soil has a particular 385 volcanic origin, and it is highly fluidizable under low lateral effort. This result is relevant for the Patagonia region because, in order to update the hazard maps, we need to consider not only the primary events but also the potential reactions from the downstream elements of the chain of events, which as in this case, may not be evident prior to the event.
What actually triggered the landslide event is not clear and not part of this research. We focused on what conditions 390 in the valley enabled the landslide to evolve into a mudflow and traveled all the way to Villa Santa Lucia. We have no information to guess how much water was available at the glacier terminal so we wanted to understand if this event was possible without a lake or large water reservoir at the glacier. Our results show that in fact the water available around the Burritos River was sufficient to transform the avalanche into a mudflow. Our numerical models and geotechnical work show that the mudflow event was possible due to pre-existing water in the saturated soil. 395 Geotechnical evidence demonstrates a high plasticity soil (CH) with a low internal friction angle associated with volcanic soils of the Holocene age. Hence, weak soils are prone to mass wasting events due to their geomechanical properties.

6
Conclusions 400 Given the complexity and the potential increase in the future of extreme event occurrence that can trigger landslides, we suggest that hazard studies should consider the structural conditions present in the area of influence of the mudflows. Soil characteristics ought to be included because they may be a crucial factor amplifying the impacts of local events triggered by hydroclimatic events. Such was the case of Villa Santa Lucia, where the water necessary to fluidize the mudflow mixture was in the soil in a volatile system that was easy to mobilize. 405 The combination of geotechnical tests and r.avaflow, which is a freely available computational software, to model mudflow are useful tools to characterize and reproduce mass wasting events, such as the mudflow occurred on December 2017 in Villa Santa Lucia in Chile. Our results present the possibility of open-source software implementation to represent mudflow events in the Patagonian Andes with a good performance. These types of studies will allow for the integration of better methodologies to enhance risk scenarios related to mudflow events in active 410 subduction zones like the Andes.