Articles | Volume 20, issue 7
Research article
03 Jul 2020
Research article |  | 03 Jul 2020

Sensitivity and identifiability of rheological parameters in debris flow modeling

Gerardo Zegers, Pablo A. Mendoza, Alex Garces, and Santiago Montserrat

Over the past decades, several numerical models have been developed to understand, simulate and predict debris flow events. Typically, these models simplify the complex interactions between water and solids using a single-phase approach and different rheological models to represent flow resistance. In this study, we perform a sensitivity analysis on the parameters of a debris flow numerical model (FLO-2D) for a suite of relevant variables (i.e., maximum flood area, maximum flow velocity, maximum height and deposit volume). Our aims are to (i) examine the degree of model overparameterization and (ii) assess the effectiveness of observational constraints to improve parameter identifiability. We use the Distributed Evaluation of Local Sensitivity Analysis (DELSA) method, which is a hybrid local–global technique. Specifically, we analyze two creeks in northern Chile (∼29 S, 70 W) that were affected by debris flows on 25 March 2015. Our results show that SD (surface detention) and β1 (a parameter related to viscosity) provide the largest sensitivities. Further, our results demonstrate that equifinality is present in FLO-2D and that the final deposited volume and maximum flood area contain considerable information to identify model parameters.

1 Introduction

In steep mountain environments, intense and localized storms can trigger the sudden movement of sediments, generating flash floods with solid volumetric concentrations up to 40 %–60 % (Takahashi1981; O'Brien and Julien1988; Calvo and Savi2009). These events, also known as debris flows, differ from water floods because – in addition to fluid stress – solid–fluid and solid–solid interactions dominate the flow motion (Takahashi1981; Iverson et al.1997). In recent years, debris flows have been recognized as a major natural hazard (Calvo and Savi2009), affecting infrastructure, economic activities and human life. For instance, debris flow events in Switzerland produced 24 fatalities and overall losses of USD 380 million during the period 1972–2007 (Hilker et al.2009). In Chile, estimated economic losses associated with the five biggest debris flow events recorded over 1980–2017 were at least USD 1.6 billion with nearly 1000 people dead or missing (Servicio Nacional de Geología y Minería2017).

Over the last decades, numerical models have emerged as a powerful tool to understand the behavior and magnitude of debris flow events, since they allow for the quantification of key variables used by engineers and decision-makers for risk management (Quan Luna et al.2011; Frey et al.2016; Calvo and Savi2009) and urban planning (Hürlimann et al.2006; Lucà et al.2014; Naef et al.2006; Arattano et al.2006). However, the application of debris flow models requires several assumptions and simplifications that make results diverge from reality at various levels (Sosio et al.2007). For example, uncertainties in terrain elevation models (e.g., satellite product and horizontal resolution), physical parameters (e.g., rheology parameters and solid concentration) and hydrological fluxes (e.g., precipitation and streamflow) used to force debris flow simulations can substantially impact relevant variables, such as flood area, sediment volumes or maximum flow depth.

The use of debris flow models for practical problems typically requires the implementation of single-phase numerical models (Rickenmann et al.2006; Naef et al.2006) that solve one-dimensional or two-dimensional Saint Venant equations, using different rheological approaches to account for frictional stress Sf. Many studies have reported good agreement between debris flow model results and post-event measurements, e.g., runout distance, flow velocity, deposit depth and flood area (D'Agostino and Tecca2006; Naef et al.2006; Sosio et al.2007; Rickenmann et al.2006; Cesca and D'Agostino2008; Lin et al.2011; Hungr1995). Nevertheless, it is recognized that, because complex debris flow dynamics change in time and space (Coussot and Meunier1996), the appropriate choice of rheological parameters is critical for a good agreement between debris flow model output and field data (Sosio et al.2007). In this context, various approaches have been adopted to characterize the sensitivity of debris flow model results to variations in model parameters, i.e., the coefficients in the model equations. For example, D'Agostino and Tecca (2006) compared FLO-2D (O'Brien and Garcia2009) simulations performed with two sets of rheological parameters and three values of the laminar coefficient K (six simulations in total), concluding that K controls the flood area and that rheological parameters control the maximum depth. Boniello et al. (2010) compared FLO-2D model results from a set of 12 back-calculated rheological parameters selected from previous studies, with another set of parameter values obtained from laboratory rheological analyses, finding a better representation of debris flow behavior with back-calculated parameters. Chow et al. (2018) conducted simulations with FLO-2D using 26 different sets of rheological parameters obtained from previous studies, combined with different values of volumetric sediment concentration Cv, specific gravity Gs and the surface detention SD, a parameter used in FLO-2D to represent flow detention. They found that the most important parameters were Cv, SD and β1, which characterizes fluid viscosity. All these studies used fixed sets of rheological parameters in their numerical experiments, and, therefore, the relative importance of such coefficients on relevant simulated variables – specifically, flow depth, flow velocity, deposit volume and flood area – remains unknown. Therefore, this paper addresses the following questions:

  1. How sensitive are debris flow model results to uncertain – and typically fixed – rheological parameters?

  2. What are the most effective post-event measurements to constrain the parameter search towards more realistic simulations?

To answer these questions, we perform a sensitivity analysis on the parameters of a numerical debris flow model and examine the effects of using post-event in situ measurements on expected parameter ranges. In particular, we analyze a debris flow event that occurred in two creeks located in the Atacama region (northern Chile; ∼29 S, 70 W) during March 2015. This event was the consequence of heavy precipitation over a 3 d period, which exceeded 60 mm at several locations (Bozkurt et al.2016), resulting in the loss of human lives and massive infrastructure damage. Therefore, our intention is to provide guidance on the choice of uncertain rheological parameters, contributing to more reliable numerical simulations for debris flow risk assessments and land use planning. The remainder of this paper is organized as follows: Sect. 2 describes the case study creeks and data; Sect. 3 describes the numerical debris flow model, sensitivity analysis and parameter search strategies; Sect. 4 presents the results and discussion; and Sect. 5 summarizes our main findings.

2 Study domain and data

We choose two ephemeral creeks located nearby in the upper Huasco River basin, Acerillas and La Mesilla (Fig. 1a), where debris flows were triggered by an extreme precipitation event on 24–26 March 2015 (Bozkurt et al.2016; Ortega et al.2019). The Huasco Valley is a semi-arid fluvial system located at the southern edge of the Atacama region, Chile. This valley is characterized by perennial rivers that only exist in the trunk valleys, while tributaries only show ephemeral streams. In these areas, heavy rainfall events may induce catastrophic debris flows and mud floods that greatly contribute to erosion (Aguilar et al.2020).

Figure 1(a) Location of the two case study creeks and reference model results. The maximum observed flood areas and modeled flow depth (reference models) are shown for (b) Acerillas creek and (c) La Mesilla creek. Elevations bands created from a satellite digital elevation model (DEM) at 12×12 m: © JAXA/METI ALOS PALSAR L1.0 2007 (Japan Aerospace Exploration Agency/Ministry of Economy, Trade and Industry Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar). Accessed through the ASF DAAC (Alaska Satellite Facility Distributed Active Archive Center) on 11 June 2017.

The Acerillas creek (15 km2 basin area; Fig. 1b) has a markedly narrow channel with almost no alluvial fan, allowing for the transportation of sediments towards the El Carmen River. Post-event measurements indicate a deposited sediment volume of 6000 m3 (Cabré et al.2020) and a maximum flood area of 37 000 m2. Conversely, the La Mesilla creek (2.5 km2 basin area; Fig. 1c) is characterized by a big alluvial fan where considerable sedimentation occurs, and post-event measurements show a deposited sediment volume of 102 000 m3 (Cabré et al.2020) and a maximum flood area of 246 500 m2. These flood areas were estimated by comparing pre- and post-event satellite Google Earth imagery. Also, a post-event topography lidar scan (acquired in February–March 2017 by the Chilean Ministry of Public Works) was available for this study. This dataset has a 1×1 m2 horizontal resolution and was post-processed in order to eliminate vegetation and buildings.

Flow discharge data at the outlet of each creek were obtained from a distributed hydrological model, HEC-HMS (Hydrologic Engineering Center Hydrologic Modeling System) version 4.2 (USACE2015). The model was configured for the entire Huasco River basin (7242 km2), upstream of the Santa Juana irrigation reservoir, as part of a debris flow mitigation project for the Chilean Ministry of Public Works. Hydrologic model simulations were forced using data from point measurements at 14 meteorological stations, spatially distributed with the inverse-distance-weighting (IDW) interpolation method (Teegavarapu et al.2006). Total rainfall records range from 20 to 76 mm, with a maximum registered intensity of 16 mm h−1. The HEC-HMS model parameters were calibrated against hourly streamflow observed at two gauge stations located in the upper part of the basin – Río El Carmen en El Corral and Río Conay en Las Lozas (Fig. 2) – obtaining a Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe1970) of 0.78 for the former and 0.64 for the latter. Although all the other gauging stations were buried or destroyed by debris flows, simulated total water volumes were similar with those captured by the Santa Juana reservoir, whose levels were low before the event. Estimated peak flow discharges for the event analyzed are 8 m3 s−1 at La Mesilla and 12.7 m3 s−1 at Acerillas.

Figure 2Calibration records for stream gauge stations (a) Río El Carmen en El Corral and (b) Río Conay en Las Lozas.


3 Methods

3.1 Debris flow model

We use the two-dimensional FLO-2D debris flow model (O'Brien et al.1993), configured at a 10 m horizontal resolution. FLO-2D is a finite difference model that simulates water or debris flows in channels or unconfined surfaces. The governing equations solved by FLO-2D are the depth-averaged continuity and momentum conservation (Eqs. 1 and 2), and the flood wave progression is controlled by topography and flow resistance (O'Brien and Julien1988). FLO-2D can also simulate debris flow rheologies using a “quadratic” rheological model that combines components associated with creep, viscous, dispersive (collisions) and turbulent stresses (O'Brien and Julien1988; O'Brien and Garcia2009; Naef et al.2006). Based on this quadratic rheology, the friction slope Sf is estimated as (Eq. 3):


where h is the local flow depth, t is time, Vx and Vy are depth-averaged velocity components along the x and y directions, g is gravitational acceleration, Sf is the friction slope, So is the bed slope, τy is the yield stress, K is a laminar strength parameter, η is the interstitial fluid dynamic viscosity, and ntd is the conventional Manning's roughness coefficient corrected by Cv (ntd=0.0538ne6.0896Cv; O'Brien and Julien1988). O'Brien et al. (1993) proposed the following empirical relationships to calculate the viscosity and yield stress as a function of the volumetric sediment concentration cv (Eqs. 4 and 5):


where α1,2 and β1,2 are experimentally defined empirical coefficients (O'Brien et al.1993; O'Brien and Garcia2009).

3.2 Sensitivity analysis

Sensitivity analysis (SA) is a powerful tool to characterize the effects of variations in input factors on environmental-model responses (Razavi and Gupta2015; Gupta and Razavi2018). When the factors of interest are the model parameters, SA helps to identify those that are redundant for the modeling purposes, contributing to a more efficient parameter search (Mendoza et al.2015). Different types of SA techniques have been proposed in the literature depending on specific objectives or even the meaning of sensitivity (see, for example, reviews by Razavi and Gupta2015, and Pianosi et al.2016).

In this work, we apply the Distributed Evaluation of Local Sensitivity Analysis (DELSA) method (Rakovec et al.2014), which is a frugal local–global hybrid technique, to identify the parameters that have the largest impact on simulated debris flow variables. Although our implementation only examines first-order sensitivities across the parameter space – as in Rakovec et al. (2014) – it should be noted that DELSA has considerable unexplored potential to characterize parameter interactions, which could be achieved by including additional terms in the prediction total variance, as suggested by Sobol and Kucherenko (2010).

First-order sensitivities are obtained using local gradients that quantify the sensitivity of a modeled output Ψ relative to individual variations of a parameter θj. The local gradients Ψθj|k are used to compute the first-order sensitivity of each parameter j at each point k of the parameter space:

(6) S 1 k j = | Ψ k j θ j | 2 1 12 θ i , max - θ i , min 2 V K ( Ψ ) ,

where VK(Ψ) is the total local variance at point k

(7) V K ( Ψ ) = j = 0 J | Ψ k j θ j | 2 1 12 θ i , max - θ i , min 2 .

The first-order sensitivity measures S1kj vary between 0 and 1, and the sum of first-order sensitivities from all parameters is equal to 1 at each sample point. In this work, parameter sampling is performed using the Latin hypercube sampling (LHS) method. LHS is a statistical method to generate an almost random sample of parameter values from a multidimensional distribution and has proven to be more efficient than other methods like Monte Carlo sampling (Olsson et al.2003; Olsson and Sandberg2002).

Local sensitivities can be analyzed in a disaggregated manner – through their cumulative frequency distribution across the parameter space – or aggregated by computing a specific statistical property, e.g., the median of all local sensitivity measures for a particular pair of parameter and target variables. We use both approaches to analyze SA results.

In this paper, we focus on the effects of debris flow model parameters on four response variables: maximum average runoff speed Vmean (m s−1), maximum average runoff height Hmean (m), maximum flood area Amax (m2) and deposited volume Voldep (m3). These response variables are calculated using the outputs from FLO-2D as



where j is the cell index; NWC is the total number of wet cells (h>0); dx and dy indicate the cell size along the x and y axis (numerical grid); and Vmax(j), Hmax(j) and Hfinal(j) are the maximum flow speed, maximum flow depth and final runoff height of the cell j, respectively.

The FLO-2D parameters considered for DELSA are those that describe the fluid rheology (Table 1): α1,2, β1,2, Cv, K, n, the SD parameter and the total volume of sediments mobilized VolT.

The detention coefficient SD is a model parameter that controls flow detention. The FLO-2D user's manual and previous studies (D'Agostino and Tecca2006) suggest that SD acts as the minimum physically plausible flow depth (i.e., flow stops if flow depth < SD). Further, D'Agostino and Tecca (2006) noted that this coefficient has a strong influence on the results, and it can be used as a surrogate of the rheology. The estimated total sediment volume mobilized VolTest by each debris flow event was calculated with the equation proposed by Chang et al. (2011) (Eq. 12).

(12) Vol Test = 0.023 A w + 0.064 A L + 13264.6 GI - 1399.2 D + 38.47 C R ,

where Aw is the watershed area, AL is the landslide area (zero in these cases), GI is the geological index (where a value of 2.5 is assumed based on a study zone report made for the Chilean Ministry of Public Works), D is the rainfall duration (D=48 h for this event) and CR is the cumulative rainfall (CR=76 mm). We obtain VolTest values of 185 000 m3 for Acerillas and 154 000 m3 for La Mesilla.

Debris flow concentration is assumed to vary with streamflow between a minimum concentration Cvmin=0.1–0.4 and a maximum concentration Cvmax at the time of peak flow, which is treated as a model parameter. To this end, we propose the following function for Cv:

(13) C v ( t ) = C v max - C v min erf Q ( t ) - Q m / Q max - Q m ϕ C v max - C v min C v max - C v min erf Q max - Q m / Q max - Q m ϕ + C v min , if Q ( t ) / Q mean 0.5 C v min , if Q ( t ) / Q mean < 0.5 ,

where ϕ is a coefficient that changes the shape of the concentration curve. ϕ and Cvmin are calculated in order to match the total volume and minimize Cvmin value:

(14) minimize C v min subject to t = 0 T C v ( t ) Q ( t ) d t = Vol T .

Rheological parameter ranges are obtained from previous debris flow studies (O'Brien and Julien1988; Boniello et al.2010; D'Agostino and Tecca2006; Sosio et al.2007; Rickenmann et al.2006; Chang et al.2011; O'Brien and Garcia2009) and are summarized in Table 1. However, additional restrictions are imposed for τy and η, with maximum values of 35 000 dyn cm−2 (dyne; g cm s−2) and 100 000 P (poise; 0.1 kg m−1 s−1), respectively (Rickenmann et al.2006). Since τy and η are functions of rheological parameters (Eqs. 5 and 4), such limits impose restrictions for α1,2, β1,2 and Cvmax that we ensure are implemented in DELSA. To this end, we develop a Python script that allows for running FLO-2D in parallel and sequentially, reducing computational cost considerably. For VolT, we assume a range of variation of ±30 % with respect to estimated values (Eq. 12).

Sosio et al. (2007)O'Brien and Julien (1988)O'Brien and Julien (1988)O'Brien and Julien (1988)D'Agostino and Tecca (2006)O'Brien and Julien (1988)Sosio et al. (2007)O'Brien and Julien (1988)Rickenmann et al. (2006)O'Brien and Garcia (2009)O'Brien and Garcia (2009)D'Agostino and Tecca (2006)Chang et al. (2011)O'Brien and Julien (1988)Rickenmann et al. (2006)O'Brien and Julien (1988)Rickenmann et al. (2006)

Table 1Values range of the model parameters. P: poise; dyn: dyne.

Download Print Version | Download XLSX

3.3 Parameter selection via constrained search

We explore the effects of parameter uncertainty on simulated debris flow variables at the two case study creeks. We also examine the utility of using reference values for specific variables to constrain the search of physically plausible parameter sets. Such values are obtained from a reference simulation conducted by Zegers (2017), who reproduced flood area and sediment volume observed during the 2015 debris flow events at Acerillas and La Mesilla creeks using FLO-2D. Zegers (2017) calibrated model parameters by contrasting results against measured flood areas, deposited volumes and flow velocity estimated from a video captured with a cell phone by a local person at Acerillas. Such a validation strategy has provided reliable results for several other creeks in the area. Parameter values used for the reference simulation are provided in Table 2. Modeled flood areas and maximum flow depth are shown in Fig. 1, and discrepancies with respect to observations can be attributed to the use of post-event topography.

Based on the reference simulation, we choose reference values of Vmean=1 m s−1 and Hmean=1.5 m at both creeks. Additionally, we estimate the reference maximum flood area using Google Earth imagery – obtaining values of 246 500 m2 for La Mesilla and 37 000 m2 for Acerillas – and use deposited sediment volumes reported by Cabré et al. (2020) as reference values, which correspond to 102 000 and 6000 m3 for La Mesilla and Acerillas, respectively.

Table 2Parameter values for reference models on Zegers (2017). P: poise; dyn: dyne.

Download Print Version | Download XLSX

4 Results and discussion

4.1 Choice of sample size

First, we seek to identify the minimum sample size Nk for which stable DELSA results can be obtained in order to minimize computational cost (Rakovec et al.2014). Therefore, we explore the effects of the choice of Nk on the cumulative frequency distributions (CDFs) of DELSA first-order sensitivity indices. Since we include Nj=9 parameters, the total number of simulations required for each case is Nt=(Nj+1)Nk. Figure 3 illustrates the sensitivity of DELSA results to variations in sample size Nk for four variables simulated by FLO-2D, with respect to parameters β1 and Cvmax, at the Acerillas creek. Since the curves obtained for Nk=500 and Nk=1000 show slight differences, departing from the CDF for Nk=100, we conclude that a sample size of Nk=500 is adequate for further analyses. The sensitivity of DELSA results to Nk was also examined at La Mesilla creek, obtaining the same conclusion regarding sample size (not shown).

Figure 3Effects of sample size Nk on the cumulative distribution function (CDF) of DELSA indices for parameters β1 (a) and Cvmax (b) at the Acerillas creek.


Figure 3 also shows that, depending on the target variable and parameter analyzed, first-order sensitivity indices can be highly heterogeneous across the parameter space. In particular, the modeled response is highly sensitive to variations of β1, with first-order sensitivities larger than 0.2 for approximately 60 % of cases and sensitivities greater than 0.5 for 20 %–40 % depending on the variable analyzed. On the other hand, the modeled variables are less sensitive to Cvmax in most cases, with DELSA indices smaller than 0.1 for approximately 70 % of the parameter sets.

4.2 Sensitivity of model responses to model parameters

Figure 4 displays the median of the full frequency distribution (obtained with Nk=500) of local first-order sensitivity indices for the two study domains: (i) Acerillas creek (ACE) and (ii) La Mesilla creek (MES). The uncertainty bands are obtained by performing bootstrapping with replacement (resampled 1000 times). In general, β1 provides the largest sensitivities for the simulated variables analyzed, which means that the fluid rheology, in particular the viscosity coefficient (η=η(β1)), is a main parameter controlling flow behavior. Moreover, the results of Acerillas seem to be more sensitive than those obtained at La Mesilla with respect to β1. This could be better explained when analyzed together with SD, the detention coefficient.

Figure 4DELSA sensitivity indices (synthetized as the median from the cumulative frequency distribution) for all parameters and model responses. Results are displayed for Acerillas (red) and La Mesilla (black) creeks, and the sampling uncertainty (bootstrapping with 1000 times resampling) is indicated by boxplots. The vertical bold line in the boxplot is the median; the body of each boxplot shows the interquartile range (Q75Q25); and the whiskers represent the sample minima and sample maxima. DELSA indices are displayed in log space for a better visualization of inter-parameter differences.


As expected, the simulated deposited volume Voldep is very sensitive to SD because this parameter controls flow detention. For the remaining simulated variables, SD also rises as an important parameter at La Mesilla, but it shows secondary importance at Acerillas, which can be explained by catchment differences. While the fluid rheology explains flow behavior at Acerillas creek (mainly sensitive to β1), depositional or detention processes – represented by SD – gain importance across the larger alluvial fan of La Mesilla creek.

Another parameter that provides large sensitivities – especially in simulated mean flow velocity Vmean – is the total sediment volume VolT, whose sensitivities have the same order of magnitude as those produced by β1. The high sensitivity of VolT on Vmean is explained because the former influences fluid rheology through Cv(t) (see Eqs. 4, 5 and 13). However, VolT does not produce large variations in the total deposited volume, which could be explained because, as deposition occurs, the flow is channelized between the deposited margins, preventing flow spreading. For example, when increasing SD values, Amax decreases as the flow is forced to stop at deeper heights. This could be a structural weakness of FLO-2D, which lacks proper representations of complex depositional and dewatering processes.

The large sensitivities in model response to variations of β1 suggest that the viscous stress (second term in Eq. 3) is the main contributor to sensitivities in simulated frictional slope. On the other hand, DELSA sensitivity indices associated with yield stress τy and Manning's roughness coefficient – the other components of the frictional slope – are of second-order importance. Even more, model results are practically insensitive to Manning's roughness coefficient.

Our results show some differences with related previous work. D'Agostino and Tecca (2006) compared FLO-2D simulations performed with two sets of rheological parameters and three values of K (i.e., six simulations), concluding that the latter controls the flood area, while rheological parameters control the maximum flow depth. On the other hand, our results indicate that the laminar coefficient K provides small sensitivities in simulated flood areas and that the most influential rheological parameter on simulated maximum flow depth is β1. Although D'Agostino and Tecca (2006) provided the first insights into FLO-2D parameter sensitivities, the number of parameters involved and the sample size were not large enough to draw robust conclusions. Recently, Chow et al. (2018) conducted simulations with FLO-2D using 26 different sets of rheological parameters obtained from 45 previous studies. They found that the most influential parameters – in order of importance – were Cv, SD and β1. These results are different from ours due to discrepancies in the parameter sampling method, the use of fixed sets of rheological parameters and their parameter ranking definition, which does not consider separate effects on key simulated variables. For example, Chow et al. (2018) ignored the effect of the total sediment volume when analyzing Cv, whereas we examine the separate effects of maximum sediment concentration and total sediment volume, obtaining that the latter is more important than the maximum sediment concentration.

4.3 Parameter uncertainty effects

Figure 5 shows the full range of variation in model responses (produced by Nk=500 parameter sets), normalized by reference values obtained from post-event measurements (Voldep and Amax) and results obtained from the reference models (for Hmean and Vmean). These are compared with the simulated ensemble that results from screening model outputs imposing five different observational constraints: (i) ±20 % reference mean flow velocity (FVEL), (ii) ±20 % reference mean flow depth (FH), (iii) ±20 % reference maximum flood area (FAREA), (iv) ±40 % reference volume deposit (FVOL), and (v) ±20 % reference maximum flood area and ±40 % reference volume deposit (FAREAVOL). To be clear, constraint (i) results from keeping all those parameter sets that provide a simulated mean flow velocity within the range 0.8–1.2 Vref. Constraints (ii)–(iv) work in a similar way for other observed variables, while constraint (v) filters all parameter sets that simultaneously provide flood areas and deposit volumes within the ranges 0.8–1.2Aref and 0.6–1.4 Volref, respectively. We assume a weaker observational constraint in the case of the deposited volume because of possible uncertainties associated with different measurements techniques.

Figure 5Effects of parametric uncertainty in normalized model responses for the original parameter sample (top panels) and alternative observational constraints – flow velocity (FVEL), flow depth (FH), flood area (FAREA) and sediment volume (FVOL) – and joint area–volume constraint (FAREAVOL). Results are displayed for Acerillas (black) and La Mesilla (red) creeks. The hatched area for Vmean, Hmean and Amax corresponds to ±20 % of their reference values, and for Voldep it is up to ±40 % of the estimated sediment volume. The vertical bold line in the boxplots is the median; the body of each boxplot shows the interquartile range (Q75Q25); and the whiskers represent the sample minima and sample maxima.


Figure 5 also shows that, in general, the effects of parametric uncertainty on simulated variables are considerable and the ensemble median can be substantially different from the reference boundaries. This is somewhat expected, since the literature provides large ranges for model parameters (see Table 1 for details). Overall, uncertainties arising from the original parameter samples (top panels) are larger in Acerillas, except for mean flow velocity. This could be explained by the larger sensitivity of flow velocity to flow rheology, while the rest of simulated variables are more sensitive to deposition processes, mainly represented by SD.

Most simulations overestimate deposited volumes and flow depth, especially at Acerillas. For flow velocity, the ensemble of parameter sets provides mixed results in both creeks, with underestimation in most cases (median values lower than the reference values); however, there are still several parameter sets that produce an overestimation of flow velocity. The results obtained for maximum area reveal differences among both creeks: in Acerillas, most parameter sets tend to overestimate the flood area, whereas most simulated values are within the expected range at La Mesilla. This could happen because Voldep VolT at Acerillas, while Voldep VolT at La Mesilla; moreover, the maximum flood area at La Mesilla is approximately 6 times larger than in Acerillas. Thus, small variations in Voldep imply important fractional changes with respect to the volume reference values at Acerillas.

Figure 5 shows that the velocity constraint FVEL does not have an impact on the rest of simulated variables. Nevertheless, the application of alternative observational constraints helps to reduce the spread of the remaining variables. For example, the height filter FH improves simulations of maximum flood area, although it does not have much effect on velocity or deposit volume. The area restriction FAREA improves simulated flow depth only at Acerillas, as most of the original ensemble members were already inside the expected reference boundaries at La Mesilla. The volume constraint FVOL reduces the uncertainty in all variables, with the smallest improvement for flow velocity. This is because the volume is directly linked to flow height and flood area, but not to flow velocity. Finally, the largest reductions in ensemble spread are obtained when parameter sets are constrained by using area and volume observations (FAREAVOL).

The maximum flood area and deposited volume are relatively easy to measure and are probably the most used post-event measurements for calibrating debris flow models (Chow et al.2018; Cesca and D'Agostino2008; D'Agostino and Tecca2006; Sosio et al.2007; Frey et al.2016; Quan Luna et al.2011).

4.4 Parameter identifiability

Figure 6 illustrates the effects of applying observational constraints, specifically flood area and deposited volume constraints, on parameter identifiability. Results show that the resampled values of α1,2, β2, Cvmax, VolT, K and n cover practically the entire original range. However, the application of observational restrictions provides substantial reductions in the ranges of β1 and SD, which mainly explain the flow rheology and depositional processes in our study areas. Further, lower parameter values are obtained in comparison to the full range, especially SD at the Acerillas creek. These results indicate that viscosity η=η(α1, β1, Cvmax) is the most restricted parameter when applying these constraints, discarding all medium–high values. On the other hand, resampled values of τy=τy(α2, β2, Cvmax) cover almost the entire original range.

Figure 6Effects of applying a flood area–volume observational constraint (FAREAVOL) on parameter identifiability. Results are displayed for Acerillas (black) and La Mesilla (red) creeks. The vertical bold line in the boxplots is the median; the body of each boxplot shows the interquartile range (Q75Q25); and the whiskers represent the sample minima and sample maxima. The grey–black diamonds represent parameter values for the reference simulations.


The reference SD value is close to the upper range in La Mesilla after applying the FAREAVOL constraint, and also much larger than the resulting maximum values filtered at Acerillas. This result is somewhat expected, since SD does not provide large model sensitivities in that domain. Similarly, the reference β1 value is within the upper range of filtered values (>Q75). However, filtered η values around the baseline model parameter result from the compensation of low values of α1 (near the minimum in the filtered range) and Cvmax (below the median of filtered values). A similar effect is observed for τy, whose filtered range results from the compensation of α2, β2 and Cvmax. In summary, different combinations of α1,2, β1,2 and Cv can generate viscosity and yield stress values that are suitable to reproduce the 2015 debris flow events in Acerillas and La Mesilla. This is a well-known problem in environmental models – referred to as the equifinality, nonuniqueness or nonidentifiability of model parameters – that has been widely discussed for more than 3 decades in the hydrology literature (e.g., Beven2006; Kelleher et al.2015) but not carefully addressed in the debris flow modeling community.

Figure 6 also shows different behavior in other parameters. For example, the reference values for K and n are in the lower range of the filtered ensembles, while the reference Vol T is in the upper body of the boxplot. Low K values produce low Sf2, the second term at the right hand of Eq. (3), representing viscous stress. This could be compensated for using larger values of SD, as in the reference model. These results demonstrate that equifinality in FLO-2D does not only involve rheological parameters and that SD could be an important parameter to correct the unrealistic model representation of rheology (D'Agostino and Tecca2006).

The main goal of this study was to characterize the sensitivity of model responses to variations in uncertain rheological parameters, using only independent information on parameter values (i.e., the situation comparable to Sobol, as proposed by Rakovec et al.2014). Hence, the only verification dataset available (flood area and sediment volume from the March 2015 event) was used to examine the identifiability of model parameters. However, the relative importance of additional observations could be assessed through the Observation-Prediction (OPR) statistic (Tiedeman and Hsieh2004), and the potential new information provided by field data (e.g., sedimentological and morphological characteristics) for a specific parameter (e.g., α1 and β1) could be quantified with the Parameter-Prediction (PPR) statistic (Tonkin et al.2007). It should be noted that, in both cases, the equation for the total local variance (Eq. 7) would be different, as additional information should be incorporated (see Appendix A in Rakovec et al.2014).

5 Summary and conclusions

We performed a sensitivity analysis on the parameters of a widely used numerical debris flow model (FLO-2D) and assessed the effects of applying observational constraints on parameter identifiability. Our study domains are two morphologically different ravines, Acerillas and La Mesilla, located in the Atacama region (∼29 S, 70 W), Chile. While Acerillas is characterized by a straight and well-defined channel with almost non-alluvial fans, La Mesilla has a big alluvial fan where deposition is prone to occur.

We found that β1 – a parameter used to estimate the fluid mixture viscosity – provides the largest sensitivities in the variables analyzed, followed by SD – a model parameter used to represent flow detention. Interestingly, the relative importance of β1 and SD depends on the study site, with the former being more important for Acerillas and the latter being more important for La Mesilla. These results suggest that, while rheological processes dominate flow behavior at Acerillas (straight channel with small alluvial fan), sedimentation and detention processes control flow in La Mesilla (big alluvial fan). Although the total mobilized sediment VolT does not have an effect on Voldep, it is important for representing flow velocity, as VolT is used to estimate Cv(t), a key parameter for fluid rheology. Finally, model results seem to be almost insensitive to Manning's roughness coefficient n, while DELSA sensitivities for the remaining parameters are of second-order importance and provide similar indices.

The comparison between the original model parameter ranges (N=500) and the ensemble resulting from applying observational restrictions shows that SD and β1 (i.e., η) are the parameters whose identifiability is mostly improved, while others practically preserve their original range. In addition, we obtain that different combinations of model parameters (including those that describe rheology) can provide very similar results, indicating that equifinality is present in FLO-2D. Our results also support the idea that single-phase rheological models lack a strong physical basis (Iverson2003), and, therefore, their determination requires expert knowledge. However, an encouraging finding is that the final deposited volume Voldep and maximum flood area Amax contain considerable information for identifying model parameters.

We obtain that SD strongly affects model results at La Mesilla, having also large effects on simulated deposited volumes at Acerillas. Moreover, this study provides evidence that SD is one of the most important parameters controlling flow behavior and could possibly surrogate rheology in the model (D'Agostino and Tecca2006). One-phase debris flow models still lack robust representations of complex process interactions during flow stopping that produce temporal and spatial changes in fluid rheology. Thus, these rheological changes have been replaced by simpler approaches (e.g., the incorporation of SD).

Future investigations should aim to improve the structure of debris flow models and hence achieve better simulations of deposition and erosion processes, stopping phases, and changing rheologies. Further, the development of computationally frugal methods to improve the understanding of parameter interactions in environmental models emerges as an attractive avenue for future research.

Data availability

Lidar topography and used streamflows are available at (Zegers2020).

Author contributions

GZ, AG, PM and SM were involved in the conceptualization. GZ, PM and SM contributed to the methodology and analysis and drafted the paper. GZ configured the model, conducted simulations, analyzed the results and made the figures. GZ, AG, PM and SM contributed to writing, review and editing.

Competing interests

The authors declare that they have no conflict of interest.


We thank Tomás Gómez and Miguel Lagos for contributing the streamflow data used in this study and Martin Mergili and Mary Hill for their thoughtful comments and suggestions for improving this article. We also thank the support of the Faculty of Physical and Mathematical Sciences, Universidad de Chile, through the Department of Civil Engineering and “Centro de Investigación, Desarrollo e Innovación de Estructuras y Materiales” (IDIEM). Finally, the authors thank the Chilean Ministry of Public Works (MOP-DOH) and the CONICYT/PIA (project no. AFB180004). Pablo A. Mendoza acknowledges the support of a Fondecyt postdoctoral grant (no. 3170079).

Financial support

This research has been supported by the Faculty of Physical and Mathematical Sciences of the Universidad de Chile through the Advanced Mining Technology Center (AMTC).

Review statement

This paper was edited by Thomas Glade and reviewed by Mary Hill and Martin Mergili.


Aguilar, G., Cabré, A., Fredes, V., and Villela, B.: Erosion after an extreme storm event in an arid fluvial system of the southern Atacama Desert: an assessment of the magnitude, return time, and conditioning factors of erosion and debris flow generation, Nat. Hazards Earth Syst. Sci., 20, 1247–1265,, 2020. a

Arattano, M., Franzi, L., and Marchi, L.: Influence of rheology on debris-flow simulation, Nat. Hazards Earth Syst. Sci., 6, 519–528,, 2006. a

Beven, K.: A manifesto for the equifinality thesis, J. Hydrol., 320, 18–36, 2006. a

Boniello, M. A., Calligaris, C., Lapasin, R., and Zini, L.: Rheological investigation and simulation of a debris-flow event in the Fella watershed, Nat. Hazards Earth Syst. Sci., 10, 989–997,, 2010. a, b

Bozkurt, D., Rondanelli, R., Garreaud, R., and Arriagada, A.: Impact of warmer eastern tropical Pacific SST on the March 2015 Atacama floods, Mon. Weather Rev., 144, 4441–4460, 2016. a, b

Cabré, A., Aguilar, G., Mather, A. E., Fredes, V., and Riquelme, R.: Tributary-junction alluvial fan response to an ENSO rainfall event in the El Huasco watershed, northern Chile, Prog. Phys. Geogr.: Earth Environ., in preparation, 2020. a, b, c

Calvo, B. and Savi, F.: A real-world application of Monte Carlo procedure for debris flow risk assessment, Comput. Geosci., 35, 967–977,, 2009. a, b, c

Cesca, M. and D'Agostino, V.: Comparison between FLO-2D and RAMMS in debris-flow modelling: a case study in the Dolomites, in: Monitoring, Simulation, Prevention and Remediation of Dense Debris Flows II, vol. I of WIT Transactions on Engineering Sciences, WIT Press, Southampton, UK, 197–206,, 2008. a, b

Chang, C.-W., Lin, P.-S., and Tsai, C.-L.: Estimation of sediment volume of debris flow caused by extreme rainfall in Taiwan, Eng. Geol., 123, 83–90,, 2011. a, b, c

Chow, C., Ramirez, J., and Keiler, M.: Application of Sensitivity Analysis for Process Model Calibration of Natural Hazards, Geosciences, 8, 218,, 2018. a, b, c, d

Coussot, P. and Meunier, M.: Recognition, classification and mechanical description of debris flows, Earth-Sci. Rev., 40, 209–227,, 1996. a

D'Agostino, V. and Tecca, P. R.: Some considerations on the application of the FLO-2D model for debris flow hazard assessment, in: Monitoring, Simulation, Prevention and Remediation of Dense and Debris Flows, vol. 1 of WIT Transactions on Ecology and the Environment, Vol. 90, WIT Press, Southampton, UK, 159–170,, 2006. a, b, c, d, e, f, g, h, i, j, k, l

Frey, H., Huggel, C., Bühler, Y., Buis, D., Burga, M. D., Choquevilca, W., Fernandez, F., García Hernández, J., Giráldez, C., Loarte, E., Masias, P., Portocarrero, C., Vicuña, L., and Walser, M.: A robust debris-flow and GLOF risk management strategy for a data-scarce catchment in Santa Teresa, Peru, Landslides, 13, 1493–1507,, 2016. a, b

Gupta, H. V. and Razavi, S.: Revisiting the Basis of Sensitivity Analysis for Dynamical Earth System Models, Water Resour. Res., 100, 1628–1633,, 2018. a

Hilker, N., Badoux, A., and Hegg, C.: The Swiss flood and landslide damage database 1972–2007, Nat. Hazards Earth Syst. Sci., 9, 913–925,, 2009. a

Hungr, O.: A model for the runout analysis of rapid flow slides, debris flows, and avalanches, Can. Geotech. J., 32, 610–623, 1995. a

Hürlimann, M., Copons, R., and Altimir, J.: Detailed debris flow hazard assessment in Andorra: A multidisciplinary approach, Geomorphology, 78, 359–372,, 2006. a

Iverson, R. M.: The debris-flow rheology myth, Debris-flow hazards Mitigation: Mechanics, Prediction, and Assessment, Millpress, Rotterdam, 2003. a

Iverson, R. M., Reid, M. E., and LaHusen, R. G.: Debris-flow mobilization from landslides, Annu. Rev. Earth Planet. Sci., 25, 85–138,, 1997. a

Kelleher, C., Wagener, T., and McGlynn, B.: Model-based analysis of the influence of catchment properties on hydrologic partitioning across five mountain headwater subcatchments, Water Resour. Res., 51, 4109–4136, 2015. a

Lin, P.-S., Lee, J., and Chang, C.: An application of the FLO-2D model to debris-flows simulation – A case study of Song-Her district in Taiwan, Ital. J. Eng. Geol. Environ., 3, 947–956,, 2011. a

Lucà, F., D'Ambrosio, D., Robustelli, G., Rongo, R., and Spataro, W.: Integrating geomorphology, statistic and numerical simulations for landslide invasion hazard scenarios mapping: An example in the Sorrento Peninsula (Italy), Comput. Geosci., 67, 163–172,, 2014. a

Mendoza, P. A., Clark, M. P., Barlage, M., Rajagopalan, B., Samaniego, L., Abramowitz, G., and Gupta, H.: Are we unnecessarily constraining the agility of complex process-based models?, Water Resour. Res., 51, 716–728, 2015. a

Naef, D., Rickenmann, D., Rutschmann, P., and McArdell, B. W.: Comparison of flow resistance relations for debris flows using a one-dimensional finite element simulation model, Nat. Hazards Earth Syst. Sci., 6, 155–165,, 2006. a, b, c, d

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models part I – A discussion of principles, J. Hydrol., 10, 282–290, 1970. a

O'Brien, J. S. and Garcia, R.: FLO-2D Reference manual, available at: (last access: 1 March 2019), 2009. a, b, c, d, e, f

O'Brien, J. S. and Julien, P. Y.: Laboratory analysis of mudflow properties, J. Hydraul. Eng., 114, 877–887, 1988. a, b, c, d, e, f, g, h, i, j, k, l

O'Brien, J. S., Julien, P. Y., and Fullerton, W. T.: Two-dimensional water flood and mudflow simulation, J. Hydraul. Eng., 119, 679–683, 1993. a, b, c

Olsson, A. M. and Sandberg, G. E.: Latin hypercube sampling for stochastic finite element analysis, J. Eng. Mech., 128, 121–125,, 2002. a

Olsson, A. M., Sandberg, G., and Dahlblom, O.: On Latin hypercube sampling for structural reliability analysis, Struct. Safe., 25, 47–68,, 2003. a

Ortega, C., Vargas, G., Rojas, M., Rutllant, J. A., Muñoz, P., Lange, C. B., Pantoja, S., Dezileau, L., and Ortlieb, L.: Extreme ENSO-driven torrential rainfalls at the southern edge of the Atacama Desert during the Late Holocene and their projection into the 21th century, Global Planet. Change, 175, 226–237, 2019. a

Pianosi, F., Beven, K., Freer, J., Hall, J. W., Rougier, J., Stephenson, D. B., and Wagener, T.: Sensitivity analysis of environmental models: A systematic review with practical workflow, Environ. Model. Softw., 79, 214–232,, 2016. a

Quan Luna, B., Blahut, J., van Westen, C. J., Sterlacchini, S., van Asch, T. W. J., and Akbas, S. O.: The application of numerical debris flow modelling for the generation of physical vulnerability curves, Nat. Hazards Earth Syst. Sci., 11, 2047–2060,, 2011. a, b

Rakovec, O., Hill, M. C., Clark, M. P., Weerts, A. H., Teuling, A. J., and Uijlenhoet, R.: Distributed Evaluation of Local Sensitivity Analysis (DELSA), with application to hydrologic models, Water Resour. Res., 50, 409–426,, 2014. a, b, c, d, e

Razavi, S. and Gupta, H. V.: What do we mean by sensitivity analysis? The need for comprehensive characterization of “global” sensitivity in Earth and Environmental systems models, Water Resour. Res., 51, 3070–3092, 2015. a, b

Rickenmann, D., Laigle, D., McArdell, B. W., and Hübl, J.: Comparison of 2D debris-flow simulation models with field events, Comput. Geosci., 10, 241–264,, 2006. a, b, c, d, e, f, g

Servicio Nacional de Geología y Minería: Principales desastres ocurridos desde 1980, Ministerio de Chile, p. 45, available at: (last access: 1 July 2020), 2017. a

Sobol, I. and Kucherenko, S.: Derivative based global sensitivity measures, Procedia – Social Behav. Sci., 2, 7745–7746, 2010. a

Sosio, R., Crosta, G. B., and Frattini, P.: Field observations, rheological testing and numerical modelling of a debris-flow event, Earth Surf. Proc. Land., 32, 290–306,, 2007. a, b, c, d, e, f, g

Takahashi, T.: Debris Flow, Annu. Rev. Fluid Mech., 13, 57–77,, 1981. a, b

Teegavarapu, R., Tuffail, M., and Ormsbee, L.: Optimal function forms for spatial interpolation of precipitation data, Environm. Inform. Arch., 4, 343–353, 2006. a

Tiedeman, C. R. and Hsieh, P. A.: Evaluation of longitudinal dispersivity estimates from simulated forced-and natural-gradient tracer tests in heterogeneous aquifers, Water Resour. Res., 40, W01512,, 2004. a

Tonkin, M. J., Tiedeman, C. R., Ely, D. M., and Hill, M. C.: OPR-PPR, a computer program for assessing data importance to model predictions using linear statistics, Tech. rep., United States Geological Survey – Nevada, Henderson, Nevada, 2007. a

USACE – US Army Corps of Engineers Hydrologic Engineering Center, Hydrologic Modelling System HEC‐HMS, User's Manual, version 4.1, July 2015. a

Zegers, G. A. M. S.: 2D modeling of debris flows. Application to the Quebrada La Mesilla, Huasco Province, in: XXIII Chilean Conference on Hydraulic Engineering, October 2017, Chile, 2017. a, b, c

Zegers, G. A. M. S.: Sensitivity and identifiability of rheological parameters in debrisflow modeling – DATA, figshare, Collection,, 2020. a

Short summary
We perform a sensitivity analysis on the parameters of a numerical debris flow model and examine the effects of using post-event measurements on two creeks in Chile. Our results demonstrate the utility of sensitivity analysis in debris flow modeling and the benefits of post-event observations on parameter identifiability. This study provides guidance on the choice of uncertain parameters, contributing to more reliable simulations for debris flow risk assessments and land use planning.
Final-revised paper