Hydrogeomorphic Impacts of Floods in a First-Order Catchment: Integrated Approach Based on Dendrogeomorphic Palaeostage Indicators, 2D Hydraulic Modelling and Sedimentological Parameters

: Floods represent frequent hazards in both low- and ﬁrst-order catchments; however, to date, the investigation of peak ﬂow discharges in the latter catchments has been omitted due to the absence of gauging stations. The quantiﬁcation of ﬂood parameters in a ﬁrst-order catchment (1.8 km 2 ) was realised in the moderate relief of NE Czechia, where the last ﬂash ﬂood event in 2014 caused considerable damage to the infrastructure. We used an integrated approach that included the dendrogeomorphic reconstruction of past ﬂood activity, hydraulic modelling of the 2014 ﬂash ﬂood parameters using a two-dimensional IBER model, and evaluation of the channel stability using sedimentological parameters. Based on 115 ﬂood scars, we identiﬁed 13 ﬂood events during the period of 1955 to 2018, with the strongest signals recorded in 2014, 2009 and 1977. The modelled peak ﬂow discharge of the last 2014 ﬂood was equal to 4.5 m 3 · s − 1 (RMSE = 0.32 m) using 26 scars as palaeostage indicators. The excess critical unit stream power was observed at only 24.2% of the reaches, representing predominantly bedrock and ﬁne sediments. Despite local damage during the last ﬂood, our results suggest relatively stable geomorphic conditions and gradual development of stream channels under discharges similar to that in 2014.


Introduction
Extreme rainfall resulting in flood events is a common phenomenon in different environments, including both mountain ranges and lowlands [1,2]. The systematic monitoring of flow stages at gauging stations, including precipitation and flow discharge prediction, is currently well applied and documented in many medium-and large-sized rivers [3][4][5]. In contrast, data from mountain steep headwater catchments and, in general, first-order catchments [6] are still poor, due to the insufficient network of stream gauges and the sporadic amount of processed documentary evidence [7]. Not only do mountain headwater streams generate sediment-laden flows with aftermaths within and along fans [8,9], but also streams and gullies of first-order catchments (up to 10 km 2 ) in moderate relief can be responsible for local damage to infrastructure. Ozturk et al. [10] analysed extraordinary flash flood events (140 mm per 2 h) in a small catchment (6 km 2 ) that resulted in damage to infrastructure and a high amount of suspended sediments (t/km 2 ) due to intense hillslope-channel coupling. In addition, Terti et al. [11] pointed to a short response time of small catchments to flash floods, thereby increasing the probability of trapping people during outdoor activities.
Despite the lack of gauging records from forested first-order catchments, several approaches exist to describe the hydrogeomorphic impacts of extremely high flow stages. The dating of flood scars on riparian vegetation using dendrogeomorphic methods [12,13] and recording their position and height above the channel bottom with a combination of hydraulic models is a well-established approach for the estimation of peak flood discharge, flow velocity, and unit stream power [14]. Using flood scars as a palaeostage indicator (PSI; maximum height of the scar above the channel bottom) allowed the interpretation of past flood events in medium-and large-scale rivers (catchment areas larger than 40 km 2 ), for example, in the Western Mediterranean [15], North America [16,17], the Carpathians [18], and the Himalayas [19]. The dendrogeomorphic response of flash floods in small catchments is generally considered lower [20] but may increase due to the presence of erodible sediments amplifying lateral bank erosion. In such conditions, while small streams may not generate flood waves as large as those of large rivers, the presence of exposed and scarred tree roots [21] may complete a relatively low number of scarred tree stems.
The alluvial streams draining first-order catchments often do not display clear relationships between channel geometry, bed substrate, unit stream power and drainage area, and their resulting form and evolution trajectory are unpredictable unless local conditions (e.g., bedrock resistance, intensity of hillslope-channel coupling processes, land use history, and presence of large instream wood) are constrained [22][23][24][25][26]. These streams are characterized by a more or less developed stepped-bed morphology with a wide range of sediment size, where individual steps controlling channel bed stability consist of interlocked boulders, bedrock outcrops or large wood pieces [27]. This implies the relative stability of their channel beds under relatively high discharges (up to floods of 20-50-year recurrence intervals) and thus only limited adjustments of the channel morphology and geometry to lower (e.g., bankfull) flows owing to the presence of generally shallow flows, particle-size interactions (hide/protrusion effect) and additional bed form resistance [25,[28][29][30][31]. The correlations between the unit stream power, sediment calibre and prevailing fluvial process may exist at the reach scale when spatially limited depositional reaches can be accompanied by local bed sediment fining and an abrupt decrease in the unit stream power of a high-magnitude flood [23]. Nevertheless, our knowledge of the direct relationships between the transport capacity of a particular flood event, bed stability and the resulting fluvial processes in first-order catchments is still somewhat limited by the lack of detailed post-flood field surveys of geomorphic consequences.
As mentioned, two-dimensional hydraulic models have been successfully applied to peak discharge reconstructions in ungauged or poorly gauged catchments of various sizes and environments [18,[32][33][34][35]. This peak discharge is mostly defined not as the single deterministic value but rather as a range of values due to uncertainties inherent in the process of the estimation [36]. Discharges predicted in such a way could be used in flood frequency analyses [18], where they serve as outliers to the measured discharges. Consequently, the predicted discharges could serve as a basis for the evaluation and mitigation of flood risk [34] or as information about the flood magnitude in the historical period [32]. Moreover, the results of hydraulic simulations can be used for the estimation of stream transport capacity and channel stability during a particular flood event at a very detailed scale. In this sense, the bed shear stress and unit stream power are relevant parameters for calculating the incipient motion of coarse bed particles and thus evaluating the stability of stepped-bed channels consisting of relatively stable cobble to boulder steps [30,37].
For our purposes, we selected a first-order catchment (1.8 km 2 ) in the moderate relief of the Eastern Sudetes (NE Czechia), where, during the last 15 years, the occurrence of several flash flood events caused substantial geomorphic imprints [38]. In particular, the last intense precipitation event (27 May 2014) resulted in a moderate flood risk in the case of medium-sized rivers (2-year recurrence interval), but a discharge of an approximately 100-year recurrence interval was estimated at ungauged small streams. The short-lived storm, with a total rainfall amount of 40 mm (locally up to 80 mm), had an intensity of between 40 and 60 mm/h. Its hydrogeomorphic response was primarily due to unprecedented rainfall intensity and large antecedent precipitation. The financial costs of the flash flood aftermaths in the affected region were calculated as approximately EUR 200,000 [39].
Our aims were to (i) create the chronology of the flash flood events in this small catchment using dendrogeomorphic approaches, (ii) estimate the parameters (peak flow discharge, flow velocity, bed shear stress, and unit stream power) of the last 2014 flash flood event using the combination of PSI and 2D hydraulic modelling, and (iii) describe the stream transport capacity and channel stability during the 2014 flash flood event based on the hydraulic simulation data and sedimentological parameters. Such a comprehensive approach may help to better quantify the flash flood parameters of ungauged streams and thus contribute to more reliable management of small streams in future extreme climate events.

Study Site
The hydrogeomorphic impacts of floods were studied on a small tributary of the Bělá River in the foreland of the Eastern Sudetes (NE part of Czechia, Central Europe; 50 • 17 N, 17 • 17 E; Figure 1). The area comprises a Proterozoic and Palaeozoic basement composed of orthogneiss, which is overlaid by glacifluvial deposits of the Saalian and Elsterian glaciations with different grain-size distributions of till sediments. These deposits are the result of deglaciation phases in proglacial areas on the ice margins and were recently covered by the Holocene sandy loam colluvium [40,41]. The study site is characterised by a temperate climate, with a mean annual precipitation between 850 and 900 mm. Most precipitation falls during the spring and summer months [42], with an occasional occurrence of extreme daily rainfall (more than 50 mm per day). Documentary evidence provides information about several flood events (e.g., 1903, 1921, 1971, 1977, 2007, 2009, and 2014), relating to both advective rains and short-term intense rainfall events that are currently responsible for intense gully incision and damage to infrastructure within the studied region [38,43].
Water 2020, 12, x FOR PEER REVIEW 3 of 19 due to unprecedented rainfall intensity and large antecedent precipitation. The financial costs of the flash flood aftermaths in the affected region were calculated as approximately EUR 200,000 [39]. Our aims were to (i) create the chronology of the flash flood events in this small catchment using dendrogeomorphic approaches, (ii) estimate the parameters (peak flow discharge, flow velocity, bed shear stress, and unit stream power) of the last 2014 flash flood event using the combination of PSI and 2D hydraulic modelling, and (iii) describe the stream transport capacity and channel stability during the 2014 flash flood event based on the hydraulic simulation data and sedimentological parameters. Such a comprehensive approach may help to better quantify the flash flood parameters of ungauged streams and thus contribute to more reliable management of small streams in future extreme climate events.

Study Site
The hydrogeomorphic impacts of floods were studied on a small tributary of the Bělá River in the foreland of the Eastern Sudetes (NE part of Czechia, Central Europe; 50°17'N, 17°17'E; Figure 1). The area comprises a Proterozoic and Palaeozoic basement composed of orthogneiss, which is overlaid by glacifluvial deposits of the Saalian and Elsterian glaciations with different grain-size distributions of till sediments. These deposits are the result of deglaciation phases in proglacial areas on the ice margins and were recently covered by the Holocene sandy loam colluvium [40,41]. The study site is characterised by a temperate climate, with a mean annual precipitation between 850 and 900 mm. Most precipitation falls during the spring and summer months [42], with an occasional occurrence of extreme daily rainfall (more than 50 mm per day). Documentary evidence provides information about several flood events (e.g., 1903, 1921, 1971, 1977, 2007, 2009, and 2014), relating to both advective rains and short-term intense rainfall events that are currently responsible for intense gully incision and damage to infrastructure within the studied region [38,43]. We focused on an approximate 700-m-long reach of an unnamed first-order stream draining a catchment area of 1.8 km 2 with a mean elevation of 417.5 m a.s.l. Fields and dwellings dominate the upper and middle catchment area, while the lower part (i.e., the study reach) is predominantly covered by a mixed forest (Figure 1c). The mean stream gradient of the study reach is 0.05 m/m, with a maximum up to 0.17 m/m. The channel cuts into the Pleistocene and Holocene deposits, with typical alternations of stepped-bed morphology and bedrock outcrops (Figure 2a) accompanied by poorly We focused on an approximate 700-m-long reach of an unnamed first-order stream draining a catchment area of 1.8 km 2 with a mean elevation of 417.5 m a.s.l. Fields and dwellings dominate the upper and middle catchment area, while the lower part (i.e., the study reach) is predominantly covered by a mixed forest (Figure 1c). The mean stream gradient of the study reach is 0.05 m/m, with a maximum up to 0.17 m/m. The channel cuts into the Pleistocene and Holocene deposits, with typical alternations of stepped-bed morphology and bedrock outcrops (Figure 2a) accompanied by poorly sorted bed sediments (Figure 2b,c). Fresh slope/bank failures (up to 50 m long and 10 m high) and generally unstable banks have caused channel widening during recent floods, resulting in the frequent occurrence of exposed tree root systems and flood scars on riparian vegetation along the entire reach (Figure 2b-d). In contrast, there is no evidence of deep channel incision due to the resistant gneiss bedrock and large interlocked boulders at several reaches within the stream. Slopes and floodplains surrounding the channel are overgrown by a mixed forest composed predominantly of Pinus sylvestris L., Picea abies (L.) Karst., Alnus glutinosa (L.) Gaertn, and Tilia cordata Mill. At 0.5 river km (r. km), the stream drains to a 50-m-long concrete trough within a viaduct under a railway. Moreover, culverts are presented under a cycle path in the lower part and under a road through a village in the middle part of the catchment.

Dendrogeomorphic Fieldwork and Analyses
First, the terrain fieldwork focused on dendrogeomorphic sampling to create the chronology of past flood events and to select the flood scars caused by the 2014 flood event. The tree sampling followed the standard dendrogeomorphic procedure for flood reconstruction [13,44]. The sampling was focused on flood scars occurring either on tree stems or exposed tree roots. In the case of tree stems, increment cores were extracted from the edge of the scar and from the undisturbed part of the tree stem using the Pressler increment borer (40 × 0.5 mm). Tree stems of a small diameter (up to 10 cm) were cut by handsaw in the position of the flood scar to gain the stem disc or wedge while approximately 2-cm-wide cross-sections were sampled from exposed and scarred living roots. Only those scars oriented against the supposed direction of the flow path were considered for sampling.

Dendrogeomorphic Fieldwork and Analyses
First, the terrain fieldwork focused on dendrogeomorphic sampling to create the chronology of past flood events and to select the flood scars caused by the 2014 flood event. The tree sampling followed the standard dendrogeomorphic procedure for flood reconstruction [13,44]. The sampling was focused on flood scars occurring either on tree stems or exposed tree roots. In the case of tree stems, increment cores were extracted from the edge of the scar and from the undisturbed part of the tree stem using the Pressler increment borer (40 × 0.5 mm). Tree stems of a small diameter (up to 10 cm) were cut by handsaw in the position of the flood scar to gain the stem disc or wedge while approximately 2-cm-wide cross-sections were sampled from exposed and scarred living roots. Only those scars oriented against the supposed direction of the flow path were considered for sampling. In addition, only scars that occurred on exposed but stabilized roots with the lowest tendency to flex during higher discharges were used to avoid the underestimation of modelled peak flow discharge [45]. All sampled trees were carefully described, and the height of the sampled scar was carefully noted, labelled with a visible object, photographed, and targeted using GPS. Reference trees growing near the study reach without any geomorphic influence on tree growth were sampled to cross-date with the disturbed samples and to eliminate false and missing rings.
Laboratory processing followed the standard dendrogeomorphic procedure [13,46]. All increment cores were glued into woody supports and-together with cross sections-dried and polished to be ready for dendrogeomorphic analysis. The tree rings of the increment cores and stem discs were counted and measured using TimeTable and PAST4 software [47], and their growth patterns were compared with the appropriate reference chronology (compiled in Arstan software [48] using a double detrending procedure) to ensure the reliability of dating. As root segments are more problematic regarding the occurrence of missing and wedging rings, the zig-zag segment tracing method [49] was applied to carefully count the years of each ring. In the case of problematic roots (i.e., dense growth increment and small roots), we used microslides cut by GLS-1 microtome and processed according to standard chemical procedures to precisely define the position of the flood scar within the root section [50,51].
In the next step, we identified the years with the occurrence of scars and onset of callus tissues ( Figure 3) and compiled the chronology of past flood events. Scars represent an unequivocal signal of flood events and are considered the most reliable growth disturbance in dendrogeomorphic flood reconstructions [44]. The event identification was based on the event-response index (I t index [52]), calculated as: where R is the number of scars in a year t and A is the total number of sampled trees living in a year t. Then, a certain event was considered when I t ≥ 10% and the number of scars ≥ 3, while a probable event was determined when 10% > I t ≥ 5% and the number of scars ≥ 2. From the whole dataset of scars, we eventually selected the scars dated to 2014 as a PSI of the May 2014 flash flood event, whose parameters were modelled.
Water 2020, 12, x FOR PEER REVIEW 5 of 19 In addition, only scars that occurred on exposed but stabilized roots with the lowest tendency to flex during higher discharges were used to avoid the underestimation of modelled peak flow discharge [45]. All sampled trees were carefully described, and the height of the sampled scar was carefully noted, labelled with a visible object, photographed, and targeted using GPS. Reference trees growing near the study reach without any geomorphic influence on tree growth were sampled to cross-date with the disturbed samples and to eliminate false and missing rings. Laboratory processing followed the standard dendrogeomorphic procedure [13,46]. All increment cores were glued into woody supports and-together with cross sections-dried and polished to be ready for dendrogeomorphic analysis. The tree rings of the increment cores and stem discs were counted and measured using TimeTable and PAST4 software [47], and their growth patterns were compared with the appropriate reference chronology (compiled in Arstan software [48] using a double detrending procedure) to ensure the reliability of dating. As root segments are more problematic regarding the occurrence of missing and wedging rings, the zig-zag segment tracing method [49] was applied to carefully count the years of each ring. In the case of problematic roots (i.e., dense growth increment and small roots), we used microslides cut by GLS-1 microtome and processed according to standard chemical procedures to precisely define the position of the flood scar within the root section [50,51].
In the next step, we identified the years with the occurrence of scars and onset of callus tissues ( Figure 3) and compiled the chronology of past flood events. Scars represent an unequivocal signal of flood events and are considered the most reliable growth disturbance in dendrogeomorphic flood reconstructions [44]. The event identification was based on the event-response index (It index [52]), calculated as: where R is the number of scars in a year t and A is the total number of sampled trees living in a year t. Then, a certain event was considered when It ≥ 10% and the number of scars ≥ 3, while a probable event was determined when 10% > It ≥ 5% and the number of scars ≥ 2. From the whole dataset of scars, we eventually selected the scars dated to 2014 as a PSI of the May 2014 flash flood event, whose parameters were modelled.

Channel Parameters and Channel Geometry
This part consisted of (i) the measurement of the longitudinal profiles (channel/floodplain geometry) to model the 2014 flood behaviour and (ii) the classification of erosive/depositional segments and the measurement of the largest clasts for the determination of channel stability during the 2014 floods: (i) The channel/floodplain geometry and the height of the 2014 scars (PSI) were measured in cross sections using a total station (GTS-212) with an accuracy of ±2.5 cm and a GNSS geodetic receiver (Trimble R2) with an accuracy of ±4 cm. All measurements were taken in the spring of 2019. To register the complexity of the channel geometry, the density of cross sections was as follows: the position of the first cross section was located directly in the position of the scars. Afterwards, five cross sections were measured, with the distance of one meter among each other, followed by two cross sections at a distance of 2.5 m and one cross section at a distance of 5 m. Next, cross sections were measured at a distance of 10 m. In the case of close distances between the scars, no overall process was applied. In total, 341 cross sections were measured, resulting in 4300 surveyed points. The GNSS geodetic receiver was only applied to place the relative position of points to the absolute location based on the local geographical reference system. A digital terrain model (DTM) in raster format was created from these cross sections with a grid size equal to 0.1 m and was used as a source for the hydraulic modelling.
(ii) To reveal the channel stability during the last flood event, the middle axes of the five largest bed particles were measured by a tape (with ±0.01 m accuracy; with a less-precise ±0.05 m accuracy only in a few cases of partially buried boulders) in 10 ± 1 m intervals along the stream's longitudinal profile, except for the reach located in the culvert (0.51-0.54 r. km). Only the particles within the bankfull channel were accounted for. Consequently, the representative mean boulder diameter MBD (mm) for each of the channel cross sections was calculated as the arithmetical mean of these five measurements [23]. We classified each cross section as erosional, stable, or depositional by the observed signs of the present stability of the adjacent channel reach. The erosional reaches indicated trends of incision together with the frequent presence of exposed roots or bedrock outcrops in the channel banks, whereas the depositional reaches were typified by the occurrence of locally widened channels with developed unvegetated bars. The stable reaches represented the transport-balanced segments without evident signs of recent incision or bed aggradation. We observed no wood obstructions in the channel that would influence the channel bed stability and sediment transport processes.

Hydraulic Modelling
This phase comprised three steps: (i) hydraulic model creation, setup and calibration; (ii) estimation of scar peak discharges (SPD) and reach peak discharge; and (iii) scenario modelling of the 2014 flood event. The two-dimensional IBER model (version 2.5.1) was applied to the hydraulic modelling of the selected reach [53]. This is an established software that was applied to the estimation of palaeoflood discharges of small [18] to large rivers [35,45]. An unstructured mesh, which comprised almost 200,000 elements with an average size of 0.3 m, was developed over the reach. As an initial condition at t = 0, the river was set dry. The flow was subcritical throughout the whole domain for all used discharges. We imposed inlet boundary conditions based on the uniform discharge and the critical depth at the outlet of the studied reach. The first PSI was located 20 m from the outlet of the reach, which allowed the model to overcome inaccuracy in the selected boundary conditions. A wet-dry threshold of 0.01 m and 2nd order roe scheme was chosen. The Courant-Friedrichs-Lewy was set to 0.9 and the mixed length turbulence model was selected. Although the erosion and accumulation of alluvium could occur during the 2014 flood, we considered the stable riverbed during the modelling [54]. The model was calibrated to a single value of water stage in the cross section where we measured the discharge (equal to 25.23 cubic litres per second), using the velocity meter. The cross section was located in the downstream part of the selected reach and the measurement took place at the end of May 2019. Based on the calibration results, a uniform value for Manning's roughness coefficient (n) equal to 0.08 was applied to the overall reach. The selection of the roughness value was based on the studies dealing with high-gradient streams [55,56]. In the second step, an iterative process was run to find the SPD [57] that produced the best fit between the scars and modelled water depths [58] and to select the reach peak discharge that produced the best RMSE value [32]. Finally, the three scenarios were established to model the characteristics (unit stream power and bed shear stress) of the 2014 flood event: Q min (discharge equal to the 1st quartile of the SPD), Q optimum (reach peak discharge which produced the optimum/lowest value of RMSE), and Q max (the 3rd quartile of the SPD). Although the main results were created for the Qoptimum, the remaining scenarios allowed the evaluation of the uncertainty inherent in palaeoflood discharge estimation [36].

Relations between the Hydraulic and Sedimentologic Parameters and Calculation of Channel Stability
To evaluate the flow competence of the 2014 flood (i.e., the ability of this flood to transport coarse bed material and destabilize the channel bed), we used the critical unit stream power ω ci -transported particle diameter D i (mm) relationship. To approximate the local conditions of the frequent occurrence of a stepped-bed morphology, relatively high channel gradients and a small catchment area, we applied the relationship developed by field observations in similar steep channels (0.06 ≤ S ≤ 0.14 m/m) draining small catchments (0.2 ≤ A ≤ 2.2 km 2 ) in Central European medium-high mountain settings [23,59]: This relationship (2) was derived from direct observations of the largest boulders that were transported by a high-magnitude flood event and from the displacements of marked particles during lower discharges (covering grain-sizes 20-400 mm) in stepped-bed streams with poorly sorted bed sediments. As D i , we substituted MBD and calculated the potential critical unit stream power ω cMBD , which will lead to the incipient motion of MBD in a given cross section. The resulting value was compared with the unit stream power ω simulated for the 2014 event, and the excess critical unit stream power ω E was expressed in the following form: This implies that ω E ≥ 1 indicates a potentially unstable bed structure consisting of coarse grains (i.e., rapids or individual step units) during the 2014 event. We simulated ω E for all three scenarios (Q min , Q optimum and Q max ).
In addition, we employed the peak bed shear stress of the 2014 event (τ b ) to find a possible relationship between the parameters of τ b , ω and MBD and to assess potential differences between the groups of cross-sections by their present stability (erosional, vertically stable and depositional) during the Q optimum scenario. Due to the non-normality of the data, we used Spearman's correlation r sp to test the potential relationships between τ b , ω and MBD. A non-parametric Kruskal-Wallis test was used to compare τ b , ω and MBD between the reaches with depositional, transport-balanced and erosional tendencies when the Z-value test with Bonferonni corrections was used to distinguish significantly different groups. We used a significance level of 0.05 for all the tested data.

Chronology of Past Flood Events and Botanical Evidence of the May 2014 Flash Flood Event
In total, we successfully sampled and dated 75 scarred individuals (20 tree stems and 55 scarred roots) with a predominance of samples from P. sylvestris (37.3%), P. abies (18.7%), and A. glutinosa (17.3%), whereas five samples (root cross sections) had to be excluded from the chronology due to uncertainties during the dating procedure (Table 1). Eventually, it was possible to date 115 scars (1.5 scars per tree). Overall, we determined 13 flood events (seven certain and six probable events) during the period of 1955 to 2018 (limited by the minimum number of 10 sampled trees; Figure 4).  Focusing on the botanical evidence of the 2014 flood event, we identified 26 scars throughout the whole study reach (6 scars on tree stems and 20 scars on exposed roots; Table 1). The mean height of the scars above the thalweg was 96.6 ± 34.6 cm. The minimum height (41 cm) was observed at 0.33 r. km at a straight channel reach, while the maximum height (156 cm) was recorded at 0.07 r. km at the failure of a concave bank. The distribution of the 2014 scars ( Figure 5) was zero at the reach between 0.35 and 0.5 r. km, where the total number of scarred trees was generally lower, and the scars were generally of older dates. In contrast, the highest abundance of the 2014 flood scars was recorded between 0.2 and 0.35 r. km (11 scars overall).   Focusing on the botanical evidence of the 2014 flood event, we identified 26 scars throughout the whole study reach (6 scars on tree stems and 20 scars on exposed roots; Table 1). The mean height of the scars above the thalweg was 96.6 ± 34.6 cm. The minimum height (41 cm) was observed at 0.33 r. km at a straight channel reach, while the maximum height (156 cm) was recorded at 0.07 r. km at the failure of a concave bank. The distribution of the 2014 scars ( Figure 5) was zero at the reach between 0.35 and 0.5 r. km, where the total number of scarred trees was generally lower, and the scars were generally of older dates. In contrast, the highest abundance of the 2014 flood scars was recorded between 0.2 and 0.35 r. km (11 scars overall).

Results of Hydraulic Modelling
The flow conditions during the simulations were subcritical in the major part of the reach (91%-94% of the total area) for all simulated flow scenarios ( Table 2). Other flow characteristics (depth, velocity, and Froude number) resulting from the raster of hydraulic modelling are described in Table  2. The characteristics of the RMSE curve of Qoptimum during the 2014 flood and the SPD characteristics are shown in Figure 6. The value of the reach peak discharge (Qoptimum) of the 2014 flood estimated by the RMSE (0.32 m) was equal to 4.5 m 3 ·s −1 . The values of the SPD showed high variability (coefficient of variability = 0.79) when the minimal SPD was 1.5 m 3 ·s −1 and the maximal SPD was equal to 16.5 m 3 ·s −1 . The values of the Qmin and Qmax scenarios were calculated from all estimated SPD and were equal to 2.63 and 9.38 m 3 ·s −1 , respectively. The deviation between the measured PSI elevations and the simulated water level for Qoptimum ranged from −0.60 to 0.63 m with a median equal to −0.05 m. Similarly, the Qmin deviations ranged from −0.30 to 0.75 m, and the median was equal to 0.08 m. For Qmax, we registered deviations from −1.07 to 0.42 m, and the median was equal to −0.31 m. The highest magnitudes of velocity (up to 11.97 m·s −1 ; Figure 7) and bed shear stress (up to 1787.10 N.m −2 ) were registered at the positions of the highest changes in riverbed topography.
We further investigated the effect of the hydraulic conditions in the positions of the PSIs using non-parametric Spearman's correlation. However, we did not observe any significant relationship between the PSI deviations and the velocity magnitudes for all three scenarios: Qoptimum (rsp = 0.19, p = 0.36), Qmin (rsp = 0.08; p = 0.71), and Qmax (rsp = 0.004; p = 0.98; Figure 8).

Results of Hydraulic Modelling
The flow conditions during the simulations were subcritical in the major part of the reach (91-94% of the total area) for all simulated flow scenarios ( Table 2). Other flow characteristics (depth, velocity, and Froude number) resulting from the raster of hydraulic modelling are described in Table 2. The characteristics of the RMSE curve of Q optimum during the 2014 flood and the SPD characteristics are shown in Figure 6. The value of the reach peak discharge (Q optimum ) of the 2014 flood estimated by the RMSE (0.32 m) was equal to 4.5 m 3 ·s −1 . The values of the SPD showed high variability (coefficient of variability = 0.79) when the minimal SPD was 1.5 m 3 ·s −1 and the maximal SPD was equal to 16.5 m 3 ·s −1 . The values of the Q min and Q max scenarios were calculated from all estimated SPD and were equal to 2.63 and 9.38 m 3 ·s −1 , respectively. The deviation between the measured PSI elevations and the simulated water level for Q optimum ranged from −0.60 to 0.63 m with a median equal to −0.05 m. Similarly, the Q min deviations ranged from −0.30 to 0.75 m, and the median was equal to 0.08 m. For Q max , we registered deviations from −1.07 to 0.42 m, and the median was equal to −0.31 m. The highest magnitudes of velocity (up to 11.97 m·s −1 ; Figure 7) and bed shear stress (up to 1787.10 N.m −2 ) were registered at the positions of the highest changes in riverbed topography.

Channel Stability and Relations between the Hydraulic and Sedimentologic Parameters
A positive significant relationship exists between the bed shear stress and the unit stream power calculated for all cross sections (n = 62) and the Qoptimum scenario during the 2014 event (rsp = 0.65, p < 0.0001) (Figure 9a). In contrast, there were no significant correlations between the parameters of τb and MBD (rsp = 0.05, p = 0.68) or ω and MBD (rsp = −0.04, p = 0.74) (Figure 9b,c). After the removal of

Channel Stability and Relations between the Hydraulic and Sedimentologic Parameters
A positive significant relationship exists between the bed shear stress and the unit stream power calculated for all cross sections (n = 62) and the Qoptimum scenario during the 2014 event (rsp = 0.65, p < 0.0001) (Figure 9a). In contrast, there were no significant correlations between the parameters of τb and MBD (rsp = 0.05, p = 0.68) or ω and MBD (rsp = −0.04, p = 0.74) (Figure 9b,c). After the removal of We further investigated the effect of the hydraulic conditions in the positions of the PSIs using non-parametric Spearman's correlation. However, we did not observe any significant relationship between the PSI deviations and the velocity magnitudes for all three scenarios: Q optimum (r sp = 0.19, p = 0.36), Q min (r sp = 0.08; p = 0.71), and Q max (r sp = 0.004; p = 0.98; Figure 8).

Channel Stability and Relations between the Hydraulic and Sedimentologic Parameters
A positive significant relationship exists between the bed shear stress and the unit stream power calculated for all cross sections (n = 62) and the Qoptimum scenario during the 2014 event (rsp = 0.65, p < 0.0001) (Figure 9a). In contrast, there were no significant correlations between the parameters of τb and MBD (rsp = 0.05, p = 0.68) or ω and MBD (rsp = −0.04, p = 0.74) (Figure 9b,c). After the removal of

Channel Stability and Relations between the Hydraulic and Sedimentologic Parameters
A positive significant relationship exists between the bed shear stress and the unit stream power calculated for all cross sections (n = 62) and the Q optimum scenario during the 2014 event (r sp = 0.65, p < 0.0001) (Figure 9a). In contrast, there were no significant correlations between the parameters of τ b and MBD (r sp = 0.05, p = 0.68) or ω and MBD (r sp = −0.04, p = 0.74) (Figure 9b,c). After the removal of 15 cross-sections assigned to deep pools or bedrock sections containing a limited number of coarse grains, the final relationship was significant for τ b and MBD (r sp = 0.33, p = 0.024) and remained insignificant for ω and MBD (r sp = 0.13, p = 0.38).  The mutual comparison of the MBD of cross sections located in depositional (n = 8), stable (n = 30) and erosional reaches (n = 24) showed significant differences between the groups (p = 0.002), when the cross sections of stable reaches indicated significantly higher values of MBD than those of locations with contemporary signs of prevailing depositional or erosional processes (Figure 10a). The analysis of the bed shear stress revealed significant differences between these groups (p = 0.012) when the cross sections of erosional reaches indicated significantly higher values of τb than those of the other groups (Figure 10b). On the other hand, we observed no significant differences among these groups in terms of ω (p = 0.095), although the erosional reaches were again characterised by somewhat higher values of ω than the other reaches (Figure 10c).
The evaluation of flow competence and the potential stability of coarse bed material in individual cross sections during the 2014 event showed that 15 of 62 cross sections (24.2%) indicated ωE ≥1 for the Qoptimum scenario ( Figure 11). These cross sections were frequently representative of bedrock reaches with erosional tendencies with relatively fine grain-size character of the limited alluvial cover. An extraordinarily high excess was observed immediately downstream of the culvert (0.56-0.57 r. km). For the Qmax scenario, half of the cross sections (31 of 62; 50.0%) were perceived as those with unstable alluvial cover, still leaving the part upstream of the culvert (0.38-0.48 r. km) as relatively stable. On the other hand, the Qmin scenario predicted only seven unstable cross sections (11.3%) with a notable excess of critical unit stream power immediately downstream of the culvert (0.56-0.57 r. km). The mutual comparison of the MBD of cross sections located in depositional (n = 8), stable (n = 30) and erosional reaches (n = 24) showed significant differences between the groups (p = 0.002), when the cross sections of stable reaches indicated significantly higher values of MBD than those of locations with contemporary signs of prevailing depositional or erosional processes (Figure 10a). The analysis of the bed shear stress revealed significant differences between these groups (p = 0.012) when the cross sections of erosional reaches indicated significantly higher values of τ b than those of the other groups (Figure 10b). On the other hand, we observed no significant differences among these groups in terms of ω (p = 0.095), although the erosional reaches were again characterised by somewhat higher values of ω than the other reaches (Figure 10c). 15 cross-sections assigned to deep pools or bedrock sections containing a limited number of coarse grains, the final relationship was significant for τb and MBD (rsp = 0.33, p = 0.024) and remained insignificant for ω and MBD (rsp = 0.13, p = 0.38). The mutual comparison of the MBD of cross sections located in depositional (n = 8), stable (n = 30) and erosional reaches (n = 24) showed significant differences between the groups (p = 0.002), when the cross sections of stable reaches indicated significantly higher values of MBD than those of locations with contemporary signs of prevailing depositional or erosional processes (Figure 10a). The analysis of the bed shear stress revealed significant differences between these groups (p = 0.012) when the cross sections of erosional reaches indicated significantly higher values of τb than those of the other groups (Figure 10b). On the other hand, we observed no significant differences among these groups in terms of ω (p = 0.095), although the erosional reaches were again characterised by somewhat higher values of ω than the other reaches (Figure 10c).
The evaluation of flow competence and the potential stability of coarse bed material in individual cross sections during the 2014 event showed that 15 of 62 cross sections (24.2%) indicated ωE ≥1 for the Qoptimum scenario ( Figure 11). These cross sections were frequently representative of bedrock reaches with erosional tendencies with relatively fine grain-size character of the limited alluvial cover. An extraordinarily high excess was observed immediately downstream of the culvert (0.56-0.57 r. km). For the Qmax scenario, half of the cross sections (31 of 62; 50.0%) were perceived as those with unstable alluvial cover, still leaving the part upstream of the culvert (0.38-0.48 r. km) as relatively stable. On the other hand, the Qmin scenario predicted only seven unstable cross sections (11.3%) with a notable excess of critical unit stream power immediately downstream of the culvert (0.56-0.57 r. km).  for the Q optimum scenario (Figure 11). These cross sections were frequently representative of bedrock reaches with erosional tendencies with relatively fine grain-size character of the limited alluvial cover. An extraordinarily high excess was observed immediately downstream of the culvert (0.56-0.57 r. km). For the Q max scenario, half of the cross sections (31 of 62; 50.0%) were perceived as those with unstable alluvial cover, still leaving the part upstream of the culvert (0.38-0.48 r. km) as relatively stable. On the other hand, the Q min scenario predicted only seven unstable cross sections (11.3%) with a notable excess of critical unit stream power immediately downstream of the culvert (0.56-0.57 r. km).

Discussion
The dating of flood events, together with the determination of the 2014 flash flood parameters in the moderate relief of Central Europe, were introduced using the combination of dendrogeomorphic methods, 2D hydraulic modelling, and sedimentological parameters. Based on these approaches, we provide a quantitative estimation of the last flash flood event and a determination of channel reaches that tended to be (un)stable. Unlike the palaeoflood reconstructions from larger catchments, we included scarred roots as a possible PSI in the first-order catchments due to the generally lower peak flow discharges and lower volume of transported sediments. If the scarred root is a part of a stable root system and does not evidence a high rate of flexibility, it can be used not only as a helpful tool to date the time of (flood) erosion [21,60,61], but also as a PSI of recent floods.

Hydrogeomorphic Response of Flash Floods in the First-Order Catchment
The dendrogeomorphic results confirm that, despite the strong geomorphic impact of the last 2014 flash flood, there is evidence of former flood events within this catchment. We recorded the hydrogeomorphic impacts of floods in years with the occurrence of debris flows and rockfalls in the surrounding mountains (e.g., in 1991, 2006 and 2010 [62,63]). In addition, several identified years coincided with documentary data about local and/or regional flooding there (e.g., in 1971, 1977, 1997 and 2009). Polách and Gába [43] described a spatially limited downpour in July 1971 after strong antecedent precipitation, resulting in local damage to small streams within this region. This situation was likely similar to that of the last 2014 flash flood, which was considered unprecedented regarding the rainfall intensity, but the hydrogeomorphic response seems to be comparable to that of the 1971 event. Moreover, the transported parts of damaged culverts within the channel and several older scars on tree stems (Figure 2d) suggest an even higher hydrogeomorphic impact in the past.

Discussion
The dating of flood events, together with the determination of the 2014 flash flood parameters in the moderate relief of Central Europe, were introduced using the combination of dendrogeomorphic methods, 2D hydraulic modelling, and sedimentological parameters. Based on these approaches, we provide a quantitative estimation of the last flash flood event and a determination of channel reaches that tended to be (un)stable. Unlike the palaeoflood reconstructions from larger catchments, we included scarred roots as a possible PSI in the first-order catchments due to the generally lower peak flow discharges and lower volume of transported sediments. If the scarred root is a part of a stable root system and does not evidence a high rate of flexibility, it can be used not only as a helpful tool to date the time of (flood) erosion [21,60,61], but also as a PSI of recent floods.

Hydrogeomorphic Response of Flash Floods in the First-Order Catchment
The dendrogeomorphic results confirm that, despite the strong geomorphic impact of the last 2014 flash flood, there is evidence of former flood events within this catchment. We recorded the hydrogeomorphic impacts of floods in years with the occurrence of debris flows and rockfalls in the surrounding mountains (e.g., in 1991, 2006 and 2010 [62,63]). In addition, several identified years coincided with documentary data about local and/or regional flooding there (e.g., in 1971, 1977, 1997 and 2009). Polách and Gába [43] described a spatially limited downpour in July 1971 after strong antecedent precipitation, resulting in local damage to small streams within this region. This situation was likely similar to that of the last 2014 flash flood, which was considered unprecedented regarding the rainfall intensity, but the hydrogeomorphic response seems to be comparable to that of the 1971 event. Moreover, the transported parts of damaged culverts within the channel and several older scars on tree stems (Figure 2d) suggest an even higher hydrogeomorphic impact in the past. Furthermore, we identified that more than 40% of scars can be dated to the period 1951-2000, pointing to slow and progressive channel bank erosion during several flood events rather than to abrupt changes in channel geometry and the fast decay of exposed roots, as is typical for mountain regions [64].
We found no significant correlation by plotting the MBD and τ b for all 62 cross sections, but a significant positive relationship existed after removal of fifteen cross sections assigned to deep pools and bedrock segments. This increase in statistical significance reflected the presence of the fine-grained character of the alluvial cover in these cross sections, the calibre of which did not necessarily correspond to the stream transport capacity (expressed by τ b ) during the examined flood event. Thus, we assumed more frequent transport of these fine particles during lower-magnitude events (e.g., bankfull or even lower flows). On the other hand, the set of remaining 47 cross sections had adjusted the calibre of the MBD to the calculated bed shear stresses. This positive relationship clearly pointed to the adjustment of the stepped-bed architecture during high-magnitude events towards a condition that provides the maximum possible bed stability, as was documented by previous flume experiments [65].
In general, the highest values of τ b and ω were calculated for cross sections with erosional tendencies when compared to those of stable or depositional cross sections, although only the parameter of τ b indicated statistically significant differences. These high values of τ b and ω typically reflected a low width-depth ratio in erosional reaches with limited inundation ability together with the concentrated energy of flood flows. The presence of alluvial pockets of relatively fine sediments in the erosional reaches with exposed bedrock (resulting in very low values of MBD) explained our observations of the coarser bed particles in the stable reaches when one may expect local coarsening in erosional reaches with the highest calculated τ b and, contrarily, fining in depositional reaches [23,66,67]. In addition, no differences in the MBD were observed between the erosional and stable reaches in the case of stepped-bed headwater streams based in flysch rocks, despite higher unit stream power being calculated for the erosional reaches [23]. These observations from small streams contradict the relation between the sediment coarsening and channel incision that was reported for gravel-bed rivers [68] by taking into consideration the specifics of steep headwater streams (e.g., mixture of alluvial and semi-alluvial reaches or strong effect of local lithology predicting patterns of sediment supply).
The parameter ω is frequently used to determine the capacity of a stream to mobilize specific bed grain sizes [37,69,70]. In our case, the excess critical unit stream power ω E identified the segments with unstable alluvial cover during the examined flood event. Only one-quarter of the cross sections indicated transport of the MBD under the Q optimum scenario, when these cross sections were often representative of bedrock reaches with limited alluvial cover consisting of relatively fine grains. Even by applying the Q max scenario, approximately 50% of the cross sections remained stable. Despite the large number of scars possibly related to intensive sediment transport and bank erosion during the 2014 flood, this event likely did not reach the critical threshold discharge to completely rework the stepped-bed character of the studied stream. Previous field measurements in steep mountain streams have perceived high-magnitude floods of up to a 50-year recurrence interval as those that mobilize most step-forming particles [71,72]. This suggests either a lower-than-expected magnitude of the 2014 flood described by local authorities or high channel bed stability of this first-order stream owing to the presence of large interlocked boulders.

Benefits and Limits of the Approach Used in the First-Order Catchment
Introducing the scarred roots as a PSI in a first-order catchment entails uncertainties similar to those in the case of scars on tree stems in larger rivers. As noted by several authors [18,36], PSIs located in straight channel reaches or on the inner side of channel bends are more suitable for peak discharge reconstructions than those of trees located on the outer side of channel bends or growing in overbank sections with dense vegetation cover. In our study, the majority of the PSIs (20) were located in the straight channel and/or inner side, reducing the overestimation of peak flow discharge. Moreover, the avoidance of root sampling in concave banks and slope failures for the PSI estimation reduces the probability of dating scars caused by bank/slope erosion. This situation may frequently occur, even during floods [73], but cannot be used as information about peak flow discharge.
A tricky situation may occur when several consecutive intense rainfalls affect the catchment. Therefore, more than one peak flow discharge may occur and thus enter as the uncertainty in hydraulic modelling. In our case, several storms generated during May 2014 could result in higher-than-average discharges [38], so we cannot exclude the possibility that several of the dated scars with lower PSI belong to some lower spring peak flows rather than to the surveyed flood of 27 May 2014. In contrast, anatomical analysis of some scarred roots revealed the position of the 2014 scar within the earlywood cells, pointing to the spring floods ( Figure 3d) and thus eliminating the possible inclusion of scars during intense summer or autumn rainfall in 2014.
The results of the peak flow estimations show a high variation in SPD and high deviation in differences between the predicted water level and the observed PSI height. We hypothesize that this could be related to the small size of the catchment, where the role of the input data (PSI, DTM and roughness coefficient) is more crucial than in a large catchment (which produces large peak discharges). Similarly, Ballesteros-Cánovas et al. [18] reported the highest deviation of the estimated peak discharge in a catchment of the lowest size, while the opposite conclusions are described by Ballesteros-Cánovas et al. [57], where the authors reported much larger uncertainties for larger catchments than for smaller catchments. Following the approach of Bohorquez et al. [32], who defined the critical value of deviation between the simulated water level and the observed height of PSI as 0.2 m, for the Q min , fifteen scars were within this range, followed by Q optimum (12 PSI) and Q max (10 PSI).
Further uncertainties associated with peak discharge reconstruction using the hydraulic model are related to the DTM accuracy [58] and model calibration. Although we tried to create an accurate DTM, especially at the position and close surroundings of the PSIs, our DTM is not error free, particularly between the cross sections where the interpolation took place. Victoriano et al. [34] found a 7% difference in the estimated peak discharge of 316 m 3 ·s −1 (catchment area: 5.72 km 2 ) when two different DTMs were compared. We call for further research to assess how the level of DTM accuracy influences the results of peak discharge reconstructions, even in small catchments. The uniform value of Manning's roughness coefficient was applied over the overall reach, and the final value was a product of the calibration exercise to single "low-flow" discharge. Indeed, differences in the calibration for the low-and high-flow discharges are well known [56] and the influence of various values of Manning's roughness coefficient on the results of the hydraulic model are described by Ballesteros-Cánovas et al. [57]. In fact, the hydraulic modelling of flood waves in large, low gradient rivers could cause serious problems with the approach used in this study (i.e., calibration to the discharge several orders of magnitude lower than simulated flood). However, in the case of small, high-gradient streams, this approach is better than the selection of the roughness value from the literature. Overall, we believe our approach reduced the uncertainty in the roughness coefficient to an acceptable level.
Although we considered the channel bed as a stable component with limited erosion and/or deposition during the 2014 flood, this does not reflect the natural condition in the modelled reach, because the majority of the riverbed is formed by poorly sorted alluvium with stripes of bedrock in a few places. This simplification could further influence the height of the PSI and the estimated discharge [14]. In addition, we did not observe any significant relationship between the velocity of the flood wave and the absolute deviation of the PSI from the modelled discharges (Q min , Q optimum , and Q max ), which is in agreement with the study of Ballesteros-Cánovas et al. [14]. We assume that the position of the PSI regarding the channel geometry (convex/concave/straight reach) is more important in small streams where the velocity is primarily changing with the stream gradient and channel width. Unfortunately, we were not able to test the influence of PSI deviations among convex, concave, and straight reaches due to the limited number of PSIs in each group.
We simulated subcritical flow over the overall reach, although supercritical flow occurred during the modelling in some parts. Flow conditions in steep channels can be transitional, with supercritical flow over the steps and subcritical flow in the pools. However, flow will never be supercritical throughout the simulated domain, because such flow would destabilize most channels [74]. This fact is supported by other studies [75,76], which described that, despite high velocity and extreme turbulence, flow in mountain streams is critical or subcritical, leading to the assumption that supercritical flow in natural streams does not exist for any extended length.

Conclusions
The estimations of peak flood discharges and the quantifications of related morphological changes and bedload transport in small streams remain a challenge, especially in the case of poorly gauged first-order catchments, despite their ability to transfer high amounts of water and sediment into low-order catchments. In this study, we presented an integrated approach based on palaeoflood reconstruction, which was performed in a first-order catchment in the moderate relief of Central Europe, to provide quantitative data about the last 2014 flash flood event and to create the chronology of former floods. The dendrogeomorphic results show the regular occurrence of flood events during at least the last 60 years, although the intensity of the last 2014 flash flood was considered to be unprecedented according to local authorities. The results of 2D hydraulic modelling are in favour of the optimum scenario of the 2014 peak flow discharge of 4.5 m 3 ·s −1 , which resulted in the formation of root and stem scars. Nevertheless, a sedimentological survey confirms the rather progressive development of the studied channel reach with limited ability to transport step-forming boulders during peak flow discharges similar to that in 2014 and thus relative bed stability during events of such magnitude. Despite notable geomorphic imprints and ongoing lateral erosion caused by the most recent flash flood events, with damage to infrastructure (especially in 2009 and 2014), we conclude that substantial geomorphic transition is practically excluded during the similar rainfall episodes that occurred within the last 60 years. The limited precision of the results of hydraulic modelling lead to the relatively high variability in possible peak flow in 2014 (from modelled Q min to Q max ). This level of accuracy is still comparable to larger rivers, but may slightly change the interpretation of suggested event magnitudes in the case of smaller catchments. Therefore, further research dealing with the quality of input data into palaeoflood reconstruction (e.g., DTM accuracy and variable channel roughness) is crucial to investigation into first-order catchments.
Funding: This research received no external funding.