Flood Hazard Management in Public Mountain Recreation Areas vs . Ungauged Fluvial Basins . Case Study of the Caldera de Taburiente National Park , Canary Islands ( Spain )

Las Angustias River is an ungauged stream in the Caldera de Taburiente National Park (Spain), where frequent intense flash-flood events occur. The aim of this research is to analyze the flood hazard at the Playa de Taburiente. Based on the limited information available (short time-series of daily precipitation), a statistical frequency analysis of 24 h rainfall was completed and the precipitation results were transformed into surface runoff. To determine if the model underestimates the flows that are generated in the basin, the dendro-geomorphological information available was used to calibrate results. The results of the HMS model were significantly lower. At this point, both the rainfall data and the rainfall-runoff model were re-analyzed to maximize the rainfall intensity values and the runoff generated (increasing the CN value for the basin). For the 1997 flood event, a 1250 m3·s−1 flood minimizes the RMSE for the disturbed tree sample; this flow value also clearly exceeds any peak flow derived from the rainfall-runoff analysis. It is only when rainfall intensity and surface runoff are maximized that the peak flows obtained approximate those associated with dendro-geomorphological data. The results highlight the difficulties of flood hazard management in ungauged torrential basins in mountain recreational areas (such as National Parks). Thus, in the absence of flow records, when considering the maximum rainfall intensity scenario may be a useful and effective tool for flood risk management.


Introduction
Flash floods are common in mountain basins, as a result of their high slopes and quasi-circular morphology, the resulting strong connectivity, and high order streams [1][2][3]. The precipitation in these basins also has an important orographic component, and as a result, is very variable in spatio-temporal terms [4,5]. This means that mountain basins are highly prone to extreme precipitation events, in both

Case Study Area
The Caldera de Taburiente National Park is on the island of La Palma, the north-westernmost island of the Canary Islands archipelago, Spain (Figure 1). Declared a national park in 1954, this area of 47 km 2 is centrally situated in the northern part of the island (28°40′-28°46′ N, 14°08′-14°13′ W (Greenwich Meridian)). The altitude ranges from 2426 m above sea level (a.s.l.) in the Roque de Los Muchachos to 430 m a.s.l. in the Barranco de Las Angustias just on the edge of the park. The main part of the park is formed by the headwaters of the Barranco de Las Angustias, a semicircular basin 8 km in diameter and 2000 m from top to bottom, resembling a great volcanic crater or caldera, although its morphology is in fact the result of the superposition of several volcanic edifices, erosive phases, and large landslides [21]. This research focuses on the central area of the park, known as 'Playa de Taburiente', which is a Y-shaped 2 km long stream reach, formed by two gorges, the 'Verduras de Alfonso' and 'Cantos de Turugumay' ephemeral rivers, and the Taburiente River, the result of their confluence (Figure 1).

Geology
The Canary Islands were formed by volcanic activity around 20 Ma ago with an approximately westward age progression, with La Palma and El Hierro the two most recent islands of the archipelago. The volcanic history of La Palma originated around 4 Ma ago when a seamount formed on the Jurassic oceanic crust and was subsequently elevated to an altitude of up to 1500 m a.s.l., where it now forms the Basal Complex [22]. Subaerial activity started around 2.0 Ma [23], when a multi- The main part of the park is formed by the headwaters of the Barranco de Las Angustias, a semi-circular basin 8 km in diameter and 2000 m from top to bottom, resembling a great volcanic crater or caldera, although its morphology is in fact the result of the superposition of several volcanic edifices, erosive phases, and large landslides [21]. This research focuses on the central area of the park, known as 'Playa de Taburiente', which is a Y-shaped 2 km long stream reach, formed by two gorges, the 'Verduras de Alfonso' and 'Cantos de Turugumay' ephemeral rivers, and the Taburiente River, the result of their confluence (Figure 1).

Geology
The Canary Islands were formed by volcanic activity around 20 Ma ago with an approximately westward age progression, with La Palma and El Hierro the two most recent islands of the archipelago. The volcanic history of La Palma originated around 4 Ma ago when a seamount formed on the Jurassic oceanic crust and was subsequently elevated to an altitude of up to 1500 m a.s.l., where it now forms the Basal Complex [22]. Subaerial activity started around 2.0 Ma [23], when a multi-stage shield volcano was formed. This central edifice experienced at least one major collapse, resulting in a landslide and subsequent rebuilding of the edifice [23,24]. Subaerial volcanic reactivation, with explosive volcanism predominating during its initial stages producing abundant volcanoclastic and phreatomagmatic materials at the base of the subaerial edifice, continued until at least 0.41 Ma. Landslides on the Taburiente and Cumbre Nueva around 0.7 Ma resulted in the formation of the Caldera de Taburiente by incision and retrogressive erosion [23].

Precipitation
The three main characteristics of the rainfall in the Caldera de Taburiente are: (i) inter-annual irregularity; (ii) seasonality; and, (iii) torrentiality. These characteristics are also evident on the E and SE slopes of other islands in the Canary archipelago [25,26]. The mean annual rainfall is 1018.6 mm, but this is irregularly distributed (i), so that the coefficient of variance is 51.7% and the standard deviation 526.8. The marked seasonality (ii) is evident as over half of all precipitation falls during the winter, but in the summer months, when the islands are affected by the Azores anticyclone, less than 1% of the annual rainfall is recorded. Finally, the torrential rainfall (iii) is one of the most significant traits as in some months several times the average monthly rainfall has been recorded; e.g., more than seven times in January 1979; or six times in February 2010. The torrential nature of the precipitation is also evidenced by the 200 mm of rainfall in a single day recorded on 44 occasions in just 34 years ( Figure 2); this amount of rain therefore represents a significant percentage of the monthly and annual totals (average 40% and 25%, respectively). stage shield volcano was formed. This central edifice experienced at least one major collapse, resulting in a landslide and subsequent rebuilding of the edifice [23,24]. Subaerial volcanic reactivation, with explosive volcanism predominating during its initial stages producing abundant volcanoclastic and phreatomagmatic materials at the base of the subaerial edifice, continued until at least 0.41 Ma. Landslides on the Taburiente and Cumbre Nueva around 0.7 Ma resulted in the formation of the Caldera de Taburiente by incision and retrogressive erosion [23].

Precipitation
The three main characteristics of the rainfall in the Caldera de Taburiente are: (i) inter-annual irregularity; (ii) seasonality; and, (iii) torrentiality. These characteristics are also evident on the E and SE slopes of other islands in the Canary archipelago [25,26]. The mean annual rainfall is 1018.6 mm, but this is irregularly distributed (i), so that the coefficient of variance is 51.7% and the standard deviation 526.8. The marked seasonality (ii) is evident as over half of all precipitation falls during the winter, but in the summer months, when the islands are affected by the Azores anticyclone, less than 1% of the annual rainfall is recorded. Finally, the torrential rainfall (iii) is one of the most significant traits as in some months several times the average monthly rainfall has been recorded; e.g., more than seven times in January 1979; or six times in February 2010. The torrential nature of the precipitation is also evidenced by the 200 mm of rainfall in a single day recorded on 44 occasions in just 34 years ( Figure 2); this amount of rain therefore represents a significant percentage of the monthly and annual totals (average 40% and 25%, respectively).

Vegetation and Land Use
Pinus canariensis forest is the most characteristic vegetation of the Caldera de Taburiente where it forms extensive forests and nowadays is one of the best preserved forests of this species in the Canary Islands. On the riverbed, there are occasional groves of willows (Salix canariensis) and some stands of the less demanding laurisilva species in the most protected shaded areas; while at higher altitudes, the plant cover is formed by different shrub species [27]. Although the Caldera de Taburiente is very sparsely populated due to its rugged topography, the research results of Garzón-Machado et al. [28] indicate that the Canary pine forest has been severely impoverished by herbivore activity, probably resulting from the introduction of goats by the Auaritas, the original inhabitants of La Palma. Other land uses include forest harvesting, hunting, and a few relatively recent rainfed crops of low economic value which share the land with grazing.

Precipitation Analysis
Data was used from two stations: a pluviometer (C106U, daily data) installed inside the Caldera de Taburiente N.P., and a continuous record pluviograph (C139E, hourly data) installed in La Palma airport, thirty km from the basin (Figure 1). Other measuring stations are installed in the Caldera, including another rain gauge situated at a lower altitude than the area of interest with a short, poor quality record series; there are also other gauges (monthly total data) that are installed at higher altitudes, but their monthly values were not considered useful for the aims of this case study.
Using data obtained from the rain gauge installed in the Caldera, a frequency distribution analysis was carried out to estimate the precipitation quantiles that are associated with different return periods. These quantile estimates were obtained from annual maximum series (SQRT-ET MAX , an extreme value distribution proposed by Etoh et al. [29], which has been widey used in Spain for annual maximum rainfall data statistical analysis; Gumbel and GEV functions) [30], and from exceedance series (Generalized Pareto-POT). The use of exceedance series to estimate or assign frequencies to extreme events is supported by many previous studies [31][32][33], enabling a better fit of the high range of the frequency function [34,35], and has been applied to both flow and precipitation estimates. These exceedance series were derived from precipitation events with a duration of 24 h, defined by an absence of precipitation recorded in the 24 h before and after this rainfall event. The precipitation threshold value (the process was carried out by the graphic method) was 25 mm rainfall. This value produced an exceedance series of 357 values above the threshold, which is significantly higher than the length of the annual maximum series (37 years).
Based on the results of the goodness-of-fit test performed (Kolmogorov-Smirnov test and Mean Square Error test), the Generalized Pareto frequency function was chosen, with parameters of fit estimated by the L-Moments method, applied to the exceedance series over the threshold. Kolmogorov-Smirnov test shows very similar values for all of the distribution functions: while MSE test shows that GP-POT distribution has the lowest error value (Table 1).
Daily rainfall distribution in the form of hourly time series has been obtained with the IDF curves currently used in Spain for the study area [36]. The precipitation time distribution throughout the event is important, as this will have a clear influence on parameters, such as initial abstractions and runoff generation, and therefore on the flows generated. These data were used to generate synthetic hyetographs using the alternating blocks method [37], for the three day duration of the event.
The attempts to maximize rainfall intensity from the available rainfall data focused on two aspects. On the one hand, we considered the influence of the rainfall event duration by using different values. Rainfall events of 72, 24, and 1 h were considered; these values are related to the duration of the maximum volume rainfall volume event in 1997 (selected for the availability of dendro-geomorphological information), duration of the event considered in the maximum rainfall frequency analysis, and duration of a local convective rainfall event. On the other hand, we considered the possible effect of orography on the total basin rainfall volume. This effect was considered because of the high basin relief, with the resulting difference in altitude between the location of rainfall gauging stations and the upper basin where the study area is located. The orographic effect was estimated using linear regression from the available rainfall data.

Rainfall-Runoff Modeling
The basin model ( Figure 3) was derived from a 5 m Digital Elevation Model (DEM). Disaggregation was based on both the spatial distribution of the physiographic factors that determine a homogenous hydrologic response (i.e., lithology, soils, cover type, hydrologic conditions, and slope) and those control sections where dendro-geomorphological data was available. The use of simulation methods for the various processes involved was based on the information available for the basin. The use of the SCS-CN loss method, which is well-developed and widely accepted in the US and other countries [38][39][40], was based on the availability of high resolution orthophotos for the area and also on the possibility of conducting fieldwork to test and adjust the values that were previously obtained by photointerpretation. The use of other methods of estimating initial losses or abstractions (e.g., Green-Ampt) was discarded due to the lack of more detailed soil composition and texture data.
The absence of rainfall and flow data associated with the same event reduces the availability of transformation methods in the HEC-HMS environment to only two options, so the runoff hydrograph was obtained using the Clark unit hydrograph [41], with estimated concentration time (tc) and storage coefficients. Time of concentration was derived using the Témez formula, which depends on the length of the main channel for each sub-basin and its slope [42]. The storage coefficient was calculated assuming that it represents 0.6 tc. The Muskingum-Cunge method was implemented for the flood wave routing [43], due to the absence of observed flow data. The Muskingum-Cunge parameters were calculated based on the flow and channel characteristics. This method involved the use of a finite difference scheme for solving the Muskingum equation, where the parameters in the Muskingum equation were determined based on the grid spacing for the finite difference scheme and the channel geometry characteristics. The channel geometry and roughness were accurately defined using LIDAR (0.5 points per m 2 ) and orthophoto data combined with GIS tools.
Although we considered the associated uncertainties an essential point, the absence of flow data to calibrate and validate the hydrological output model flows makes it impossible to estimate this in each hydrological calculation, and its overall effect as the methodology unfolds. However, the sensitivity of the HEC-HMS model to SCS-CN value variations was considered (CN increase in 10% and 25%) in fluvial sub-basins upstream of the "Playa de Taburiente".

Flood Discharge Estimation Using Manning's Equation
In paleoflood reconstruction, event dating and frequency estimation are essential, and also a description of the magnitude of discharges, volumes, and water circulation times [44]. Estimates of peak discharges of past flood events in ungauged streams, such as the Taburiente river basin, are usually obtained using different types of paleostage indicators (PSI), hydraulic formulae, and models. PSIs can be high-water marks (HWM) from slack water deposits, floating wood jams or dendroevidence, such as tree stem scars. Hydraulic methods range from very simple formulae (Manning, Chèzy, Froude, Darcy-Weisbach...) to very complex mathematical-numerical models (simultaneous simulations from one-dimensional to three-dimensional and from monophasic to biphasic).
For the reconstruction of peak discharges in the Taburiente streams the height of scar dendroevidence and Manning's equation were used. These were chosen in preference to more complex methods, such as two-dimensional (2D) or three-dimensional (3D) numerical models since: (i) there are several sources of uncertainty, such as the mobile riverbed, that only allow estimation of the order of magnitude of peak flows; (ii) the flow regime during floods is mainly one-dimensional and generally uniform; and, (iii) the geometric data does not allow for the adjustment of very complex numerical models, taking into account turbidity, eddy structures, etc. Sometimes it is more reasonable to use simpler methods than more complex techniques without the necessary data sources.

Flood Discharge Estimation Using Manning's Equation
In paleoflood reconstruction, event dating and frequency estimation are essential, and also a description of the magnitude of discharges, volumes, and water circulation times [44]. Estimates of peak discharges of past flood events in ungauged streams, such as the Taburiente river basin, are usually obtained using different types of paleostage indicators (PSI), hydraulic formulae, and models. PSIs can be high-water marks (HWM) from slack water deposits, floating wood jams or dendro-evidence, such as tree stem scars. Hydraulic methods range from very simple formulae (Manning, Chèzy, Froude, Darcy-Weisbach...) to very complex mathematical-numerical models (simultaneous simulations from one-dimensional to three-dimensional and from monophasic to biphasic).
For the reconstruction of peak discharges in the Taburiente streams the height of scar dendro-evidence and Manning's equation were used. These were chosen in preference to more complex methods, such as two-dimensional (2D) or three-dimensional (3D) numerical models since: (i) there are several sources of uncertainty, such as the mobile riverbed, that only allow estimation of the order of magnitude of peak flows; (ii) the flow regime during floods is mainly one-dimensional and generally uniform; and, (iii) the geometric data does not allow for the adjustment of very complex numerical models, taking into account turbidity, eddy structures, etc. Sometimes it is more reasonable to use simpler methods than more complex techniques without the necessary data sources.
Manning's equation, or, more accurately, the Strickler-Manning equation, is widely used (see Díez-Herrero et al. [45], and many other authors), to relate flow velocity to channel morphometry parameters with a roughness coefficient listed for different bed materials: where: v is average streamflow velocity; R is hydraulic radius (A/P); S is energy grade line slope (coincident with the bed); and, n is roughness coefficient. The channel morphometry parameters come from a digital elevation model with 1 × 1 m spatial resolution, from LiDAR images. The transverse cross-sections were placed where trees with dendro-evidence have been dated and assigned to a specific event; and, obtained using ArcGIS 10.3 tools (ESRI Geosystems, Redlands, CA, USA). The roughness values for the riverbed and margins were assigned using ortophoto reclassification and fieldwork, using visual guide classes [46]. Finally, the computations were performed using special macro tools implemented in a Microsoft Excel spreadsheet.

2D Hydraulic Modeling
The hydraulic modelling was carried out with numerical simulation using a finite differences model (St Venant 2D equations), using the Iber hydraulic simulation free software [47], with the following characteristics: -Peak flows associated with precipitation recorded during 11-13 January 1997 event, obtained from the hydro-meteorological model.
-Peak flow range (300-2500 m 3 ·s −1 ), to minimize errors (RMSE) when compared with the vegetation disturbance data (tree scars), by comparing the water sheet depths with the height of the disturbances. The aim here was to define the liquid flows needed to minimize this error.
The topographical information relating to the channel geometry of the Caldera de Taburiente examined in this case study was included in a 1 m resolution DEM covering the whole area where disturbed trees were found and the reaches immediately upstream in the Cantos de Turugumay and Verduras de Alfonso ravines. The altimetric accuracy of the data in this DEM is shown by the RMSE data, with an estimated value of z ≤ 0.5 m.
To set the boundary conditions, two independent inflow points were identified, corresponding to the Cantos de Turugumay and Verduras de Alfonso ravines; and, a single outflow point was identified in the channel of the Taburiente, approx. 500 m from the area of the disturbed trees. A dry scenario was also considered as an initial condition for the model since the function of the river is ephemeral, with an estimated basin concentration time of 4.5 h to its outlet into the sea. Finally, the roughness of the terrain (expressed by the Manning n coefficient) was derived from the photointerpretation of the satellite orthoimages with 0.5 m spatial resolution, and field work. These delimited a series of zones with homogeneous relief roughness, grain size, and appearance, and vegetation density, which were then classified following the relevant tables and graphic guides [46]. Once the hydraulic model was developed, the most important results obtained from model simulation include data related to depth, absolute elevation level, and velocity field variables. But, just as with the results of the rainfall-runoff model, the absence of flow gauge stations in the basin did not allow for the calibration of the model, and, therefore, the validation of results.

Dendrogeomorphological Data Sources
Dendro-geomorphological evidences consisted of 266 samples collected from 60 disturbed trees of P. canariensis (i.e., those presenting scars, exposed roots, resprouting or apical bud loss, and dead trees) following Díez-Herrero et al. [48]. In some cases, trees showed several wounds corresponding to different events, but more often only one impact signal per tree was observed. Other nearby trees with no visible injuries were also sampled, as scars may have already closed if they are old or not very wide. Furthermore, 16 undisturbed trees were sampled to provide a reference chronology.
63 wounds were dated, in some cases, more than one for each tree, with very different percentages in the two cohorts and in total. Especially relevant were the injuries occurring in 1962 and in 1997, both presenting a large number of replications (48% and 22%, respectively, of the total number of trees).Wounds dated to 1993 and 2003 are represented in 23% and 20% of the younger trees. Over the last decade injuries have been dated to 2001, 2003, 2007, and 2009, which together account for 18%.
Other anomalies in the samples were also found with evidence of injuries but were not reliably dated [49].

Maximum Precipitation Results
The results obtained, from precipitation frequency analysis shown in Figure 4 and Table 1, show how the GEV distribution function systematically underestimates the values of the precipitation quantiles when compared with SQRT-ET MAX and GP-POT distribution functions. On the other hand, the SQRT-ET MAX function applied to the maximum series yields results higher than those obtained by the GP function applied to the exceedance series. From goodness-of-fit analysis (Table 1), we concluded that the best overall results could be assigned to General Pareto-POT distribution function.
In the analysis of rainfall data associated with the 1997 flash flood event, the first attempt considered the greatest rainfall volume for that year (11-13 January 1997), which recorded 191 mm rainfall in 72 h. However, the 1997 rainfall record at meteorological station C106U shows that the rainfall event for 19-20 January reached a value of 156 mm, and therefore greater rainfall intensity ( Figure 2). In addition, this second event suggests the possibility that this volume of precipitation, although recorded in two days, was associated with heavy rains occurring in a time interval of only one hour, divided between these two days. This would be the most favorable situation for generating a high magnitude flash flood.
On the other hand, the influence of a rainfall-elevation gradient in the basin has also been considered, and for this, rainfall that was recorded at three meteorological stations located on the island of La Palma has been observed: Station C139E (La Palma airport), located at 33 m a.s.l.; C106U (Caldera de Taburiente) at 820 m a.s.l.; and, ESICA3800000238750A at 568 m a.s.l. These stations share records for years 2005 to 2013. Within this time period, four high rainfall events were selected, and a linear regression was fitted to the daily rainfall data (Figure 4). The degree of fit (R 2 ) for linear regression varies between 0.74 and 0.99. For the best adjustment obtained (R 2 = 0.99), a linear rainfall-elevation gradient of +0.071 mm m −1 was estimated, with +0.079 mm m −1 being the mean value for the four events.
Applying the best-fit rainfall-elevation, it was estimated that a hypothetical rainfall event lasting 1 h (between 19-20 January 1997) could reach mean values of 236 mm ("Verduras de Alfonso" gorge-W1580), 234 mm ("Cantos de Turugumay" gorge-W1530), and 184 mm ("Playa de Taburiente" sub-basin-W1570). All of these represent a significant increase in the possible rainfall intensities that were associated with triggering the flash flood analyzed here. Finally, a situation was envisaged in which the return period rainfalls considered were also produced in a time interval of 1 h, or even less, in view of other historical events. Table 2 shows hourly rainfall intensity for all rainfall event models considered.    SQRT-ETmax, Gumbel and GEV functions were applied to annual maximum data series, while the Generated Pareto (GP) function has been applied to exceedances over threshold data series (POT technique). The linear regression fit for rainfall vs. altitude analysis from rainfall data available is shown at the bottom of the table.

Flood Discharges
When considering the three precipitation distribution scenarios over the 72 and 24 h duration of the event mentioned above allows for a study of the basin behavior in each case, especially of aspects that are related to the ground infiltration capacity in each scenario, and therefore to the surface runoff that produces variations in the flow rates generated. These flow rates ( Figure 5) show how the three scenarios that are considered represent situations, which may be more or less favorable to flash flood generation. In all cases, the right skewed hyetograph produces the highest peak flows, as usual.
The different hyetograph models used obtain high variability in their associated peak flows (Table 3), ranging from tens to hundreds of cubic meters. These differences are linked in the first instance to the different precipitation return periods considered. In view of the results obtained, the variability linked to peak flow values cannot be related only to differences in rainfall return periods (which may explain approximately 50% of the observed variability), highlighting the influence of the other variables that are considered on the results. These variables would therefore be responsible for the remaining variability.
The results also highlight the importance of rainfall intensity (ranging from 72 h to 1 h), which by itself can increase 2-3 times the peak flow rate generated for both the 1997 flash flood event and the considered return periods, a much higher increase than the one associated with the variation in the CN value, which may be linked to a 35%-50% variation (when considering a 25% increase in the CN values). Shortening the duration from 24 h to 1 h in the hyetograph associated with the different return periods considered obtains an increase in peak flow values close to 50%. In addition, the effect of the CN parameter is conditioned by the rainfall intensity value, showing a positive correlation. Table 3. Rainfall-runoff model (HEC-HMS) peak flows (m 3 s −1 ) at Taburiente fluvial basin location points.  Directly related to rainfall intensity, the rainfall-elevation gradient derived from the precipitation record analysis at three stations on La Palma increased the values of the generated peak flows by up to 50%. When considering all these factors (Table 3), may mean a four to five-fold increase in the 1997 flash flood peak flow.

Rainfall-Runoff Event
On the other hand, using the available dendro-geomorphological data for the study area [49], the peak flow rate that is associated with each case was then estimated using the Manning´s formula. The results obtained for the data set for the 1997 flood event show high variability, which is reduced if only six cases are used, shown grouped (Table 4 and Figure 6). Independently of whether all the data are used or only the data for the grouped cases, the mean data value (~1800 m 3 s −1 ) is double the maximum peak flow value obtained by processing the rainfall data and implementing it in a hydrological model. The first results obtained from a preliminary 2D hydraulic model for the study area obtain a peak flow value approaching 1250 m 3 s −1 , minimizing the RMSE coefficient when comparing flow depth and the height of the dendro-geomorphological evidence. This would be an intermediate value between the maximum value obtained in the hydrological model (HEC-HMS), and the average peak flow value obtained using Manning's formula.   Directly related to rainfall intensity, the rainfall-elevation gradient derived from the precipitation record analysis at three stations on La Palma increased the values of the generated peak flows by up to 50%. When considering all these factors (Table 3), may mean a four to five-fold increase in the 1997 flash flood peak flow.
On the other hand, using the available dendro-geomorphological data for the study area [49], the peak flow rate that is associated with each case was then estimated using the Manning´s formula. The results obtained for the data set for the 1997 flood event show high variability, which is reduced if only six cases are used, shown grouped (Table 4 and Figure 6). Independently of whether all the data are used or only the data for the grouped cases, the mean data value (~1800 m 3 s −1 ) is double the maximum peak flow value obtained by processing the rainfall data and implementing it in a hydrological model. The first results obtained from a preliminary 2D hydraulic model for the study area obtain a peak flow value approaching 1250 m 3 s −1 , minimizing the RMSE coefficient when comparing flow depth and the height of the dendro-geomorphological evidence. This would be an intermediate value between the maximum value obtained in the hydrological model (HEC-HMS), and the average peak flow value obtained using Manning's formula.

Comparing Results
The first and most important point to highlight when comparing the results obtained using hydrological-hydraulic models is the significant difference between the flows obtained from the rainfall-runoff transformation process of the statistical approximation (~50 m 3 s −1 , for the 11-13 January, rainfall recorded data), and those that were obtained taking the dendro-geomorphological data into account (~1250-1800 m 3 s −1 ). This latter flow value clearly exceeds the return periods considered (maximum 500 years, Table 3) in the rainfall-runoff hydrological model developed for the study area, and the hydrological models where rainfall intensity, CN, and rainfall-altitude gradient are considered. The main outcome of the results shown in Table 3 is variability, which obviously does not facilitate decision making. Depending on the data used and how they are processed, the results show a high standard deviation.
In addition, these flow differences mean a distinct behavioral change in this gorge. As shown in Figure 7, it functions either as a braided channel with emerged bars (~50 m 3 s −1 ) or as a single channel occupying the whole river bed (~1250 m 3 s −1 ). The difference in channel behavior is less evident (Figure 7) if we compare the flows that minimize RMSE value versus dendro-geomorphological data height (~1250 m 3 s −1 ), and those that are associated with the 500-years return period (510 m 3 s −1 , from 500-years return period rainfall according to the GP-POT distribution function). In these cases, a single channel occupying the whole river bed is the main river configuration. Figure 8 shows the clear differences in flow-depth data in this channel for each peak flow considered. study area, and the hydrological models where rainfall intensity, CN, and rainfall-altitude gradient are considered. The main outcome of the results shown in Table 3 is variability, which obviously does not facilitate decision making. Depending on the data used and how they are processed, the results show a high standard deviation.
In addition, these flow differences mean a distinct behavioral change in this gorge. As shown in Figure 7, it functions either as a braided channel with emerged bars (~50 m 3 s −1 ) or as a single channel occupying the whole river bed (~1250 m 3 s −1 ). The difference in channel behavior is less evident (Figure 7) if we compare the flows that minimize RMSE value versus dendro-geomorphological data height (~1250 m 3 s −1 ), and those that are associated with the 500-years return period (510 m 3 s −1 , from 500-years return period rainfall according to the GP-POT distribution function). In these cases, a single channel occupying the whole river bed is the main river configuration. Figure 8 shows the clear differences in flow-depth data in this channel for each peak flow considered.  Therefore, the peak flow results obtained using different types of information show significant differences when considering dendro-geomorphological evidence versus rainfall data (processed to maximize precipitation intensity); and, these differences become extreme if peak flows obtained using dendro-geomorphological evidence are compared with those obtained from a classic statistical maximum frequency analysis of 24 h rainfall events. Therefore, the peak flow results obtained using different types of information show significant differences when considering dendro-geomorphological evidence versus rainfall data (processed to maximize precipitation intensity); and, these differences become extreme if peak flows obtained using dendro-geomorphological evidence are compared with those obtained from a classic statistical maximum frequency analysis of 24 h rainfall events.

Uncertain Sources
Hydrological models, such as HEC-HMS, are inherently uncertain, due both to the random nature of floods and to the uncertainty component of epistemic origin (due to limited data and knowledge; [50]). In this context, model calibration and validation enable an understanding and mitigation of the model uncertainty, by iterative variation of parameter values within a predefined range. However, this proved impossible in this case due to the lack of observed flow data. The above limitation is mitigated to some extent by an in-depth knowledge of the basin physiography and hydrological response. Flash floods in mountain basins may trigger the mobilization of varying amounts of solid load (i.e., sediment and woody debris), modifying the channel dynamics and morphology [51]. In this respect, the HEC-HMS version here obtains hydrographs that only represent clear water flow as output data. This means that discharges calculated in this way may be underestimated, and therefore cannot be used to assess risks [52], since estimated hazard parameters (i.e., depth and velocity) obtain values lower than they should be.
Ruiz-Villanueva et al. [51] present the differences in the precipitation curve recorded in a small mountain basin from the effect of altitude and slope orientation, concluding that these factors play an essential role in the basin precipitation distribution, and that the curves generated (positive with altitude) have a direct influence, which explains total precipitation values and gauged flow rates in the basin.
On the other hand, applying frequency distribution functions to extrapolate the maximum precipitation values associated with extreme quantiles (with very low occurrence probability) always generates a degree of uncertainty in the results obtained. This uncertainty (epistemic uncertainty; [53]) is the result of the insufficient length of the data series used as the basis of these predictions; so that the insufficiently representative data series, loss of high magnitude records, measurement errors, undetected anomalous values, etc., may lead to significant variations in the values that are associated with the lowest frequency quantiles. Using threshold exceedance series, thus incrementing the length of the data series and enabling better fit of the high range of the frequency function [34,35], was the method chosen in an attempt to limit the uncertainties linked to these statistical analyses.

Uncertain Sources
Hydrological models, such as HEC-HMS, are inherently uncertain, due both to the random nature of floods and to the uncertainty component of epistemic origin (due to limited data and knowledge; [50]). In this context, model calibration and validation enable an understanding and mitigation of the model uncertainty, by iterative variation of parameter values within a predefined range. However, this proved impossible in this case due to the lack of observed flow data. The above limitation is mitigated to some extent by an in-depth knowledge of the basin physiography and hydrological response. Flash floods in mountain basins may trigger the mobilization of varying amounts of solid load (i.e., sediment and woody debris), modifying the channel dynamics and morphology [51]. In this respect, the HEC-HMS version here obtains hydrographs that only represent clear water flow as output data. This means that discharges calculated in this way may be underestimated, and therefore cannot be used to assess risks [52], since estimated hazard parameters (i.e., depth and velocity) obtain values lower than they should be.
Ruiz-Villanueva et al. [51] present the differences in the precipitation curve recorded in a small mountain basin from the effect of altitude and slope orientation, concluding that these factors play an essential role in the basin precipitation distribution, and that the curves generated (positive with altitude) have a direct influence, which explains total precipitation values and gauged flow rates in the basin.
On the other hand, applying frequency distribution functions to extrapolate the maximum precipitation values associated with extreme quantiles (with very low occurrence probability) always generates a degree of uncertainty in the results obtained. This uncertainty (epistemic uncertainty; [53]) is the result of the insufficient length of the data series used as the basis of these predictions; so that the insufficiently representative data series, loss of high magnitude records, measurement errors, undetected anomalous values, etc., may lead to significant variations in the values that are associated with the lowest frequency quantiles. Using threshold exceedance series, thus incrementing the length of the data series and enabling better fit of the high range of the frequency function [34,35], was the method chosen in an attempt to limit the uncertainties linked to these statistical analyses.

Discussion
One of possible reasons for the general non-coincidence between the event dates documented or in the systematic record of intense precipitations, and those deduced from the dendrochronological analysis of samples is the non-linear nature of the rainfall-runoff process and its effect on riparian vegetation.
Thus, the most intense precipitation events in the meteorological station do not necessarily produce the highest flood flow rates. The spatial-temporal variability between precipitation and flow rates is influenced by the local or general character of the maximum precipitation intensities, so that the hydrological response may or may not be limited to one sector of the basin; and, by the ground moisture content previous to the most intense precipitation event, conditioned by the preceding precipitation (previous three days, month, or even hydrological year), which strictly controls the runoff threshold, and therefore the runoff coefficient. Thus, the possible combination of the factors mentioned above could justify the considerable discrepancy found in the model results.
For this reason, peak flow assessment in ungauged basins has been approached from different points of view. The most widely used method is probably regionalization [54][55][56], where statistical methods are used to transfer the hydrological characteristics of gauged basins located in the same area as the ungauged basin in question. The use of images and data collected by different satellite sensors has also been considered [57,58], and has mainly proved effective for large rivers (even taking into account the limitations presented by the use of satellite images [59]). Other authors, such as Ballesteros-Cánovas et al. [60][61][62], show that dendro-geomorphological information is useful for estimating these peak flows. The method proposed by Diakakis [63], uses the fluvial morphometry for peak flow assessment in ungauged basins. Another interesting approach use Bayesian statistic (Monte Carlo analysis) for hydrological modeling [20], allowing for an uncertainty estimation of peak flow results, which improve information for manage makers.
Based on the available information for the study area, this analysis follows the Ballesteros-Cánovas proposal, using dendro-geomorphological data for peak flow assessment in ungauged basins. However, here, this dendro-geomorphological data is used first to calibrate the results obtained from the hydrological models, which are based on statistical processing of available rainfall information.
Our aim was to see if this processing of rainfall data could be considered as a new way to estimate peak flows in ungauged basins, where all of the approaches mentioned above do not offer optimum results or are not applicable. We must consider that these approaches to peak flow assessment in ungauged basins require various indispensable information sources: gauged basins with similar hydrological characteristics; sufficiently high-quality satellite images, where the vegetation does not limit interpretation; and, dendro-geomorphological data associated with the flood event occurrence; or a compilation of locations damaged by flood events.
Another approach to using rainfall data in ungauged basin analysis is atmospheric event modeling [64], although the results are not conclusive. This method obtains the space-time rainfall distribution, which is an advantage for hydrological modeling, but offers low spatial resolution for small basin analysis, and may underestimate rainfall volume, which are disadvantages of this approach.
Early warning systems (EWS) may also be considered in flood hazard mitigation in ungauged basins. An alert is activated when a rainfall event exceeds a pre-defined threshold. However, EWS present weaknesses in this study area: absence of automatic rainfall gauging stations in the upper part of the basin, and a very limited basin lag time (about 1 h).
It is important to remember here that the study area is located in the exceptional environment of a National Park with an important influx of visitors. Earlier studies by the USGS [14][15][16] opted for regionalization methods that are based on regional equations obtained by multiple regressions of data from gauged basins. These papers, dating back to the 1980s and 1990s, are prior to the expansion and generalization of the use of paleostage indicators. The analysis presented in this paper could therefore be considered as reactivating and updating studies of peak flows in ungauged basins located in exceptional environments that are used by the general public, such as National Parks.

Flooded Areas in Caldera de Taburiente N.P. Past Flood Events
Combining the use of rainfall data and dendro-geomorphological evidence with detailed hydrological and hydraulic models enabled us to describe the hydrodynamics of the 1997 flash flood event analyzed here. This is of vital importance in flood hazard assessment and for generating the relevant maps to determine the zones that are affected, or unaffected, by flash floods (Figure 7). For example, taking into account the particular characteristics of the study area (a National Park used by a large number of hikers), delimiting risk areas where visitors may suffer loss of stability due to flow [65,66], will be very important for the park managers (even more in an area with 65,000 visitors each year and previous human dies, isolated visitors on the campsite, and trekking paths damages by flood events). The results of the hydraulic models (Figure 9a) show considerable differences in terms of high risk areas. Thus, in the 50 m 3 s −1 peak flow model, the high-risk areas in the "Playa de Taburiente" are marginal, limited to the most highly developed channels within a braided channel pattern. In the 510 m 3 s −1 peak flow model, the high-risk areas are much more extensive, with all of the channels and lower topographic development bars included in the high-risk zone. Finally, in the 1250 m 3 s −1 peak flow model, the high-risk zones represent about 85% of the total area, defining a main central channel with important secondary channels on each side.
Geosciences 2018, 8, 6 18 of 23 Combining the use of rainfall data and dendro-geomorphological evidence with detailed hydrological and hydraulic models enabled us to describe the hydrodynamics of the 1997 flash flood event analyzed here. This is of vital importance in flood hazard assessment and for generating the relevant maps to determine the zones that are affected, or unaffected, by flash floods (Figure 7). For example, taking into account the particular characteristics of the study area (a National Park used by a large number of hikers), delimiting risk areas where visitors may suffer loss of stability due to flow [65,66], will be very important for the park managers (even more in an area with 65,000 visitors each year and previous human dies, isolated visitors on the campsite, and trekking paths damages by flood events). The results of the hydraulic models (Figure 9a) show considerable differences in terms of high risk areas. Thus, in the 50 m 3 s −1 peak flow model, the high-risk areas in the "Playa de Taburiente" are marginal, limited to the most highly developed channels within a braided channel pattern. In the 510 m 3 s −1 peak flow model, the high-risk areas are much more extensive, with all of the channels and lower topographic development bars included in the high-risk zone. Finally, in the 1250 m 3 s −1 peak flow model, the high-risk zones represent about 85% of the total area, defining a main central channel with important secondary channels on each side.  It has been shown that when meteorological information alone is used, even when it is processed to maximize both precipitation intensity and volume in a theoretical storm event, the peak flows obtained are still far lower than when dendro-geomorphological data is included in the analysis. This is the case for both the hydrodynamic aspects (peak flow, water depth, flow velocity) and the related flood risk. Thus, improved spatial probability of the areas affected when flash floods occur in the zone can be obtained by including dendro-geomorphological data in the analysis, taking into account the inherent uncertainties in the limited data available to compare and contrast the results.
Nevertheless, processing rainfall data to maximize both the intensity and precipitation volume brings the associated peak flow values closer to those that are obtained using dendro-geomorphological information. This is especially true if the results that are offered by this technique are compared with those associated with classic rainfall frequency analysis of maximum precipitation in 24 h. Perhaps the uncertainty associated with this approach to peak flow assessment in ungauged basins does not allow for the correct link between frequency and event. However, using the highest intensity rainfall event calculated for the area can provide the approximate area of the flood-prone areas, thus determining the flood hazard zone. The availability of this simple information (flood hazard zones) can be very useful in management and decision-making in natural environments, such as National Parks (Figure 9b).
The discussion above highlights the dilemma of whether or not these studies of this type are feasible and the results reliable. Even taking this into account, another and perhaps more important question must be asked: When systematic meteorological-hydrological data is not available or very limited, should the study of river dynamics and the associated risks in basins be ruled out? We believe that the results of this study show that the answer here is no, it should not be ruled out. Even with all of the uncertainties raised, the combined use of all available sources of flood event data enables an approximation of the fluvial dynamics and associated flood risks to be obtained. More importantly, if the answer to this question is yes, then it should be ruled out, as this research is not feasible, then most river basins worldwide would not be studied, due to the very limited statistically consistent hydrological data available.

Conclusions
The combined and successive use of systematic meteorological and dendro-geomorphological data sources with hydrological analysis, rainfall-runoff models, and hydraulic models has enabled improved flow magnitude assessment for the 1997 flash-flood event in the Playa de Taburiente.
The main findings and contributions of this paper are: -Work undertaken in areas with very scarce availability of systematic information, and disparity in time-space data resolution. Thus, the proposed methodology can be applied to other areas with similar characteristics, a high proportion of drainage basins in mountain areas worldwide (especially outside Europe and the US). -Different techniques and methods (hydrological, hydraulic, dendro-geomorphological ...) have been combined to make best use of the scarce data available. -Most attention has been paid to the rainfall data, processed first in a classic statistical frequency analysis of 24 h rainfall data. This rainfall data was subsequently re-processed (duration, altitudinal gradient, hyetograph bias) to generate maximum intensity and volume events. -Is has been noted the results differences derived from a classical rainfall-runoff process modelling (limited by ungauged fluvial basin site place), and those obtained through the incorporation into the analysis of re-worked rainfall data (looking for maximum rainfall intensity events). - The incorporation of indirect flood dendro-evidence (FDEs) has been included in peak flow assessment, and to calibrate the results obtained previously from statistical analysis of the rainfall data.
-The peak flow results obtained from the rainfall data processed to maximize intensity are lower than those obtained using FDE, but enable an approximate definition of the flooded area, and therefore the flood hazard area. -These results show the high variability in estimated discharge and water-stage values in study areas where data is scarce, and how this variability, when transferred to risk assessment, can generate important errors and uncertainties.
This paper is a new example of the feasibility of using dendro-geomorphological data as PSI, and of its evident potential usefulness in peak flow assessment in both ungauged and gauged basins. In addition, the study proposes an alternative method to assess peak flows and flood hazard zones in areas where the information available is limited to scarce and inconveniently located rainfall data. Looking for the most favorable events (in terms of rainfall intensity and volume) for flash-flood generation, has allowed for us to estimate peak flows and obtain flood hazard distribution. This information, in areas with very scarce data availability, can be useful for decision-making and management in the context of public natural environments such as National Parks.