Different scenarios of an arch-dam breach and their impact on the time-space evolution of flood waves are analysed using numerical modelling. As the accidents involving this type of dam are among the most catastrophic ones, the 108 m in height Paltinu arch-dam, Romania, was chosen as a case study due to its problems in the past. Three dam breach magnitudes and two inflow hydrographs for the worst-case scenario of Normal Operating Pool elevation in the reservoir were chosen as variable parameters, to assess their influence on the dam break wave characteristics and downstream flooded areas. The flood was routed along the 18 km reach of the Doftana River down to the confluence with the Prahova River. A 2D numerical model was set up with the help of HEC-RAS software, which was also used to analyse the resulting hazard maps under a GIS environment. Comparison of inundation boundary, maximum depths and velocities, as well as the arrival time at control sections allow for conclusions to be drawn. These predictive results of shape, magnitude, and time to peak of the flood waves are essential for flood risk management to obtain the risk maps, estimated damage costs, and possible affected areas.

  • The study investigates through 2D numerical modeling the consequences of a possible arch-dam failure in Romania.

  • The study proposes a new methodology to assess the influence of different dam breach scenarios on the time-space evolution of flood waves.

  • The proposed methodology could be applied to other similar flood studies.

  • The research shows the advantages of using the proposed numerical method.

Dams have been constructed worldwide for diverse purposes, such as water supply, flood protection, hydropower, irrigation, and recreational activities (International Rivers 2007). According to the International Commission on Large Dams (ICOLD), there are over 58,000 dams out of a total of 800,000 in the world, which are considered ‘large’ in the World Register of Dams (ICOLD 2022; Ferrari et al. 2023). To qualify as ‘large,’ a dam must have a height of at least 15 m, or a height between 10 and 15 m and an impounded volume of more than 3 million m3 (ICOLD 2023).

Dams can be classified based on their building material and shape. The traditional classifications include: (i) embankment and rock-filled dams, which are made of earth and/or rock and account for more than 78% of all dams worldwide, (ii) gravity dams, which are made of concrete, stone, or other masonry, and (iii) arch and buttress dams, which are made of concrete (Stematiu & Ionescu 1999; Deangeli et al. 2009; Tanchev 2014).

However, the history of dam construction has been marked by dam accidents, which can result in catastrophic consequences. Since the 12th century, around 2,000 dam failures have occurred (Singh 1996). During the 20th century, 10% of these failures resulted in the loss of at least 8,000 lives and caused material damage worth millions of dollars (Singh 1996; Luino et al. 2014).

At the failure of a dam, a massive amount of water is suddenly released through the breach, leading to flash flooding characterised by a high discharge of short duration. The resulting wave can sweep away buildings, bridges, roads, and railways, among other things (Singh 1996).

More than 2,000 dams have been built in Romania on the main rivers to meet various water needs, and about 241 of these are considered large by ICOLD (ICOLD 2022). Most of the dams in Romania were constructed between 1950 and 1990 and have already been operational for 30–70 years. Because many of them are located within 100 km of the major seismic area of Vrancea, they must be constantly monitored. Belci Dam on the Tazlău River in 1991 and Cornățel Dam on the Vedea River in 2003 were the only notable dam break accidents in Romania (Filotti 2015).

To minimise the loss of life and property damage caused by dam breaks it is crucial to use hydraulic numerical modelling to determine the downstream flooded areas, flow depth, wave velocity, flood intensity, and the peak time of propagation for a given accident scenario (Singh 1996; Nistoran Gogoașe et al. 2016). Over time, numerical models have evolved from 1D approximation to 3D, whereas graphical interfaces for data management and visualisation of results have been improved and integrated into a GIS environment. The two-dimensional approximation significantly decreases the computational complexity as compared to the 3D approach, data storage requirements, and computation time needed to study the behaviour of real flood waves. Various hydraulic modelling packages integrate the 2D Shallow Water Flow Equations, such as TUFLOW, Mike 21, Telemac 2D, Delft3D Flow, SOBEK, InfoWorks, HEC-RAS, LISTFLOOD-FP, DIVAST, and JFLOW. Of these, the free HEC-RAS 2D has been chosen to be used for the present dam break study, since this software has been intensively developed in the last few years by the USACE (2022), it is widely used and it has a very good graphical interface with visualisation tools for displaying the hydraulic variables in a GIS environment. As a result of numerical simulations, flood mapping helps visualise the spatial and time distribution of flood hazards and risks for a given flooding event (of either natural or anthropic cause). This is an indispensable tool for flood risk management, helping to plan emergency actions, assess potential losses, reduce existing risks and prevent the occurrence of new ones, at the same time preparing the population for such catastrophic events (Tsai et al. 2019; Maranzoni et al. 2022).

Engineers and scientists face the challenge of making accurate estimates of dam-breach peak discharges from limited data sets. Therefore, statistical approaches, such as probabilistic methodologies, have been proposed to address this issue (Moglen et al. 2019; Bello et al. 2022). These allow for a more comprehensive understanding of the reservoir-flooded area system behaviour and the associated risks. However, they are very complex and require expertise to develop and interpret. Even though deterministic models cannot capture the inherent uncertainties and variability of natural systems, they rely on fixed equations to predict the behaviour of the system being studied and provide a simplified and straightforward approach to dam breach modelling.

The primary reasons for the failure of dams are diverse. These include design errors relating to flooding events, such as overtopping (which accounts for 34% of cases) (Adamo et al.2020b), issues with foundation rock stability and bank failures (30% of cases), and internal or underlying dam seepage (20%). Other reasons include malfunction or failure of spillway structures (10%), landslides (Adamo et al. 2020a, 2020b), structural failures, upstream dam failure in case of cascades, operation or maintenance errors, earthquakes, and sabotage. Failure modes can occur through overtopping, piping or seepage, cracking, sliding, and overturning (Fiedler 2016).

The breach formation time depends on the type of dam material and mode of failure. In embankment dams, the rupture process is relatively slow and can take from several hours to days, while for concrete dams, the failure can occur within a matter of minutes. Of all the aspects of a dam break, the estimation of the characteristics and parameters of the breach is typically associated with the highest degree of uncertainty (Ahmadi & Yamamoto 2021).

Among all types of dams, arch dams have the shortest breach progression time, which can result in significant human casualties and material damage (USBR 2019). The most well-known case of such an occurrence is the Malpasset arch-dam in France, which experienced a failure in 1959 due to geological issues in the left abutment. The dam collapsed by overturning, leading to the loss of 421 lives and an estimated 68 million US dollars in damages.

Accurate estimations of the location, shape, and dimensions of a dam breach and the time it takes to develop are crucial for flood risk assessment studies. The type and construction material of the dam primarily influence these parameters. The water discharge through or over the breach is calculated using equations similar to those for flow through an orifice or over a weir, with the head and opening area increasing in time according to a specified mathematical relationship.

As a result of the potential dangers associated with dam failures, scientists have developed physical models to analyse the failure mechanisms and breach shapes of various types of dams, including arch dams (Morris et al. 2007). In the case of narrow canyon arch dams that are exposed to earthquakes, physical model tests have shown the development of diagonal cracks that run parallel to the side slope abutments and horizontal cracks, leading to trapezoidal-shaped breaches (USBR 2002).

When estimating overtopping, the peak outflow magnitude is directly affected by the estimated dimensional discharge coefficient, C, of the breach (associated with a weir flow). Even though the exact value of (where is the non-dimensional weir discharge coefficient) cannot be determined, according to USACE (2014), for arch dams being overflown, the recommended values of C should be in the range of 1.1–1.8 m0.5/s. Additionally, the average breach width of arch dams is typically in the range (0.8–1) L (with L the top length of the dam at the deck), its slope being equal to or lower than the valley wall slope, and a failure time of 10 min or less. Many researchers have developed peak outflow equations based on data from historic dam failures as a function of breach/weir head, water volume in the reservoir, the height of the dam, and other factors (USACE 2014).

The present study examines the computed results of a possible failure at the Paltinu arch-dam located on the Doftana River in Romania, by exploring several scenarios. The primary objective of the study is to analyse the impact of various dam breach assumptions on the resulting flood wave and its routing downstream. Computations are performed with the HEC-RAS software version 6.3.1 (USACE 2022) which incorporates a 2D simplifying approach. Its practicality was demonstrated in previous works (Kumar et al. 2017; Albu et al. 2020; Gogoașe Nistoran et al. 2023). Geometric and hydrologic data are processed and results are computed under GIS, which makes it easier to process the flood maps (Albano et al. 2019).

The research is a 2D modelling investigation for a possible arch-dam failure in Romania, to assess the coupled influence of different dam breach scenarios and hydrologic events on the time-space evolution of flood waves in the inundated areas. The paper is organised as follows: Section 2 presents the site, data, equations and the numerical method used by the HEC-RAS software, employed to reproduce the dam break accident with multiple scenarios. In Section 3, simulation results are analysed comparatively, in order to formulate the Conclusions in Section 4.

Site and data

A case study for the analysis is the Doftana River, which is a 51 km-long tributary of the Prahova River, flowing in the Carpathians of the Curvature region of central Romania. Doftana River watershed has a 419 km2 area, with a difference in elevation of 1,543 m and a mean annual flow rate of 4.8 m3/s at its confluence (Armaș 1999) (Figure 1). The main gauging station along the Doftana River is at Teșila, where the mean annual flow is 4.59 m3/s (measured between 1950 and 1998) (Armaș 1999) and the 100-year flood is 435 m3/s (FRMP Buzău-Ialomița 2022).
Figure 1

(a) Location of the Doftana watershed, Paltinu reservoir, and Teșila Gauging Station. (b) Reservoir synthetic inflow hydrographs for the1,000-year and 100-year floods.

Figure 1

(a) Location of the Doftana watershed, Paltinu reservoir, and Teșila Gauging Station. (b) Reservoir synthetic inflow hydrographs for the1,000-year and 100-year floods.

Close modal
Along the Doftana River and some of its 17 tributaries, 13 micro- and one hydropower plants were developed (Figure 2). Of these, seven are located upstream and six downstream of the main reservoir, Paltinu, which is embanked in a narrow valley by a concrete arch-dam.
Figure 2

Location of the hydropower plants along the Doftana River. The most important is Paltinu HPP downstream the largest reservoir with the same name.

Figure 2

Location of the hydropower plants along the Doftana River. The most important is Paltinu HPP downstream the largest reservoir with the same name.

Close modal
The reservoir has multiple functions, such as flood control, supply with drinking water for Câmpina-Ploiești-Brazi towns, hydropower, irrigation and leisure activities. Its characteristic features and volume-elevation curve are given in Figure 3 (Asman & Boboc 2021). The drained watershed area in the cross-section of the dam is 334 km2.
Figure 3

(a) Volume-elevation curve and (b) characteristic features of the Paltinu reservoir.

Figure 3

(a) Volume-elevation curve and (b) characteristic features of the Paltinu reservoir.

Close modal
Paltinu is a double-curvature arch-dam, with a perimeter joint (Figure 4(a)), which was built between 1966 and 1971, at a distance of about 18 km upstream of the confluence of Doftana River with Prahova. This constructive solution was chosen because of the geo-morphological conditions (faulted and brittle rock on the left wall slope, the bare rock unit consisting of sandstone with interlayered marl shales and clays, specific to a landslide-prone area). The characteristics of the dam are given in Figure 3(b)). The dam is equipped with two types of appurtenant structures for flow release: a morning glory spillway (Figure 4(b)) and two bottom outlets positioned by design at the elevations of 574 and 625 m.a.s.l., each discharging 49 and 180 m3/s at Normal Operating Pool level (NOP). The NOP elevation (or Full Reservoir Level – FRL) is the maximum level to which water can be stored in the reservoir during ordinary operating conditions, without spillway discharge or through sluiceways (US Society on Dams 2023). The dam is classified in the special importance category, which requires continuous surveillance.
Figure 4

(a) Paltinu arch-dam and consolidation works and (b) the morning glory spillway.

Figure 4

(a) Paltinu arch-dam and consolidation works and (b) the morning glory spillway.

Close modal

During the first filling of the Paltinu reservoir, in 1974, when the water elevation reached 642.85 m.a.s.l., an incident occurred that required the reservoir to be emptied: seepage and abnormal displacements were noticed at the left abutment. Subsequent works followed such as abutment consolidations, increased evacuation discharge through the spillways, improved joint sealings, and an improved dam safety surveillance monitoring system. Some of these works took place with the reservoir partially full. The consolidation works ended in 1986 and a second reservoir filling was performed up to the NOP level in 1995 (Asman & Boboc 2021).

All these additional consolidation works were also necessary because Paltinu dam is located near the first two most important Romanian earthquake epicentres in seismic activity: Vrancea (at 86 km NE), with a 50 year return period of earthquakes of 8-degree intensity (on Medvedev-Sponhouer-Karnik-64 Scale) and Câmpulung – Făgăraș (at 10–15 km W), with an intensity of 9 on the same scale and a return period of 100 years (NWARW 2012).

The design of the dam in the selected cross-section was carried out for a 1,000-year flood with a peak flow of 718 m3/s – and was verified for a 10,000-year flood with a peak flow of 920 m3/s, under natural regime flow (Table 1). The installed discharge of the hydropower is 15 m3/s, producing through the two Francis turbines an installed hydropower of 10.2 MW and a mean of 30.2 GWh/year.

Table 1

Characteristic features of the Paltinu dam

Crt no.Characteristic featureValue (meas. unit)
Dam height 108 m 
Length at deck 460 m.a.s.l. 
Deck elevation 652 m.a.s.l. 
Width at deck 6 m 
Width at foundation 22 m 
Surface area at NOP 182 ha 
Design peak discharge 0.1% 718 m3/s 
Verification peak discharge 0.01% 920 m3/s 
Terrain elevation downstream the dam 555 m.a.s.l. 
Crt no.Characteristic featureValue (meas. unit)
Dam height 108 m 
Length at deck 460 m.a.s.l. 
Deck elevation 652 m.a.s.l. 
Width at deck 6 m 
Width at foundation 22 m 
Surface area at NOP 182 ha 
Design peak discharge 0.1% 718 m3/s 
Verification peak discharge 0.01% 920 m3/s 
Terrain elevation downstream the dam 555 m.a.s.l. 

Mathematical model and numerical scheme

With the aim of preventing loss of life and damage to property caused by a dam break, it is crucial to accurately and quickly estimate the flood wave routing. However, predicting such events can be difficult due to their complex nature and their evolution in both time and space. Therefore, a mathematical model is necessary to study how these flows progress downstream.

The 2D Shallow Water Equations (SWEs) were proposed by Barré de Saint-Venant (1871) to model channel flows and express the physical conservation principles of mass and momentum. They are commonly used to simulate the free surface dynamics of water in hydraulic and environmental studies such as flow in channels, rivers, estuaries and coastal areas, flood – forecasting, and marine waves, under the assumptions that the height scale of the flow is much smaller than its length scale and the pressure is hydrostatic. This nonlinear system of partial differential equations is based on the conservation of mass and linear momentum, and describes how the depth or the elevation of the water surface and the average velocities in both directions x and y change, over time.

The complete set of 2D SWE can be further simplified for small-scale analyses, by neglecting the Coriolis and wind forces. These approximations provide the following system of continuity and full-momentum equations in the differential form (USACE 2022):
formula
(1)
In Equation (1), t is the time, is the water depth, is the bottom elevation (considered fixed, with no erosion/deposition), is the water surface elevation relative to the datum, u and v are the depth-averaged velocity components in the x and y Cartesian directions respectively, q is a source/sink flux term, is the horizontal eddy viscosity coefficient, and are the horizontal bed shear stresses, given as a function of the drag coefficient, and the modulus of the velocity vector, . The horizontal bed shear stresses can be written as:
formula
(2)
The drag coefficient may be expressed in terms of the Manning coefficient [L−1/3T], n, and hydraulic radius, R, through
formula
(3)
or in terms of the nonlinear bottom friction coefficient, , hydraulic radius, Manning coefficient and velocity magnitude, through
formula
(4)

Further simplifications of the full-momentum shallow water Equation (1) can be made by neglecting the inertial and diffusion terms. Thus the reduced diffusive and kinematic wave representations are obtained, which cannot reproduce the dynamics of dam break flows.

Over the last few decades, there has been a significant interest in investigating the efficiency and precision of numerical techniques to solve the SWE (Toro & Garcia-Navarro 2007). Different numerical methods, such as finite difference, finite element or finite volume (FV), based on diverse numerical grids, such as Cartesian or boundary fitted, either structured or unstructured, can be used to obtain a solution to these equations. All of these have advantages and disadvantages in the context of inundation modelling, but not all of them are equally effective in capturing the complex dynamics of the system. Routing a dam break flood wave downstream is one of the most difficult unsteady flow problems to solve (USACE 2022), particularly because of the irregular topography and hydraulic properties, wet and dry cells, initial and variable boundary conditions, time step, all leading to computational instabilities. Notable results were obtained with 2D SWE models based on approximate Riemann solvers, such as those developed by Roe (Roe 1981), Harten (Harten et al. 1983) and Toro (Toro 2001, 2009) and high-resolution Godunov-type methods (Harten 1997), Anastasiou & Chan 1997). Accurate simulations were performed for transcritical flows, shock waves and moving wet–dry fronts (Ginting & Mundani 2019). Most of the schemes have to comply with the Courant–Frederichs–Lewy (CFL) < 1 of the time step for stability. The introduction of the Large Time Step (LTS) approach by LeVeque (LeVeque 1998) led to the emergence of remarkably efficient computational solvers for the SWEs (Xing 2017; Xu et al. 2022).

HEC-RAS software has implemented two methods to solve the SWEs, that differ in how they treat the acceleration terms: the Eulerian Method (SWE-EM) and the Eulerian-Lagrangian Method (SWE-ELM). The former is easier to compute in the space and time-fixed frame of reference, whereas the latter is more expensive to compute in its movable frame of reference along a flow path, and it allows for larger time steps. Both solvers use a combination of finite difference and FV methods on an unstructured polygonal mesh with subgrid bathymetry. For the current dam break case study, the SWE-ELM method was used to compute subcritical, supercritical, and transcritical flow regimes over irregular topography (USACE 2022).

FV schemes are a class of numerical methods that are particularly well-suited for solving the SWEs (LeVeque 2004). These schemes have several advantages, including their ability to handle nonlinear terms and discontinuities, and their ability to operate with irregular geometries and boundary conditions (Kurganov 2018). In particular, FV schemes can accurately capture transcritical flows and shock-type discontinuities, which are sudden changes in flow properties that occur over very short distances or time (Causon et al. 1999).

The HEC-RAS solver for the SWE is based on an implicit FV approach. This solution method can use larger time steps compared to explicit methods, allowing an enhanced level of stability. A notable feature of the technique is the robust management of wetting and drying within 2D cells, from entirely dry to a sudden inflow of water.

To improve the efficiency of numerical simulations that solve the SWEs, non-uniform meshes can be used. The mesh representation can be unstructured (irregular connectivity) or structured (regular connectivity), raster-based, flexible and allow for different levels of resolution across the domain, for which the topography is usually acquired with LiDAR technology. The approach is particularly useful for modelling complex terrain features, where high-resolution grids can be used, while coarser grids are sufficient in other areas, reducing the overall computational cost and runtime (Ferrari et al. 2023).

Regarding the HEC-RAS mesh structure, the software provides flexibility in utilising either unstructured or structured computational meshes. While its primary design centres around unstructured meshes, it's also capable of handling structured meshes. The software assumes the orthogonality of the cells, thus enabling better computational efficiency. As a result, the computational cells can adopt various shapes and sizes, including triangles, squares, rectangles, and even polygons with a maximum of eight sides per element. The outer boundary of the computational mesh is established using a polygon.

A particular feature of each computational cell and its corresponding face in HEC-RAS is that it relies on the characteristics of the underlying terrain model, known as ‘high-resolution subgrid model’ (Casulli 2008). The software uses a dedicated pre-processor to thoroughly establish a relationship between water surface elevation and detailed geometric and hydraulic property tables of the cells and cell faces with the purpose to reduce the overall computation time. This relationship is established using the terrain data at the grid-cell resolution within each individual cell. Thus, a cell can accommodate the precise water volume and be just partially wetted in accordance with the provided water surface elevation and based on the grid-cell resolution. Each computational cell face is treated as a cross-section, and subsequently subjected to the pre-processing stage, resulting in the creation of comprehensive hydraulic property tables that detail the connection between elevation and various attributes such as wetted perimeter, area, roughness, and others. The flow across the interfaces of the cells is reliant on this comprehensive data. Therefore, larger computational cells may be used while substantially retaining most of the complexity of the terrain that dictates the flow dynamics, with the advantage of significantly accelerated run times.

Numerical model setup

Two main types of data are required for the flood modelling study: a topo-bathymetric survey or digital terrain model to describe the geometry of the terrain, and hydrological data to specify the boundary conditions. In recent years, flood modelling has been performed using GIS software, which allows for the integration of multiple results layers on different maps or on digital terrain models.

To represent the spatial data, a freely available digital elevation model with a spatial resolution of 25 m was downloaded from the Copernicus Land Monitoring service (EU-DEM 2016). The river bathymetry was not used, as the main channel is very narrow and the mean annual river discharge is very low (4.5 m3/s) in comparison with the dam break peak flood (tens of thousands of m3/s). The reservoir was modelled as a storage area characterised by an elevation-capacity curve.

The 2D computation area was delineated along the Doftana River downstream to its confluence with the Prahova River, largely enough to accommodate the flood wave (Figure 5(b)). Two average values of the 2D Manning roughness coefficient were delimited on the map: one of 0.07 s/m1/3 for the vegetated (agro-forestry) areas, and another one of 0.11 s/m1/3 for the inhabited (medium intensity developed) areas (Figure 5(b)). These values were chosen according to the specifications given by Papaioannou et al. (2018), based on CORINE land cover data for this area (Copernicus 2018).
Figure 5

(a) Dam breach shape for the considered three cases; (b) computation domain with reservoir storage area, dam, 2D mesh and different roughness areas for the inhabited zones.

Figure 5

(a) Dam breach shape for the considered three cases; (b) computation domain with reservoir storage area, dam, 2D mesh and different roughness areas for the inhabited zones.

Close modal

Water initially present in the Doftana River's main channel or in other small reservoirs of the downstream hydropower plants was neglected. The reservoir's tail was defined as the boundary condition at the upstream end of the system, while the confluence of the Prahova River was defined as the boundary condition at the downstream end. Additionally, the dam was incorporated as a connection between the storage area of the reservoir and the downstream 2D computation zone, based on its constructive features. Since the interest of the study is related to the inundated areas downstream of the dam, the simplified level pool routing of the flood wave is used for the reservoir (Tsai et al. 2019), instead of a dynamic wave routing. Therefore, in the modelled Paltinu reservoir the water surface remains horizontal at all time steps during the computations (Ionescu & Gogoașe 2019).

The model's boundary conditions include the upstream and downstream conditions embodying the inflow hydrographs at the tail of the Paltinu reservoir (as shown in Figure 1(b)) and a uniform flow applied along a cross-section at the downstream confluence, respectively. This downstream condition is determined by the terrain slope of around 1% characteristic of this sub-Carpathian region. The inflow hydrographs considered are those of the design and check floods, having the characteristics outlined in Table 2. Other inflows from tributaries of the Doftana River were also disregarded since their discharges are insignificant compared to those of the flood waves. The initial conditions for the model are the water elevation in the reservoir at the designed NOP Level and water elevation in the river valley below the ground level.

Table 2

Characteristics of the Paltinu reservoir inflow hydrographs

CaseQT–year (m3/s)Return periodTpeak – time to peak (h)T total duration (h)
718 1,000-year 10 48 
435 100-year 10 48 
CaseQT–year (m3/s)Return periodTpeak – time to peak (h)T total duration (h)
718 1,000-year 10 48 
435 100-year 10 48 

In this particular study, the mesh generated by the HEC-RAS software to model the potentially flooded area was composed of 78,612 cells. Most of these cells were square-shaped covering surfaces of 25 m × 25 m. To ensure numerical stability and prevent instabilities that can arise in rapidly developing phenomena, relatively short computation times between 1 and 0.1 s were chosen, depending on simulation cases.

Typically, numerical models are calibrated using experimental data. However, in the context of dam break simulations, such data can only be obtained either in a laboratory or from real accidents after they occur. Consequently, software developers often use benchmarking tests, which involve comparing the performance of the numerical model with established standards or known results from previous studies, as studies for scenarios of hypothetical accidents cannot be calibrated. The HEC-RAS software covers verification, validation, and benchmarking tests for 2D simulations (USACE 2018a, 2018b), featuring various examples of flooding, including the Malpasset dam break case.

Scenarios and modelling process

In line with the findings obtained from physical modelling studies (Morris et al. 2007; Aureli et al. 2021), the failure of the Paltinu dam was assumed to be triggered by the water level overtopping the dam deck elevation (652 m.a.s.l.) for all investigated scenarios. The breach is considered to develop into its maximum area in approximately 6 min, assuming a trapezoidal shape with lateral slopes matching the abutment slopes: a right slope of 1:5 and a left slope of 1:2 (see Figure 5(a)). To model the dam breach, the case studies presented in Table 3 were selected, according to the guidelines provided by USACE (2014) and NRDWR (2010). A synthesis of the modelling process is presented in the flow chart from Figure 6.
Table 3

Breach characteristics for the case studies

CasesPeak inflow at reservoir tail (m3/s)Initial NOP level in reservoir (m.a.s.l.)Breach width at base (m)Final breach elevation (m.a.s.l.)Breach formation time (min)Breach area (m2)
Q100-year = 435 649 70 627 3,925 
Q100-year = 435 649 140 623 6,990 
Q100-year = 435 649 160 621 8,310 
Q1,000-year = 718 649 70 627 3,925 
Q1,000-year = 718 649 140 623 6,990 
Q1,000-year = 718 649 160 621 8,310 
CasesPeak inflow at reservoir tail (m3/s)Initial NOP level in reservoir (m.a.s.l.)Breach width at base (m)Final breach elevation (m.a.s.l.)Breach formation time (min)Breach area (m2)
Q100-year = 435 649 70 627 3,925 
Q100-year = 435 649 140 623 6,990 
Q100-year = 435 649 160 621 8,310 
Q1,000-year = 718 649 70 627 3,925 
Q1,000-year = 718 649 140 623 6,990 
Q1,000-year = 718 649 160 621 8,310 
Figure 6

The flow chart of the modelling process.

Figure 6

The flow chart of the modelling process.

Close modal

Different decreasing time steps were used in order to avoid computation instabilities and errors. Computed hydraulic parameters were saved in the output database file every 10 min.

Table 4 shows the most important computed features of the flood waves downstream of the dam and at the confluence of Doftana with the Prahova River. Flow attenuation was calculated as the ratio of the peak discharges at the dam breach and at the downstream river confluence, respectively.

Table 4

Comparative results of the flood wave characteristics for the considered simulation cases

Characteristics CasesArrival time from dam break at the downstream confluence (min)Peak Q and time at the dam/at the confluence, (m3/s), (hh:mm)Maximum width of the inundation boundary near confluence (m)Peak flow attenuation (%)Max. velocity V DS the dam (m/s)Max. flood depth, h DS the dam (m)Max. flood intensity v·h DS of the dam (m2/s)
C1 60 18,896 (9:50)/7,820 (10:50) 410 41 21 14 298 
C2 50 31,844 (9:50)/12,073 (10:40) 470 38 23 18 416 
C3 40 50,129 (9:50)/15,781 (10:30) 520 31 25.6 26.2 540 
C4 60 19,947 (8:20)/7,713 (9:20) 440 39 22.1 14.4 318 
C5 50 33,410 (8:20)/11,735 (9:10) 533 35 23.4 18.6 438 
C6 40 55,584 (8:20)// 16,872 (9:00) 615 29 25.7 26.7 550 
Characteristics CasesArrival time from dam break at the downstream confluence (min)Peak Q and time at the dam/at the confluence, (m3/s), (hh:mm)Maximum width of the inundation boundary near confluence (m)Peak flow attenuation (%)Max. velocity V DS the dam (m/s)Max. flood depth, h DS the dam (m)Max. flood intensity v·h DS of the dam (m2/s)
C1 60 18,896 (9:50)/7,820 (10:50) 410 41 21 14 298 
C2 50 31,844 (9:50)/12,073 (10:40) 470 38 23 18 416 
C3 40 50,129 (9:50)/15,781 (10:30) 520 31 25.6 26.2 540 
C4 60 19,947 (8:20)/7,713 (9:20) 440 39 22.1 14.4 318 
C5 50 33,410 (8:20)/11,735 (9:10) 533 35 23.4 18.6 438 
C6 40 55,584 (8:20)// 16,872 (9:00) 615 29 25.7 26.7 550 

Figure 7 depicts the breach discharge hydrographs at the dam and the attenuated ones at the downstream confluence of Doftana with the Prahova River. The breach hydrographs have almost identical shapes, with slightly higher peaks for the C4–C6 cases than for the corresponding C1–C3 cases. As seen from the computed results in Table 4, the largest peak flow of 55,584 m3/s belongs to case C6, corresponding to the 1,000-year flood and the largest breach, whereas the lowest one of 15,781 m3/s belongs to case C3, corresponding to the 100-year flood and the smallest breach. A similar trend is observed for the downstream hydrographs. The corresponding peak flow attenuation and associated travel times of the peak flow between the dam and the downstream confluence are practically the same (for cases C1 and C4, C2 and C5, C3 and C6): 60, 50 and 40 min, for all three increasing breach cases, irrespective of the magnitude of the inflow hydrograph at the reservoir's tail (Figure 7 and Table 4). Therefore, the reservoir inflow hydrograph magnitude seems not to have an important influence on the travel time and flow attenuation of the dam break flood wave.
Figure 7

Flood routing between the dam and the downstream confluence with the Prahova River, for simulation cases: (a) C1–C3 (Q100-year) and (b) C4–C6 (Q1,000–year).

Figure 7

Flood routing between the dam and the downstream confluence with the Prahova River, for simulation cases: (a) C1–C3 (Q100-year) and (b) C4–C6 (Q1,000–year).

Close modal

However, the shape of the downstream flow hydrograph is different in the two case sets: the time to peak for the downstream flood waves is larger for the highest flood hydrograph, which means the flood waves last longer in cases C4–C6 than in cases C1–C3, expecting to produce more damage in this extra time. If multiple control sections are drawn along the river valley at different distances from the dam, in order to investigate the flood wave hydrographs characteristics, the travel time of the wavefront/peak can be obtained. These results can be very useful for predicting available evacuation times for the residents of the area, and the operation of sirens located on the valley slopes.

Since the dam breach is triggered by the water stage reaching the elevation of the deck elevation (652 m.a.s.l.) in the reservoir, the flood wave starts 90 min earlier for the 1,000-year flood (cases C4–C6) than for the 100-year one (cases C1–C3) (Table 4). This happens because it needs more time for the reservoir to fill up between the NOP level and the dam deck level at a smaller incoming discharge. The reservoir stage hydrographs in Figure 8 show a sharp decrease between the dam deck elevation towards the breach crest elevation after the dam failure 90 min later for the 100-year flood This means that the higher the flood, the higher the risk of dam breach due to overtopping. Even if no scenarios were considered for different water surface elevations in the reservoir, it seems obvious that the risk of dam overtopping increases when the initial reservoir level is higher. Also, one may see that the shape of the reservoir emptying curve in Figure 8 depends only on the dam breach magnitude, not on the incoming hydrograph.
Figure 8

Time variation of the reservoir water surface elevation for the six simulation cases.

Figure 8

Time variation of the reservoir water surface elevation for the six simulation cases.

Close modal

Through the functionalities of its integrated GIS interface, RAS Mapper, HEC-RAS offers the capability to generate maps and rasters (that can be exported to any GIS software), of the flood hazard indicators such as flow depth (h), flow velocity (), flood intensity (), inundated area, flooding duration, flooding arrival time, return period/scenario probability, released water volume, etc. Flood hazards can be quantified or classified into classes using qualitative attributes such as low, moderate, medium, high, and extreme, according to different ranges of the values of these indicators (Marangoni et al. 2022).

Figures 9 and 10 present comparatively the maps of maximum depth and velocity for all simulation cases C1 ÷ C6. They represent the maximum values in a computation cell, over the entire simulation period. These map layers, which can be overlapped onto the road and building maps, are indispensable in the development of flood risk management plans. From the colour intensity/legend one may see that the water depth and velocities along the flooded Doftana Valley are the largest in the worst-case scenario, C6. Also, velocity increases in the narrow areas of the river valley and is larger in the upstream part, where the discharges and the slope of the valley are greater.
Figure 9

Maximum depths in the flooded area for computation cases C1–C6.

Figure 9

Maximum depths in the flooded area for computation cases C1–C6.

Close modal
Figure 10

Maximum velocities in the flooded area for computation cases C1–C6.

Figure 10

Maximum velocities in the flooded area for computation cases C1–C6.

Close modal

Several rural and urban communities exist along the Doftana River Valley, downstream of the Paltinu Dam, such as Seciuri, Lunca Mare, Plaiu Câmpinei, Brebu and Câmpina town, together with a number of plants and factories. Therefore, the dam break flood wave may cause important social and economic damage with an increased level of risk, due to the population growth over time in this area.

Figure 11 shows two inundation boundaries (maximum waterfront limit): for the worst-case (C6) and the best-case (C1) scenarios in two different areas along the river. As the Doftana River Valley is quite narrow, the horizontal difference between the two boundary lines is less than 100 meters, yet enough to flood a few more houses. These boundaries are important for the insurance companies, not only for the development of flood management plans.
Figure 11

Inundation boundaries for simulation cases C1 (red) and C6 (yellow): (a) detail of an industrial plant and (b) a bridge, near Câmpina.

Figure 11

Inundation boundaries for simulation cases C1 (red) and C6 (yellow): (a) detail of an industrial plant and (b) a bridge, near Câmpina.

Close modal
Figure 12 depicts the computed velocity vector magnitudes superimposed over the water depth layer and the satellite image. Knowing the direction the flow comes from helps in the design of flood protection barriers and defence systems.
Figure 12

Velocity vectors over the water depth layer, showing direction of flow in the flooded areas.

Figure 12

Velocity vectors over the water depth layer, showing direction of flow in the flooded areas.

Close modal
Figure 13 shows the map of maximum flood intensity, defined as the product between the depth and velocity magnitude (in m2/s). One may see this flood parameter has much higher values (up to 200 m2/s) downstream, close to the dam, where the flood wave peak flows are the largest and the river valley is the narrowest. Since the travel time of the flood wave is also the shortest, these are the most important areas to protect. Even though the risk of a dam break is very low, people and insurance agencies should be aware of the potential risk per se when building a house in these areas.
Figure 13

Maximum flood intensity, defined as the product of flow depth and velocity (m2/s).

Figure 13

Maximum flood intensity, defined as the product of flow depth and velocity (m2/s).

Close modal

The findings of this study provide a basis for the future development of a methodology to estimate the spatial and temporal evolution of flood waves generated by the breach of an arch-dam, aiming at supporting the improvement of flood risk management.

Dam failures produce devastating floods with catastrophic consequences upon the downstream areas, social and economic losses. Computer simulations provide a useful instrument for predicting the flood hazard, aiming at mitigating the risk of such events on the terrain and inhabited areas. Numerical results are useful for the development of emergency planning, flood risk management strategies and protection measures to improve the resilience of the communities. However, flooding dynamics depend on the chosen scenario and initial/boundary data, from the magnitude of the inflow, and initial level/volume in the reservoir, to dam failure mechanism, breach characteristics, numerical method, mesh features, roughness parameters, etc. In comparison with other constructive types of dams, such as embankment, gravity or buttress dams, the quickest ruptures triggering the highest intensity events occur for the arch dams.

This is why the proposed paper analyses the impact of different breach scenarios and hydrologic events upon the main characteristics of the disastrous flood waves released downstream an arch-dam. As a case study the Paltinu dam in Romania was chosen, since it had some safety problems in the past, at the joint between the abutments and the rock.

In this study, 2D numerical simulations were performed with the HEC-RAS software, version 6.3.1. The research covered 18.4 km of the inhabited Doftana River Valley, from the Paltinu reservoir down to its junction with the Prahova River, near Campina. The simulations focused on three different breach scenarios for two extreme hydrologic events, the 100-year and 1,000-year floods. These events were used as input hydrographs at the reservoir tail, and the simplified level pool was used to route the flood wave through the reservoir. The analysis compared several flood hazard indicators of the computed dam break waves and inundation maps, such as the water depth, velocity, and flood intensity, peak discharge, arrival time, time to peak, in the control sections and the reservoir water surface variation in time.

According to the values presented in Table 4 and from the plots shown in Figures 7 and 8 it can be inferred that the reservoir inflow hydrograph has a minimal impact on the maximum flood wave discharge resulting from the dam break, as long as the reservoir water stage at the time of failure and the characteristics of the breach remain unchanged. This effect can be explained by a small ratio of the inflow to the reservoir volume (Ferrari et al. 2023), as in the case of Paltinu Dam. Also, the reservoir inflow hydrograph magnitude does not have an important influence on the travel time and flow attenuation of the dam break flood wave. In all three corresponding scenarios of the flood event, the dam break occurs 90 min earlier compared to the flood event. This happens because of the time required for the reservoir to fill up from Normal Operating Pool to deck level, considered as the breach-triggering event. Therefore, in cases when the reservoir water level is the dam breach-triggering event, the fuller the reservoir, the sooner the dam breaks and the quicker the flood wave travels downstream. So the emergency actions have to be prompter and the people have to be evacuated faster.

As anticipated, the breach area has the most significant influence on the maximum discharge at the dam. The shape of the reservoir emptying curve in Figure 8 depends mainly only on the dam breach magnitude. As a result, for each increase in breach area considered, the arrival time of the dam break wave at the downstream confluence with the Prahova River decreases by 10 min. Consequently, the dam breach size is a very important parameter that should be included in defining dam failure scenarios for such simulations.

One may see from Figure 7 that the shapes of the breach and routed hydrographs are very similar and have nearly identical peak discharge values for the corresponding breach scenarios in both the 100-year and 1,000-year flood events. This similarity is due to the matching shape and magnitude of the dam break hydrographs that enter the downstream 2D computational area.

The values of water depth, velocity and flood intensity in the inundated region, computed in the present study are extremely high (Figures 9, 10, 12 and 13), with most of the downstream areas experiencing total destruction. Romanian regulations (NRDWR 2010) specify simultaneous limits for total destruction the inundation intensity of 7 m2/s and velocity of 1 m/s.

By overlaying the computed inundation boundaries on maps or satellite images of inhabited areas as shown in Figure 11, vulnerability and risk maps can be developed for flood risk management and prevention of life loss and property damage. For more accurate results, it is recommended to use similar studies with Lidar topographic data and Corine land cover information (Albu et al. 2020; Urzică et al. 2021).

The results obtained in the present study allow the definition of useful scenarios to enhance and facilitate the effective management of flood risk. This could lead to more efficient, adapted and successful strategies to minimise the impacts of flooding on communities, infrastructure, and the environment.

All relevant data are included in the paper or its Supplementary Information.

The authors declare there is no conflict.

Adamo
N.
,
Al-Ansari
N.
,
Sissakian
V.
,
Laue
J.
&
Knutsson
S.
2020a
Dam safety: General considerations
.
Journal of Earth Sciences and Geotechnical Engineering
10
(
6
),
1
20
.
Adamo
N.
,
Al-Ansari
N.
,
Sissakian
V.
,
Laue
J.
&
Knutsson
S.
2020b
Dam safety and overtopping
.
Journal of Earth Sciences and Geotechnical Engineering
10
(
6
),
41
78
.
Albano
R.
,
Mancusi
L.
,
Adamowski
J.
,
Cantisani
A.
&
Sole
A.
2019
A GIS tool for mapping dam break flood hazards in Italy
.
International Journal of Geoinformation
8
(
250
),
1
23
.
Anastasiou
K.
&
Chan
C. T.
1997
Solution of the 2D shallow water equations using the finite volume method on unstructured triangular meshes
.
International Journal for Numerical Methods in Fluids
24
(
11
),
1225
1245
.
Armaș
I.
1999
Doftana River Watershed (Bazinul Hidrografic Doftana)
.
Editura Enciclopedică
.
Asman
I.
&
Boboc
S.
2021
Paltinu dam – surveillance during time (in Romanian)
.
Hidrotehnica Journal
64
,
32
44
.
Barré de Saint Venant
A. J.-C.
1871
Théorie du mouvement non-permanent des eaux, avec application aux crues des riviéres et á l'introduction des marées dans leur lit
.
Comptes Rendus de L'Acad’Emie des Sciences
73
,
147
154
.
Casulli
V.
2008
A high-resolution wetting and drying algorithm for free-surface hydrodynamics
.
International Journal for Numerical Methods in Fluids
60
(
4
),
391
408
.
https://doi.org/10.1002/fld.1896
.
Causon
D. M.
,
Mingham
C. G.
&
Ingram
D. M.
1999
Advances in calculation methods for supercritical flow in spillway channels
.
Journal of Hydraulic Engineering
125
(
10
),
1039
1050
.
Copernicus Corine Land Cover (CLC)
2018
Deangeli
C.
,
Giani
G. P.
,
Chiaia
B.
&
Fantilli
A. P.
2009
Dam failures, Chapter 1
.
WIT Transactions on State of the Art in Science and Engineering
6
.
Available from: https://www.witpress.com/Secure/elibrary/papers/9781845641429/9781845641429001FU1.pdf?ref=patrick.net.eeeguide.com. Components of Hydroelectric Power Plant, General terminology, Full Reservoir level (F.R.L.). https://www.eeeguide.com/components-of-hydroelectric-power-plant/ (accessed 18 August 2023)
.
EU-DEM
2016
version 1.1, European Environmental Agency. EU Copernicus programme. Available from: https://land.copernicus.eu/imagery-in-situ/eu-dem/eu-dem-v1.1?tab=metadata.
Fiedler
W.
,
2016
Managing dam safety risks related to hydraulic structures
. In:
Hydraulic Structures and Water System Management
(
Crookston
B.
&
Tullis
B.
eds.).
6th IAHR International Symposium on Hydraulic Structures, Portland, OR, 27-30 June. doi:10.15142/T3720628160853
.
Filotti
A.
2015
Atitudini față de pericolele naturale provocate de ape (Human attitudes related to the dangers induced by extreme hydrologic events)
. In
Noema, XIV
.
Academia Română
.
FRMP Flood Risk Management Plan – Buzău-Ialomița
2022
(in Romanian, Planul de Management al Riscului la Inundații Buzău-Ialomița). Available from: https://inundatii.ro/wp-content/uploads/2022/10/PMRI_actualizat_ciclul-II_ABA-Buzau-Ialomita_versiune-preliminara-1.pdf.
Gogoașe Nistoran
D. E.
,
Ionescu
C.-S.
&
Simionescu
Ș.-M
.
2023
Influence of dam break scenarios on flood wave characteristics. Case study-Paltinu Dam, Romania
.
IOP Conference Series: Earth and Environemntal Science
1136
,
012031
.
Harten
A.
,
Lax
P.
&
van Leer
B.
1983
On upstream differencing and Godunov-type schemes for hyperbolic conservation laws
.
SIAM Rev.
1983
(
25
),
35
61
.
ICOLD, International Commision of Large Dams
2022
.
ICOLD, International Commision of Large Dams
2023
.
International Rivers, Damming Statistics
2007
Available from: https://archive.internationalrivers.org/damming-statistics (accessed 25 March 2023)
.
Ionescu
C. S.
&
Gogoașe Nistoran
D. E.
2019
Influence of reservoir shape upon the choice of Hydraulic vs. Hydrologic reservoir routing method
.
E3S Web Conf.
85
,
07001
.
EENVIRO 2018 – Sustainable Solutions for Energy and Environment. https://doi.org/10.1051/e3sconf/20198507001
.
Kumar
S.
,
Jaswal
A.
,
Pandey
A.
&
Sharma
N.
2017
Literature review of dam break studies and inundation mapping using hydraulic models and GIS
.
International Research Journal of Engineering and Technology
4–5
,
55
61
.
LeVeque
R. J.
2004
Finite Volume Methods for Hyperbolic Problems
.
Cambridge University Press
,
UK
.
ISBN 0-521-81087-6
Luino
F.
,
Tossatti
G.
&
Bonaria
V.
2014
Dam failures in the 20th century: Nearly 1000 avoidable victims in Italy alone
.
Journal of Environmental Science and Engineering A3
3–1
,
19
31
.
Maranzoni
A.
,
D'Oria
M.
&
Rizzo
C.
2022
Quantitative flood hazard assessment methods: a review
.
Journal of Flood Risk Management
e12855
.
https://doi.org/10.1111/jfr3.12855
.
Morris
M. W.
,
Hassan
M. A. A. M.
&
Vaskinn
K. A.
2007
Breach formation: Field test and laboratory experiments
.
Journal of Hydraulic Research
45
(
S1
),
9
17
.
Nistoran Gogoașe
D. E.
,
Popovici
D. A.
,
Savin
B.
,
Armaș
I.
,
2016
GIS for Dam-break flooding: Study area: Bicaz-Izvorul Muntelui (Romania)
. In:
Space and Time Visualization
(
Boștenaru Dan
M.
&
Crăciun
C.
eds.).
Springer, Cham
, pp.
253
283
.
NRDWR Departament of Natural Resources Division of Water Resources
2010
Guidelines for dam Breach Analysis
.
February 10
.
NWARW – National Waters Administration ‘Romanian Waters’
2012
Project for special surveillance of Paltinu dam (Proiect special de urmărire a barajului Paltinu)
.
Papaioannou
G.
,
Efstratiadis
A.
,
Vasiliades
L.
,
Loukas
A.
,
Papalexiou
S. M.
,
Koukouvinos
A.
,
Tsoukalas
I.
&
Kossieris
P.
2018
An operational method for flood directive implementation in ungauged urban areas
.
Hydrology
5
,
24
.
doi:10.3390/hydrology5020024
.
Singh
V. P.
1996
Dam Breach Modeling Technology
.
Springer Science + Business, B. V
,
Dordrecht.
Stematiu
D.
&
Ionescu
Ș
.
1999
Siguranță și risc în construcții hidrotehnice (Safety and risk in hydrotechnic constructions), Editura Didactică și Pedagogică, București
.
Tanchev
L.
2014
Dams and Appurtenant Hydraulic Structures
, 2nd edn.
CRC
,
Boca Raton, FL
.
Toro
E. F.
2001
Shock-Capturing Methods for Free-Surface Shallow Flow
.
John Wiley
,
Chichester, UK
.
Toro
E. F.
2009
Riemann Sovlers and Numerical Methods for Fluid Dynamics
, 3rd edn.
A practical Introduction, Springer-Verlag
,
Berlin Heidelberg
.
ISBN 978-3-540-25202-3
.
Toro
E. F.
&
Garcia-Navarro
P.
2007
Godunov-type methods of free-surface shallow flows. A review
.
J. of Hydraulic Research
45
(
6
),
737
751
.
Tsai
C. W.
,
Jing-Jie
Y.
&
Chi-Hao
H.
2019
Development of probabilistic inundation mapping for dam failure induced floods
.
Stochastic Environmental Research and Risk Assessment
33
,
91
110
.
https://doi.org/10.1007/s00477-018-1636-8(01)
.
Urzică
A.
,
Mihu-Pintilie
A.
,
Stoleriu
C. C.
,
Câmpianu
C. I.
,
Huț̦anu
E.
,
Pricop
C. I.
&
Grozavu
A.
2021
Using 2D HEC-RAS modeling and embankment dam break scenario for assessing the flood control capacity of a multi-Reservoir System (NE Romania)
.
Water
13
,
57
.
USACE
2014
HEC Engineering Center, Using HEC-RAS for Dam break studies, TD-39
.
USACE
2018a
Hydraulic Engineering Center, HEC-RAS Verification and validation tests. Available from: https://www.hec.usace.army.mil/software/hec-ras/documentation/RD-52_HEC-RAS%20Verification%20and%20Validation.pdf.
USACE
2018b
Hydraulic Engineering Center, HEC-RAS, Benchmarking of the HEC-RAS 2D hydraulic modeling capabilities. Available from: https://www.hec.usace.army.mil/software/hec-ras/documentation/RD-52_HEC-RAS%20Verification%20and%20Validation.pdf
.
USACE
2022
US Army Corps of EngineersHydraulic Engineering Center – River Analysis System, HEC-RAS, Hydraulic Reference Manual. Available from: https://www.hec.usace.army.mil/confluence/rasdocs/ras1dtechref/latest
USBR, USACE
2002
Investigation of the Failure Modes of Concrete Dams – Physical Model Tests Dam Safety Office
.
USBR, USACE
2019
Risk Analysis for Concrete Arch Dams
.
(Chapter E-4), Best practices in Dam and Levee safety Risk Analysis
US Society on Dams
2023
Glossary. Available from: https://www.ussdams.org/resource-center/glossary/#N, Normal Reservoir Level, (accessed 18 August 2023).
Xing
Y.
2017
Numerical Methods for the Nonlinear Shallow Water Equations (Chapter 13)
. In:
Handbook of Numerical Analysis
, Vol.
18
,
(Abgrall, R. & Shu, C-W., eds.). Elsevier, Amsterdam
, pp.
361
384
Xu
R.
,
Alistair
G. L.
&
Xu
B. B.
2022
Application of large time step TVD high order scheme to shallow water equations
.
Atmosphere
13
(
11
),
1856
.
https://doi.org/10.3390/atmos13111856
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY-NC-ND 4.0), which permits copying and redistribution for non-commercial purposes with no derivatives, provided the original work is properly cited (http://creativecommons.org/licenses/by-nc-nd/4.0/).