There are no hydrometric data for the Barcarena watershed. For this reason, peak discharges for the return periods (R_T) of 2, 5, 10, 20, 50 and 100 years were obtained using designed hyetographs and a rainfall-runoff model performed in HEC-HMS software (Figure 2). The designed hyetographs for the considered return periods were developed from the IDF (intensity-duration-frequency) curves of the São Julião do Tojal (SJT in Figure 1), a rain gauge located about 15 km NE to the study area. The IDF curves were built from the parameters available in Brandão et al. (2001), which was established from a period with 34 years of data (equation 1):

where i(d,Rt) is the rainfall intensity i (mm) associated with a given d duration (min) and Rt return period (years); a and b are parameters resulting from the adjustment between rainfall intensities and durations (associated with a given return period) by the least-squares method.

The hyetographs were designed using the alternating block method and are composed by blocks of 10 minutes. The alternating block hyetograph (Chow et al., 1988) is one of the methods based on the IDF curves, as the Chicago hyetograph (Keifer and Chu, 1957). These methods produce the highest peak flows, resulting in critical flood situations and reliable design values, but are not intended to represent real event patterns (Grimaldi et al., 2012; Sordo-Ward et al., 2014; Na and Yoo, 2018; Chimene and Campos, 2020). The duration of the estimated storms equals the time of concentration of the Barcarena watershed. The Soil Conservation Service curve number (SCS-CN) for AMC III (antecedent moisture conditions -saturated soils) was chosen as the loss method and lag time as the transform method in HEC-HMS.

The peak discharges were obtained for the Barcarena and Massamá streams. For the Barcarena stream, peak discharges were estimated upstream and downstream of the confluence with the Massamá stream. Later, these values were inserted in HEC-RAS software for hydraulic modelling (Figure 2).

The quality of a Digital Surface Model (DSM) is affected by the spatial distribution, density and accuracy of the elevation data, the spatial resolution of the model and the interpolation methods (Hutchinson and Gallant, 2000; Zhang et al., 2016; Habib et al., 2020). The first version of the DSM for the study area was built in GIS (ArcMap software) using elevation data (contour lines and elevation points) at 1:2000 scale provided by the Oeiras City Hall. However, and despite the detailed scale, there is no elevation data along the stream channels, except when crossed by the contour lines, creating flat areas upstream and downstream this point, separated by a 2-m step. This will be an obvious misrepresentation of the channel’s morphology and of its longitudinal variation in the DSM, leading to consequent errors in the floodable areas, water depths and flow velocities. These weaknesses can be eliminated by creating a “hydrologically correct” DSM from the Topo to Raster tool of ArcMap. This allows the insertion of the location of the thalwegs in the DSM, i.e., the line of the lowest elevation in the stream bed, ensuring the transversal and longitudinal variations in elevation along the stream channels. Buildings were also included in the DSM as polygons with higher elevations than the surrounding areas to estimate their effects on the floodable areas, water depths and flow velocities. The survey and vectorization of the buildings were performed using orthophotos (resolution of 10 cm) provided by the Oeiras City Hall. The final DSM includes contour lines, elevation points, buildings and the position of the thalwegs (Figure 2). The correctness and accuracy of the DSM is essential for the quality of the hydraulic modelling results (Horritt and Bates, 2002; Werner, 2004; Merz et al., 2007; Lastra et al., 2008).

Hydraulic modelling was performed through the interconnection between HEC-RAS and HEC-GeoRAS for ArcMap in three steps: pre-processing, modelling and post-processing (Geravand et al., 2020). In the first step, stream centerlines, bank lines, flow path centerlines, cross-sections (XS cut lines) and the Manning’s roughness coefficients (Chow et al., 1988) for channels and floodplains were defined with HEC-GeoRAS using the DSM, according to the orientations of USACE (2009). Cross-sections were drawn perpendicularly to the streamlines and separated by 10 m, in average. Cross-sections were placed closer in confluences, bends or natural/artificial narrowing of the stream channels to capture changes in their geometry and morphology.

The input data generated in HEC-GeoRAS and the peak discharges obtained in HEC-HMS (Figure 2) were imported to HEC-RAS for modelling (second step). Peak discharges were used to determine steady flow water surface profiles for the 2-, 5-, 10-, 20-, 50- and 100-year R_T along the study area. In steady-state condition, water surface elevation and flow velocity were calculated by HEC-RAS for each discrete cross-section by solving continuity, energy and roughness (Manning) equations (Cook and Merwade, 2009). Mixed flow regime and the boundary condition normal depth slope were chosen for the HEC-RAS computation. There is only one infrastructure in the study area that is relevant to be included in hydraulic modelling: a hydraulic passage/culvert in the Massamá stream close to the confluence with the Barcarena stream (Figure 1).

The HEC-RAS results were then exported back to ArcMap (third step), where floodable areas, water depths and flow velocities were spatialized using HEC-GeoRAS (Figure 2). This ArcMap extension can interpret the velocities calculated at each cross-section and interpolate velocities for the entire flood perimeter using the DSM, producing quasi-2D flow velocity grids. Water depths and flow velocities are presented in a raster data structure with a resolution of 0.5 m. For each considered return period, HEC-GeoRAS creates raster files with water depth and flow velocity values for each pixel. Raster structure also allows the determination of the spatial distribution of water depths and flow velocities, to identify maximum values or to stablish spatial relationships between the two flow parameters. Facing the absence of peak flood discharges for the past events that would allow an approximate validation of flood extents for the considered return periods, the February 2008 event and the respective return period were taken as references. Resorting to the IDF curves of the SJT rain gauge, the rainfall that triggered this event was estimated in 20-year R_T for a duration equal to the time of concentration of the Barcarena watershed (Leal et al., 2019). Thus, the flood extents obtained for the 20-year R_T were confronted with the observation of flood marks, testimonies from the affected populations, newspapers’ information, television reports and photographs of the 2008 event to assess the quality of the hydraulic modelling results. The calibration of these results was performed through an iterative process, adjusting the DSM and the Manning’s roughness coefficients.

Pearson linear correlations were determined for the Barcarena and Massamá streams to measure the strength of the association between the longitudinal values of stream bed slope, channel width, flood width (i.e., width of flood extent for an event with a certain return period), water depth and flow velocity. Cross sections every 5 m were drawn perpendicularly to the flow direction to define the longitudinal variations in the channel widths and the flood widths for the events with minimum and maximum R_T considered: 2- and 100-year R_T, respectively. Longitudinal water depths and flow velocities were also determined every 5 meters (coincident with the width cross sections) and for the lines with the highest values during 2- and 100-year R_T events. The correlation coefficients and their significance were obtained using SPSS software.

The floodable areas for the 2-, 5-, 10-, 20-, 50- and 100-year R_T in the study area can be seen in Figure 1. The estimated peak discharges for the Barcarena (upstream and downstream of the confluence) and the Massamá streams are presented in Table 1. The small size of the Barcarena watershed and the location of the study area in the middle section of the watershed implies that peak discharges are less than 200 m³/s. The shape of the valley bottom and the natural and anthropic features of the stream channels, namely the existence of retaining walls, prevent wide floodable areas (Figure 1 and Table 1). Floodwaters remains confined to the stream channel in some sectors, not overflowing to the floodplain, even during high- magnitude events, which favors a high-energy environment with high flow velocities. The results also point to the interruption of the floodable area on the Massamá stream, where it is culverted under a road, causing the flow accumulation upstream of the hydraulic passage (Figure 1).

Floodable areas for a 2-year R_T event represent just over half (51%) of those referring to a 100-year recurrence (Table 1). The floodable area of a 50-year R_T event represents 93% of the floodable area for a 100-year R_T event. What seems to be a small increase in area (1980 m²) for an event of twice the magnitude (50- to 100-year R_T) represents an important increase in the height and volume of floodwaters (m³), expressed by the increment of water depth.

As expected, the highest water depths coincide with the stream channels and, more specifically, with the thalwegs. As the magnitude increases, there are increases in water depths and a consequent progressive occupation of the floodplains (Figure 3). The average water depth increases between 0.98 m (2-year R_T) and 1.54 m (100-year R_T). Five classes were defined to understand how the water depths are distributed in the floodable areas and how they vary in events with different return periods: ≤ 1 m, 1-2 m, 2-3 m, 3-4 m and > 4 m (Table 2). The percentage of area occupied by the first class (≤ 1 m) ranges between 52% for 2 years and 38% for the 20-year R_T. Its importance decreases until the 20-year R_T and rises again for the 50- and 100-year R_T. This is explained by the overflow of the Barcarena stream in several sections during events with higher return periods. Water depths are naturally lower in the floodplains than in the stream channel, which justifies the growing importance of the lower classes. In the second class (1-2 m), the occupied area decreases with a decrease in frequency, ranging between 42% for the 2-year R_T and 25% for the 100-year R_T. Between 2 and 3 m, the occupied area increases up to the 10-year R_T (22%) and decreases progressively at the highest magnitudes. In the last two classes (3-4 m and > 4 m), the relevance grows with the increase of the return period. However, the higher class is only represented in the 20-, 50- and 100-year R_T, while the second higher class is not represented in the estimated results for the 2-year R_T (Table 2).

As with the water depth, the highest flow velocities occur along the stream channels (Figure 4), in the sections where the bed slope is more pronounced. The greater roughness justifies the much lower velocities in the floodplains. The flow velocity also tends to be higher downstream and lower upstream of narrower sections, whether they are natural or artificial. This happens close to the hydraulic passage on the Massamá stream (Figures 3 and 4). There is flow accumulation upstream from this infrastructure, widening the floodable areas, while downstream there is a considerable increment in flow velocities. The visible damage in the Massamá stream bed (red circle in Figure 5) demonstrates the high velocities of the floodwaters and the destructive capacity of this tributary until the confluence with the Barcarena stream. On the other hand, flow velocity decreases when the channel becomes wider or there is an overflow of the watercourse into the floodplains. Flow velocity losses also occur upstream of the confluences. Reductions are identified in the Barcarena stream, upstream of the confluence with the Massamá stream, most notable for the higher return periods (Figure 4). This can be explained by the turbulent flow, typical of the convergence between two watercourses during flash floods.

The average flow velocities increase with increasing return periods, ranging between 1.55 m/s for the 2-year and 2.12 m/s for the 100-year R_T flood. The five classes (≤ 1 m/s, 1-3 m/s, 3-5 m/s, 5-7 m/s and > 7 m/s) show a more linear pattern when compared to water depths. Thus, higher velocities are more relevant in events of greater magnitude (Table 3). The only exception is the 1-3 m/s class, whose importance grows up to the 10-year R_T, reaching more than 50% of the total, and decreases thereafter.

Spatial correlations were determined between water depth and flow velocity in the floodable areas for the 2-, 5-, 10-, 20-, 50- and 100-year R_T using GIS. Correlation coefficients (R) range between 0.39 for the 5-year R_T and 0.50 for the 100-year R_T (Table 4). Despite the values are not particularly relevant, which means that the spatial relationships between water depth and flow velocity are not directly proportional, there is an increasing (linear) trend with return periods from 10- to 100-year R_T. The increasing discharges lead to the expansion of floodable areas and the progressive occupation of the banks by floodwaters. Both water depths and velocities tend to be low in these areas, unlike what happens along the stream channels, where longitudinal changes in channel geometry/morphology or flow turbulence can affect the flow behavior.

Additionally, two samples of 100 random pixels with a minimum distance of 5 m between them were selected in GIS. The first sample uses water depths and flow velocities of a 2-year R_T event and the second uses values of a 100-year R_T event (Figure 6). The correlation coefficients are like those presented in Table 4: 0.35 for the 2-year and 0.47 for the 100-year R_T (both significant at the 0.01 level).

By representing the water depth in a curve of increasing values and associating their respective flow velocities, i.e., obtained for the same 100 pixels, two distinct trends in flow velocity in both return periods are identifiable (Figure 7). There is a noticeable increasing trend in flow velocity during the first phase (light orange), following the growing line of water depth. In a second phase (dark orange), the flow velocity tends to decrease as the water depth increases. There is a moment when high water depths are no longer associated with high velocities. There is even a propensity for flow velocity to decrease when the highest water depths are reached. It is confirmed that there is no direct relationship between water depth and velocities. The values of flow velocity can differ by more than 4 m/s for similar water depths in the 2-year R_T (Figure 7a) and the 100-year R_T (Figure 7b).

The opposite exercise represents flow velocities in a curve of increasing values, associating their respective water depths (Figure 8). The results are similar, revealing a first phase (light blue) with an increasing trend of water depth as the flow velocity increases and a second phase (dark blue) with an inverse behavior. When high velocities are achieved it is less likely that water depths will also be high. Still, for similar velocities the results show that water depths can range by more than 2 m in the 2-year R_T (Figure 8a) and more than 4 m in the 100-year R_T (Figure 8b).

It is also important to analyze what occurs in the same pixel in all return periods to assess whether the increase in magnitude has consequences on the association between water depth and flow velocity. Of the 100 pixels in the random sample for the highest return period (100-year R_T), four represent different types of flow behaviors along the Barcarena stream channel (Figure 9). The location of these four points can be found in Figure 6. The constant increase in water depth and flow velocity with the increase in R_T occurs in 60% of the pixels, configuring the most frequent behavior (represented in Figure 9a). The increment in magnitude implies gradual increases in both flow components. However, this type of relationship rarely occurs for very high values. Figure 9b shows the typical flow behavior upstream of confluences and obstacles, in which velocity tends to decrease with increasing water depth. In this case, the highest flow velocity is reached for the 5-year R_T, decreasing continuously thereafter, while the water depth is constantly rising. Figure 9c shows an increase in the water depth up to the 50-year R_T and an abrupt loss in the 100-year R_T, while the flow velocity has a sharp rise. This occurs immediately upstream of an overbank flooding situation, in which the velocities achieved in the previous dozens of meters are substantially higher for the 100-year R_T when compared to the lower return periods (visible in Figure 4). The high flow velocity reached consequently leads to a shallower water depth at this point (upstream of the overbank flooding). In the latter type, both components grow progressively, but flow velocity starts to decrease after the 20-year R_T as a result of high values of water depth (Figure 9d), causing the overbank flooding (Figure 3). In short, Figure 9a illustrates the most common flow behavior, Figure 9b shows the behavior upstream of a confluence or obstacles, Figure 9c demonstrates the relationship between water depth and velocity immediately upstream of an overbank flooding situation, and Figure 9d represents the flow behavior when overbank flooding occurs.

The maximum values of water depth range between 2.77 m and 5.01 m (Table 5). The highest maximum value was obtained for the 50-year R_T and not for the 100-year R_T floods, although the difference is not significant (8 cm). This can be explained by the larger extension of the floodable areas for the 100-year R_T and by the way it affects the distribution of water depths. The maximum values of flow velocity reach 6.41 m/s for the 2-year and 8.52 m/s for the 100-year R_T (Table 6).

The location of the maximum values of water depth and flow velocity is different and dependent on the magnitude of the flash flood events (Figure 10). On the other hand, increases in the return period can change the location of the maximum water depth and the maximum flow velocity. This can be explained by the progressive flooding of the floodplain with increasing discharges and the consequent greater flood extents, causing changes in the behavior of the flow parameters.

Two locations were identified for the maximum values of water depth: one for the 2- to 20-year R_T and another for the 50- and 100-year R_T (Figure 10a). These are separated by about 120 meters. Three locations were obtained for the maximum velocities: 1) 2-year R_T; 2) 5 to 20-year R_T; and 3) 50 and 100-year R_T (Figure 10b). Except for the 2-year R_T, in which the maximum flow velocity occurs downstream of the hydraulic passage on the Massamá stream (Figure 10b), the maximum values for the remaining recurrences are close to each other in the Barcarena stream. The erosive capacity of both streams in these sectors during flash floods is shown in Figure 10. The maximum values of flow velocity for the 5-, 10- and 20-year R_T are located about 10 meters upstream of the maximum values for the 50- and 100-year R_T, representing the most critical section regarding flow velocity (Figure 10b). This occurs in a narrow section of the channel and downstream of a bend apex, confirming that the highest velocities in a curved bend can be found further downstream with increasing discharges, as mentioned by Harrison et al. (2015). On the other hand, this critical section is coincident with the lowest roughness segment, where the Barcarena stream is confined between stone walls, and is located downstream of the highest bed slope value in this watercourse: 9.9% (Figure 11).

The maximum water depth for the 2- to 20-year R_T (Figure 10a) is located upstream of the highest flow velocity section (Figure 10b). For the 50- and 100-year R_T, this point was found in a section where the floodable area is much wider when compared to what happens upstream and downstream (Figure 10a). Both points are in low bed slope sections: < 2% (Figure 11). In addition to the maximum values not being spatially coincident, the maximum water depths do not correspond to particularly high values of flow velocity and vice versa (Tables 5 and 6). On the other hand, the corresponding flow velocity is higher for events with a 20-year R_T (3.55 m/s) than for the 50 (2.37 m/s) or the 100-year R_T (2.83 m/s) due to the change in the location of the maximum values of water depth (Figure 10a and Table 5). Even so, all these values are far from the maximum flow velocities estimated for any of the return periods (Table 6), which also occurs for the water depths (Table 5).

Values of bed slope, roughness, channel width, flood width, water depth and flow velocity along the stream channels were obtained to understand their longitudinal variations in the Barcarena and Massamá streams. The results of the hydraulic modelling showed that, while the highest values of water depth follow the thalwegs, this is not always the case regarding the highest values of flow velocity. So, the lateral position of these lines is not always coincident along the channels.

Correlation coefficients were determined using the longitudinal values of bed slope, channel width, flood width, water depth and flow velocity for the 2- and 100-year R_T in the Barcarena and Massamá streams (Table 7). Unlike bed slope and channel width, which do not change with increasing return periods, flood width varies when the magnitude of an event increases. Thus, flood widths can only be associated with water depths and flow velocities with the same return period. Bed slope tends to be negatively correlated to water depth and positively correlated to flow velocity, which means that higher values of slope lead to lower depths and higher velocities. Stronger correlations were found for Massamá stream when compared to the Barcarena stream. Channel width presents a similar behavior to the bed slope (higher widths correspond to lower depths and higher velocities), although with more neglectable values, especially in the Barcarena stream. Flood width has a distinct behavior, with some noteworthy aspects: 1) higher flood widths tend to be more associated with higher depths and lower velocities, which is understandable because wide flood extents are caused by high water depths, resulting in decreases of longitudinal flow velocities; 2) the correlation coefficients are always higher for a 100-year R_T event when compared to a 2-year R_T event, meaning that the longitudinal variation of water depth and flow velocity are more dependent on the flood width during high-magnitude events; and 3) the correlation coefficient with the longitudinal flow velocity in the Barcarena stream for a 100-year R_T event is the highest presented in Table 7 (-0.68).

The influence of bed slope in the longitudinal variations (along the thalwegs) of water depths and flow velocities can be confirmed in Figure 11. Higher depths are spatially related to low values of slope, while the maximum velocities are coincident with high slope and low roughness sections. Nevertheless, other factors like channel width, flood width/overbank flooding or the presence of obstacles can affect longitudinal water depths and flow velocities. Despite the relationship between depth and velocity is not particularly strong when considering the floodable areas (Table 4), this becomes more relevant along the streamlines. There are strong negative relationships between longitudinal water depths and flow velocities, with the two flow components moving in opposite directions. The correlation coefficients (R) for the Barcarena stream reach -0.78 and -0.66 for the 2- and 100-year R_T events, respectively. More relevant values are achieved for the Massamá stream (-0.83 and -0.87) for the same R_T events because this is a smaller watercourse with limited overbank flooding, leading to a more direct relationship between the two flow components. The smaller flood width/extent also explains the higher correlation coefficients during a 2-year R_T event when compared to a 100-year R_T event along the Barcarena stream. If floodwaters are confined between the walls of the stream channel, there are fewer factors that interfere with their longitudinal flow behavior.

The authors would like to thank the Oeiras City Hall for providing the elevation data. This work was funded by the Research Unit (UIDB/00295/2020 and UIDP/00295/2020). Miguel Leal and Pedro Pinto Santos were funded by the FCT-Portuguese Foundation for Science and Technology, I.P. through the grant SFRH/BD/96632/2013 and the contract CEEIND/ 00268/2017, respectively. Eusébio Reis was financed by national funds through FCT, under the framework of the project BeSafeSlide-Landslide Early Warning soft technology prototype to improve community resilience and adaptation to environmental change (PTDC/ GES-AMB/30052/2017).

The data that support the findings of this study are available upon request to the authors. The primary data belong to the Oeiras City Hall and, therefore, cannot be made publicly available.