Improving Flood Maps in Ungauged Fluvial Basins with Dendrogeomorphological Data. An Example from the Caldera de Taburiente National Park (Canary Islands, Spain)

Flash floods represent one of the more usual natural hazards in mountain basins, and, combined with the lack of reliable flow data and the recreational use of the drainage basin by tourists and hikers, there is a significant risk of catastrophe. Here, we present a dendro-geomorphological reconstruction of a past flash flood event in the Caldera de Taburiente N.P. (Canary Islands, Spain), an ungauged drainage basin in the SW side of the volcanic island of La Palma. We couple two-dimensional hydraulic modelling in a highly-resolved topographic environment (LiDAR data) with (1) peak flow data for various Tyear return periods from an uncalibrated hydrological model and (2) a data set of scars on trees, to investigate the magnitude of a 1997 dated flash-flood. From the results, flood hazards and associated risks would be clearly underestimated by using only the unique available hydrological data (a rainfall gauge station downstream of the study area). Hydraulic models using scars data show a higher flood hazard scenario, improving the flood hazard map by using all available flood evidence. Moreover, all this will allow for better implementation of appropriate adaptation policies by National Park managers, and therefore the mitigation of future disasters.

Improving Flood Maps in Ungauged Fluvial Basins with Dendrogeomorphological Data.An Example from the Caldera de Taburiente National Park (Canary Islands, Spain)

Introduction
Historical flood data and the reconstruction of past floods have previously been used for extreme flood events and flood hazard assessment [1][2][3].The reconstruction and mapping of past floods in ungauged rivers would lead to a more complete understanding of fluvial dynamics and magnitudes, especially since reconstruction mapping places more emphasis on ascertaining the greatest flood magnitudes along a stretch of river.The reconstruction of past floods makes it possible to map hazard areas related to a given flood magnitude and, in some cases, to a known frequency, which will improve the flood risk assessment for the area.Flood risk analysis and management is conditioned by data availability.Statistical analysis (FFA; flood frequency analysis) of annual maximum peak flow can be used to get flow quantile estimations when the time series are statistically significant.When flow data are not representative, rainfall-runoff models can be used if there are enough precipitation data.These hydrological models provide output data, which are useful to determine the flood hazard area using hydraulic models [4].However, there is frequently an absence of both flow and precipitation data for mountain basins, because either the data simply do not exist or what is available has a low spatio-temporal significance.As a result, indirect flood characterization methods are increasingly being used as an alternative [5].In this context, different palaeostage indicators (PSI) and high water marks (HWM) have been used to characterize past floods, forming the basis of palaeohydrology [6] and post-event flood analysis [7,8].Indirect data sources can be grouped as (i) geological/geomorphological; (ii) documentary; and (iii) botanical.
Geological and geomorphological data sources include sediment deposits and erosion lines or marks, enabling the reconstruction of peak flow data and associated parameters by using empirical formulas or mathematical models [9].Documentary methods are based on the documentation available in historical archives.Permanent high-water flood level markers on buildings can be used as additional data sources.These enable the reconstruction of the flooded area and therefore of the flood depths reached by integrating indirect evidence data of floods into hydraulic models.The levels inferred from this data can then be transformed into flows and integrated into the systematic records as non-systematic data, thereby reducing frequency analysis uncertainty [2,10].Finally, indirect evidence of flood occurrence can also be recorded from riparian vegetation, using a dendro-geomorphological approach [11], in which these data are used in a similar approach to that used with documentary data.It also reduces uncertainty in hazard and risk analysis, since roughness and other sensitive parameters can be calibrated [12].
As part of the latter approach, the use of the information contained in tree rings has been shown to be an important source of evidence of past flood events.Ballesteros-Cánovas et al. [13] and Benito and Díez-Herrero [14] provide an extensive review of the main contributions and future prospects of tree ring data analysis in flood record assessment.In addition to dating procedures from tree ring information [15][16][17][18][19][20], they have been used successfully for the estimation of flood magnitude in combination with paleo-hydraulic techniques [21][22][23].Dendro-geomorphological data has also proved very useful for flood reconstruction in mountain areas (where data availability is scarcer), which, as a rule, are more prone to flash-floods [24][25][26], increasing the flood hazard of those events.
Therefore, the main objective of this paper is to evaluate and analyse the differences and possible improvements in the estimation of flood hazard zones through flood maps obtained from the use of dendro-geomorphological information.The singular characteristics of the study area (a national park with an average annual influx of 65,000 visitors) make it necessary to precisely delimitate the risk levels in those areas in which people, goods, and services are exposed to flood hazards.The results are expected to provide a sounder basis for defining flood hazard zones in an area characterized by intense tourist activity [27].To achieve our objective, we have endeavored to reconstruct a 1997 dated flood event from both the limited rainfall and flow information available for the area [28] and the available dendro-geomorphological information [29].The results obtained through the implementation of both data sets within bidimensional hydraulic models will allow us to define different flood hydraulic parameters, which will play the role of a classification criterion for flood hazard level spatial distribution in the study area.Finally, validity and uncertainties in the results are theoretically considered, as well as the hydraulic model's sensitivity to value ranging in some control variables.

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 situated centrally 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 a.s.l.(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 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 [30].This research focuses on the central area of the park, known as 'Playa de Taburiente', a Y-shaped 2 km long (gravel bed braided type) river reach, formed by two gorges, the 'Verduras de Alfonso' and 'Cantos de Turugumay' ephemeral rivers, and the Taburiente River as the result of their confluence (Figure 1).

Greenwich Meridian
).The altitude ranges from 2426 m a.s.l.(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 [30].This research focuses on the central area of the park, known as 'Playa de Taburiente', a Y-shaped 2 km long (gravel bed braided type) river reach, formed by two gorges, the 'Verduras de Alfonso' and 'Cantos de Turugumay' ephemeral rivers, and the Taburiente River as the result of their confluence (Figure 1).

Dendrogeomorphological Data Sources
Dendro-geomorphological fieldwork consisted of collecting samples of damaged trees (i.e., those presenting scars, exposed roots, resprouting, or apical bud loss and dead trees; Figure 2), other nearby trees with no visible injuries (as scars may have already closed if they are old or not very wide), and undamaged trees of P. canariensis.In some cases, trees showed several wounds corresponding to different events, but more often only one impact signal per tree was observed.Finally, we collected 233 samples from 60 disturbed trees and 37 samples from 16 undisturbed trees to provide a reference chronology [29].

Dendrogeomorphological Data Sources
Dendro-geomorphological fieldwork consisted of collecting samples of damaged trees (i.e., those presenting scars, exposed roots, resprouting, or apical bud loss and dead trees; Figure 2), other nearby trees with no visible injuries (as scars may have already closed if they are old or not very wide), and undamaged trees of P. canariensis.In some cases, trees showed several wounds corresponding to different events, but more often only one impact signal per tree was observed.Finally, we collected 233 samples from 60 disturbed trees and 37 samples from 16 undisturbed trees to provide a reference chronology [29].Two or more increment cores (up to ten in some cases) were taken with an increment borer near the tree base if the scar was close to the ground or, more often, at breast height (approximately 1.30 m above the ground).At least one core was taken in the non-affected tissue, and one or more cores were extracted successively closer to the edge of the scar, until one core intersected the scar edge or showed a distinctive growth pattern to identify the year when the scar was produced [17].In addition, some wedges were taken from the overgrowing callus with a handsaw, and whole sections were cut from two dead trees.Each sampled tree was recorded in a spatial database with additional information including (i) scar size and other evidence of disturbance; (ii) tree height, tree diameter at breast height (DBH), and other parameters (Figure 3); and (iii) sketches, drawings, and photographs.Two or more increment cores (up to ten in some cases) were taken with an increment borer near the tree base if the scar was close to the ground or, more often, at breast height (approximately 1.30 m above the ground).At least one core was taken in the non-affected tissue, and one or more cores were extracted successively closer to the edge of the scar, until one core intersected the scar edge or showed a distinctive growth pattern to identify the year when the scar was produced [17].In addition, some wedges were taken from the overgrowing callus with a handsaw, and whole sections were cut from two dead trees.Each sampled tree was recorded in a spatial database with additional information including (i) scar size and other evidence of disturbance; (ii) tree height, tree diameter at breast height (DBH), and other parameters (Figure 3); and (iii) sketches, drawings, and photographs.Once the samples had been dried and prepared, rings were measured using the Lintab measuring system and associated Tsap software [31].For the cross sections or wedges, a minimum of four radius measurements per sample was made.Tree ring series were crossdated using classical methods, including visual, graphic, and statistical techniques [32], and checked using Cofecha software [33].
After successfully crossdating the series within each tree and within each cluster of trees (according to their location and age), injuries were dated in most cases, checking the series from affected and non-affected tissues.Finally, synchrony between the tree ring series of damaged trees and the reference chronology have ensured the accuracy/quality of the wound dating.Growth releases and suppressions appearing synchronously in each group of trees were detected using the Jolts program [34] and were used as other indicators of the injury dates.
63 wounds were dated, in some cases more than once for each tree.Especially relevant were the injuries occurring in 1962 and in 1997 (Figure 3), both presenting a large number of replications (48% and 22%, respectively, of the total number of trees).Wounds dated to 1993 and 2003 were, respectively, found 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 [29].

Hydrological Data
There is no flow data available for past flood events in the basin; it is therefore impossible to apply flood frequency analysis (FFA) using discharges.As a result, the flow data for the "Playa de Taburiente" reach came from a hydrological flood analysis that was carried out using rainfall-runoff models [28].The available rainfall information was statistically treated to determine the precipitation volumes associated with different return periods and durations; the authors also considered an altitudinal gradient for the intensity and precipitation volumes, due to the unfavourable location of the rainfall gauging stations (all located in the middle-lower part of the basin, without information Dendro-geomorphological data of interest for hydraulic assessment (B).Z s , scar topographic elevation; Z w , water topographic elevation during event; Z g , present day ground elevation; h s , scar height; h w , water depth during flood event regarding present day topography; D, height difference between scar and water surface (deviation).
Once the samples had been dried and prepared, rings were measured using the Lintab measuring system and associated Tsap software [31].For the cross sections or wedges, a minimum of four radius measurements per sample was made.Tree ring series were crossdated using classical methods, including visual, graphic, and statistical techniques [32], and checked using Cofecha software [33].
After successfully crossdating the series within each tree and within each cluster of trees (according to their location and age), injuries were dated in most cases, checking the series from affected and non-affected tissues.Finally, synchrony between the tree ring series of damaged trees and the reference chronology have ensured the accuracy/quality of the wound dating.Growth releases and suppressions appearing synchronously in each group of trees were detected using the Jolts program [34] and were used as other indicators of the injury dates.
63 wounds were dated, in some cases more than once for each tree.Especially relevant were the injuries occurring in 1962 and in 1997 (Figure 3), both presenting a large number of replications (48% and 22%, respectively, of the total number of trees).Wounds dated to 1993 and 2003 were, respectively, found 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 [29].

Hydrological Data
There is no flow data available for past flood events in the basin; it is therefore impossible to apply flood frequency analysis (FFA) using discharges.As a result, the flow data for the "Playa de Taburiente" reach came from a hydrological flood analysis that was carried out using rainfall-runoff models [28].The available rainfall information was statistically treated to determine the precipitation volumes associated with different return periods and durations; the authors also considered an altitudinal gradient for the intensity and precipitation volumes, due to the unfavourable location of the rainfall gauging stations (all located in the middle-lower part of the basin, without information about rainfall in the upper part where the study area is located).However, as the authors themselves indicate, the absence of flow data prevents any type of hydrological model calibration, beyond estimating their sensitivity to some parameters such as the Soil Conservation Service-Curve Number (SCS-CN) or the duration and intensity of rainfall events.The hydrological model was developed using HEC-HMS software.The use of simulation methods for the various processes involved was based on the information available for the basin.So, the loss method used was SCS-CN; the precipitation-runoff transformation method used in the basin model was 'Clark unit hydrograph', and the flood wave routing was implemented using the Muskingum-Cunge method.Due to the absence of any other type of flow information for the study area, those presented by Garrote et al. [28] will be used for their implementation in the hydraulic models to be developed.

Hydraulic Modelling
The hydraulic modelling was carried out with numerical simulation using a finite differences model (St Venant 2D equations) and the Iber hydraulic simulation free software [35], whose characteristics are outlined in the following paragraphs.
The boundary conditions related to inflows were restricted to the peak flow, since they are the most reliable flow data available for the study area.Thus, two flow groups were modelled: -Peak flows associated with precipitations recorded during the 11-13 January 1997 event, obtained from Garrote et al. [28].-Peak flow range (300-2500 m 3 s −1 ), seeking to minimize errors (RMSE) when they are related with the dendro-geomorphological 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 was included in a 1 m × 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 and downstream in the Taburiente River.The altimetric accuracy of the data has an estimated value of z ≤ 0.2 m.
A theoretical second representation of the relief was generated in an attempt to limit possible differences between the date of the LIDAR topography data collection (2009) and channel geometry at the date of the event analysed (1997).The aim of this secondary DEM was to assess possible variations in the bed topography due to the mobile bed characteristics of the study area and the time gap between LiDAR data and flash flood event.
This second topography (referred to below as v1997) presents surface modifications of the riverbed (compared with the original, v2009) within the distribution area of disturbed trees.The geomorphological criterion used to include these modifications was as follows: during the initial growth period, the trees must have been outside the flood prone area; therefore, the topographic elevation of the tree base can be used as the upper elevation limit of the floodplain and riverbed at that point.To generate the v1997 DEM, we combined two topographic surfaces: the current topography (v2009) and the best adjust of a 3D surface derived from the topographic position of the disturbed trees.Conservation of the highest elevation between these two surfaces was the criterion used for the surface implementation.Based on that criterion, the resulting combined surface retained the current morphology of the slopes outside the channel, while, in the streambed and banks, the surface was related to the elevation of the dendro-geomorphological data.
Two three-dimensional meshes were generated from these topographic models, formed by 1, 3, and 5 m non-structured elements, depending on their location.The 1 m elements were used for the areas with disturbed trees, with 3 m elements for the upstream and downstream areas and 5 m elements for the valley slopes.
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 selected in the channel of the Taburiente, approx.500 m downstream from the area of the disturbed trees.
Initial conditions for the hydraulic model were defined by a dry scenario, since in natural conditions the river is ephemeral.The whole basin has a concentration time of 4.5 h, but this time is reduced to 1.5 h for the "Playa de Taburiente" outlet.
The last boundary condition of the model to note is the surface roughness of the terrain, expressed by the Manning n coefficient.This variable was derived from the photointerpretation of 0.5-m spatial resolution aerial orthoimages and with field work.Different areas with homogeneous relief roughness and vegetation density were defined, and the corresponding Manning's n value (Figure 4) was assigned to them following the most common tables and graphic guides [36,37].Due to the flash flood event date, other approaches to generating surface roughness maps using satellite images [38] were not feasible due to the low spatial resolution of the images (the channel width was covered by only 3-5 pixels).Two surface roughness terrain models were developed in order to analyse the model sensitivity to Manning roughness coefficients, due to the absence of high resolution orthoimages (pixel size lower than 0.5 m) for the 1997 flash-flood event.Orthoimage dates used for sensitivity analysis were 2002 (the date closest in time to the flood event) and 2013 (the most current orthoimage for the study area).This sensitivity analysis allows us to check the possible influence of surface roughness variations over time in peak flow generation.Initial conditions for the hydraulic model were defined by a dry scenario, since in natural conditions the river is ephemeral.The whole basin has a concentration time of 4.5 h, but this time is reduced to 1.5 h for the "Playa de Taburiente" outlet.
The last boundary condition of the model to note is the surface roughness of the terrain, expressed by the Manning n coefficient.This variable was derived from the photointerpretation of 0.5-m spatial resolution aerial orthoimages and with field work.Different areas with homogeneous relief roughness and vegetation density were defined, and the corresponding Manning's n value (Figure 4) was assigned to them following the most common tables and graphic guides [36,37].Due to the flash flood event date, other approaches to generating surface roughness maps using satellite images [38] were not feasible due to the low spatial resolution of the images (the channel width was covered by only 3-5 pixels).Two surface roughness terrain models were developed in order to analyse the model sensitivity to Manning roughness coefficients, due to the absence of high resolution orthoimages (pixel size lower than 0.5 m) for the 1997 flash-flood event.Orthoimage dates used for sensitivity analysis were 2002 (the date closest in time to the flood event) and 2013 (the most current orthoimage for the study area).This sensitivity analysis allows us to check the possible influence of surface roughness variations over time in peak flow generation.Flow depth, specific flow, absolute elevation level, and velocity field variables were obtained from the hydraulic simulations.Hydraulic model results were implemented in a GIS (ArcGIS 10.3.ESRI Geosystems) for further analysis, among them flow depth, flood prone area, and flow velocity differences in space between different peak flow hydraulic models, or the product of flow depth by flow velocity.The flow depth by flow velocity product is linked to the flood hazard for people by loss of stability due to flow [39,40], and it becomes extreme when the product result is higher than 2 m 2 s −1 .Flow depth, specific flow, absolute elevation level, and velocity field variables were obtained from the hydraulic simulations.Hydraulic model results were implemented in a GIS (ArcGIS 10.3. ESRI Geosystems) for further analysis, among them flow depth, flood prone area, and flow velocity differences in space between different peak flow hydraulic models, or the product of flow depth by flow velocity.The flow depth by flow velocity product is linked to the flood hazard for people by loss of stability due to flow [39,40], and it becomes extreme when the product result is higher than 2 m 2 s −1 .

1997 Flood Reconstruction
From the data sources and using the abovementioned methodological approach, we have reconstructed the main 1997 flood characteristics, both the event peak flows and the floodable and hazard maps, and obtained partial conclusions.
The first attempt at hydraulic simulation was carried out from the flows obtained in the hydrological model that used the precipitations registered at the station C106U for the three days between 11 and 13 January 1997.It soon became clear that the Q pm (peak flow value minimizing the RMSE) exceeded the Q r-r (peak flow range value derived by rainfall-runoff transformation [28] for this event).The RMSE value refers to the fit error between the heights of the scars (dendro-geomorphological evidence, FDEs) and the water surface height associated with flood modelling in those locations.To determine this peak flow value, a simulation was then run using the estimated 500-year return period discharge value as the second peak flow value.The results of the modelling continued to show a high RMSE, and the variable flow depth obtained was always below the variable scar height.From these results, simulation using different increasing peak flow values was developed until a lower limit value of the RMSE was obtained.This value was established (Figure 5A) as a peak flow of 1235 m 3 •s −1 , and it does not appear to be driven by topographic and surface roughness variables, as the sensitivity analysis results show (Figure 5).

1997 Flood Reconstruction
From the data sources and using the abovementioned methodological approach, we have reconstructed the main 1997 flood characteristics, both the event peak flows and the floodable and hazard maps, and obtained partial conclusions.
The first attempt at hydraulic simulation was carried out from the flows obtained in the hydrological model that used the precipitations registered at the station C106U for the three days between 11 and 13 January 1997.It soon became clear that the Qpm (peak flow value minimizing the RMSE) exceeded the Qr-r (peak flow range value derived by rainfall-runoff transformation [28] for this event).The RMSE value refers to the fit error between the heights of the scars (dendrogeomorphological evidence, FDEs) and the water surface height associated with flood modelling in those locations.To determine this peak flow value, a simulation was then run using the estimated 500-year return period discharge value as the second peak flow value.The results of the modelling continued to show a high RMSE, and the variable flow depth obtained was always below the variable scar height.From these results, simulation using different increasing peak flow values was developed until a lower limit value of the RMSE was obtained.This value was established (Figure 5A) as a peak flow of 1235 m 3 •s −1 , and it does not appear to be driven by topographic and surface roughness variables, as the sensitivity analysis results show (Figure 5).The notable difference in the peak flow values leads to two possible approaches that will be discussed below: the first, that the rainfall record for the event analysed is not correct, either because of measurement errors or because the spatial heterogeneity of the precipitation meant that the amount recorded by the pluviograph was not representative of the rainfall in the upper part of the Caldera; and the second, the role of solid (sediment and wood) flow in the flash floods that occurred The notable difference in the peak flow values leads to two possible approaches that will be discussed below: the first, that the rainfall record for the event analysed is not correct, either because of measurement errors or because the spatial heterogeneity of the precipitation meant that the amount recorded by the pluviograph was not representative of the rainfall in the upper part of the Caldera; and the second, the role of solid (sediment and wood) flow in the flash floods that occurred at the study site.
Another point to highlight in these results is that the hydrological model peak flow values [28] are always higher for the "Cantos de Turugumay" gorge, even considering that the rainfall intensities are similar in the two gorges.However, best fit between hydraulic simulation and dendro-geomorphological data require the peak flow to be greater in the "Verduras de Alfonso" gorge.This would indicate that the rainfall intensities cannot be considered similar in these two basins, because the sizes of the precipitation convective cells are smaller than the dimensions of basins.
Furthermore, the hydraulic model shows that the peak flows associated with the "Verduras de Alfonso" gorge need to be higher than those flowing along the "Cantos de Turugumay" ravine.For the latter, the necessary peak flow (400 m 3 s −1 ) to minimize the differences with respect to the dendro-geomorphological data is similar to that obtained by Garrote et al. [28] for a 200-year return period rainfall (one-hour duration hyetograph).However, taking into account these flows in the "Cantos de Turugumay" ravine, the necessary peak flows in "Verduras de Alfonso" ravine (between 675 and 835 m 3 s −1 ) clearly exceed all those obtained from all hydrological models considered.This situation is similar to that previously noted by Ballesteros-Cánovas et al. [41] in the Kullu district (Western Indian Himalayas), which could be related to the highly localized nature of extreme rainfall events in the headwater catchment of the Caldera de Taburiente N.P.
The limited flooded area linked to the 2D hydraulic modelling peak flows (derived from the rainfall-runoff model without rainfall-altitude gradient applied) is one of the most interesting items of the results obtained for the 1997 flash flood event (Figure 6(A1)).This can be seen from a simple visual analysis of the flash flood prone areas.Thus, for the v2009 topographical model, only 6 of the 11 trees give a depth value, and this is a low value, so that the RMSE obtained is around 1 m compared with the height of the identified disturbances.In the 1997 model, the number and RMSE are the same as in the previous model, although it is not the same set of trees which provide the depth values.
In addition, these flow differences mean a distinct behavioural change in this gorge.As shown in Figure 6(A1), 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 (~1200 m 3 s −1 ).The difference in channel behaviour is less evident (Figure 6(A2)) if we compare both the peak flow obtained using dendro-geomorphological data height (1235 m 3 s −1 ) and that associated with the 500-year return period (510 m 3 s −1 ).In these cases, a single channel occupying the whole river bed is the main river shape.However, significant flow depth differences can be observed between the hydraulic models, which range from 0.25-1 m in the area of "Playa de Taburiente".Water depth differences (Figure 6B) for the entire area covered by the hydraulic model show an average value of 1.2 m.However, this number is clearly affected by the values calculated for the narrower sections of both gullies.For the flow velocity data, the differences obtained between the hydraulic models with 1235 m 3 s −1 and 510 m 3 s −1 peak flows have an average value of 2.6 m/s, although in this case, as with the water depth, the average value is affected by the flow velocities reached in those narrower sections of the gullies, with the average flow velocity in the area of "Playa de Taburiente" being around 2.2 m s −1 .This water depth and flow velocity variability presents us with a situation in which the flow energy increases noticeably even though the area affected by the flood does not present significant topological differences, or not to the point of modifying the channel plan.Thus, it can be observed that the flow velocity in the 510 m 3 s −1 peak flow hydraulic model has a clear distribution of high velocity zones in the channels on a braided riverbed, and low speed zones on the bars of that braided pattern, which are also the low water depth areas.This flow velocity distribution is much less clear in the 1235 m 3 s −1 peak flow model.Here, although velocities are still higher in the channels, velocities on the bars are also high, showing smoother transitions in the velocity distribution.
The last point is of great importance in flood risk assessment and for generating the relevant mapping to determine the zones affected, or unaffected, by flash floods (Figures 6(C1-C3) and 7).For example, taking into account the particular characteristics of the study area (a National Park with a large influx of hikers, over 60.000 people per year), the delimitation of risk areas where people could suffer loss of stability due to flow [39,40] will be of great importance for Park management.The results from the hydraulic models (Figures 6(C1-C3) and 8) show large differences in relation to the high-risk area.Thus, in the 50 m 3 s −1 peak flow model, the high-risk areas in the "Playa de Taburiente" This water depth and flow velocity variability presents us with a situation in which the flow energy increases noticeably even though the area affected by the flood does not present significant topological differences, or not to the point of modifying the channel plan.Thus, it can be observed that the flow velocity in the 510 m 3 s −1 peak flow hydraulic model has a clear distribution of high velocity zones in the channels on a braided riverbed, and low speed zones on the bars of that braided pattern, which are also the low water depth areas.This flow velocity distribution is much less clear in the 1235 m 3 s −1 peak flow model.Here, although velocities are still higher in the channels, velocities on the bars are also high, showing smoother transitions in the velocity distribution.
The last point is of great importance in flood risk assessment and for generating the relevant mapping to determine the zones affected, or unaffected, by flash floods (Figures 6(C1-C3) and 7).For example, taking into account the particular characteristics of the study area (a National Park with a large influx of hikers, over 60.000 people per year), the delimitation of risk areas where people could suffer loss of stability due to flow [39,40] will be of great importance for Park management.The results from the hydraulic models (Figures 6(C1-C3) and 8) show large differences in relation to the high-risk area.Thus, in the 50 m 3 s −1 peak flow model, the high-risk areas in the "Playa de Taburiente" are small and are limited to the most developed river channels within a braided channel pattern.In the 510 m 3 s −1 peak flow model, the high-risk zones are much larger, with all of the channels and the lower bars included in the high-risk zone.Finally, in the 1235 m 3 s −1 peak flow model, the high-risk zones cover about 85% of the bottom valley area, defining a main central channel with considerable secondary channels on each side.are small and are limited to the most developed river channels within a braided channel pattern.In the 510 m 3 s −1 peak flow model, the high-risk zones are much larger, with all of the channels and the lower bars included in the high-risk zone.Finally, in the 1235 m 3 s −1 peak flow model, the high-risk zones cover about 85% of the bottom valley area, defining a main central channel with considerable secondary channels on each side.The variability in the results must be due to the different sources of information used.It has been shown that when only meteorological information is used, even when it is processed with the objective of maximizing both the intensity and the volume of precipitation in a theoretical storm event, the peak flows obtained are still far from those derived when dendro-geomorphological data is included in the analysis.Therefore, the use of indirect flood characterization data (in this case dendro-geomorphological data) has revealed that the magnitude of the 1997 flash-flood would have been underestimated if the assessment was only based on the limited meteorological information available, as previously noted [41].This is the case for both the hydrodynamic parameters (peak flow, water depth, and flow velocity) and the related flood hazard.Thus, an improvement in the spatial probability of the areas affected when flash floods occur in the zone can be established due to inclusion in the analysis of dendro-geomorphological data (with potentially major impacts on Disaster Risk Management-DRM) and to the design of mitigation and adaptation plans [41], taking into account all the uncertainties inherent in the limited data available to compare and contrast the results.The variability in the results must be due to the different sources of information used.It has been shown that when only meteorological information is used, even when it is processed with the objective of maximizing both the intensity and the volume of precipitation in a theoretical storm event, the peak flows obtained are still far from those derived when dendro-geomorphological data is included in the analysis.Therefore, the use of indirect flood characterization data (in this case dendro-geomorphological data) has revealed that the magnitude of the 1997 flash-flood would have been underestimated if the assessment was only based on the limited meteorological information available, as previously noted [41].This is the case for both the hydrodynamic parameters (peak flow, water depth, and flow velocity) and the related flood hazard.Thus, an improvement in the spatial probability of the areas affected when flash floods occur in the zone can be established due to inclusion in the analysis of dendro-geomorphological data (with potentially major impacts on Disaster Risk Management-DRM) and to the design of mitigation and adaptation plans [41], taking into account all the uncertainties inherent in the limited data available to compare and contrast the results.

Uncertainty Sources of Flood Hazard Results
Our data sources and analysis methods have identified several uncertainties from different origins and ranges that can be grouped into four sets: alluvial movable river bed, topography, spatial distribution of disturbed trees, and scar height and dating data.
The movable alluvial river bed, which is linked to processes other than solid load erosion, transport, and deposition during flash flood occurrence, can be considered the first source of error and uncertainty.This solid load modifies the flow hydrodynamics by raising the height of the water sheet and incrementing the density and fluid viscosity, which also increase the flow hydrodynamics' transport capacity.Other processes can be related to the increment of solid load: forming pools, narrowing the effective section, concentrating the flow, and initiating the sudden release of pooled water.As regards the detritic solid load, the floating woody load promotes similar processes, leading to major hydrodynamic and turbulent changes that may be responsible for the generation and position of the dendro-geomorphological evidence itself [42,43].
Therefore, channel erosion-sedimentation processes may lead to topographical variations, which may be significant in the results obtained from hydraulic modelling [44,45], which is why most previous studies have been located over stable bed topography channel reaches [46,47].As the "Playa de Taburiente" river channel was found to have high sediment mobility, a theoretical v1997 DEM was defined with the aim of assessing the influence of this mobile riverbed on the results obtained.However, only slight variations were observed in flow depth values for the different hydraulic models.It was therefore concluded that the possible changes in the river bed topography associated with incision or deposition processes (linked to the liquid and solid flow volumes required to cause the tree disturbances observed) do not affect the results of the analysis.
Moreover, a geomorphological phenomenon that counteracts the uncertainty related to the movable river bed has been observed.It is as if the solid load transport in this reach were to some extent in dynamic equilibrium (there appears to be a compensatory phenomenon), so that the wet

Uncertainty Sources of Flood Hazard Results
Our data sources and analysis methods have identified several uncertainties from different origins and ranges that can be grouped into four sets: alluvial movable river bed, topography, spatial distribution of disturbed trees, and scar height and dating data.
The movable alluvial river bed, which is linked to processes other than solid load erosion, transport, and deposition during flash flood occurrence, can be considered the first source of error and uncertainty.This solid load modifies the flow hydrodynamics by raising the height of the water sheet and incrementing the density and fluid viscosity, which also increase the flow hydrodynamics' transport capacity.Other processes can be related to the increment of solid load: forming pools, narrowing the effective section, concentrating the flow, and initiating the sudden release of pooled water.As regards the detritic solid load, the floating woody load promotes similar processes, leading to major hydrodynamic and turbulent changes that may be responsible for the generation and position of the dendro-geomorphological evidence itself [42,43].
Therefore, channel erosion-sedimentation processes may lead to topographical variations, which may be significant in the results obtained from hydraulic modelling [44,45], which is why most previous studies have been located over stable bed topography channel reaches [46,47].As the "Playa de Taburiente" river channel was found to have high sediment mobility, a theoretical v1997 DEM was defined with the aim of assessing the influence of this mobile riverbed on the results obtained.However, only slight variations were observed in flow depth values for the different hydraulic models.It was therefore concluded that the possible changes in the river bed topography associated with incision or deposition processes (linked to the liquid and solid flow volumes required to cause the tree disturbances observed) do not affect the results of the analysis.
Moreover, a geomorphological phenomenon that counteracts the uncertainty related to the movable river bed has been observed.It is as if the solid load transport in this reach were to some extent in dynamic equilibrium (there appears to be a compensatory phenomenon), so that the wet section area remains constant over time.This process implies that the position of the erosive and sedimentary forms changes in the cross section over time, but not the channel wet section area.Additionally, this effect is not limited to the cross-section scale, but also operated at the "Playa de Taburiente" scale.This means that the dendro-geomorphological evidence (such as tree scars), although they cannot be used in this case study for a precise flow estimate, give an idea of the order of magnitude of the flow volumes (water plus sediment), over time, for a given section.
Topographic data precision is another source of uncertainty to be considered, as it plays a key role in all hydraulic modelling and controls how representative the results actually are [48][49][50].If the DEM used in the hydraulic model is not able to represent the terrain morphology, the reliability of the results obtained is questionable [51].Our model is adapted to the conclusions by Casas et al. [50], in the use of LIDAR data and 3D mesh element size for the topographic representation of the terrain.As with topographic data, the spatial resolution of the different variables used in hydrological and hydraulic modelling is another source of uncertainty [52,53].Although it has been impossible to eliminate the uncertainty related to this phenomenon, attempts have been made to limit it as far as possible by using the topographic information with the highest resolution available for the area (1 m pixel size), in both the hydraulic analyses and later in the topographic analysis.
Another source of uncertainty may be related to changes in surface roughness over time, which could be caused by the existence of a movable riverbed, although it can be observed from the aerial orthoimages used (dated in 2002 and 2013) that spatially and temporally most of the channel is occupied by boulder bars.Orthoimages taken at earlier dates, particularly in 2002, show small variations in the extent of the willow forests near the National Park services centre.However, the hydraulic models derived from the two different surface roughness maps do not show significant variations in flow depth values (Figure 5A), so the uncertainty related to changes in the surface roughness due to movable bed may be considered negligible.
Another source of uncertainty which cannot be ignored refers to the spatial distribution and location of the disturbed trees, which were located through a topographic survey.However, due to the abrupt relief of the area, any minimal displacement or uncertainty in the real location of the disturbed trees may produce significant differences in its topographic elevation.Therefore, these minimal displacements may similarly produce significant differences in hydraulic results.To limit this uncertainty when analysing the match between the height of the disturbances and of the water sheet, a 3 × 3 m buffer was taken into account around the position of the disturbed tree (Figure 9), and its 9 elevation values were used.The results obtained lead to the following conclusion: although the associated RMSE values (Figure 5A) vary according to the topographical values used inside the buffer (minimum, maximum, and tree location values), these do not modify the peak flow value associated with the lower RMSE.
Finally, the uncertainty produced by the dendro-chronological methods must also be considered.Up to 63 different wounds were dated in the study area corresponding to 8 event years, but the dendro-chronological methods may possibly not have recorded all the flash floods that occurred [29].There is no guarantee that all flood events had been dendro-geomorphologically dated [24,54] because, for example, the injuries may be masked by new tree tissues formed after the flash flood occurrence if the wound was old or not extensive [55], or the damaged trees may have disappeared, swept away by flash floods after uprooting [29].With the aim of reducing the uncertainties linked to the scar height, Ballesteros-Cánovas et al. [46] used a size-gradation classification of scars to estimate peak discharge in ungauged mountain streams.Due to the low number of reliable disturbed trees observed for each flash-flood event, the size-gradation classification is not available in our study area, so the sampling strategy followed was exclusively focused on tree scars in the flow direction of flash floods [56,57].
The validity of the results obtained and their application are determined by these previously noted problems and uncertainties.These uncertainties are not exclusive to the case study of this reach in the Caldera de Taburiente National Park; on the contrary, they are common to all gorges and ravines in the Canary Islands (with limited or no gauging), to most of the Macaronesian archipelagos, and to many mountain basins in the Iberian Peninsula, especially ungauged reaches [27,57].It must be taken into account that the lack of gauge data or the possibility of obtaining quantiles from calibrated rainfall-runoff models are the basis for most of the uncertainty sources considered.For this reason, and in spite of its limited statistical reliability, the information on flood flow magnitude derived from dendro-geomorphological evidence is justified [41] as the only evidence objectively available and therefore useful for flood risk and hazard assessment.be taken into account that the lack of gauge data or the possibility of obtaining quantiles from calibrated rainfall-runoff models are the basis for most of the uncertainty sources considered.For this reason, and in spite of its limited statistical reliability, the information on flood flow magnitude derived from dendro-geomorphological evidence is justified [41] as the only evidence objectively available and therefore useful for flood risk and hazard assessment. .Sketch of the buffer approach followed for the sensitivity analysis of the peak flow estimation results against the location of dendro-geomorphological evidence.Red lines show the elevation variability inside the buffer area.

Conclusions
The use of dendro-geomorphological data sources with two dimensional hydraulic models has enabled the detection, dating, and improved flow magnitude assessment of the 1997 flash-flood event in the Playa de Taburiente.It has been noted that the results derived from a classical rainfall-runoff modelling process (limited by the fluvial basin's ungauged location) differed from those obtained through the incorporation into the analysis of indirect flood characterization evidence (such as flood dendro-geomorphological evidence-FDE) for the 1997 flash flood event in the Caldera de Taburiente National Park; these differences are significant in terms of estimated peak flows and respective flood-prone areas.The incorporation of indirect flood dendro-evidence (FDEs) has meant an increase in both flood-prone zones and areas where people may suffer loss of stability (based on .Sketch of the buffer approach followed for the sensitivity analysis of the peak flow estimation results against the location of dendro-geomorphological evidence.Red lines show the elevation variability inside the buffer area.

Conclusions
The use of dendro-geomorphological data sources with two dimensional hydraulic models has enabled the detection, dating, and improved flow magnitude assessment of the 1997 flash-flood event in the Playa de Taburiente.It has been noted that the results derived from a classical rainfall-runoff modelling process (limited by the fluvial basin's ungauged location) differed from those obtained through the incorporation into the analysis of indirect flood characterization evidence (such as flood dendro-geomorphological evidence-FDE) for the 1997 flash flood event in the Caldera de Taburiente National Park; these differences are significant in terms of estimated peak flows and respective flood-prone areas.The incorporation of indirect flood dendro-evidence (FDEs) has meant an increase in both flood-prone zones and areas where people may suffer loss of stability (based on flash-flood hydraulic variables, water depth, and flow velocity).Moreover, the inclusions of FDE data into the hydraulic model show a clear change in channel cross section and plant form: from a multi-channel braided river with multiple islands to a single main channel with a few well-developed secondary channels.In all, the results allow for more reliable flood hazard maps to be made and point to the necessary use of all available direct and indirect flood evidence.
These results show the high variability in the estimation of discharge and water-stage values in study areas where data are scarce, and how this variability, when transferred to flood risk assessment, can generate major errors and uncertainties.These results show the relative understanding we possess about flood hazards when we only apply the more traditional methods (such as hydrologic-hydraulic modelling), and that this methodology has a clear potential to be replicated in other fluvial contexts, where it could enhance our knowledge of flood magnitude and frequency.
Furthermore, the use of indirect flood evidence (such as dendro-geomorphological data) improves flood hazard management instruments, especially in areas with scarce available hydrologic data.In particular, in areas such as the one in the present study (mountain National Parks with both, ungauged basins and a high influx of visitors), the use of this indirect flood evidence might improve the human safety in the event of floods.

Figure 1 .
Figure 1.Location map of the "Playa de Taburiente" area and tributaries catchments, within the limits of the Caldera de Taburiente National Park, as well as its main elements (rain gauges, National Park boundary) of interest.

Figure 1 .
Figure 1.Location map of the "Playa de Taburiente" area and tributaries catchments, within the limits of the Caldera de Taburiente National Park, as well as its main elements (rain gauges, National Park boundary) of interest.

Figure 2 .
Figure 2. Mosaic images of some of the damaged trees during the 1997 flash-flood event in "Playa de Taburiente".

Figure 2 .
Figure 2. Mosaic images of some of the damaged trees during the 1997 flash-flood event in "Playa de Taburiente".

Figure 3 .
Figure 3. Map location of "Playa de Taburiente" dendro-geomorphological evidence (A); x-and yaxis units are meters, according to European Terrestrial Reference System 1989 (ETRS89).Dendrogeomorphological data of interest for hydraulic assessment (B).Zs, scar topographic elevation; Zw, water topographic elevation during event; Zg, present day ground elevation; hs, scar height; hw, water depth during flood event regarding present day topography; D, height difference between scar and water surface (deviation).

Figure 3 .
Figure 3. Map location of "Playa de Taburiente" dendro-geomorphological evidence (A); x-and y-axis units are meters, according to European Terrestrial Reference System 1989 (ETRS89).Dendro-geomorphological data of interest for hydraulic assessment (B).Z s , scar topographic elevation; Z w , water topographic elevation during event; Z g , present day ground elevation; h s , scar height; h w , water depth during flood event regarding present day topography; D, height difference between scar and water surface (deviation).

Figure 5 .
Figure 5. (A) Flow depth value for different modelled flows vs. the heights of the dendrogeomorphological evidence: results from the v2009 DEM (left) and v1997 DEM (right); (B) flow depth value difference on tree location due to surface roughness change sensitivity analysis.

Figure 5 .
Figure 5. (A) Flow depth value for different modelled flows vs. the heights of the dendro-geomorphological evidence: results from the v2009 DEM (left) and v1997 DEM (right); (B) flow depth value difference on tree location due to surface roughness change sensitivity analysis.

Figure 7 .
Figure 7. Flood stage graphic near the Information Centre of the Caldera de Taburiente National Park.Shown (A) flood stage associated to recorded rainfall during 1997 event (49.1 m 3 s −1 ; dashed-dot line), flood stage for RMSEmin relative to water surface-scar evidence height (1235 m 3 s −1 ; solid line), and other modelled flows (dashed lines).(B) Location of the cross section (white line) associated with flood stage values, and Caldera de Taburiente National Park Information Center.

Figure 7 .
Figure 7. Flood stage graphic near the Information Centre of the Caldera de Taburiente National Park.Shown (A) flood stage associated to recorded rainfall during 1997 event (49.1 m 3 s −1 ; dashed-dot line), flood stage for RMSE min relative to water surface-scar evidence height (1235 m 3 s −1 ; solid line), and other modelled flows (dashed lines).(B) Location of the cross section (white line) associated with flood stage values, and Caldera de Taburiente National Park Information Center.

Figure 8 .
Figure 8. Extreme hazard areas (EHA, where: v × d > 2 m 2 s −1 ) for people by loss of stability due to flow; each colour area is associated with a different peak flow hydraulic model.X-and y-axis units are meters, according to European Terrestrial Reference System 1989 (ETRS89)

Figure 8 .
Figure 8. Extreme hazard areas (EHA, where: v × d > 2 m 2 s −1 ) for people by loss of stability due to flow; each colour area is associated with a different peak flow hydraulic model.X-and y-axis units are meters, according to European Terrestrial Reference System 1989 (ETRS89)

Figure 9
Figure 9. Sketch of the buffer approach followed for the sensitivity analysis of the peak flow estimation results against the location of dendro-geomorphological evidence.Red lines show the elevation variability inside the buffer area.

Figure 9
Figure 9. Sketch of the buffer approach followed for the sensitivity analysis of the peak flow estimation results against the location of dendro-geomorphological evidence.Red lines show the elevation variability inside the buffer area.