Can the discharge of a hyperconcentrated flow be estimated from paleoflood evidence?

Many flood events involving water and sediments have been characterized using classic hydraulics principles, assuming the existence of critical flow and many other simplifications. In this paper, hyperconcentrated flow discharge was evaluated by using paleoflood reconstructions (based on paleostage indicators [PSI]) combined with a detailed hydraulic analysis of the critical flow assumption. The exact location where this condition occurred was established by iteratively determining the corresponding cross section, so that specific energy is at a minimum. In addition, all of the factors and parameters involved in the process were assessed, especially those related to the momentum equation, existing shear stresses in the wetted perimeter, and nonhydrostatic and hydrostatic pressure distributions. The superelevation of the hyperconcentrated flow, due to the flow elevation curvature, was also estimated and calibrated with the PSI. The estimated peak discharge was established once the iterative process was unable to improve the fit between the simulated depth and the depth observed from the PSI. The methodological approach proposed here can be applied to other higher‐gradient mountainous torrents with a similar geomorphic configuration to the one studied in this paper. Likewise, results have been derived with fewer uncertainties than those obtained from standard hydraulic approaches, whose simplifying assumptions have not been considered.


Introduction
[2] When studying torrential flows in steep mountainous catchments, one of the most important tasks is to correctly identify the process that occurred in the catchment to avoid errors in characterization. In this context, different terminologies and classifications are cited in the literature to describe flows involving water and sediments: hyperconcentrated flows, mudflows, debris flows, etc. [Pierson and Scott, 1985;Siviglia and Cantelli, 2005;Alexandrov et al., 2007;Desilets et al., 2008;Berzi and Jenkins, 2008;Sosio and Crosta, 2009;Iverson et al., 2010]. One of the greatest difficulties when characterizing floods in this geographic context is that several concurrent processes may be involved (i.e., debris flow, hyperconcentrated flow, and clear water flow), with different characteristics such as rheology or number of phases involved. The correct identification of the flows (i.e., Newtonian versus plastic) has long been regarded as the key to successful interpretation, modeling, and prediction of sediment-laden flow behavior. Hyperconcentrated flows were originally defined as those with a solid fraction (mainly clay minerals) equivalent to 2-4 weight percent, which provides the flow with a finite yield strength [Beverage and Culbertson, 1964]. Subsequently, Coussot and Meunier [1996] extended this rank from 1% to 25%. Other authors paid attention to the onset of a shear strength that for the specific case of noncohesive sand-grade particles occurs at 20%-60% solid volume fraction [Pierson and Costa, 1987;Costa, 1988;Pierson, 2005a]. The definition of debris flow is mainly focused on the existence of a fully nonlinear relationship between shear rate and strain [Coussot and Meunier, 1996;Coussot et al., 1998;Iverson, 1997Iverson, , 2003, where its solid fraction is ranked between 50% and 90% [Coussot and Meunier, 1996].
[3] Additionally, the distinction between water floods, hyperconcentrated, and debris flow depends on the rheological properties, the number of phases involved (multiphased flows), and deposit structure [Pierson and Costa, 1987;Coussot and Meunier, 1996;Svendsen et al., 2003;Lavigne and Suwa, 2004;Pierson 2005a;Hessel, 2006;Van Maren et al., 2009;Bisantino et al., 2010]. Costa and Jarrett [1981] distinguished between water flows and debris flows, and to a small degree, hyperconcentrated flows. They also differentiated between them based on sediment characteristics remaining in the channel. Hungr [2005] made a distinction between debris flow and hyperconcentrated flow based on the observed or potential peak discharge of a given event and Pierson [2005b] drew attention to the transitional behavior between water flow and debris flow.
[4] This transitional behavior means that different phases can be distinguished in the flow and so the rheology of these processes can be established between true Newtonian fluid, which is characteristic of normal streamflow, and non-Newtonian or plastic fluid that, among other things, determines the dynamic behavior of debris flow [Smith, 1986;Pierson and Costa, 1987;Sohn et al., 1999;Sosio and Crosta, 2009]. A Newtonian flow is one in which the stress state at any given point is proportional to the strain rate at that point. As a result, turbulence is the most important support process for sediment in the flow [Pierson and Scott, 1985;Costa, 1988;Wan and Wang, 1994;Winterwerp, 2001;Baas and Best, 2008]. In contrast, non-Newtonian flows are the result of fluids in which viscosity changes with the applied strain rate. The rate of shear is therefore not proportional to the corresponding stress that determines viscoplastic behavior. As a result, turbulence can be almost negligible and the most important sediment-supporting processes are buoyancy, dispersive stress, structural support, and cohesion and so solids and water move together as a single viscoplastic body [Costa, 1984;Coussot and Meunier, 1996;Hessel, 2006;McCoy et al., 2010].
[5] In debris flows, the whole mass undergoes very considerable and approximately continuous deformation, because there is little difference between the velocity of the two elements (water and sediment). However, in clean water or hyperconcentrated flows the mean velocity of the coarsest solid particles, which are mobilized by a bed load transport mechanism, differs significantly from water and fines in suspension [Xu, 1999;Svendsen et al., 2003;Shu and Fei, 2008]. In relation to the deposit structure, debris flows can be composed of sediment of all sizes ranging from boulders up to a few meters wide to clay-sized particles with negligible grading or internal layer structures [Sohn et al., 1999]. On the contrary, deposits left by hyperconcentrated flows tend to be poorly sorted and weakly stratified [Costa, 1984[Costa, , 1988].
[6] The ambiguous definition of hyperconcentrated flow, especially with reference to its rheological behavior, has allowed many flood events involving water and sediments to be characterized using classic hydraulics principles. The critical-depth method analysis [Chow, 1959] is recommended for high-energy channel conditions [Jarrett, 1987;Trieste and Jarrett, 1987;Grant, 1997;Jarrett and England, 2002], and it is increasingly being used to estimate flood discharges [House and Pearthree, 1995;Tickler, 1997;House and Baker, 2001;Rico et al., 2001;Gaume et al., 2004;Roca et al., 2009], as it depends on the geometry of just a single specified critical cross section, (i.e., flow surface width and cross-sectional flow area), which allows better estimations to be obtained for this kind of channel. The critical-depth method has frequently been used to determine paleodischarges in steep-gradient mountain torrents using a number of simplifying assumptions. This may have led to uncertainties in estimates of 615%, provided reach selection and modeling are made properly. In addition, for clear water flow or with modest sediment transport, air entrainment and flow resistance may lead to additional unknown uncertainties [Jarrett and England, 2002]. The Barranco de Arás flood (NE Spain) on 7 August 1996, which claimed 87 casualties, is one case of the use of the critical-depth method in the Iberian Peninsula to obtain flood discharges [Benito et al., 1998;Alcoverro et al., 1999]. This flood was very destructive, mainly because 37 check dams along the stream were destroyed, releasing most of the sediments they contained. As a result, the solid fraction of the flood was significantly increased, especially in the lower part of the catchment.
[7] The objective of this paper is to determine the discharge of a flow event that took place on 17 December 1997 in the Arroyo Cabrera, Spanish Central System. This event was described as a hyperconcentrated flow, because during a postevent flood inspection accomplished a few days after the occurrence, some deposits were observed in which sedimentation was poorly sorted with some stratification and imbrication. To carry out the estimation, the critical-depth method was used, combined with the use of PSI left by the event. The equations governing the physical processes linked to this hydraulic regime are complex. In fact, it seems impossible to obtain explicit equations in which some of the variables involved can be solved, unless important simplifications are made (i.e., the location of the control section just where the hydraulic drop or narrow constrictions in channel geometry are located; flow density equal to 1 T m À3 ; the angle between the thalweg and the horizontal is neglected ; assumption of hydrostatic pressure distribution in the hydraulic fall; no consideration of shear stresses and flow superelevation). Such simplifications were avoided using an iterative method. Additionally, the methodology described above also aims to be versatile enough so it can be used in other basins with a hydrogeomorphic response similar to the one studied in this paper.

Study Area and Physical Setting
[8] The Arroyo Cabrera, which is subject to torrential floods, debris flows, and hyperconcentrated flows, is a torrential tributary of the Alberche River in the Tagus Basin (Central Iberian Peninsula, Spain). It is formed from the confluence of several streams ( Figure 1). Its catchment extends over 15.5 km 2 and is roughly teardrop in shape. Watershed relief is 1224 m (from 735 m above sea level [asl] to 1959 m asl), and the main channel is 5500 m long, with an average slope of 21.6%. The bedrock consists of Upper Paleozoic granitoid and superficial quaternary formations, made up of conglomerates, gravels, sands, and silts, covering the slopes, the valley floor, and depressions known as navas. The thin soils developed on weathered profiles of granitic materials are mainly entisols, fluvisols, and inceptisols, with hydrologic soil group B predominating [Soil Conservation Service, 1972].
[9] Several factors determine the torrential behavior of the Arroyo Cabrera catchment. First, the occurrence of intense and/or large rainstorms due to orographic factors, combined with thin, poorly developed soils. In addition, the headwater has an average slope of 50%, although on some channel hillsides it can reach 100%. As a result, the runoff threshold and concentration time are decreased in relation to catchments having lower slopes. In addition, the steep slope of the main channel allows fast flood wave translation almost without basin storage. Other factors to take into account are its high drainage density, the almost circular morphology pattern, the existence of patches of impervious substratum, the occurrence of heavy rainfall events due to orographic factors, and fine grained and poorly developed soils.
[10] On the night of 17 December 1997, the most important catastrophic flow event in recent history took place in the Arroyo Cabrera. It was a complex multiphased event, having several types of geomorphic phenomena (i.e., shallow landslide, debris avalanche, hyperconcentrated flow, clear water flow), which operated simultaneously. The heavy rainfall event triggered a shallow landslide in the upper part of the catchment ($1800 m asl). During the months of November and December of 1997, the cumulated rainfall on this region, according to the available data was 821 mm recorded in the daily rain gage located in Los Serranillos (1230 m asl), 617 mm recorded in the daily rain gage located in El Burguillo (750 m asl), and 489 mm in La Adrada (720 m asl).
[11] In the abovementioned rain gages, the measured rainfall during the studied event was 141, 87, and 87 mm. These daily rain gages are located 25, 12, and 7 km, respectively, from the study area. Although they are not very far away from the study site, they are in the lower ranges of altitude in the Sierra de Gredos, so the values recorded are not representative of the characteristics of the orographic component of precipitation. The trigger ID-threshold of I ¼ 127.93 Â D À0.7077 (where I is the rainfall intensity and D is the rainfall duration), as well as critical rainfall intensity and cumulative rainfall, as a function of storm duration, have been also established by using indirect methods [Ruiz-Villanueva et al., 2011]. The boulder and gravel mass of the landslide tongue then evolved into a debris avalanche on the steeper slopes, supplying sediment load to the flood wave transforming it into a hyperconcentrated flow. Downstream, where the longitudinal profile slope becomes less steep and the channel wider, much of the sediment load transported by the flow was deposited. The hyperconcentrated flow continued losing solid load downstream until it was transformed into a water-dominated flow with floating woody debris and fine-grained sediment in the lower part of the catchment ( Figure 2). The postevent visit enabled recognizing, mapping, and interpreting the landforms and deposits left by the event. Among the erosive landforms, stand out channel erosion including the undercutting of outer margins of small meander bends. The study of the deposits consisted in the recognition and interpretation of the sedimentary structures and granulometric sequences. In this respect, poorly sorted deposits with some stratification and boulder imbrication could be recognized (see Figure 3).

Materials and Methods
[12] The methodology used to achieve the objectives stated in this research was based on the development of four methodological steps. The order of implementation of these steps, the required input data and the resulting output data are specified in Figure 4.

Suitability of the Critical-Depth Method
[13] The major requirement for using the critical-depth method is the existence of cross sections where critical flow is assumed to have occurred due to local drops or narrow constrictions in channel geometry, and provided that the subcritical regime can be demonstrated upstream. The channel must also be made of hard rock and paleostage indicators (PSI) or high-water marks (HWM) must be present to enable past flood levels to be estimated [Jarrett and Tomlison, 2000]. A thorough field survey was carried out on the entire Arroyo Cabrera to search for reaches that meet the characteristics described above. A topographic map of the chosen reach was produced using a combination of total station and GPS. From this map, detailed channel cross sections were defined with an average separation of 1 m. The cross sections were labeled using numbers between 157 (located on the hydraulic drop) and 200, which is the most upstream cross section. The bed of the studied reach is composed of granitoid bedrock, which satisfies a basic hydraulic criteria for discharge estimation, guaranteeing that changes in the channel geometry are minimal after a given flood event.
[14] As a preliminary assessment and to estimate the existence of critical flow, the HEC-RAS-1-D hydraulic model [U.S. Army Corps of Engineers (USACE), 2008] was used, however, its numerical results have not been used in this  paper. As flow data was not available, the whole reach was modeled considering a discharge range of between 25 and 1200 m 3 s À1 .

Discharge Estimation for the Observed Psi 3.2.1. Determining the Critical Depth in Each Cross Section
[15] Specific energy, H 0 , is given by: where y is the flow depth at a given cross section, is the angle between the thalweg and the horizontal, V m is the mean channel velocity, g is the acceleration of gravity, A is the area of the wetted cross section, Q is the discharge, and is the Coriolis coefficient. This coefficient takes a value of 1.45, instead of 1 as is usually used in a hydraulic calculation for clear water channels, following the recommendations of different authors [Shiono and Knight, 1991;Ervine et al., 1993;Bousmar and Zech, 1999;Sleiti and Kapat, 2008]. Normally, hydraulic models ignore the parameter , because for small slopes the cosine of this angle is practically equal to 1. However, for high-slope streams, as studied in this paper, this parameter had to be determined to ensure a high level of estimate certainty.
[16] The objective was to find a minimum energy flow state, associated with the critical regime, so that the PSI left by the flood could be used to take advantage of the existing stage-discharge relationship. To find this energy state the mathematical minimum was found by equating the derivative to zero, [17] As a result, [18] Because depth variation in a given cross section is the flow surface width (B), then [20] To obtain the critical depth, an iterative process was implemented applying (4). In this expression is known from topography, whereas is assumed constant in all cross sections.

Location of Critical Flow Section
[21] The approach used to pinpoint the location of the critical regime was to apply the momentum principle ( Figure 5). It was implemented starting from the cross section, where a hydraulic drop is located (hereafter defined as the final cross section, XS number ¼ 157) to a cross section i situated upstream ( Figure 5), where F x is the resultant force acting on the control volume in the direction of flow, is the flow density, v x is the velocity in direction x, and is the Boussinesq coefficient.
[22] To apply the momentum principle, the x-axis was considered perpendicular to the cross section with the hydraulic drop. The other forces and velocities are projected in this direction.
[23] The goal was to determine which cross section satisfies equation (5), knowing that in one of them the minimum energy condition is reached. Applying the momentum principle to a control volume between the final cross section (XS number ¼ 157) and the others (XS number ¼ 158 to XS number ¼ 200) located upstream, in which the minimum energy state is also determined, equation (5) becomes: where F Xfinal is the resultant force acting on the final cross section, is the angle between the bed channel and the vertical existing immediately downstream from the final cross section, F Xi is the resultant force acting on a cross-section i located upstream of the final cross section, is the slope in cross-section I, F zo are the existing shear stresses in the wetted perimeter, F w is the weight of flow, and is the angle between the bed channel and the vertical existing immediately upstream of the final cross section.
[24] If the curvature of the bed channel is also considered, (6) can be complemented using the angle ', corresponding to the deflection that occurs in the cross section just upstream of the hydraulic drop. The momentum equation will then be, [25] All of the angles cited above vary from one cross section to another and all of them were obtained from the detailed topographic survey carried out at the study site.

Pressure Distribution
[26] To estimate the pressure distribution, the water surface elevation curvature was determined. This plays a major role in the distribution of pressures [Chow, 1959]. The streamline radius of curvature was extracted as follows : because a PSI is available in the vicinity of the hydraulic drop, a typical S2 backwater surface profile associated with high-gradient and decreasing depths could be drawn. The bed channel was also defined precisely from the available topography. First, the tangent circle of the streamline situated at middepth of a given cross section was drawn (see Figure 8). The radius of curvature, R c , depends on depth, but estimating R c on the flow surface is complex, as the shape of the flow surface downstream from the hydraulic drop is unknown. It is more precise to determine the R c of the streamlines located close to the bed channel. Here the tangent circle corresponding to the streamline closest to the bed channel was also drawn (see Figure 9). A linear variation of the streamlines with depth was assumed. In this way, the R c of the flow surface was obtained and so the equation governing the pressure distribution could be defined.
[27] The application of the momentum principle meant a nonhydrostatic pressure distribution in the immediate vicinity of the final crosssection had to be considered, due to the abrupt change in slope. This distribution could not be obtained immediately, as it depends directly on the flow velocity, which is unknown.
where P f is the nonhydrostatic pressure distribution in the final cross section and z is the depth measured from the bed channel.
[28] Similarly, the resultant force, F X (see equation (5)), cannot be calculated directly because the geometry of the cross section is irregular, so that the pressure at each depth is applied in different areas. For this reason, the wetted area was divided into horizontal bands and the pressure increment at a given depth and in relation to the flow surface was considered on each one. This means that the pressure distribution was calculated by applying a finite difference method.
[29] In contrast, in the cross sections located upstream of the final cross-section curvature is much less, so that a hydrostatic pressure distribution could be assumed. The resultant pressures in each crosssection were therefore obtained from the formula: where P i is the resulting force of the hydrostatic pressure distribution in a cross-section i situated upstream of the final cross section, A is the area of the wetted section, and Z G is the center of gravity of the wetted section in relation to the water surface. Z G has to be calculated in each cross section for variable levels as the flow surface is unknown: 150 Z G values were determined using this process.

Estimation of Shear Stresses
[30] The resultant shear stresses, F z , between cross-section i and the final cross section was obtained from the expression, where Z o is the boundary shear stress and Pe is the wetted perimeter. F z, can be also obtained using the formula where U is the shear velocity, R h is the hydraulic radius, and S f is the friction slope. Similarly, U is related to Z o by the expression, [31] Substituting obtains [32] Shear stress on the sides of the channel was taken as one third of that corresponding to the bed channel, a usual approach in sediment transport analysis [Nouh and Townsend, 1979;Bathurst, 1985;Dietrich and Whiting, 1989;Biron et al., 1998]. The average shear stress in each cross section was calculated with weighting based on the wetted perimeter between the final cross section and those located upstream. Additionally, the density considered for the flow was 1.7 T m À3 [Beverage and Culbertson, 1964;Costa, 1984].

Iterative Resolution Scheme
[33] With the methodological approach used, it is impossible to solve the implicit equations entailed directly, because the discharge and therefore the velocity are unknown. Most of the parameters integrated in the equations are therefore also unknown. This meant that a trial and error process had to be used involving a large number of calculations and iterations. In this iterative process, many discharges were assumed, so that the variables integrated in the equations could be characterized. Many of them also had to be determined by iterations and so iterative processes had to be run within this trial and error process. To be more precise, it was used to apply the minimum specific-energy condition, so that variations depending on the depth of the cross-section areas and width flow surface were analyzed. It aimed to calibrate the moment at which equation (4) was satisfied. The same approach also had to be used for the implementation of the momentum principles as the wetted areas could not be established, as they depend on both discharge and channel geometry.
[34] The method used was completed with the application of the momentum principle and with the correction of the flow surface due to lateral superelevation. To do this, the results obtained were compared with the PSI to check the fit between the simulated and observed values, requiring the determination of the flow surface. The iterative process was repeated as many times as required to achieve optimum fit. Lateral superelevation was estimated using the PSI available and the formulas of Grashof and Woodward [Chow, 1959]. Grashof's formula, which assumes a logarithmic flow surface, is given by where R O and R I are the inner and outer radii of the bend, respectively.
[35] Woodward's formula assumes the velocity to be zero at the banks and to reach maximum value at the center of the bend with the velocity distribution varying in between according to a parabolic curve.
where V max is the maximum channel velocity and b is the flow surface width.

W12535
[36] Finally, the average value of the lateral water surface elevation was obtained by considering both methods. The iterative process finished when the optimal fit between the simulated depth, including consideration of the effect of superelevation, and the depth defined by the PSI was achieved.

Suitability of the Critical-Depth Method
[37] At the study site, at least one reach includes cross sections with appropriate morphological characteristics for critical flow to occur. There is an area where the longitudinal profile slope changes abruptly reaching a value of 78%, which can be considered a hydraulic drop. In contrast, upstream of this zone the slope is reduced considerably (see Figure 6). This suggests the existence of a control section in which the Froude number is equal to unity, provided upstream the flow regime can be assumed to be subcritical.
Similarly, in this reach there is evidence of PSI on both riverbanks, consisting of an erosional form with a sigmoid profile, which also supports this assertion (see Figure 6). The results derived from the application of the 1-D hydraulic model also make the assumption of the existence of critical flow possible. In fact, this hydraulic regime occurs over a range of discharges (Figure 7). In this respect, it is also worth noting that these findings are in agreement with those obtained by other authors [e.g., Jarrett and England, 2002].

Estimation of Discharge of the 17 December 1997 Event
[38] A peak discharge of 1080 m 3 s À1 was obtained from our analysis. Table 1 shows superelevation data obtained from the Grashof and Woodward formulae. If the average value of both is taken as representative, then the difference between the simulated and observed depths is only 7.2 cm. The location of the control section where the critical flow took place required the prior determination of the critical   (Table 2), as well as the corresponding wetted areas and critical velocities which are estimated immediately once critical depths have been calculated.
[39] Applying the momentum principle meant considering a nonhydrostatic distribution of pressures in the final cross section where the water surface elevation curvature is important. In addition, at this point the transition from subcritical to supercritical regime has already occurred. Figures 8 and 9 show the radii of curvature associated with the streamlines situated both in the middle of the cross section and closest to the bed channel. As a linear variation was assumed between radius of curvature and depth based on the expression, R c ¼ 1:74 Â y þ 4:45, then the average radius of curvature of the flow surface is 16.05 m.
[40] As a result of applying equation (9), the nonhydrostatic pressure distribution is shown in Table 3 and Figure  10. For comparison, the values and resulting graph of the hydrostatic pressure distribution for the final cross section are also included. The resultant F x on the control volume, as a result of the nonhydrostatic pressure distribution is included in Table 4. As the geometry of the final cross section is irregular, the values calculated correspond to different horizontal bands that were then combined. Table 5 shows the hydrostatic distribution established upstream of the final cross section. Table 6 includes the shear stress values for the final cross section. In addition, the resultants of those cross sections located upstream of the hydraulic drop are shown in Table 7. Finally, by iterating equation (7) and considering the estimated critical depths (see Table 2), the location of the control section was determined at a point situated between cross-sections 169 and 170 (12-13 m upstream from the hydraulic drop), and a discharge of 1080 m 3 s À1 was obtained.

Discussion
[41] We have used the critical-depth method, ignoring the simplifications usually used in the standard hydraulic approaches, to characterize a hyperconcentrated flow event. Thus, uncertainties in estimates can be minimized, especially if the control section is located with precision. Hydrologic assessments in ungaged basins, in which flow is highly torrential and where hyperconcentrated flows are common, have often been carried out using critical-depth methods that establish a stage-discharge relationship. So, it has been proved that hyperconcentrated flows exhibit similar characteristics such as water flows (i.e., hydraulic jump, subcritical, and supercritical regimes, instability, etc.) [Coussot and Meunier, 1996]. Therefore, when critical flow can be established, to be able to carry out estimates of flow discharges in ungaged basins requires both the geometry of a given selected cross section and PSI information as final input data. The use of this methodology for basins/channels having highly torrential flows normally assumes the existence of cross sections where clear water flow occurs [Benito et al., 1998;Alcoverro et al., 1999], as well as control sections that exactly coincide with local drops or narrow constrictions in channel geometry. This means that discharges calculated in this way can either underestimate discharge or not adequately reflect the risks associated with these events because the solid load is not taken into account and the method is implemented in the wrong locations.
[42] Most authors who have applied the critical-depth method have assumed numerous simplifications, which could significantly increase the degree of uncertainty in the estimates. These simplifications were minimized in our work, by using an iterative methodological approach that   allowed us to reduce uncertainty in the estimates. The formulae sometimes used [e.g., Daugherty et al., 1989, p. 361] took into account an unrealistic geometric configuration of the channel (i.e., rectangular or regular), and so their use in mountainous streams is not recommended due to their complex shape. On the other hand, critical flow does not appear only in the cross section that marks the change in the geometry of the channel, it actually occurs a little way upstream and its location is given by the cross section that satisfies the momentum principle.
[43] The accurate location of the critical flow section plays a major role in minimizing uncertainty. Usually, it is because the torrential streams having a complex and variable geometry both along the longitudinal profile and the cross section. Furthermore, in geomorphologic contexts similar to that studied here, the methodological approach used so far assumes that critical flow occurs where local drops or narrow constrictions in channel geometry are present [Benito et al., 1998], which is inaccurate, as critical flow actually begins slightly further upstream. Therefore, this is one of the main sources of uncertainty arising from the use of standard hydraulic approaches. In fact, in our study site discharges vary between 692 and 4790 m 3 s À1 depending on the cross section chosen and assuming that in each one of them the critical flow condition takes place ( Figure 11). Moreover, there is an additional source of uncertainty, as standard approaches simplify the geometry of the cross sections to get estimates. With regard to the   abovementioned, in our paper location of the control section has been determined accurately. To this purpose, we have applied the basic equations of hydraulics, having considered all the parameters involved as well as the depths associated to the PSI. Thus, the error was minimized so that the flow obtained is the one with less uncertainty, due to any change in channel geometry (both in the cross section and in the longitudinal profile) implies important variations in the estimates.
[44] Another simplification is the assumption of a flow density equal to 1 T m À3 , which is not representative of a typical sediment-laden flow occurring in a torrential mountainous catchment. Furthermore, the angle between the thalweg and the horizontal is frequently neglected. This is because normally the longitudinal slope of the river is low enough for the cosine of this angle to be close to unity. However, in high slope streams such as the one studied here, this parameter must be considered and included in the calculations to obtain the most reliable results. A hydrostatic pressure distribution is not obtained in the hydraulic fall either. On the contrary, because of the existence of a flow elevation curvature a nonhydrostatic pressure distribution occurs. This pressure distribution depends on the velocity magnitude and flow direction (i.e., the flow direction in the free hydraulic fall is not the same as the critical depth because they are not parallel). Therefore, these two factors must be taken into account to obtain the resulting pressure. Shear stresses were also considered, although their influence is not decisive in the results due to the short distance involved. Finally, to optimize the fit between the critical depth simulated and observed in the PSI, the fact that the flow surface elevation is not a horizontal plane must be considered. In fact, it is curved due to the centrifugal acceleration forced by the changes in the river geometry.
[45] To maximize the level of certainty of the flow estimate, the sediment load was considered taking into account all the equations involved in the process. However, to solve these equations some assumptions had to be made, which are explained in detail in section 3. Perhaps the most questionable point was to assume a fixed value for the flow density (1.7 T m À3 ). However, this assumption does not have a decisive impact on the methodological approach proposed. Flow density appears in the momentum equation (see equation (7)), where the terms F Xfinal and F Xi include density in their calculation. Therefore, this parameter is present in almost all of the terms of the equality, so that density does not decisively affect both the location of the control section and the flow estimated. As a result of the choice of the Coriolis coefficient the uncertainty is negligible as well. The reason is because it only affects the differences in the term velocity between the two cross sections and the velocity between them varies little due to cross-sections are next enough. With regard to the Boussinesq coefficient, uncertainty is also not significant because this coefficient only appears in the right term of the momentum equation (see equation (7)).
[46] The discharge calculated, 1080 m 3 s À1 , is only partially related to the rainfall characteristics and more on the availability of sediment contributed upstream, and because the sediment source and supply is variable with time. This estimate exceeds the discharge that could be produced by the PMP (probable maximum precipitation) [Cudworth,   1989] of the catchment. The feneration of a discharge of this magnitude would require almost instantaneous precipitation values close to 1900 mm uniformly distributed across the whole catchment. Therefore, this estimate is a result of the sediment load that was translated together with the water flow. It also shows that in physiographic contexts similar to that presented in this paper, flow estimates from the consideration of critical regime and clear water flow, which is an approach used in the literature, do not adequately reflect the characteristics of flows with low annual exceedance probability.
[47] Finally, it is worth outlining the value of combining hydraulic approaches with geomorphological approaches, given the limited systematic record of flood events in mountainous basins. In this geomorphic context, the availability of precise topography or hydrometeorological data is rare, which makes reconstruction of events based solely on hydraulic criteria difficult. In this regard, geomorphology allows us to infer the areal extent, the rheologic behavior, and eventually the frequency of events, as well as to calibrate hydraulic models from PSI or HWM [Ballesteros et al., 2011]. Therefore, the uncertainty of the models is substantially reduced and as a result, risk analysis is more credible. Even when there is insufficient data to assess hazard from a hydraulic methodology, coupling flood mapping (based on geomorphologic criteria) with land use can help to approximate vulnerability and risk if frequency can be estimated.

Conclusions
[48] We assessed hyperconcentrated flow discharge by using an iterative methodological approach complemented by the use of PSI. Additionally, a geomorphologic model was carried out. It enabled recognition of processes linked to the studied event and, therefore, a determination of which rheology is the one that best represents the 1997 event. The critical-depth method can be implemented beyond the realm of the basin studied, provided that in the assessed reach rock is outcropped and PSI evidence is available. The discharge estimated here, 1080 m 3 s À1 , reflects combined water and sediment volumes, which is common in high-gradient streams of mountainous basins. Concerning the above, this is one of the major limitations of standard approaches that assume clear water flow, which is not representative of highly torrential flows. In addition, other simplifications assumed in the standard hydraulic approaches, such as the use of basic formulae based on the existence of a stagedischarge relationship were not considered, so that uncertainty in the results was minimized. The use of such simplistic formulae can lead to incorrect flow estimates and, as a result, to errors when analyzing the hazards of flood events that occur in high-gradient basins. This major source of uncertainty in discharge estimates may be minimized through the implementation of an iterative methodology. In the frame of the EU Flood Directive, our work may contribute to a better analysis and management of flood risk, especially in Mediterranean countries where the majority of casualties from flooding occur in high-gradient basins and as a result of processes similar to the one described in this paper.