Using High-Density LiDAR Data and 2D Streamflow Hydraulic Modeling to Improve Urban Flood Hazard Maps: A HEC-RAS Multi-Scenario Approach

Abstract

The ability to extract streamflow hydraulic settings using geoinformatic techniques, especially in high populated territories like urban and peri-urban areas, is an important aspect of any disaster management plan and flood mitigation effort. 1D and 2D hydraulic models, generated based on DEMs with high accuracy (e.g., Light Detection and Ranging (LiDAR)) and processed in geographic information systems (GIS) modeling software (e.g., HEC-RAS), can improve urban flood hazard maps. In this study, we present a small-scale conceptual approach using HEC-RAS multi-scenario methodology based on remote sensing (RS), LiDAR data, and 2D hydraulic modeling for the urban and peri-urban area of Bacău City (Bistriţa River, NE Romania). In order to test the flood mitigation capacity of Bacău 1 reservoir (rB1) and Bacău 2 reservoir (rB2), four 2D streamflow hydraulic scenarios (s1–s4) based on average discharge and calculated discharge (s1–s4) data for rB1 spillway gate (Sw1) and for its hydro-power plant (H-pp) were computed. Compared with the large-scale flood hazard data provided by regional authorities, the 2D HEC-RAS multi-scenario provided a more realistic perspective about the possible flood threats in the study area and has shown to be a valuable asset in the improvement process of the official flood hazard maps.

1. Introduction

1.1. State-of-the-Art

In the last decades, with climate change and global warming, the associated natural disasters have reached more disastrous and catastrophic scales [1,2,3,4,5,6]. Many natural disaster patterns have diversified and modified under the pressure of climate change [7,8,9]. Among them, the intensification of the hydrological cycle has made an unprecedented impact on the magnitude, spatial extend, duration and frequency of hydro-meteorological disaster events [10,11,12]. Supporting this statement, many scientific publications support the fact that the occurrence of climate-related disasters (e.g., mainly floods, severe storms, cyclones, typhoons, droughts) have significantly increased under the abrupt changes in hydrological climatic conditions, the ecosystems resilience to transitional (wet–dry), and other related disturbances [13,14,15]. According to the Emergency Events Database (EM-DAT), in the last 50 years, climate related disasters have increased trifold between 1970 and 2017, from an average of 90 events in the period 1970–1989, to an average of more than 300 events in the period that followed after the year 2000 [16,17]. Among them, flooding phenomena are the most widespread, frequent, and costly natural disaster for the human societies [18,19], and is the most common natural hazard and the third most damaging hazard globally after storms and earthquakes [20,21].

In this context, North East of Romania is no exception. Crossed by the Siret and Prut rivers, which together constitute the biggest river basins in the country according to their surface, this territory is the most vulnerable area in terms of global climate change and modification of climate patterns. Modifications in the thermal, pluviometric and hydrological regimes have been noticed [22]. In the last 30 years, historical flood events were recorded in this area for the entire territory of Romania [23]. Floods have become a constant threat (especially after the year 2000) as, every two years, a major flood event was documented (e.g., 2004, 2005, 2006, 2008, and 2010) [24,25,26,27].

Besides climate change phenomenon implications, the high frequency of floods events is also determined by a whole range of socio-economic factors such as: Land use practices (e.g., changing floodplain functionalities), the development of socio-economic activities in flood-prone areas, changing living standards (e.g., urbanization—construction in high flood risk areas, urban sprawl), and poor regulations reflected by the incapacity of the responsible authorities to implement efficient planning policies [28,29,30,31,32]. The anthropic intervention on land has profoundly modified the natural standard behavior, and the changes on land use and land management are affecting, in particular, the river hydrology that determines the flood hazard [33]. This has ended up being a continuous increased trend in terms of vulnerability and exposure to disasters and flood hazard [34,35]. However, according to [36], the population has, voluntarily and under the pressure of modern society, exposed the environment to the possibility of flooding.

In Romania, there is a lack of proper flood defense measures against floods within many river basins. In general, the current flood defense infrastructure and hydro-technical constructions were built-up in the 1970s and very few improvements have been made in the years after. In the same time, over the past decades, the river regularization measures, the intensification of deforestation in the Carpathian Mountains, and housing development on the floodplains have changed the flood regime in the Siret and Prut river basins and their sub-basins. This has also the case for one of the most modified watercourses in the Siret river basin: Bistriţa River. Ten reservoirs, including canals, dams, collectors, transfer flows, protection works, and banks have been developed on the main course of the Bistriţa River, in order to mitigate the flood effects, and for the production of electricity using hydro-power plants [22,26,27]. Although no significant flood event has occurred on the Bistriţa River due to this considerable protection, the recent changes in pluviometric and hydrological regimes due to climate change, anthropic intervention on land, and reservoir clogging prove to be real challenges in the future flood hazard assessment of this area [22].

According to [37], human societies have always tried to reduce flood impact and have always sought protection from natural disasters by settling in safe areas or by building defense infrastructures. Throughout human history, there has been a constant preoccupation to understand, assess, and predict flood events and their impact [38]. Many ways were introduced in order to deal with flooding phenomena. In this context, aspects like streamflow hydraulic modeling, flood hazard assessment, and flood risk management have received increasing attention in the second half of the 20th century and has become a fundamental issue in the beginning of 21st century. Flood risk management can include: Hazard assessments, exposure assessments, vulnerability assessments, and risk assessments [39,40]. Flood hazard assessment is the most important step in the development of effective flood risk management; the flood hazard maps are used for danger estimation and for taking preventive conservation measures for risk mitigation. In this way, valuable information for different decision makers regarding spatial planning, and the design of infrastructure and emergency response preparation, become available.

In order to conduct flood hazard assessments and propose flood risk management measures, it is important to know what kind of comprehensive tools are available to be used in the analysis of flood susceptibility [41]. Generally, in flood hazard assessment, field-observed data from gauging stations are used. The existence of long records of accurate river discharge measurements is requested. However, these observations can have several uncertainties and the data are not always available [39,42]. Earth Observation (EO) datasets (e.g., space borne, aerial images, and satellite images) together with geographic information systems (GIS) can be used to determine the extent of flood areas and for the production of flood hazard and risk maps [43]. Combined with the use of remote sensing (RS) and GIS, EO datasets provide a favorable environment for relevant information processing in order to obtain the spatial extent of flood hazard areas and flood mapping [44,45,46,47]. Even if the benefits are high, this technique offers valuable information only for specific flood events, and this aspect is perceived as a disadvantage because future flood details and impact cannot be investigated [42,45].

Another available tool used in the determination of flood hazard areas is streamflow hydraulic modeling. The usage of this technique witnessed a considerable development in the recent years. As flood events have increased significantly over the last decades under the current unpredictable climate change behavior and human activities, the development of flood inundation models have become a necessity [39]. Streamflow hydraulic modeling is a useful and efficient tool in the simulation of flood events, identification of vulnerable areas, and estimation of spatially-distributed variables such as flow velocity or depth [48]. Since 1970, considerable efforts have been made in order to improve the capabilities and functionalities of this technique. Recent computational advances in hydraulic modeling offer new opportunities to support decision-making and adaptation [49]. The applicability of streamflow models have proved efficient in flood hazard and risk mapping studies [50,51], real time flood forecasting [52], and in remodeling past flood events [53]. Also, the availability of Digital Elevation Models (DEMs) [54] based on high-density Light Detection and Ranging (LiDAR) data have considerably improved the accuracy of flood parameters [47]. Depending on the selected mathematical method and scale, the complexity of models can be classified in simple interpolation methods (1D) and spatially detailed models which solve the water equations in 2D [55]. The use of a 1D simple model is recommended for channels, as long as the water remains in the steep slopes of the channel and no significant water spread is expected. 1D modeling is the perfect solution in cases of flood propagation along the main river [56,57]. Two-dimensional models are recommended in cases when the water is expected to overtop and the flow direction may change, spreading across a large area [42,58]. However, there are many works which make comparisons between 1D and 2D models [56,58,59,60].

1.2. Case Study: Aims and Objectives

In this study, we developed a method for flood vulnerability assessment under real (average discharge) and mathematical (calculated discharge) hydrological data based on HEC-RAS, LiDAR data, and 2D hydraulic modeling. Different from other studies [23,24,25,26,27,47], which provide flood hazard maps based on hydrological data calibrated at river basin scale [23,24,26,27], we developed for the first time in the study area flood hazard maps adapted to local environment settings and calculated discharge. In this context, four 2D streamflow hydraulic scenarios (s1–s4) were computed in order to test the flood mitigation capacity of hydro-technical constructions located downstream on the Bistriţa River (North East Romania) (Figure 1a). The scenarios were based on different discharge releases (average discharge; calculated discharge: s1–s4) from the Bacău 1 reservoir (rB1), also known as the Lilieci storage reservoir, and from its hydro-power plant (H-pp) (Figure 1b). The average and calculated discharges were correlated with the official operating regulations of rB1. The results provided a more realistic perspective about the possible flood threats located downstream: rB1, an area which overlays with the urban and peri-urban area of Bacău City; and rB2, with the Bacău 2 reservoir (Figure 1c).

2. Methodology

2.1. Selected Site for 2D Streamflow Modeling

The study area selected for 2D streamflow hydraulic modeling and for urban food hazard assessment covers 29.73 km² in the urban and peri-urban area of Bacău City (Bacău County, North East Romania) (Figure 1a). The town has a population of 196,883 inhabitants according to the 2016 census, making it the 12th largest city in Romania. The considered Bistriţa River reach extends 2 km upstream from the confluence of the Bistriţa and Siret rivers (Figure 1b,c). Due to more than 10 large reservoirs and other complex hydro-technical works (e.g., dams, water gates, hydro-power plants, channels, headraces, flood protection dikes) which equip the downstream sector of the Bistriţa watershed, the predictive flood models are difficult to achieve through classical methods. According to [22], even if the Bacău City was included on the settlements list potentially affected by floods in North East Romania, the flood hazard maps provided by the authorities would not highlight this state. For this reason, in the last 10 years, new neighborhood areas have been built-up in the floodplain of Bistriţa River or in the proximity of the rB1 and rB2 reservoirs.

2.2. Data Acquisition

Figure 2 summarizes the workflow chart followed in this study with the hydrological, LiDAR, and built-up data obtaining process, and the key steps for HEC-RAS 2D streamflow and flood modeling.

2.2.1. Development of DEM

The DEM data used in this study comprehends airborne LiDAR technology. LiDAR DEMs display the terrain with a high degree of horizontal and vertical accuracy, a feature that is essential to terrain-related applications such as streamflow and flood modeling [54,61,62]. The data was obtained from the National Administration Romanian Waters (NARW)–Siret Water Basin Administration (SWBA) and consisted of raw ground point elevation data, at a spatial density between 0.5 and 3 m. More than 7.833 × 10⁶ points covering almost 100 km² was considered in the DEM computation (Figure 2a). The final product (DEM) was obtained by means of natural neighbor interpolation technique [63]. The natural neighbor interpolation technique was preferred as it finds the closest subset of input samples to a query point and applies weights to them based on proportionate areas in order to interpolate a value [64]. To capture as many topographical details as possible, the DEM was computed at 0.5 m spatial resolution (Figure 2a).

As the considered study area overlaps a highly urbanized space, crossed by one of the most developed hydro-technical rivers in the country, an overview regarding the current situation of the urban built-up areas was necessary. Initially the open source database offered by OpenStreetMap [65] was consulted. As the available database proved incomplete, the last step in built-up data acquisition consisted of the manual digitalization of the missing built-up areas using a high resolution orthophotos digital image (0.5 m resolution), dating from 2012. In this way, the minimum and maximum polygonal areas of the improved built-up vector was situated between 10 m² and 19,534 m² (Figure 2b). Furthermore, in order to assure the hydraulic streamflow accuracy, the obtained vector was integrated into the DEM. The new data was rasterized and then added to the already obtained LiDAR DEM, taking into consideration an average height of 10 m. The built-up areas affected by the inundation were identified using the built-up vector obtained in the previous steps. All residential and industrial built-up areas that intersected or shared a common spatial extend with the inundation extend, were considered impacted (Figure 2b).

2.2.2. Hydrological Data

The hydrological data used in this study was extracted from the Lilieci reservoir (rB1 and H-pp) official operating rules dated from 2012 and offered by SWBA (Table 1). According to Romanian regulations STAS—4273/61 and STAS 4273/82, the Lilieci reservoir (rB1) is classified after the height of the dam and its storage capacity volume, in the third importance category—having the flow sizing rate Q2% (50-years recurrence interval) = 510 m³/s, and Q0.5% (500-years recurrence interval) verification flow rate of 1670 m³/s. Nevertheless, the rB1 spillway gates (Sw1) and levees were designed to transit flow rates that correspond the second importance category—Q1% (100-years recurrence interval) = 1140 m³/s and Q0.1% (1000-years recurrence interval) = 2100 m³/s. The H-pp is located at approximately 1726 m from the Sw1, on a specially designed channel and it is equipped with two Kaplan turbines with a 180 m³/s maximum flow rate.

2.3. 2D Hydraulic Modelling

In relation to watercourses (e.g., streams, rivers, channels), the hydraulic models describe water movement through space in three directions [66]: A downstream direction along the river channel, a lateral direction (e.g., whenever the water begins to spill out overland), and a vertical direction which practically defines the height of the flood. Standard flood modeling practices include 1D modeling (upstream to downstream direction), 2D modeling (downstream and lateral directions), hybrid 1D–2D modeling and 3D numerical models, along with hydrograph design, specified ground roughness, and accurate digital elevation data [39,67]. For the present study, a 2D approach was adopted for a number of reasons (Figure 2c):

• The main purpose of the study was to develop overtopping multi-scenarios in order to define the areas with risk of flooding in a complex urban environment; where it cannot be assumed that all flow will be parallel to the main river and where a higher hydrodynamic accuracy is required [68] (Figure 2d,e).

• Considering that the water will not propagate in the vertical direction, the vertical water velocities will be negligible compared to the horizontal velocities and also the kinetic energy losses due to the vertical movement. That is why this study has considered a 2D model to develop its conclusions.

• The highly accurate representation of the river’s bathymetry and floodplain topography, which is represented by a high resolution digital elevation’s model.

• Two-dimensional flood propagation modeling offers additional information regarding some characteristics of the flood, such as flow velocity and water trend propagation.

The 2D functionality builds flood models in a more accurate way, considering the flow variability in time, and both spatial dimensions, the x and y, making it more suitable for case studies that consider wide floodplains, urban environments, dam/levee breach situations, where the water is expected to spread over an open, unconstrained area, in multiple directions [68]. The release of the HEC-RAS software (5.0 version), in 2016, integrated 2D unsteady flow capabilities, offering the possibility to analyse water propagation over a predefined surface, which is found in the form of a digital elevations model [69,70]. HEC-RAS is a well-known software, capable of modeling a flood inundation event [39,42]. The scenarios proposed by this study were simulated using the open source HEC-RAS software (5.0.3 version), developed by U.S. Army Corps of Engineers (USAGE). Starting with HEC-RAS software version 5.0, two-dimensional unsteady water flow modeling can be performed. The program 2D flow modeling algorithm solves either shallow water equations, also called bidimensional Saint Venant equations (Equation (1)), or the 2D diffusion wave equations (Equations (2) and (3)).

(1)$\frac{\partial \zeta }{\partial t}+\frac{\partial p}{\partial x}+\frac{\partial q}{\partial y}=0$

(2)$\frac{\partial p}{\partial t}+\frac{\partial }{\partial x}\left(\frac{p^2}{h}\right)+\frac{\partial }{\partial y}\left(\frac{pq}{h}\right)=-\frac{n^{2}pg\sqrt{p^2+q^2}}{h^2}-gh\frac{\partial \zeta }{\partial x}+pf+\frac{\partial }{\rho \partial x}\left(h{\tau }_{xx}\right)+\frac{\partial }{\rho \partial y}\left(h{\tau }_{xy}\right)$

(3)$\frac{\partial q}{\partial t}+\frac{\partial }{\partial y}\left(\frac{q^2}{h}\right)+\frac{\partial }{\partial x}\left(\frac{pq}{h}\right)=-\frac{n^{2}qg\sqrt{p^2+q^2}}{h^2}-gh\frac{\partial \zeta }{\partial y}+qf+\frac{\partial }{\rho \partial y}\left(h{\tau }_{yy}\right)+\frac{\partial }{\rho \partial x}\left(h{\tau }_{xy}\right)$

where, h is the water depth (m), p and q are the specific flow in the x and y directions (m²s⁻¹), $\zeta$ is the surface elevation (m), g is the acceleration due to gravity (ms⁻²), n is the Manning’s Roughness coefficient (m^−1/3 s), $\rho$ is the water density (kg m⁻³), ${\tau }_{xx}$, ${\tau }_{yy}$, and ${\tau }_{xy}$ are the components of the effective shear stress and f is the Coriolis (s⁻¹). When the diffusive wave is selected, the inertial terms of the momentum equations are neglected: $\partial p/\partial t+\partial /\partial x\left(p^2/h\right)+\partial /\partial y\left(pq/h\right)=0$ (Equation (2)); $\partial q/\partial t+\partial /\partial y\left(q^2/h\right)+\partial /\partial x\left(pq/h\right)=0$ (Equation (3)) [42,45,71].

Due to their faster computational time and greater stability properties, the 2D diffusion wave equations were preferred in the present study [70]. In this case, the floodplain flow was approximated as a two-dimensional diffusion wave where water can flow in any direction based on the defined topography and the resistance to flow determined by the type of land use [72,73]. A 2D model consists of a 2D computational grid which groups/discretizes the river and its adjacent areas into a collection of individual and connected cells, called grid cells (2D flow cells), that are used to characterize the underlying topography [56,74]. The 2D flow area is the region (the boundary selected by the user) where the water flow will be modeled.

A close polygonal mesh with a computation point spacing of 5 m was built in order to define the spacing between the computational grid-cells. This managed to capture all the LiDAR terrain characteristics. In this way, more than 2,000,000 grid cells were generated; therefore, no break line was needed to be introduced (e.g., top of the levees, along the main channels). The computational mesh controlled the water flow throughout the 2D flow designed area and from cell to cell. Considering the underlying terrain and the generated computational mesh, the software developed detailed elevation–volume relationships and detailed hydraulic property curves for each cell face (elevation vs. wetted perimeter, area, and roughness) [70]. In this way, the water elevation values were calculated for each centroid of the grid. The Manning’s Roughness coefficients, specific to each landcover class was set according to [75] (Figure 3).

After the geometric data is created and the roughness coefficient is set, the next step consists of establishing the boundary conditions. This means that flood data along the boundaries of the 2D flow area is assigned. Flow hydrographs and normal depth boundary conditions representing the average riverbed slope are used in order to bring water into the 2D flow area. Four flood scenarios were simulated and, for each scenario, two hydrographs and two normal depths boundary conditions were introduced. The hydrographs were 24 h long with values recorded each hour. The upstream boundary conditions (the hydrographs representing imposed flow condition) were located in two specific locations downstream of Lilieci dam: At the hydropower plant of Lilieci reservoir (on the controlled channel), and on the natural river bed (where water from the reservoir is released). The imposed water level boundary condition was located at the downstream extremes of the channel and at the end of storage lake control channel (in the limits of DEM). An energy slope of 10⁻⁴ mm⁻¹ was used.

In order to ensure the stability of the model, the time steps were estimated according to the Courant–Friedrichs–Lewy condition [70] (Equation (4)).

(4)$C=\frac{V\Delta T}{\Delta x}\le 1.0 \left(\mathrm{with} C_{max}=3.0\right)\mathrm{or}\Delta T\le \frac{\Delta x}{V} \left(\mathrm{with} C=1.0\right)$

where, C is the Courant number, V is the flood wave velocity (m/s), ΔT the computational time step (s), and Δx the average cell size (m) [67]. According to Equation (4), a time step of 10 s was selected. To run a 24 h simulation at a time step of 10 s, the 2D models took between ~3 h (for Q_avg. scenario generated based on average discharge) and ~12 h (for each s1–s4). The flood inundation (depth), flood velocity, and other related results were obtained for each hour according to the hydrograph output interval.

2.4. Streamflow and Flood Simulation Accuracy

For accuracy purposes, a 2D water flow test simulation was computed. The data needed for this flow accuracy assessment test consisted of average discharge values at the Lilieci reservoir (rB1) for the year 2012 and, respectively, an orthophoto digital image (Figure 4a) and a Landsat 5 satellite image (Figure 4b) showing the water extend in the area acquired in the same date. The results of this simulation showed a close agreement between the water extend obtained from the model computation and the available orthophoto digital image (Figure 4c). The performance assessment was realized by comparing these three sets of spatial data. A difference of 5% (16 ha) was recorded between the orthophoto digital image water extend and the 2D average discharge simulation result. A second comparison was carried out using the Landsat 5 satellite image, acquired on 26 August 2012. The water extend were extracted according to [23] and compared with the 2D average discharge simulation result. In this case a difference of 9% (29 ha) was observed. Taking into consideration the presence of vegetation which could have masked the presence of water in the digital orthophoto images and the coarse resolution of Landsat scene, the accuracy of the flood extent assessment was judged satisfactory.

2.5. Multi-Scenario Development

The flood scenarios were based on hydrological data extracted from the Lilieci reservoir (rB1) official operating rules. Upstream boundary conditions for all scenarios were developed on a simulation of water flow over a period of 24 h, at a temporal resolution of 1 h. The average discharge scenario, spatially displays the main channel at a water release of Q = 34 m³/s located at the Lilieci hydro-power plant (H-pp) and a Q = 2.8 m³/s flow located at the reservoir main gates (Sw1), where the water is released directly in the main river bed (in order to assure the river servitude discharge). Scenario 1 (s1) takes into consideration the average discharge at H-pp of the Lilieci reservoir (Q = 34 m³/s) and the maximum water release in the case of the full opening of one gate (Sw1) of the same storage lake (Q = 490 m³/s). The maximum values were used in the simulation for a period of 1 h, whereas for the remaining 23 h, the flow rates were set to their average values (Q = 2.8 m³/s). Scenario 2 (s2) considers water flow modeling of Q = 34 m³/s at the H-pp and Q = 980 m³/s at the Sw1 when two gates are fully opened. The distribution of time periods was the same as for s1. For scenario 3 (s3), the same parameters were used as for s1 and s2, with the exception that a water flow of Q = 1470 m³/s is assumed when all 3 gates of the Sw1 are opened. The last scenario 4 (s4) simulated the water flow when all 4 gates of the Sw1 are fully open at a discharge rate of Q = 1960 m³/s. The other parameters remained unchanged. In Table 2 are the simulation parameters for the 4 modeled scenarios.

3. Results

3.1. Flood Pattern

The 2D modeling was performed using the HEC-RAS 5.0.3 version which offered various output possibilities in terms of detailed animation and mapping of flood characteristics within the RAS mapper feature. Following the calculation of the 4 scenarios given in Table 1, individual layers of the flooded areas along Bistriţa River and Bacău City regarding the flood extent, velocity, water surface elevation and depth, were exported.

3.1.1. Flood Extent

The flood extent layer consists of an Inundation Boundary shapefile vector layer capturing the areas affected by the flood during the whole 24 h simulation period (Figure 5). According to s1, an area of 6.27 km² within the Bistriţa floodplain, the inflow area (fan-delta) of rB2 and downstream sector of rB2 reservoir, are potentially affected by floods. Also, six buildings located within the s1 flood extent are potentially affected. In the s2 simulation, the total flood extent occupies 8.66 km² and 173 buildings in the peri-urban area of Bacău City which are located in the flood zone. In the third scenario (s3), the water release from Sw1 and H-pp spills over protection dikes and floods the built-up area of Bacău City. The total flood extent increases to 11.46 km² (built-up area—1.09 km²) and 2461 of buildings (e.g., block with flats, houses and yards, other residential buildings, industrial units) are potentially affected by the inundation waters. The last scenario (s4) indicates the most catastrophic situation where total flood extent is 12.82 km², from which 1.5 km² are built-up area, and 3780 residential buildings are vulnerable to the flood hazard (Table 3). Overall, in the first two scenarios (s1 and s2), the flood extent affected only the built-up areas within the Bistriţa floodplain, and in the last two scenarios (s3 and s4) the flood extent occupied the peri-urban and urban areas.

3.1.2. Flood Velocity

Another important result of the 2D modeling process is flood velocity. This is computed by recording the maximum velocity for each cell in the computational mesh regardless the time when it was measured in the whole 24 h simulation period (Figure 6). According to all computed scenarios (s1–s4), the 0.01–1 m/s velocity class has the highest frequency (88.3%), which corresponds with the average velocity recorded downstream rB2 (0.6 m/s). In the first two scenarios (s1 and s2), all vulnerable buildings are potentially affected by the 0.01–1 m/s flood velocity class, except three buildings located in front of Sw1 or along the natural watercourse of Bistriţa River. In the case of s3, only 9% of total affected buildings are threatened by 2 m/s flood velocity, and in case of s4, only 14% of total built-up area is potentially affected by high velocity of water (1–5 m/s). However, due to floodplain roughness and hydro-technical works within study area, the high flood velocities (>5 m/s) are registered only at the Sw1 and Sw2 gates, and through the narrow sections between rB1 and rB2 river sector (Table 4). Overall, in each computed scenario, the flood velocity is not the main threat for built-up area, but the river morphology (e.g., river bed, banks, alluvial deposits), and the hydro-technical works conservation stage (e.g., dams, channel, bridges) can alter through time. For these reasons, the maintenance of the flood mitigation equipment within urban and peri-urban areas must be a priority for the competent authorities.

3.1.3. Flood Depth

The maximum flood depth maps were generated by the 2D model by taking into consideration the maximum depth for each cell, no matter the time when that maximum depth was registered during the whole 24 h simulation period (Figure 7). According to s1 and s2, 95.5% of vulnerable buildings are potentially affected by a flood depth which does not exceed 1 m depth. In case of s3, when water levels exceed the height of protection dikes and inundates the urban and peri-urban area of Bacău City, 1174 buildings will be affected by floods which do not exceed 1 m in depth; 39.4% of the vulnerable built-up areas are potentially affected by floods with a depth between 1–2 m, and 12.7% by floods with a depth between 2–3 m. In the same scenario (s3), only 4 buildings are located in areas with water that exceeds 3 m. In case of s4, the most catastrophic scenario taken into account in this study, 34.7% of buildings would be affected by floods with <1 m depth, 30.9% by floods with depths between 1–2 m, 25.9% by floods with depths between 2–3 m, 7.9% by floods with depths between 3–4 m, and the rest of 0.4% (18 buildings) by floods that exceed >5 m depth (Table 5).

3.2. Flood Hazard Assessment

Usually, the flood hazard assessment is based on quantifiable variables like flood extent, water velocity or water depth, and indicate the vulnerability of built-up areas to hydrological events with possible destructive impact. In this study, we provided the flood hazard assessment using only the flood extent and flood depth resulting from the four HEC-RAS 2D scenarios (s1–s4). The flood velocity was not taken into account because is a constant variable in all four flood scenarios (88.3% from total built-up area are affected by the 0.01–1 m/s velocity class). To generate the flood hazard categories, the water depth for each flood extent was classified according to the Japanese criteria of the Ministry of Land Infrastructure and Transport (MLIT) [42]. The criteria suggest five flood hazard categories: H1—very low hazard (flood depth < 0.5 m); H2—low hazard (flood depth between 0.5–1 m); H3—medium hazard (flood depth between 1–2 m); H4—high hazard (flood depth between 2–5 m); H5—extreme hazard (flood depth > 5 m) (Figure 8). The flood hazard classification methodology and hazard description based on flood depth is detailed in Table 6.

According to s1, even if all 5 hazard classes are encountered, all vulnerable buildings are situated in the very low class of hazard (H1). In case of s2, 86.7% of total vulnerable buildings are very low affected by flood (H1), 8.75% are low affected by flood (H2), and the rest of 4.6% are located in the medium exposed area to flood hazards (H3). In s3 situation, the number of potentially affected buildings is significantly increased: H1—674 buildings; H2—500 buildings; H3—971 buildings; and H4—316 buildings. Also, this is the moment when the floods will affect the urban area of Bacău City. According to s4, more the 34.3% of total flood extent are located in the high (H4) and extreme (H5) hazard classes, the number of potentially affected buildings exceed the all NARW scenarios based on recurrence flood interval probabilities: H1—694 buildings; H2—618 buildings; H3—1171 buildings; H4—1279 buildings; and H5—18 buildings (Table 7; Figure 9a,b).

4. Discussion

Floods can cause a tremendous amount of damage, especially in the highly regulated river sectors overlaid with highly urbanized areas along the floodplain. For this reason, the multiple viewpoints, methods, assumptions, and future possibilities adapted to new trends (social, economic, natural) and determined by different factors, play an important role in the flood hazard management and establishing flood-vulnerable areas [39]. Even in a highly regulated river system, flood hazards exist; properties and lives are in danger. A flash flood, dense precipitation, or an error in the discharge flow at the gates of the reservoirs can turn into a catastrophic flood event. For this reason, a good preparation regarding this topic is always welcomed. As remote sensing (RS) techniques continue to improve and the availability of data increases, more RS data will be integrated and used in flood modeling. In this context, according to [54], 2D LiDAR DEM based flood simulations provide the best results. The horizontal, vertical accuracy is better on flat terrain (such as wide floodplains) and the integration of vegetation and building heights offer the perfect support for flood hydraulic simulations and accurate water flow propagation.

Due to the unavailability of LiDAR surface elevation points, we were not able to perform the simulation on a Digital Surface Model (DSM with vegetation and building integration). We computed the DEM and manually integrated the building heights. Manually integration of building heights (in the case of data unavailability) can be a solution. Prior to this step, we wanted to check if any open source satellite imagery could be used in order to extract the built-up areas. Initially, a Landsat 8 satellite image dated 11 August 2018 was downloaded using the US Geological Survey (USGS) website [76], pre-processed and processed—including radiometric and atmospheric calibration, maximum likelihood classification (MLC), classification accuracy assessment—in order to extract the areas containing artificial surfaces, residential, commercial, industrial, and transportation infrastructure. A full methodology of the previous mentioned steps can be consulted in the works of [77]. In this way a full perspective regarding the current situation of the urban development throughout the floodplain can be obtained. However, even though this method is practical and can be used in areas where no built-up data is available, the results proved to be at a coarse resolution and impracticable for the level of study we proposed. Therefore we decided to extract the built-up areas at a higher resolution. We manually digitized the missing built-up areas using an orthophotos digital image, dating from 2012. The association of open source imagery (Landsat, Sentinel) and LiDAR data offered unsuitable results.

The 2D flow models and LiDAR DEMs prove to be the perfect combination when it comes to the flood hazard assessment and accurate inundation delineation. Due to its accurate representation of the complex hydraulic conditions that can be found in floodplains (e.g., channels, confluences, bridges, water reservoirs, roads, etc.), 2D hydraulic models can capture the hydraulic behavior of the river in a more accurate way. Water propagation, extent, velocity, and elevation can be reproduced under real and hypothetical flow data and boundary conditions input. In this way we can simulate for different purposes the past flood events, using registered hydrographs at gauging stations, or we can imagine a theoretical hypothesis and, in this way, improve flood hazard and risk maps. This method is the best option in case of flood map creation due to possible reservoir dam failure at various dam discharges. On the other hand, the unavailability of high resolution DEMs can end up underestimating the flood hazard and its dangers. A precise representation of reality is mandatory for accurate flood mapping process.

Our study was focused on a highly regulated and well-equipped hydro-technical river sector which overlays the urban area of Bacău. This proved to be an impediment in the validation process of the 2D hydraulic modeling, as no notable flood event was recorded here in the recent years. Controlled water discharges from the reservoirs located upstream protect the areas downstream. For this reason, a considerable development of built-up areas and urban sprawl was noticed downstream, on both sides of the riverbed. Having only the average official discharges for the year 2012 available at the H-pp and Sw1 to model, we still managed to model the flood hazard in the area based on the official operating rules of the dam Sw1 and taking into consideration the reservoir maximum water retention capacity (Figure 10a).

Furthermore, we managed to compare the 2D flood extent results with the official available hazard and risk maps available in [78], which was realized under the Directive 2007/60/EC of the European Parliament and Council [79,80] on the assessment and management of flood risks (Figure 10b). The comparison results showed similarities even when the input material for their computation were at different scales and resolution. According to the National Administration “Romanian Waters”, the hazard and risk maps for Bistriţa River were realized for the entire catchment at a spatial resolution between 5 and 10 m, while our analysis used a local 0.5 m input DEM where all the topographical details and hydraulic conditions were taken into consideration. Additionally, in a highly regulated river sector where several reservoirs are present, such as in our study area, the flood hazard is governed by their operating rules and management. Any mistake in the manipulation of water discharge can prove catastrophic for the areas downstream.

5. Conclusions

The 2D hydraulic modeling using HEC-RAS (v 5.0.3) and high-density LiDAR data applied for streamflow simulation within the urban and peri-urban area of Bacău City (North East Romania), produced sufficiently accurate information regarding flood hazard vulnerability. Based on multi-scenario development using different discharge releases from the rB1 (average discharge; calculated discharge: s1–s4) and its hydro-power plant (H-pp), the following concluding remarks can be summarized:

• Combining 2D hydraulic modeling with RS LiDAR data and local hydrological parameters proved to be an efficient method in order to improve the streamflow and spatial distribution of water over a densely urbanized and well-equipped hydro-technical area, in case of a flood event triggered by different discharges from the rB1 and its hydro-power plant (H-pp).

• The multi-scenario approach allowed the testing of flow capacity throughout the rB1 downstream floodplain, determining in this way, the discharge value which can cause a flood event. s1 and s2 scenarios kept the flood extent inside the river levees, meaning that a discharge up to 1160 m³/s for a period of one hour can be handled by the hydro-technical works downstream. The flood only starts to spread over the urban and peri-urban areas of Bacău City at a discharge rate over 1160 m³/s (s3 and s4).

• Scenario 2 (s2) allows a discharge of 6,636,960 m³ of water, in a 24 h period (according to the input hydrograph) and captures the flood extent in the case where almost all the contained water in rB1 (Lilieci reservoir) is discharged. At a water surface elevation of 178.64 m and taking into consideration the existing LiDAR DEM, the reservoir stores 6,770,000 m³ of water.

Overall, the 2D hydraulic modeling multi-scenario results can be exported into a set of flood hazard parameters such as flood depth, flood extent, flood velocity, or water surface elevation; and can answer real questions regarding the flood hazard threat at local level. Developing streamflow scenarios on a small-scale level is a very important aspect for any flood mitigation effort, especially in the urban areas located along main rivers, because large-scale analysis (river basin analysis) understates flood risk perception. Due to the adaptability of the 2D streamflow hydraulic model proposed in this study, the method can become a valuable asset in flood mitigation. However, a LiDAR DEM based 2D flood simulation is essential for every urbanized environment in the context of climate change and modern society development pressure and trends.

Author Contributions

Conceptualization, A.M.-P. and C.I.C.; methodology, A.M.-P., C.I.C. and M.N.P.; software, A.M.-P., C.I.C., C.C.S. and M.N.P.; validation, A.M.-P., C.I.C., C.C.S. and M.N.P.; formal analysis, A.M.-P. and C.I.C.; investigation, A.M.-P., C.I.C., C.C.S. and L.E.P.; resources, A.M.-P. and L.E.P.; data curation, A.M.-P. and C.I.C.; writing—original draft preparation, A.M.-P. and C.I.C.; writing—review and editing, A.M.-P.; supervision, A.M.-P.; project administration, A.M.-P.; funding acquisition, A.M.-P.

Funding

This work was funded by the Ministry of Research and Innovation of Romania within Program 1—Development of the national RD system, Subprogram 1.2—Institutional Performance—RDI excellence funding projects, Contract no.34PFE/19.10.2018 (beneficiary: A.M.-P.).

Acknowledgments

The authors would like to express their gratitude to the employees of the Romanian Waters Agency Bucharest, Siret Water Basin Administration Bacău, who kindly provided a significant part of the LiDAR data used in the present study. All data was processed in the Geoarchaeology Laboratory (Coordinator: A.M.-P.) of Institute for Interdisciplinary Research, Science Research Department, “Alexandru Ioan Cuza” University of Iaşi (UAIC), Romania. Our thanks go to the all anonymous reviewers, who helped us in improving the manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The founding sponsors had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, and in the decision to publish the results.

Abbreviations

The following abbreviations are used in this manuscript:

DEM’s	Digital Elevation Models
DSM	Digital Surface Models
EM-DAT	Emergency Events Database
EO	Earth Observation
GIS	Geographic Information System
HEC-RAS	Hydrologic Engineering Canter’s—River Analysis System
H-pp	Hydro-power plant
LiDAR	Light Detection and Ranging
MLC	Maximum Likelihood Classification
MLIT	Ministry of Land Infrastructure and Transport (Japan)
NARW	National Administration “Romanian Waters”
rB1	Bacău 1 reservoir (Lilieci storage reservoir)
rB2	Bacău 2 reservoir
RS	Remote Sensing
s1–s4	2D streamflow hydraulic scenarios used in this study
SWBA	Siret Water Basin Administration
Sw1	Bacău 1 dam spillway gates
Sw2	Bacău 2 dam spillway gates
USACE	U.S. Army Corps of Engineers
USGS	US Geological Survey

Figures and Tables

Figure 1. (a) Geographic location of the study area in North East Romania; (b) Elevation map and, (c) Digital Elevation Model (DEM)-derived slope map of the urban and peri-urban areas of Bacău City. Data collected by the National Administration Romanian Waters (NARW) in 2014; DEM derived from ground point’s elevation data (0.5 point density) collected using Light Detection and Ranging (LiDAR) techniques from aircraft. The abbreviations within the (b,c) maps show the locations of the study site boundaries used in the 2D hydraulic modeling for the Bistriţa River, Bacău 1 (rB1) and Bacău 2 (rB2) reservoirs, Bacău 1 spillway (Sw1) and Bacău 2 spillway (Sw2) gates, and the Bacău 1 hydro-power plant (H-pp). All spatial data was defined and projected in the Romanian national projection (STEREO 70).

Figure 2. Workflow chart of 2D streamflow and flood hydraulic modeling process using hydrological data (Sw1) and LiDAR ground point elevation (0.5–3 m spatial density) within the urban and peri-urban area of Bacău City: (a) Generate the LiDAR based DEM with 0.5 m spatial resolution based on LiDAR ground point elevation; (b) manual digitalization of buildings using a high resolution orthophotos and built-up data integration in DEM; (c) generate the four 2D multi-scenario (s1–s4) using HEC-RAS v.5.0.3, where Qavg. was generated based on average discharge at Sw1 for 2D streamflow accuracy, and s1–s4 was computed using calculated discharge at Sw1 for 2D flood multi-scenario development; (d) generate the flood pattern for each computed scenario and export the individual layers (e.g., flood extent, flood velocity, flood depth); (e) flood hazard assessment based on the built-up data and flood depth classification according to the Ministry of Land Infrastructure and Transport (MLIT) [42].

Figure 3. Landcover map within the study area and Manning’s Roughness coefficient n (m−1/3 s) value for each class according to [75].

Figure 4. Streamflow 2D accuracy: (a) orthophoto (digital image) for 2012, (b) Landsat 5 (satellite image) for 2012, and (c) results of 2D average discharge simulation result. In all images was highlighted the water extend within inflow area (fan-delta) of the rB2 reservoir for comparison.

Figure 5. Flood extent derived from the HEC-RAS 2D multi-scenario generated based on LiDAR data and average discharge; s1–s4 was computed using calculated discharge (see Table 2).

Figure 6. Flood velocity derived from the HEC-RAS 2D multi-scenario generated based on LiDAR data and average discharge; s1–s4 was computed using calculated discharge (see Table 2).

Figure 7. Flood depth derived from HEC-RAS 2D multi-scenario generated based on LiDAR data and average discharge; s1–s4 was computed using calculated discharge (see Table 2).

Figure 8. Flood hazard maps based on flood depth classification according to the MLIT [42] (see Table 5). Data derived from HEC-RAS 2D multi-scenario generated based on LiDAR data and average discharge; s1–s4 was computed using calculated discharge (see Table 2).

Figure 9. The distributions of flood hazard classes within (a) flood extent area, and (b) built-up area, according to MLIT [42] criteria (see Table 6). The statistics apply for all 2D multi-scenarios (s1–s4) in each chart.

Figure 10. Comparison between flood extent maps within the urban and peri-urban area of Bacău City: (a) 2D HEC-RAS multi-scenario based on average discharge and calculated discharge (s1–s4) developed in these study, and (b) official flood extent maps (NARW) [78] based on large-scale hydrological and spatial data, calculated for 10% (10-year), 5% (20-year), 1% (100-year), and 0.1% (1000-year) recurrence intervals.

Table 1

The volume of water contained (million m³) in the rB1 reservoir according to specific water surface elevations (m) and related scenarios (s1–s4) developed in the study.

Water surface elevation (m)	172.14	177.86	178.14	178.64	179.14	179.21	179.64	179.8	180.14
Volume of water contained (million m3)	0	4.8	5.41	6.77	8.18	8.4	9.66	10.1	11.18
Related scenario		s1		s2		s3		s4

Table 2

Streamflow features and flood multi-scenario developed in this study.

2D Streamflow Scenario Code	No. of Open Gates (Sw1 1 and H-pp 2)	No. Hours/Open Gates (Sw1 1)	No. Hours/2D Simulation	Q 3 (m3/s) (Sw1 1 + H-pp 2)	Total Q 3 (m3/s)	Q 4 (%) from Sw1 1
Avg. discharge	Servitude discharge and H-ppavg	24 h	24 h	2.8 + 34	36.8	-
Scenario 1.–s1	1 and H-ppmax	1 h	24 h	490 + 34	524	25
Scenario 2.–s2	2 and H-ppmax	1 h	24 h	980 + 34	1014	50
Scenario 3.–s3	3 and H-ppmax	1 h	24 h	1470 + 34	1504	75
Scenario 4.–s4	4 and H-ppmax	1 h	24 h	1960 + 34	1994	100

¹ Sw1: Bacău 1 spillway gates; ² H-pp: Hydro-power plant; ³ Q: Discharge (m³/s); ⁴ Q: Percentage of total spillway gates discharge capacity.

Table 3

Flood extent area, built-up area, and number of buildings potentially affected by floods computed for each scenario (s1–s4).

Flood Extent	s1	s2	s3	s4
Total flood extent (km2)	6.27	8.66	11.46	12.82
Built-up area (km2)	0.001	0.04	1.09	1.5
Number of affected buildings	6	173	2461	3780

Table 4

Number of buildings potentially affected by floods vs. flood velocity (m/s) computed for each scenario (s1–s4).

Flood Velocity (m/s)	s1	s2	s3	s4
<1	6	170	2260	3238
1–2		3	184	353
2–3			17	189
3–4
4–5

Table 5

Number of buildings potentially affected by floods vs. flood depth (m) computed for each scenario (s1–s4).

Flood Depth (m)	s1	s2	s3	s4
<1.0	6	165	1174	1312
1.0–2.0		8	971	1171
2.0–3.0			312	980
3.0–4.0			4	299
4.0–5.0				18

Table 6

Flood hazard classification based on water depth according to the MLIT [42].

Flood Hazard	Flood Depth (m)	Hazard Classes	Hazard Description
H1	<0.5	Very low	Flood does not pose hazard to people and on-foot evacuation is not difficult.
H2	0.5–1	Low	Flood water poses hazard for infants and on-foot evacuation of adults becomes difficult; evacuation becomes more complicated.
H3	1–2	Medium	Flood depth can drown people; people may be safe inside their homes.
H4	2–5	High	People are exposed to flood hazard even inside their homes and evacuate towards the roof of their homes is suggested.
H5	>5	Extreme	Built-up structures like homes may get covered by the flood; people may get drowned even if they evacuate towards the roof of their homes.

Table 7

Number of buildings potentially affected by floods vs. flood hazard classes computed for each scenario (s1–s4).

Flood Hazard Classes 1	s1	s2	s3	s4
H1 (Very low)	6	150	674	694
H2 (Low)		15	500	618
H3 (Medium)		8	971	1171
H4 (High)			316	1279
H5 (Extreme)				18

¹ Flood hazard classes according to [42] (see Table 6).