Research article 03 Jul 2020
Research article  03 Jul 2020
Sensitivity and identifiability of rheological parameters in debris flow modeling
 ^{1}Advanced Mining Technology Center, Universidad de Chile, Santiago, Chile
 ^{2}Department of Geoscience, University of Calgary, Calgary, Canada
 ^{3}Department of Civil Engineering, Universidad de Chile, Santiago, Chile
 ^{1}Advanced Mining Technology Center, Universidad de Chile, Santiago, Chile
 ^{2}Department of Geoscience, University of Calgary, Calgary, Canada
 ^{3}Department of Civil Engineering, Universidad de Chile, Santiago, Chile
Correspondence: Santiago Montserrat (santiago.montserrat@amtc.cl)
Hide author detailsCorrespondence: Santiago Montserrat (santiago.montserrat@amtc.cl)
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 singlephase 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 (FLO2D) 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 FLO2D and that the final deposited volume and maximum flood area contain considerable information to identify model parameters.
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 % (Takahashi, 1981; O'Brien and Julien, 1988; Calvo and Savi, 2009). 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 (Takahashi, 1981; Iverson et al., 1997). In recent years, debris flows have been recognized as a major natural hazard (Calvo and Savi, 2009), 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ía, 2017).
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 decisionmakers for risk management (Quan Luna et al., 2011; Frey et al., 2016; Calvo and Savi, 2009) 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 singlephase numerical models (Rickenmann et al., 2006; Naef et al., 2006) that solve onedimensional or twodimensional Saint Venant equations, using different rheological approaches to account for frictional stress S_{f}. Many studies have reported good agreement between debris flow model results and postevent measurements, e.g., runout distance, flow velocity, deposit depth and flood area (D'Agostino and Tecca, 2006; Naef et al., 2006; Sosio et al., 2007; Rickenmann et al., 2006; Cesca and D'Agostino, 2008; Lin et al., 2011; Hungr, 1995). Nevertheless, it is recognized that, because complex debris flow dynamics change in time and space (Coussot and Meunier, 1996), 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 FLO2D (O'Brien and Garcia, 2009) 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 FLO2D model results from a set of 12 backcalculated 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 backcalculated parameters. Chow et al. (2018) conducted simulations with FLO2D using 26 different sets of rheological parameters obtained from previous studies, combined with different values of volumetric sediment concentration C_{v}, specific gravity G_{s} and the surface detention SD, a parameter used in FLO2D to represent flow detention. They found that the most important parameters were C_{v}, 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:

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

What are the most effective postevent 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 postevent 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.
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 semiarid 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).
The Acerillas creek (15 km^{2} 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. Postevent measurements indicate a deposited sediment volume of 6000 m^{3} (Cabré et al., 2020) and a maximum flood area of 37 000 m^{2}. Conversely, the La Mesilla creek (2.5 km^{2} basin area; Fig. 1c) is characterized by a big alluvial fan where considerable sedimentation occurs, and postevent measurements show a deposited sediment volume of 102 000 m^{3} (Cabré et al., 2020) and a maximum flood area of 246 500 m^{2}. These flood areas were estimated by comparing pre and postevent satellite Google Earth imagery. Also, a postevent 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 m^{2} horizontal resolution and was postprocessed in order to eliminate vegetation and buildings.
Flow discharge data at the outlet of each creek were obtained from a distributed hydrological model, HECHMS (Hydrologic Engineering Center Hydrologic Modeling System) version 4.2 (USACE, 2015). The model was configured for the entire Huasco River basin (7242 km^{2}), 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 inversedistanceweighting (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 HECHMS 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 Sutcliffe, 1970) 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 m^{3} s^{−1} at La Mesilla and 12.7 m^{3} s^{−1} at Acerillas.
3.1 Debris flow model
We use the twodimensional FLO2D debris flow model (O'Brien et al., 1993), configured at a 10 m horizontal resolution. FLO2D is a finite difference model that simulates water or debris flows in channels or unconfined surfaces. The governing equations solved by FLO2D are the depthaveraged continuity and momentum conservation (Eqs. 1 and 2), and the flood wave progression is controlled by topography and flow resistance (O'Brien and Julien, 1988). FLO2D 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 Julien, 1988; O'Brien and Garcia, 2009; Naef et al., 2006). Based on this quadratic rheology, the friction slope S_{f} is estimated as (Eq. 3):
where h is the local flow depth, t is time, V_{x} and V_{y} are depthaveraged velocity components along the x and y directions, g is gravitational acceleration, S_{f} is the friction slope, S_{o} is the bed slope, τ_{y} is the yield stress, K is a laminar strength parameter, η is the interstitial fluid dynamic viscosity, and n_{td} is the conventional Manning's roughness coefficient corrected by C_{v} (${n}_{\mathrm{td}}=\mathrm{0.0538}n{e}^{\mathrm{6.0896}{C}_{\mathrm{v}}}$; O'Brien and Julien, 1988). 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 c_{v} (Eqs. 4 and 5):
where α_{1,2} and β_{1,2} are experimentally defined empirical coefficients (O'Brien et al., 1993; O'Brien and Garcia, 2009).
3.2 Sensitivity analysis
Sensitivity analysis (SA) is a powerful tool to characterize the effects of variations in input factors on environmentalmodel responses (Razavi and Gupta, 2015; Gupta and Razavi, 2018). 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 Gupta, 2015, 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 firstorder 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).
Firstorder 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 $\frac{\partial \mathrm{\Psi}}{\partial {\mathit{\theta}}_{j}}{}_{k}$ are used to compute the firstorder sensitivity of each parameter j at each point k of the parameter space:
where V_{K}(Ψ) is the total local variance at point k
The firstorder sensitivity measures $S{\mathrm{1}}_{k}^{j}$ vary between 0 and 1, and the sum of firstorder 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 Sandberg, 2002).
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 V_{mean} (m s^{−1}), maximum average runoff height H_{mean} (m), maximum flood area A_{max} (m^{2}) and deposited volume Vol_{dep} (m^{3}). These response variables are calculated using the outputs from FLO2D as
where j is the cell index; N_{WC} 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 V_{max}(j), H_{max}(j) and H_{final}(j) are the maximum flow speed, maximum flow depth and final runoff height of the cell j, respectively.
The FLO2D parameters considered for DELSA are those that describe the fluid rheology (Table 1): α_{1,2}, β_{1,2}, C_{v}, K, n, the SD parameter and the total volume of sediments mobilized Vol_{T}.
The detention coefficient SD is a model parameter that controls flow detention. The FLO2D user's manual and previous studies (D'Agostino and Tecca, 2006) 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 Vol_{Test} by each debris flow event was calculated with the equation proposed by Chang et al. (2011) (Eq. 12).
where A_{w} is the watershed area, A_{L} 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 C_{R} is the cumulative rainfall (C_{R}=76 mm). We obtain Vol_{Test} values of 185 000 m^{3} for Acerillas and 154 000 m^{3} for La Mesilla.
Debris flow concentration is assumed to vary with streamflow between a minimum concentration ${C}_{{\mathrm{v}}_{\mathrm{min}}}=\mathrm{0.1}$–0.4 and a maximum concentration ${C}_{{\mathrm{v}}_{\mathrm{max}}}$ at the time of peak flow, which is treated as a model parameter. To this end, we propose the following function for C_{v}:
where ϕ is a coefficient that changes the shape of the concentration curve. ϕ and ${C}_{{\mathrm{v}}_{\mathrm{min}}}$ are calculated in order to match the total volume and minimize ${C}_{{\mathrm{v}}_{\mathrm{min}}}$ value:
Rheological parameter ranges are obtained from previous debris flow studies (O'Brien and Julien, 1988; Boniello et al., 2010; D'Agostino and Tecca, 2006; Sosio et al., 2007; Rickenmann et al., 2006; Chang et al., 2011; O'Brien and Garcia, 2009) 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 ${C}_{{\mathrm{v}}_{\mathrm{max}}}$ that we ensure are implemented in DELSA. To this end, we develop a Python script that allows for running FLO2D in parallel and sequentially, reducing computational cost considerably. For Vol_{T}, 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)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 FLO2D. 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 postevent topography.
Based on the reference simulation, we choose reference values of V_{mean}=1 m s^{−1} and H_{mean}=1.5 m at both creeks. Additionally, we estimate the reference maximum flood area using Google Earth imagery – obtaining values of 246 500 m^{2} for La Mesilla and 37 000 m^{2} for Acerillas – and use deposited sediment volumes reported by Cabré et al. (2020) as reference values, which correspond to 102 000 and 6000 m^{3} for La Mesilla and Acerillas, respectively.
4.1 Choice of sample size
First, we seek to identify the minimum sample size N_{k} 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 N_{k} on the cumulative frequency distributions (CDFs) of DELSA firstorder sensitivity indices. Since we include N_{j}=9 parameters, the total number of simulations required for each case is ${N}_{\mathrm{t}}=({N}_{j}+\mathrm{1}){N}_{k}$. Figure 3 illustrates the sensitivity of DELSA results to variations in sample size N_{k} for four variables simulated by FLO2D, with respect to parameters β_{1} and ${C}_{{\mathrm{v}}_{\mathrm{max}}}$, at the Acerillas creek. Since the curves obtained for N_{k}=500 and N_{k}=1000 show slight differences, departing from the CDF for N_{k}=100, we conclude that a sample size of N_{k}=500 is adequate for further analyses. The sensitivity of DELSA results to N_{k} was also examined at La Mesilla creek, obtaining the same conclusion regarding sample size (not shown).
Figure 3 also shows that, depending on the target variable and parameter analyzed, firstorder sensitivity indices can be highly heterogeneous across the parameter space. In particular, the modeled response is highly sensitive to variations of β_{1}, with firstorder 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 ${C}_{{\mathrm{v}}_{\mathrm{max}}}$ 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 N_{k}=500) of local firstorder 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.
As expected, the simulated deposited volume Vol_{dep} 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 V_{mean} – is the total sediment volume Vol_{T}, whose sensitivities have the same order of magnitude as those produced by β_{1}. The high sensitivity of Vol_{T} on V_{mean} is explained because the former influences fluid rheology through C_{v}(t) (see Eqs. 4, 5 and 13). However, Vol_{T} 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, A_{max} decreases as the flow is forced to stop at deeper heights. This could be a structural weakness of FLO2D, 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 secondorder 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 FLO2D 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 FLO2D 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 FLO2D using 26 different sets of rheological parameters obtained from 45 previous studies. They found that the most influential parameters – in order of importance – were C_{v}, 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 C_{v}, 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 N_{k}=500 parameter sets), normalized by reference values obtained from postevent measurements (Vol_{dep} and A_{max}) and results obtained from the reference models (for H_{mean} and V_{mean}). 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 V_{ref}. 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.2 A_{ref} and 0.6–1.4 Vol_{ref}, respectively. We assume a weaker observational constraint in the case of the deposited volume because of possible uncertainties associated with different measurements techniques.
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 Vol_{dep} ≪ Vol_{T} at Acerillas, while Vol_{dep} ∼ Vol_{T} at La Mesilla; moreover, the maximum flood area at La Mesilla is approximately 6 times larger than in Acerillas. Thus, small variations in Vol_{dep} 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 postevent measurements for calibrating debris flow models (Chow et al., 2018; Cesca and D'Agostino, 2008; D'Agostino and Tecca, 2006; 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}, ${C}_{{\mathrm{v}}_{\mathrm{max}}}$, Vol_{T}, 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}, ${C}_{{\mathrm{v}}_{\mathrm{max}}}$) 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}, ${C}_{{\mathrm{v}}_{\mathrm{max}}})$ cover almost the entire original range.
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 (>Q_{75}). 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 ${C}_{{\mathrm{v}}_{\mathrm{max}}}$ (below the median of filtered values). A similar effect is observed for τ_{y}, whose filtered range results from the compensation of α_{2}, β_{2} and ${C}_{{\mathrm{v}}_{\mathrm{max}}}$. In summary, different combinations of α_{1,2}, β_{1,2} and C_{v} 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 wellknown 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., Beven, 2006; 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 ${S}_{{\mathrm{f}}_{\mathrm{2}}}$, 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 FLO2D 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 Tecca, 2006).
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 ObservationPrediction (OPR) statistic (Tiedeman and Hsieh, 2004), 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 ParameterPrediction (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).
We performed a sensitivity analysis on the parameters of a widely used numerical debris flow model (FLO2D) 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 welldefined channel with almost nonalluvial 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 Vol_{T} does not have an effect on Vol_{dep}, it is important for representing flow velocity, as Vol_{T} is used to estimate C_{v}(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 secondorder 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 FLO2D. Our results also support the idea that singlephase rheological models lack a strong physical basis (Iverson, 2003), and, therefore, their determination requires expert knowledge. However, an encouraging finding is that the final deposited volume Vol_{dep} and maximum flood area A_{max} 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 Tecca, 2006). Onephase 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.
Lidar topography and used streamflows are available at https://doi.org/10.6084/m9.figshare.c.5034968.v2 (Zegers, 2020).
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.
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 (MOPDOH) and the CONICYT/PIA (project no. AFB180004). Pablo A. Mendoza acknowledges the support of a Fondecyt postdoctoral grant (no. 3170079).
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).
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, https://doi.org/10.5194/nhess2012472020, 2020. a
Arattano, M., Franzi, L., and Marchi, L.: Influence of rheology on debrisflow simulation, Nat. Hazards Earth Syst. Sci., 6, 519–528, https://doi.org/10.5194/nhess65192006, 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 debrisflow event in the Fella watershed, Nat. Hazards Earth Syst. Sci., 10, 989–997, https://doi.org/10.5194/nhess109892010, 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.: Tributaryjunction 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 realworld application of Monte Carlo procedure for debris flow risk assessment, Comput. Geosci., 35, 967–977, https://doi.org/10.1016/j.cageo.2008.04.002, 2009. a, b, c
Cesca, M. and D'Agostino, V.: Comparison between FLO2D and RAMMS in debrisflow 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, https://doi.org/10.2495/DEB080201, 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, https://doi.org/10.1016/j.enggeo.2011.07.004, 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, https://doi.org/10.3390/geosciences8060218, 2018. a, b, c, d
Coussot, P. and Meunier, M.: Recognition, classification and mechanical description of debris flows, EarthSci. Rev., 40, 209–227, https://doi.org/10.1016/00128252(95)000658, 1996. a
D'Agostino, V. and Tecca, P. R.: Some considerations on the application of the FLO2D 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, https://doi.org/10.2495/DEB060161, 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 debrisflow and GLOF risk management strategy for a datascarce catchment in Santa Teresa, Peru, Landslides, 13, 1493–1507, https://doi.org/10.1007/s103460150669z, 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, https://doi.org/10.1029/2018WR022668, 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, https://doi.org/10.5194/nhess99132009, 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, https://doi.org/10.1016/j.geomorph.2006.02.003, 2006. a
Iverson, R. M.: The debrisflow rheology myth, Debrisflow hazards Mitigation: Mechanics, Prediction, and Assessment, Millpress, Rotterdam, 2003. a
Iverson, R. M., Reid, M. E., and LaHusen, R. G.: Debrisflow mobilization from landslides, Annu. Rev. Earth Planet. Sci., 25, 85–138, https://doi.org/10.1146/annurev.earth.25.1.85, 1997. a
Kelleher, C., Wagener, T., and McGlynn, B.: Modelbased 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 FLO2D model to debrisflows simulation – A case study of SongHer district in Taiwan, Ital. J. Eng. Geol. Environ., 3, 947–956, https://doi.org/10.4408/IJEGE.201103.B103, 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, https://doi.org/10.1016/j.cageo.2014.01.006, 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 processbased 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 onedimensional finite element simulation model, Nat. Hazards Earth Syst. Sci., 6, 155–165, https://doi.org/10.5194/nhess61552006, 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.: FLO2D Reference manual, available at: https://www.flo2d.com/ (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.: Twodimensional 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, https://doi.org/10.1061/(ASCE)07339399(2002)128:1(121), 2002. a
Olsson, A. M., Sandberg, G., and Dahlblom, O.: On Latin hypercube sampling for structural reliability analysis, Struct. Safe., 25, 47–68, https://doi.org/10.1016/S01674730(02)000395, 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 ENSOdriven 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, https://doi.org/10.1016/j.envsoft.2016.02.008, 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, https://doi.org/10.5194/nhess1120472011, 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, https://doi.org/10.1002/2013WR014063, 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 debrisflow simulation models with field events, Comput. Geosci., 10, 241–264, https://doi.org/10.1007/s1059600590213, 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: http://sitiohistorico.sernageomin.cl/pdf/presentacionesgeo/PrimerCatastroNacionalDesastresNaturales.pdf (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 debrisflow event, Earth Surf. Proc. Land., 32, 290–306, https://doi.org/10.1002/esp.1391, 2007. a, b, c, d, e, f, g
Takahashi, T.: Debris Flow, Annu. Rev. Fluid Mech., 13, 57–77, https://doi.org/10.1146/annurev.fl.13.010181.000421, 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 forcedand naturalgradient tracer tests in heterogeneous aquifers, Water Resour. Res., 40, W01512, https://doi.org/10.1029/2003WR002401, 2004. a
Tonkin, M. J., Tiedeman, C. R., Ely, D. M., and Hill, M. C.: OPRPPR, 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, https://doi.org/10.6084/m9.figshare.c.5034968.v2, 2020. a