Back-Analysis of the Abbadia San Salvatore (Mt. Amiata, Italy) Debris Flow of 27–28 July 2019: An Integrated Multidisciplinary Approach to a Challenging Case Study

: On 27–28 July 2019, in a catchment of the Mt. Amiata area (Italy), an extreme rainfall induced a debris ﬂow, which caused a channeled erosive process just upstream of the Abbadia San Salvatore village, the obstruction of a culvert at the entrance to the urban area, and the subsequent ﬂooding of the village. In this paper, we present the back analysis of this event. The complexity of this case study is due to several peculiar characteristics, but above all, to the clogging of the culvert, a phenomenon difﬁcult to simulate numerically. The methodology used for the reconstruction of the event is based on a multidisciplinary approach. A geological ﬁeld investigation was carried out to characterize the catchment and assess the availability of debris. Then, a cascade of numerical models was employed to reconstruct the debris ﬂow: the FLO-2D software was used to model the runoff along the hydrographic network while the mobile-bed debris ﬂow TRENT2D model, available through the WEEZARD system, was used to quantify both the erosion and deposition processes that occurred during the event. To simulate the culvert clogging, a novel modelling procedure was developed and applied. Despite the challenging framework, the results, in terms of debris volume, erosion rates, deposition area, and timing of the culvert obstruction, agree reasonably well with the observed data. It is worth noticing that these results were obtained mainly using parameters set a priori , namely calibrated on a physical basis. This proves that the proposed methodology is robust and effective, with good predictive capability. Therefore, it may be considered, according to the European Union (EU) Flood Directive, an “appropriate practice and the best available technology that does not imply excessive costs” to support predictive hazard mapping of situations as the one here considered.


Introduction
Debris flows are fluxes composed of a mixture of water and sediments with variable solid concentrations that commonly occur in steep drainage paths of mountain catchments [1][2][3]. The main triggering factor is represented by high-intensity rainfall events [4], which can lead to the formation of such flows [2]. Debris flows can show both highly impulsive features, when localized phenomena such as slope or clogged section failures occur, and smoother, although unsteady, behaviors when the source of sediments is mainly bed erosion. In these last cases, the features of the flow can be related to the hydrological forcing.
In any case, these types of flows have a high destructive potential due to the high momentum involved (with velocities indicatively in the 1-20 m/s range [5]), the erosive capability and associated deposition of debris material. For these reasons, they are one of the most dangerous hydrogeological phenomena in urbanized mountainous areas, often causing fatalities and extensive damage [6,7]. In addition, climate change contributes to an The Risola catchment study area is completely included in the units of the Pleistocene Mt. Amiata volcanic complex [30,31] (Figure 2). For most of its course, the Risola flows within the Quaranta Formation, except for the summit section, where the Bellaria Figure 1. The area between the top of Mt. Amiata and Abbadia San Salvatore. Symbology-Black line: catchment of the Fosso Acquagialla-Fosso Ermeta, Fosso Fonte Risola, and Torrente Vivo creeks; light blue and orange lines: main natural and buried hydrographic network, respectively; light blue dots: Tuscany Region rain gauges (AS: Abbadia San Salvatore-TOS07000001; LV: Laghetto Verde-TOS11000114; VA: Vetta Amiata-TOS11000115; [28]); red star: inflow point of the debris flow simulation within the Risola catchment; red dot: the beginning of the culvert in via Fosso Canali; turquoise polygons: "Laghetto Verde" (west), and "Laghetto Muraglione" (east). Background: orthophotos year 2016, modified from [29]. Coordinate system: Gauss-Boaga, west zone.
The urban area of Abbadia San Salvatore includes two artificial reservoirs "Laghetto Verde" (the western turquoise polygon in Figure 1) and "Lago Muraglione" (the eastern turquoise polygon in Figure 1), which collect part of the water coming from the Fosso Acquagialla-Fosso Ermeta and the Risola creeks, respectively. Moreover, these creeks have been buried at the entrance to the town. For the Risola, the culvert begins at the intersection with Fosso Canali street, at an altitude of about 845 m a.s.l.
The land cover, excluding urban areas, is characterized by the presence of a dense forest consisting of chestnut trees up to an altitude of 900 m a.s.l. and beech forest at higher altitudes.
The Risola catchment study area is completely included in the units of the Pleistocene Mt. Amiata volcanic complex [30,31] (Figure 2). For most of its course, the Risola flows within the Quaranta Formation, except for the summit section, where the Bellaria Formation crops out. The Quaranta Formation is made up of two members, both with essentially trachidacitic compositions, but with different textural characteristics [30]: the Marroneto member, characterized by alternating layers from a few centimetres to several meters thick with varying frequency of volcanic glass containing phenocrysts and xenoliths, and the overlying Leccio member, consisting of blocky lava flows. The Bellaria Formation consists of reddish-grey trachytic-trachidacitic lavas with rounded mafic magmatic enclaves [30,32].
Formation crops out. The Quaranta Formation is made up of two members, both with essentially trachidacitic compositions, but with different textural characteristics [30]: the Marroneto member, characterized by alternating layers from a few centimetres to several meters thick with varying frequency of volcanic glass containing phenocrysts and xenoliths, and the overlying Leccio member, consisting of blocky lava flows. The Bellaria Formation consists of reddish-grey trachytic-trachidacitic lavas with rounded mafic magmatic enclaves [30,32]. To the north of the Risola catchment, the Ermeta Member of the Ermeta-Macinaie Formation, crops out: it is a reddish to brownish, sub-porphyritic aphanitic lava flow with rare clinopyroxene and plagioclase microphenocrysts [30].
Finally, to the east of the study area, the Northern Apennines Ligurian units crops out, which are represented by the Santa Fiora Formation and the Basalts. The Santa Fiora Formation [30] consists of grey-brown argillites and calcilutites. The age of each member and formation is reported in the legend of Figure 2.
The volcanic formations are affected by widespread saprolite weathering processes [33][34][35] of spatially variable intensity, particularly marked in the case of the Quaranta Formation. Weathering generally originates from discontinuities and progressively extends to the entire mass, determining different degrees of decomposition [35]: (i) fractured bedrock; (ii) residual rock blocks (corestones) that may either be immersed in a loose sandy matrix or produce vast ground accumulations after removal of the loose fraction by runoff; (iii) sandy soils. Types (ii) and (iii) cover almost the whole study area.
These processes are important as they can affect slope stability and failure modes [36,37], as well as control the availability of unconsolidated material on slopes and along the hydrographic network that can be mobilized during debris flows. In particular, the wide distribution of saprolite allows us to state that the Risola catchment has an unlimited solid availability [38]. To the north of the Risola catchment, the Ermeta Member of the Ermeta-Macinaie Formation, crops out: it is a reddish to brownish, sub-porphyritic aphanitic lava flow with rare clinopyroxene and plagioclase microphenocrysts [30].
Finally, to the east of the study area, the Northern Apennines Ligurian units crops out, which are represented by the Santa Fiora Formation and the Basalts. The Santa Fiora Formation [30] consists of grey-brown argillites and calcilutites. The age of each member and formation is reported in the legend of Figure 2.
The volcanic formations are affected by widespread saprolite weathering processes [33][34][35] of spatially variable intensity, particularly marked in the case of the Quaranta Formation. Weathering generally originates from discontinuities and progressively extends to the entire mass, determining different degrees of decomposition [35]: (i) fractured bedrock; (ii) residual rock blocks (corestones) that may either be immersed in a loose sandy matrix or produce vast ground accumulations after removal of the loose fraction by runoff; (iii) sandy soils. Types (ii) and (iii) cover almost the whole study area.
These processes are important as they can affect slope stability and failure modes [36,37], as well as control the availability of unconsolidated material on slopes and along the hydrographic network that can be mobilized during debris flows. In particular, the wide distribution of saprolite allows us to state that the Risola catchment has an unlimited solid availability [38].
The spatial variation in the bedrock weathering degree also leads to variations in the slope steepness along the Risola course, with the development of areas where either erosive phenomena or deposition of solid load prevail. The course of the Risola follows an abrupt change of direction at an altitude of about 900 m a.s.l., passing from NW-SE to SW-NE (red star in Figure 1), as the effect of the interference with a km-scale SW-NE tectonic lineament ( Figure 2). Upstream of this area, for a few hundred meters, the watercourse gradient is very low, implying the deposition of most of the coarse solid material (Figure 3a). Figure 3b refers to the lower part of the Risola upstream the culvert, near Fosso Canali: here, the course average slope is about 7 • , and both in the riverbed and along the banks, high volumes of debris material, with boulders of size up to 0.1-1 m, are observed.
The spatial variation in the bedrock weathering degree also leads to variations in the slope steepness along the Risola course, with the development of areas where either erosive phenomena or deposition of solid load prevail. The course of the Risola follows an abrupt change of direction at an altitude of about 900 m a.s.l., passing from NW-SE to SW-NE (red star in Figure 1), as the effect of the interference with a km-scale SW-NE tectonic lineament ( Figure 2). Upstream of this area, for a few hundred meters, the watercourse gradient is very low, implying the deposition of most of the coarse solid material ( Figure  3a). Figure 3b refers to the lower part of the Risola upstream the culvert, near Fosso Canali: here, the course average slope is about 7°, and both in the riverbed and along the banks, high volumes of debris material, with boulders of size up to 0.1-1 m, are observed. Morphometric analysis of the basin, whose main data are summarized in Table 1, shows that the Melton index of the Risola catchment is 0.48, approximately corresponding to the threshold value of 0.5, beyond which the type of transport expected is that related to debris flows [39]. The propensity for the development of intense solid transport phenomena during heavy rainfall is associated with the large availability of potentially mobilizable sediment.

The 27-28 July 2019 Event
On the night between 27-28 July 2019, an intense meteorological event affected the entire territory of the Tuscany Region. In the Mt. Amiata area, starting at 22:00 on 27 July, the AS rain gauge (Figure 1) recorded rainfall of 210 mm/4 h, with maximum intensity peaks of 45 mm/15 min, never recorded before during the summer months [40].
The processing of the available data for rainfall in the Mt. Amiata area [28] with a 15′time scan shows that the highest intensity occurred between 22:45 and 23:15 h (Figure 4). In this short period, the rain gauges AS, LV and VA recorded 123, 86, 85 mm of cumulative rainfall, respectively.
The estimated return periods associated with the maximum values recorded in the different durations (1-3-6-12-24 h) are well beyond 200 years [40]. The maximum rainfall values over 12 and 24 h are greater than the recorded historical maximum by more than Morphometric analysis of the basin, whose main data are summarized in Table 1, shows that the Melton index of the Risola catchment is 0.48, approximately corresponding to the threshold value of 0.5, beyond which the type of transport expected is that related to debris flows [39]. The propensity for the development of intense solid transport phenomena during heavy rainfall is associated with the large availability of potentially mobilizable sediment. On the night between 27-28 July 2019, an intense meteorological event affected the entire territory of the Tuscany Region. In the Mt. Amiata area, starting at 22:00 on 27 July, the AS rain gauge (Figure 1) recorded rainfall of 210 mm/4 h, with maximum intensity peaks of 45 mm/15 min, never recorded before during the summer months [40].
The processing of the available data for rainfall in the Mt. Amiata area [28] with a 15 -time scan shows that the highest intensity occurred between 22:45 and 23:15 h (Figure 4). In this short period, the rain gauges AS, LV and VA recorded 123, 86, 85 mm of cumulative rainfall, respectively.
The estimated return periods associated with the maximum values recorded in the different durations (1-3-6-12-24 h) are well beyond 200 years [40]. The maximum rainfall values over 12 and 24 h are greater than the recorded historical maximum by more than 60%, which are, to all intents and purposes, absolute outliers, reflecting an extreme weather phenomenon, at least for the available documentation.
The very heavy rain caused widespread flooding in the urban area of Abbadia San Salvatore, but the situation of greatest damage occurred in the area of Via Fosso Canali (Figure 5b-e), where areas of both erosion (the asphalt was eroded) and accumulation of coarse debris, mixed with wood, developed. deposited in the urban area, as well as the debris present in the riverbed a allowed us to identify the type of phenomenon that occurred: the prevale gravel material, the occurrence of metric-sized boulders and the negligible cohesive material suggest that the phenomenon was essentially a flow of w granular material. Finally, the presence of a flat river stretch upstream of the red star in the indication that the clogging of the culvert occurred essentially in conjun

Methodologies
This Section describes the modelling chain, the specific models chosen for each chain steps and the approach developed for the simulation of culvert clogging that will be used in the back-analysis of the 27-28 July 2019 event. Nevertheless, the geological and geomorphological analysis, already described in Section 2 for convenience, must be considered as the first step of the whole methodology.  The flooding of the urban area has been due to the occlusion of the Risola culvert (dimensions 1.40 m wide × 1.00 m high), caused by a clogging phenomenon of the transported material (mainly debris with some wood). From the testimonies collected among the inhabitants of the area (video [41] and interviews), the clogging of the Risola culvert occurred between 23:30 and 23:45, i.e., almost coeval with the most intense phase of the rainfall. The event affected houses, commercial activities, and cars, causing EUR millions of damages [42].
Further quantitative information about the event, useful for back-analysis of the phenomena, was obtained from the scars caused by the impacts of solid material on the barks of some trees near the creek bed upstream of the culvert. They indicate that the maximum level of the flow above the actual stream bed was approximately 0.5 m, but in particularly narrow sections measuring (<5 m) 1.5 m were reached. The visual analysis of the material deposited in the urban area, as well as the debris present in the riverbed after the event, allowed us to identify the type of phenomenon that occurred: the prevalence of sandy-gravel material, the occurrence of metric-sized boulders and the negligible amount of fine cohesive material suggest that the phenomenon was essentially a flow of water with loose granular material.
Finally, the presence of a flat river stretch upstream of the red star in Figure 1, and the indication that the clogging of the culvert occurred essentially in conjunction with the rainfall peak, suggest that the phenomenon can be categorized among the debris flows with hydrological forcing.

Methodologies
This section describes the modelling chain, the specific models chosen for each chain steps and the approach developed for the simulation of culvert clogging that will be used in the back-analysis of the 27-28 July 2019 event. Nevertheless, the geological and geomorphological analysis, already described in Section 2 for convenience, must be considered as the first step of the whole methodology.
As described in Section 2.1, the event can be classified among the debris flows with hydrological forcing. In this case, the formation of the debris flow can be schematized as in Figure 6 (dashed rectangles): surface runoff forms in the upper part of the basin where sediment transport is negligible; then, a stretch of debris flow formation follows downstream, characterized by a flow of water and sediments deriving prevalently from a bed erosion; this is followed by a final stretch, where some further erosion may occur, but the main process is the deposition of the material transported by the flow [43].
In many cases, as in the present one, the detailed studies focus only on the area around the location of settlements subject to possible hazards, and usually this area corresponds to the deposition zone (rightmost dashed rectangle, Figure 6). In such conditions, it is possible to use an approach based on a cascade of models (filled rectangles, Figure 6): a hydrological model, closed in a suitable section close to the beginning of the study area, can be used to estimate a liquid hydrograph (Section 3.1); the debris flow formation transept is not described with a model, but we limit to evaluate the solid and mixture discharges as a function of time in the input section of the debris-flow model by using a suitable "amplification factor" (Section 3.2); finally, the simulation of the mixture flow in the last part of the basin can be obtained with a debris-flow model (Section 3.3). However, the Risola debris flow presents the clogging of a closed section (culvert), a phenomenon that debris flow models do not commonly consider. Therefore, we had to develop a specific procedure to face this aspect; a procedure which is presented in Section 3.4.
The input dataset used for the hydrological modelling was implemented employing the ArcGIS 10.7™ software (ESRI, Redlands, CA), while the one used for the debris flow modelling was implemented with both ArcGIS 10.7™ and the WEEZARD system [44].
Geosciences 2022, 12, x FOR PEER REVIEW 8 of 2 a suitable "amplification factor" (Section 3.2); finally, the simulation of the mixture flow in the last part of the basin can be obtained with a debris-flow model (Section 3.3). How ever, the Risola debris flow presents the clogging of a closed section (culvert), a phenom enon that debris flow models do not commonly consider. Therefore, we had to develop specific procedure to face this aspect; a procedure which is presented in Section 3.4. The input dataset used for the hydrological modelling was implemented employin the ArcGIS 10.7™ software (ESRI, Redlands, CA), while the one used for the debris flow modelling was implemented with both ArcGIS 10.7™ and the WEEZARD system [44].

Hydrological Modelling: The FLO-2D Model
Following the conceptual scheme of Figure 6, the first step to be considered is th hydrological modelling of the event. As the Risola is an ungauged catchment, a phys cally-based rainfall-runoff model was necessary to determine the liquid hydrograph i the required section. For this purpose, we used the FLO-2D model [12,13], a hydrologic hydraulic model that simulates the overland flow (over slopes and in channelized reaches fed by the rain (source term) and drained by infiltration (sink term). FLO-2D is a wel established model and it is widely used for hydraulic modelling of debris flows, as we as for hydrological modelling, as in this study, e.g. [45][46][47][48].
The underlying mathematical approach to flow propagation is based on the dynami wave approximation of the mass and momentum equations. The relevant differentia equations are integrated over a regular Cartesian grid defining the topography of the com putational domain. FLO-2D model has three possible different approaches to describe the infiltratio [13]. The formulation used in this study case is the Green-Ampt (G-A) equation [49], physically-based equation based on the soil porous media characteristics [50], which i widely used in hydrologic studies as it has a relatively low computational demand whil providing reasonably accurate physical mechanisms [51][52][53][54]. When the maximum infi tration capacity is exceeded, surface flow is generated. Moreover, FLO-2D has the capa bility of simulating the presence of impermeable layers by assuming a limit to the infiltra tion depth.
In the following section, we describe the physical quantities (parameters) require by FLO-2D to set up a hydrological simulation with the G-A approach.
Data Input of the Hydrological Modelling

Hydrological Modelling: The FLO-2D Model
Following the conceptual scheme of Figure 6, the first step to be considered is the hydrological modelling of the event. As the Risola is an ungauged catchment, a physicallybased rainfall-runoff model was necessary to determine the liquid hydrograph in the required section. For this purpose, we used the FLO-2D model [12,13], a hydrologichydraulic model that simulates the overland flow (over slopes and in channelized reaches) fed by the rain (source term) and drained by infiltration (sink term). FLO-2D is a wellestablished model and it is widely used for hydraulic modelling of debris flows, as well as for hydrological modelling, as in this study, e.g. [45][46][47][48].
The underlying mathematical approach to flow propagation is based on the dynamic wave approximation of the mass and momentum equations. The relevant differential equations are integrated over a regular Cartesian grid defining the topography of the computational domain. FLO-2D model has three possible different approaches to describe the infiltration [13]. The formulation used in this study case is the Green-Ampt (G-A) equation [49], a physicallybased equation based on the soil porous media characteristics [50], which is widely used in hydrologic studies as it has a relatively low computational demand while providing reasonably accurate physical mechanisms [51][52][53][54]. When the maximum infiltration capacity is exceeded, surface flow is generated. Moreover, FLO-2D has the capability of simulating the presence of impermeable layers by assuming a limit to the infiltration depth.
In the following section, we describe the physical quantities (parameters) required by FLO-2D to set up a hydrological simulation with the G-A approach.

Data Input of the Hydrological Modelling
The following parameters are required independently of the infiltration method used.
• DEM: to obtain the domain for the case study, two Digital Elevation Models (DEMs) were merged, one based on LiDAR acquisition with a cell size of 1 m, and one with a cell size of 10 m, since the first one did not cover the Mt. Amiata summit area [55]. • n: the Manning roughness [s/m 1/3 ] was spatially assigned reclassifying the land use map [56], according to the value proposed by [57,58]. • T: this parameter expresses the minimum flow depth [m] before runoff volume is exchanged with adjacent grid elements and represents the volume of water stored in small depressions (puddles) that do not become part of the overland runoff or infiltration. • Rain: rainfall simulation was obtained by spatializing the rainfall data (mm/15') from the three available rain gauges ( Figure 1) using the Inverse Distance Weight (IDW) method [59].
The G-A equation, based on Darcy's law, is: where ∆t is the computational time step, ∆F is the change in infiltration (mm) over the computational time step, F(t) is total infiltration at time t. The parameters of these equations are: K is the saturated hydraulic conductivity (mm/h) and finally γ = (PSIF + Head)·DTHETA where PSIF is capillary suction (mm), Head is the incremental rainfall for the time step plus flow depth on the grid element (mm) and DTHETA is volumetric soil moisture deficit (dimensionless) function of the effective porosity of the soil. Finally, the following spatially distributed parameters, affecting the runoff calculation, must be assigned.

•
Impervious areas (%): the portion of a grid element that is impervious to infiltration (e.g., rock outcrops or urban areas). • Initial abstraction (mm): the initial loss of rainfall that precedes infiltration and excess rainfall-runoff (e.g., vegetation interception).

The Debris-Flow Mixture Hydrograph
The mixture hydrograph to be used as boundary condition in the inflow section of the debris flow model, is calculated based on the liquid hydrograph obtained with the hydrological model. The relation, proposed by [2], is the following: where Q mix (t) is the mixture discharge, Q hydrol (t) is the liquid discharge and c is a reference concentration, which can be estimated through the following relation: where i f is the average bed slope near the inflow section, ∆ = (ρ s − ρ w )/ρ w is the relative reduced density of the solid material, and φ d is the dynamic friction angle of the material.

Debris Flow Modelling: The TRENT2D Model
Since the solid part of the debris flow we are interested in is characterized by a loose coarse material (see the end of Section 2.1) and the flow occurs both on a fixed and mobile bed, the model we have chosen is the TRENT2D model. TRENT2D (an acronym for Transport in Rapidly Evolving, Natural Torrent), represents one of the most advanced two-dimensional models for the study of granular debris flows and intense solid transport phenomena. In this Section we describe the main mathematical features of the model, addressing the reader to Appendix A and the original paper for more specific mathematical aspects. In addition, in Section 3.3.1, we highlight the characteristics of the WEEZARD system used for the management of the input data, launching the simulations, and for the results analysis.
The TRENT2D uses a two-phase isokinetic description of the mixture that composes the flow, where the solid and liquid components are considered and described, from the point of view of the equations of motion, in a separate but interacting way. A second characterizing element concerns the description of the bed, over which the debris flow flows. In the original formulation of the TRENT2D, the movement of the mixture occurs on a mobile bed, i.e., on a bed that can be subject to erosion and deposition during the motion. However, a more recent version of the model [60] also allows the consideration of cases where the mixture flows on a non-erodible bed. In the latter case, the model provides for the possibility of depositing material above the non-erodible bed (thus becoming, at least temporarily erodible), but prevents the bed from reaching, due to erosion, elevation below the level of the non-erodible layer. A third aspect, intimately related to the description of the bed, concerns the composition of the mixture: as it is a two-phase model, the average volumetric concentration of sediment on the vertical is variable in space and time and is related to the local hydrodynamic conditions. The formulation used is as follows: where c b is the concentration in the bed (equal to 1 − p, where p is the porosity of the material) considered constant, β is a transport parameter to be determined based on the flow and sediment characteristics, → u is the depth-averaged velocity vector of the mixture, g is the acceleration of gravity and, finally, h is the flow depth (difference in elevation between the free surface of the flow and the bed). The variability in concentration is related to the exchange of sediment and water that occurs between the bed and the flow during erosion and deposition phenomena. It is worth noting that with a two-phase approach, the mixture cannot stop as a whole, but only the solid phase (and the liquid part necessary for the saturation of the material) stops, while the liquid phase flows away. A fourth aspect concerns the formulation of the global law of resistance to motion. It is given by the resistance of the liquid phase plus the resistance of the solid phase, which, in the flows, is predominant and is mainly related to collisional phenomena. The law used in the TRENT2D model is that proposed by [61] and considers only the collisional contribution: where → τ 0 is the bed shear stress, ρ w is the water density, a = 0.32 according to [61], ∆ = (ρ s − ρ w )/ρ w is the relative reduced density of the solid material, φ d is the dynamic friction angle of the material, is the submergence, where d is a characteristic diameter of the transported material and h is the local mixture depth.
In the case of a fixed bed, concentration is no longer connected to the local flow conditions (namely, Equation (4) does not hold), and it becomes a simple advected quantity. As for the bed shear stress, a classical Strickler relation was used, namely: where K s is the Strickler's coefficient.
WEEZARD (WEbgis modElling and haZard Assessment for mountain flows: An integRated system in clouD) is a system that allows the management of all phases of data preparation, execution of simulations, analysis of results and production of hazard maps for a range of phenomena that may affect mountain environments (debris flows, intense solid transport phenomena, mudflows, dense snow avalanches). The key difference from a simple management interface that normally operates locally on a computer is that it is based on an innovative web services system [44] accessed through a common network browser [62] (Figure 7). preparation, execution of simulations, analysis of results and production of hazard maps for a range of phenomena that may affect mountain environments (debris flows, intense solid transport phenomena, mudflows, dense snow avalanches). The key difference from a simple management interface that normally operates locally on a computer is that it is based on an innovative web services system [44] accessed through a common network browser [62] (Figure 7). This approach has several advantages over a local system, first the fact that it is not necessary to have a high-performance computer to be able to do even computationally heavy simulations as a shared network server is used. It is just the sharing of the resource of calculation the key aspect that allows the use of sophisticated instruments of simulation at low cost, as demanded by the European directive [10]. A second element that characterizes this system is its ease of use. Indeed, WEEZARD is a word that recalls another English word, wizard, that in the computer science language indicates a procedure, generally incorporated into a more complex application, that allows the user to carry out determined operations (usually articulated) through a series of simple steps in succession. Having wizards to easily perform a complicated set of operations, such as those required to simulate natural hazards in a mountain environment, is a fundamental prerequisite for transforming an inherently complex system into an effective, efficient, and economical tool for professional use.
In extreme synthesis, the process that, from the acquisition of basic data, leads to the result of the elaboration of a hazard map, can be divided into a series of logical steps (Figure 8), independent of the type of natural hazard one wants to simulate. For each of these steps, a series of functionalities are available that allow us to easily perform a whole series of operations related to that logical step (each functionality is illustrated through video tutorials [63]) and that greatly facilitate the learning and use of the system. This approach has several advantages over a local system, first the fact that it is not necessary to have a high-performance computer to be able to do even computationally heavy simulations as a shared network server is used. It is just the sharing of the resource of calculation the key aspect that allows the use of sophisticated instruments of simulation at low cost, as demanded by the European directive [10]. A second element that characterizes this system is its ease of use. Indeed, WEEZARD is a word that recalls another English word, wizard, that in the computer science language indicates a procedure, generally incorporated into a more complex application, that allows the user to carry out determined operations (usually articulated) through a series of simple steps in succession. Having wizards to easily perform a complicated set of operations, such as those required to simulate natural hazards in a mountain environment, is a fundamental prerequisite for transforming an inherently complex system into an effective, efficient, and economical tool for professional use.
In extreme synthesis, the process that, from the acquisition of basic data, leads to the result of the elaboration of a hazard map, can be divided into a series of logical steps (Figure 8), independent of the type of natural hazard one wants to simulate. For each of these steps, a series of functionalities are available that allow us to easily perform a whole series of operations related to that logical step (each functionality is illustrated through video tutorials [63]) and that greatly facilitate the learning and use of the system.

A Novel Approach for the Simulation of Culvert Clogging
As a result of the presence of large sediments and boulders in the flow (see Section 2), we assumed that clogging occurs once the level of the debris flow reaches the upper level of the culvert. Therefore, we distinguish the following: (i) a pre-clogging phase, lasting from the beginning of the event up to the instant in which the debris flow level, near the culvert inlet, exceeds the value of the culvert coverage level (conventional clogging time); and (ii) a post-clogging phase, characterized by an obstructed section near the culvert inlet, lasting from the conventional clogging time up to the end of the phenomenon. In the following sections we present the strategy and the input data required by the TRENT2D model to implement the two phases.
It is worth noting that, since the a priori exact clogging time is unknown, a complete simulation with the modified computational domain (as specified further on) must be performed and, a posteriori, the conventional clogging time can be identified.

Strategy and Data Input of the Pre-Clogging Phase
In the pre-clogging phase, to describe the motion in the culvert, its closed section must be transformed, in the computational domain, into an open channel with walls up to the ground level in which a free-surface flow can occur.
The parameters used in the pre-clogging phase, and how they were evaluated are described below.

•
Computational domain: the first element to be defined is the computational domain constituted by a square Cartesian grid defining the altimetric base of the study area. Moreover, non-erodible areas and erodible areas with the associated depth of the erodible material can be specified.

•
Inflow conditions: the inflow condition must be imposed in the section where the debris flow is expected to enter the computational domain. By using the WEEZARD system, the mixture hydrograph is evaluated automatically, given and the inflow section as described in Section 3.2.
The following parameters can be chosen a priori, since they are closely related to measurable properties, and kept constant in all the simulations.

•
Debris parameters: include the values of the sediment bed concentration , the dynamic friction angle , and the reduced relative sediment density Δ.
On the other hand, the following parameters can be estimated as follows.

A Novel Approach for the Simulation of Culvert Clogging
As a result of the presence of large sediments and boulders in the flow (see Section 2), we assumed that clogging occurs once the level of the debris flow reaches the upper level of the culvert. Therefore, we distinguish the following: (i) a pre-clogging phase, lasting from the beginning of the event up to the instant in which the debris flow level, near the culvert inlet, exceeds the value of the culvert coverage level (conventional clogging time); and (ii) a post-clogging phase, characterized by an obstructed section near the culvert inlet, lasting from the conventional clogging time up to the end of the phenomenon. In the following sections we present the strategy and the input data required by the TRENT2D model to implement the two phases.
It is worth noting that, since the a priori exact clogging time is unknown, a complete simulation with the modified computational domain (as specified further on) must be performed and, a posteriori, the conventional clogging time can be identified.

Strategy and Data Input of the Pre-Clogging Phase
In the pre-clogging phase, to describe the motion in the culvert, its closed section must be transformed, in the computational domain, into an open channel with walls up to the ground level in which a free-surface flow can occur.
The parameters used in the pre-clogging phase, and how they were evaluated are described below.

•
Computational domain: the first element to be defined is the computational domain constituted by a square Cartesian grid defining the altimetric base of the study area. Moreover, non-erodible areas and erodible areas with the associated depth of the erodible material can be specified. • Inflow conditions: the inflow condition must be imposed in the section where the debris flow is expected to enter the computational domain. By using the WEEZARD system, the mixture hydrograph is evaluated automatically, given Q hydrol and the inflow section as described in Section 3.2.
The following parameters can be chosen a priori, since they are closely related to measurable properties, and kept constant in all the simulations.

•
Debris parameters: include the values of the sediment bed concentration c b , the dynamic friction angle φ d , and the reduced relative sediment density ∆.
On the other hand, the following parameters can be estimated as follows.
• Y : can be evaluated by using an average value of the transported grain size and a reference depth along the flow. • β : is determined by assuming a local uniform flow condition in the input section with the resistance law according to [61]: where c is given by Equation (3). β is calculated directly by the WEEZARD system based on the other input parameters.
• K s : the Gauckler-Strickler roughness coefficient can be assigned based on the land use categories [56], according to [57,58].

Strategy and Data Input of the Post-Clogging Phase
Starting from the conventional clogging time, the mixture stops flowing into the culvert and the entrance section becomes an impermeable wall. To describe this phase, the altimetry of the computational domain must be changed accordingly.

•
Instead of the artificial channel bed, the original ground level must be restored along the culvert path.

•
In the channel upstream of the culvert, from the inlet section to the culvert initial section, the original DEM must be changed according to the bed modifications obtained in the pre-clogging phase. The resulting DEM can be used as the initial condition for this phase.
The values of all the parameters used in the pre-clogging phase can also be maintained in the post-clogging phase.

Results
In this section, we present the results of the back-analysis of the Risola debris flow obtained from the application of the methodology described in the previous sections. To allow the reproducibility of the results, for each step of the modelling we also provide the values of the employed parameters and, when available, a comparison of the results with measured/estimated data.

Hydrological Modelling: Parameter Values and Results
The parameters used for the hydrological modelling are herein described.
• DEM: the two DEMs with different spatial resolutions were mosaicked and resampled to a common cell size of 5 m, obtaining a good compromise between the representation of the hydrographic network and the computational time necessary to run the simulation. • TOL: The minimum value of 0.0012 m was used, as recommended in the FLO-2D manual [13]. Using higher values would result in larger quantities of water that do not become part of the overland runoff or infiltration. A situation that, due to the characteristics of the study area, would not be justified from a hydrological point of view. • Rain: to apply the IDW method for the entire study area considering the data from the AS, LV, and VA rain gauges (Figure 1) only, an auxiliary gauge was placed south of the southern border of the Risola catchment, allowing us to extend the rainfall maps to the Risola catchment area. This was achieved by placing a barrier feature between the auxiliary gauge and the existing rain gauges without altering their rainfall data. A barrier is a polyline dataset used as a breakline that limits the searching region for input sample points: only those points located on the same side of the barrier are considered. Thus, a representation of the rainstorm distributed in time and space was obtained.
• K: a specific survey was conducted through constant-head measurements with an Aardvark-type permeameter [64] to determine the value of field saturated hydraulic conductivity. Thirteen borehole tests were performed within the sandy soil above the BEL2 and QRT volcanic formations, investigating depths between 0.30-1.50 m. Given the homogenous nature of the unconsolidated material in the study area, the average value (140 mm/h) was uniformly assigned throughout the catchment. • DTHETA: in the five days preceding the event, the three rain gauges close to the study area ( Figure 1) recorded no rainfall accumulation, and considering a longer period (i.e., 60 days), only 28 mm were recorded. This extended drought period and the permeable nature of the unconsolidated material indicate that the antecedent moisture conditions before the event were dry. Thus, the assumed DTHETA value for modelling is 0.35, as also specified in the FLO-2D manual [13]. Only in the last stretch of the Risola up to the culvert (~1 km) can the soil be considered saturated. In fact, some ephemeral springs with a flow rate of only a few l/s (negligible if compared to the discharge developed during the event) bring water into the stream bed. • PSIF: the obtained value of this parameter is equal to 30 mm. This is in agreement with the values proposed by [65].

•
Given the geological characteristics of the study area (Section 2), no impermeable boundary was imposed at the bottom of the unconsolidated material of the catchment.
Finally, the following parameters were spatially assigned by reclassifying the land use map [56] using the described criteria.
• Impervious areas (%): as described in Section 2, the bedrock of the study area is affected by widespread weathering. For this reason, infiltration is set to zero for the fractured bedrock portion (type (i) in Section 2;~15% of the total catchment area) even though no data are available on the rock mass fracturing degree and the related infiltration patterns. It may be noted that this assumption is probably reasonable for the study area only in the case of short and intense rainstorms, such as the one analyzed here. • Initial abstraction (mm): it represents the initial loss of rainfall that precedes infiltration and excess rainfall-runoff (e.g., vegetation interception). This parameter was assigned using the values proposed in the FLO-2D manual [13].
The input dataset used for the hydrological modelling is summarized in Table 2.
The resulting hydrograph is presented in Figure 9 (green line) and it is referred to where the closure section of the basin coincides with the inflow section of the debris flow model, indicated with a red star in Figure 1.
• : a specific survey was conducted through constant-head measurements with an Aardvark-type permeameter [64] to determine the value of field saturated hydrauli conductivity. Thirteen borehole tests were performed within the sandy soil above th BEL2 and QRT volcanic formations, investigating depths between 0.30-1.50 m. Given the homogenous nature of the unconsolidated material in the study area, the averag value (140 mm/h) was uniformly assigned throughout the catchment.
• : in the five days preceding the event, the three rain gauges close to the study area ( Figure 1) recorded no rainfall accumulation, and considering a longer period (i.e., 60 days), only 28 mm were recorded. This extended drought period and the per meable nature of the unconsolidated material indicate that the antecedent moistur conditions before the event were dry. Thus, the assumed value for model ling is 0.35, as also specified in the FLO-2D manual [13]. Only in the last stretch of th Risola up to the culvert (~1 km) can the soil be considered saturated. In fact, som ephemeral springs with a flow rate of only a few l/s (negligible if compared to th discharge developed during the event) bring water into the stream bed. • : the obtained value of this parameter is equal to 30 mm. This is in agreemen with the values proposed by [65].
• Given the geological characteristics of the study area (Section 2), no impermeabl boundary was imposed at the bottom of the unconsolidated material of the catch ment.
Finally, the following parameters were spatially assigned by reclassifying the land use map [56] using the described criteria.
• Impervious areas (%): as described in Section 2, the bedrock of the study area is affected by widespread weathering. For this reason, infiltration is set to zero for the fractured bedrock portion (type (i) in Section 2; ⁓15% of the total catchment area) even though no data are available on the rock mass fracturing degree and the related infiltration patterns. It may be noted that this assumption is probably reasonable for the study area only in the case of short and intense rainstorms, such as the one analyzed here • Initial abstraction (mm): it represents the initial loss of rainfall that precedes infiltra tion and excess rainfall-runoff (e.g., vegetation interception). This parameter was as signed using the values proposed in the FLO-2D manual [13].
The input dataset used for the hydrological modelling is summarized in Table 2.
The resulting hydrograph is presented in Figure 9 (green line) and it is referred t where the closure section of the basin coincides with the inflow section of the debris flow model, indicated with a red star in Figure 1.

The Debris Flow Hydrograph
The mixture discharge (Figure 9) was obtained using Equations (2) and (3) as reported in Section 3.2. The parameter values employed are described below.

•
Debris parameters: the values of c b , φ d , and ∆ parameters used in this study are those generally accepted in debris flow dynamics [61] and are reported in Table 3. This assumption stems from the fact that the material involved does not have any peculiar geotechnical characteristics that would necessitate changing the values of these three parameters.  Figure 1, where an abrupt change in direction and slope of the creek occurs (see Section 2). Table 3. WEEZARD input parameters for the pre-clogging simulation.

Pre-Clogging Debris Flow Modelling: Parameter Values and Results
The parameters used for the pre-clogging modelling are reported in the following list with a description of the reason for their choice.

•
Computational domain: an available DEM [55] with a cell size of 1 m × 1 m, based on LiDAR acquisition, was used to produce the domain, which is delimited with a red dashed rectangle in Figure 10a. The resulting morphology of the inflow section ( Figure 10b) is referred to as the red star in Figure 10a. To represent the culvert width with a reasonable approximation and considering that its width is 1.4 m, the original DEM was resampled to a cell size of 0.5 m × 0.5 m. In this way, the channel is represented by using three cells, namely with a width of 1.5 m, and the deviation of 0.1 m was assumed to be negligible for the modelling results. In addition, as the base of the culvert is located 1.5 m below the road level, the resampled DEM was lowered by 1.5 m along the path of the culvert. In the analyzed stretch, some anthropogenic structures are present, thus erodible, and non-erodible zones have been defined accordingly. One of these structures is a small building preceded by a short (~25 m) non-erodible artificial portion of the bottom of the creek, and a wall on the right bank of the stream (orange circle in Figure 10a). In addition, three check dams are located along the Risola: two of them are placed immediately downstream of the building described previously (see the green circle in Figure 10a), and the last check dam is located just upstream of the culvert inlet section. In this last area, two walls bound the creek before the entrance to the covered portion (yellow circle in Figure 10a). These structures, which are correctly represented by the DEM used for the modelling, are treated as non-erodible areas with null erodible thickness. On the other hand, the rest of the torrent is treated as an erodible area. • Y: to choose the optimal value for this parameter, we performed tests with values of 15, 10 and five. Results show that for values of 15 and 10 the flow depths are too shallow compared to those observed. At the same time, the erosions were excessive, generating a higher concentrated flow and, consequently, culvert clogging time earlier than observed. On the other hand, for Y = 5, both the flow depths along the stream and the culvert clogging time were found to be consistent with the observations. • β : the value of the transport parameter used for the numerical modelling was obtained considering the submergence as Y = 5, and considering the values of other parameters appearing in Equation (7) as reported in Table 3.
located just upstream of the culvert inlet section. In this last area, two walls bound the creek before the entrance to the covered portion (yellow circle in Figure 10a). These structures, which are correctly represented by the DEM used for the modelling, are treated as non-erodible areas with null erodible thickness. On the other hand, the rest of the torrent is treated as an erodible area.

•
: to choose the optimal value for this parameter, we performed tests with values of 15, 10 and five. Results show that for values of 15 and 10 the flow depths are too shallow compared to those observed. At the same time, the erosions were excessive, generating a higher concentrated flow and, consequently, culvert clogging time earlier than observed. On the other hand, for = 5, both the flow depths along the stream and the culvert clogging time were found to be consistent with the observations. • : the value of the transport parameter used for the numerical modelling was obtained considering the submergence as = 5, and considering the values of other parameters appearing in Equation (7) as reported in Table 3. In this phase of the back-analysis, the physical quantities we used to validate the numerical results was the clogging time, which occurred between 23:30 and 23:45, as described in Section 2.1. The results of the pre-clogging simulation, performed with the parameters in Table 3, show that the clogging time occurs at about 23:30 (i.e., ~3200 s in Figure 9), as depicted in Figure 11, thus in good agreement with what happened. In this phase of the back-analysis, the physical quantities we used to validate the numerical results was the clogging time, which occurred between 23:30 and 23:45, as described in Section 2.1. The results of the pre-clogging simulation, performed with the parameters in Table 3, show that the clogging time occurs at about 23:30 (i.e.,~3200 s in Figure 9), as depicted in Figure 11, thus in good agreement with what happened.

Debris Flow Back-Analysis: Post-Clogging Phase Results
The parameters used for this phase are the same as the pre-clogging scenario. quantities used during this phase to validate the results were the distribution of ero and deposited material along the Risola, and in the built-up area, as well as the maxim free surface levels reached by the flow.
In Figure 12a, the map of the morphological modification obtained using the clogging modelling is compared with the map of the debris flow deposition area obta by interpreting the archive photos of the event ( Figure 5). Overall, there is a good ag ment between the observed and modelled debris flow deposition area. The modelled ues of the deposit height in the built-up area result to be, on average, equal to about 0.5 m, while in some points they reach 1 m, in good agreement with the post-event ph (Figure 5b-e). The more accurate modelling results, in terms of both the flooded area maximum deposit heights, are related to the area marked in red (Figure 12a), where co material with an average thickness of 0.5 m was deposited (Figure 5b-e). On the o hand, in some of the areas characterized by thin deposit thicknesses and predomina fine material, delimited in blue (Figure 12a), the modelling results are less accura comparison with the real flooded area was not possible due to limited data being avai (see Section 2.1). However, the comparison between the numerical flooded area (i.e area where the mixture has flowed) and the estimated deposition area shows a g agreement (Figure 12b). In particular, it can be noticed that the area of the fine sedim deposition (in blue) was affected by the debris flow phenomena.

Debris Flow Back-Analysis: Post-Clogging Phase Results
The parameters used for this phase are the same as the pre-clogging scenario. The quantities used during this phase to validate the results were the distribution of eroded and deposited material along the Risola, and in the built-up area, as well as the maximum free surface levels reached by the flow.
In Figure 12a, the map of the morphological modification obtained using the postclogging modelling is compared with the map of the debris flow deposition area obtained by interpreting the archive photos of the event ( Figure 5). Overall, there is a good agreement between the observed and modelled debris flow deposition area. The modelled values of the deposit height in the built-up area result to be, on average, equal to about 0.1-0.5 m, while in some points they reach 1 m, in good agreement with the post-event photos (Figure 5b-e). The more accurate modelling results, in terms of both the flooded area and maximum deposit heights, are related to the area marked in red (Figure 12a), where coarse material with an average thickness of 0.5 m was deposited (Figure 5b-e). On the other hand, in some of the areas characterized by thin deposit thicknesses and predominantly fine material, delimited in blue (Figure 12a), the modelling results are less accurate. A comparison with the real flooded area was not possible due to limited data being available (see Section 2.1). However, the comparison between the numerical flooded area (i.e., the area where the mixture has flowed) and the estimated deposition area shows a good agreement (Figure 12b). In particular, it can be noticed that the area of the fine sediment deposition (in blue) was affected by the debris flow phenomena.
The areas where the asphalt cover was damaged were considered as erodible zones for the modelling. Simulation results indicate widespread erosion in these areas, in agreement with field observations. On the other hand, regarding the erosions that occurred along the creek shaft, no direct comparisons were possible with what actually happened during the event. However, the values obtained (maximum values of about 2 m and an average of 0.35 m) are consistent with the known evidence.
The results obtained indicate that the debris flow generated in the Risola eroded and mobilized a debris volume of 3000 m 3 , of which approximately 1500 m 3 partially buried the built-up areas, in line with the estimated value (Section 2.1). Part of this material reached the "Lago Muraglione" reservoir, contributing to its overflow, which in turn caused the flooding of its downstream area.
Finally, since the obtained results were considered satisfactory in terms of the timing of clogging, maximum heights reached by the flow, flooding extension and the volumes deposited, no further sensitivity analyses were performed. fine material, delimited in blue (Figure 12a), the modelling results are less acc comparison with the real flooded area was not possible due to limited data being a (see Section 2.1). However, the comparison between the numerical flooded area area where the mixture has flowed) and the estimated deposition area shows agreement (Figure 12b). In particular, it can be noticed that the area of the fine s deposition (in blue) was affected by the debris flow phenomena.

Discussion
In this paper, we have presented: (i) a multidisciplinary approach for the study of a debris flow, in which the clogging of a culvert plays a fundamental role in the dynamics of the occurrence; and (ii) the application of this methodology to a back-analysis of a case study. Results show that the event can be reproduced reasonably well, thus allowing the validation of the approach. In the following sections, we analyze critically what we think are the main issues of the used methodology in relation to the obtained results.

Geological and Geomorphological Analysis
The Mt. Amiata, where our case study is located, does not present geological and geomorphological features that are characteristic of the Alpine regions [66], namely with the widespread presence of Quaternary deposits of moraine and/or fluvioglacial origin that makes these areas prone to debris flows. Considering exclusively the existing cartographic data (Figure 2), the area presents a low susceptibility to these kinds of natural hazards due to the apparent lack of source sediments. In the Risola catchment, there are only small areas of eluvial-colluvial deposits, and two landslides, apparently connected to the hydrographic network. Despite this situation, the 2019 occurrence demonstrated that a source of debris had to be present. A field survey revealed the widespread availability in the Risola catchment of loose solid material caused by saprolitic weathering of the volcanic bedrock, both in the hydrographic network and on the slopes, as explained in detail in Section 2. Surveys aimed at the geological and geomorphological characterization of an area are therefore fundamental prerequisites for a back analysis but are equally important in predictive hazard mapping activities and in planning of hydrogeological risk reduction [67]. This is also emphasized by the ongoing climate change characterized by more frequent extreme rainstorms causing mass movement in areas usually not affected by such events.

The Choice of the Models
As highlighted by other studies, e.g., [68][69][70], the use of a cascade of models is an effective way to simulate the connection between rainfall, hydrological processes, transport of sediments and bed modifications in several debris flow phenomena characterized by a clear hydrologic forcing (the so-called runoff-generated debris flows, e.g., [71]). However, the choice of advanced, physically based models as cascade elements is somewhat compulsory when: (i) the available data does not allow a detailed calibration; (ii) we want to make predictive analyses. Since in this work the former applied, we choose to use the models described in Section 3 as both models use physically based parameters that can be determined a priori or with minor calibration.
As for the hydrological model, the FLO-2D with the Green-Ampt description of the infiltration appeared to be a suitable choice that utilized the knowledge of the basin derived from the field survey.
As for the debris flow model, the TRENT2D, based on a two-phase, mobile-bed approach, appeared to be the most appropriate in comparison with other commercial codes employing a single-phase approach, which is unable to catch some of the most important processes observed in the analyzed event, namely a granular debris flow. For example, the FLO-2D debris flow model uses a fixed-bed description that is unable to describe both erosions and depositions occurring during the flow, but the variation in the bed elevation can occur only when the single-phase fluid stops (final deposition), occurring when the bed shear stress is lower than a threshold (Bingham fluid [13]). The RAMMS model uses an entrainment function that changes the flowing volume, but in the calculations, the bed is kept fixed [72]. In any case, in single-phase models the mass eroded or deposited has a composition (represented by the solid volumetric concentration) that is different with respect to the flowing mass: the solid concentration of the flowing mass is always lower than the bed concentration [2,73]. For this reason, while a two-phase model is able to deal with the change in the concentration of the flowing mass and with the consequence of this change (for example, the bed shear stress is a function of the concentration), the single-phase models with constant density simply cannot.
Results of the back-analysis confirm that the choice of the models we have performed was reasonable and that this choice should ensure adequate accuracy in predictive situations.

The Culvert Clogging and Its Numerical Simulation
Although the clogging of a culvert (or other hydraulic structures) by debris flow is a widespread phenomenon affecting many mountainous urban areas, few case studies focused on its numerical modelling systematically. Moreover, the most popular and widely used hydraulic models, such as FLO-2D [13] or HEC-RAS [74], allow the clogging of hydraulic structures (e.g., bridges, culverts, and storm drain systems) but only by establishing a priori both the percentage and the time of blockage (see e.g., [75]). TELEMAC-2D software [76] also allows hydraulic modelling considering the presence of culverts, however, this model is used for the simulation of water floods in a fluvial environment (see, e.g., [77]) and is not able to simulate a blockage caused by debris flows or intense solid transport phenomena. In addition, these models are single-phase and fixed-bed models.
At the current state of the art, neither two-phase, mobile-bed models can manage directly the debris flow clogging of a closed section. Therefore, our choice to define a practical methodology was compulsory.
Indeed, we do not know when the clogging actually happens, the possible dependence on the grain size and the culvert area, and if a pressurized flow occurs before the clogging. A thorough experimental study is therefore necessary to shed light on this issue but, in the meantime, a reasonable practical approach is desirable.
Despite the simplicity of our assumptions (Section 3.4), results are consistent with reality, in terms of the timing of clogging as highlighted in the pre-clogging phase results (Section 4.3), as well as in terms of non-excessive erosion and deposition processes obtained at the end of the post-clogging simulation (Section 4.4).
We are aware that the successful application of the novel procedure in a specific case does not allow us to assert that it will work in any case. Pending further possible confirmation of the already planned experimental study, we are confident that this novel approach can be applied in hazard mapping activities.

Comparison of Simulated and Observed Data
Comparison of simulated and field data is never easy in the case of debris flows. One of the main issues is the lack of accurate post-event filed measurement. The reason is that these events are linked to dangerous situations, damage and even fatalities and therefore, since measurements are not a priority in these circumstances, they are often not promptly performed in undisturbed conditions. Moreover, many times measurements are not direct, and many others are only estimates, whose uncertainty cannot be determined.
As explained in Section 2.1, the mapping of the flooded area is subject to some inevitable uncertainty. The planimetric accuracy is therefore limited and any point comparison with simulated data has no meaning in the present case. In any case, it is possible to notice that in some areas the comparison is better while in others the difference is more significant. In particular, in the area with deposition of coarse debris and an average thickness of 0.5 m (delimited in red in Figure 12a), the results obtained with the modelling are consistent with what was observed, both in terms of thickness and flooded area. On the other hand, in the area of fine sediment deposition with an average thickness of 0.1 m (marked in blue), the results of the deposition height are less accurate, with calculated deposits lower than the estimated ones. It is worth noting that in the red area, the dimension of the sediments was significant, while the blue area was characterized by very fine sediments. This difference justifies the different performance of the used debris flow model: since the composition of the transported material is assumed to be constant, the model has been tuned to simulate the coarse material characterizing the bed in the upper part of the computational domain. Therefore, it cannot be so accurate in describing zones with very fine material.
Last but not least, in the present study, we have also used the temporal datum to make the comparison. This is a very lucky situation since such data is not commonly available.

Conclusions
After this study, thanks to the results obtained by modelling the debris flow event that occurred in the Risola catchment on 27-28 July 2019, it is possible to state that:

•
A geological-geomorphological model based on the integration of detailed mapping and specific engineering geological field investigations is a fundamental prerequisite for further accurate numerical debris flow modelling.

•
The overall methodological approach used in this study, namely a cascade of models composed of a hydrological model for the estimation of the hydraulic forcing, estimation of the mixture flow rate with the amplification factor, and two-phase modelling of the debris flow in the final reach, provides reasonable outputs for the simulation of debris flows with hydrological forcing.

•
The modelling assumptions introduced to overcome some of the limitations of the model used, namely the definition of the conventional occlusion time and the subdivision of the simulation into pre-and post-clogging, proved to be effective.

•
The two-phase approach for debris flow present in the TRENT2D model allows an effective and straightforward description of all the phenomena involving granular debris flows in challenging situations. • The WEEZARD system allowed us to handle complex and sophisticated simulations simply, achieving a methodology in line with the European Floods Directive.

•
An additional important result should be emphasized. As the back-analysis was implemented mostly using a priori fixed parameters (as they are physically based) and the parameters determined by calibration (Y and consequently β) are very close to values that can be reasonably a priori estimated, the modelling approach implemented also has a good predictive capability. Indeed, as the observed local mixture depth values range along the stream between 0.30 m and 1.50 m, and assuming an average diameter of the transported material of 0.20 m, the Y parameter can reasonably be in the range 1.5-7.5. This results in an average value of 4.5, which very close to the value of five used in the modelling.
Therefore, we are confident that the implemented modelling approach may be reasonably employed for predictive analysis, such as for hazard mapping and hydraulic design, as well as verification of mitigation structures, and that it may be considered in accordance with the EU Flood Directive as an "appropriate practice and the best available technology that does not imply excessive costs". Funding: This research was partially funded by CARITRO Foundation-Cassa di Risparmio di Trento e Rovereto within the project "Progetto WEEZARD: un sistema integrato di modellazione matematica a servizio della sicurezza nei confronti di pericoli idrogeologici in ambiente montano"; and it was also funded by Dipartimento di Scienze Fisiche, della Terra e dell'Ambiente within the project "Accordo di collaborazione scientifica tra Unione Comuni Amiata-Val d'Orcia e Università di Siena" (Project leader: L.D.), under grant agreement n. 6100. Data Availability Statement: All data generated or analysed for this study are included in this article, and the related datasets are available from the corresponding author upon reasonable request.
Referring to Figure A1, where the reference system, the variables used, and the infinitesimal control volume (in the plane) used for the derivation of the equations are reported, the system of differential equations of the TRENT2D model, considering the assumptions listed above, turns out to be the following: The first two equations express the conservation of the mass of the mixture (sum of the mass conservation equation of both phases) and the solid mass, respectively. Note that the densities of the phases are simplified in the equations as constants. The second two equations express the conservation of momentum of the mixture in the two directions of the plane divided by the density of the liquid phase ρ w . They are derived from the sum of the conservation of momentum of each phase and therefore the interphase forces, being equal and opposite in the two equations, elide and disappear from the equations. The isokinetic hypothesis also allows the momentum of the mixture, the sum of the momenta of the two phases, to be expressed as: where ∆ = (ρ s − ρ w )/ρ w . Finally, note that the term ρ w (1 + c∆)gh 2 /2 represents the thrust of the mixture on the vertical surfaces of the control volume (term P in Figure A1) while ρ w (1 + c∆)gh∂z b /∂x and ρ w (1 + c∆)gh∂z b /∂y, where z b is the height of the bed with respect to a horizontal reference plane, and expresses the thrusts of the bed on the control volume (term p b in Figure A1) in the x and y directions, respectively.
The first two equations express the conservation of the mass of the mixture (sum of the mass conservation equation of both phases) and the solid mass, respectively. Note that the densities of the phases are simplified in the equations as constants. The second two equations express the conservation of momentum of the mixture in the two directions of the plane divided by the density of the liquid phase . They are derived from the sum of the conservation of momentum of each phase and therefore the interphase forces, being equal and opposite in the two equations, elide and disappear from the equations. The isokinetic hypothesis also allows the momentum of the mixture, the sum of the momenta of the two phases, to be expressed as: (1 − ) ⃗ℎ + ⃗ℎ = (1 + Δ) ⃗ℎ (A2) where Δ = ( − )/ . Finally, note that the term (1 + Δ) ℎ /2 represents the thrust of the mixture on the vertical surfaces of the control volume (term P in Figure A1) while (1 + Δ) ℎ / and (1 + Δ) ℎ / , where is the height of the bed with respect to a horizontal reference plane, and expresses the thrusts of the bed on the control volume (term in Figure A1) in the x and y directions, respectively. As for the expression of bed tangential stress and concentration as a function of local hydrodynamic variables, they are expressed by Equations (1) and (2), respectively, in the case of an erodible zone. When the thickness of the erodible material is null, the zone is non-erodible and the bed tangential stress is expressed by Equation (3), while the concentration becomes advected by the flow and does not require a closure relation. More details regarding the derivation of the equations can be found in [43,60,78]. As for the expression of bed tangential stress and concentration as a function of local hydrodynamic variables, they are expressed by Equations (1) and (2), respectively, in the case of an erodible zone. When the thickness of the erodible material is null, the zone is non-erodible and the bed tangential stress is expressed by Equation (3), while the concentration becomes advected by the flow and does not require a closure relation. More details regarding the derivation of the equations can be found in [43,60,78].