Extreme flood reconstruction by using the 3DNet platform for hydrological modelling

Flood flow simulations (FFS), aimed at flood forecasting or at assessing the impact of flood protection measures, represent an important application of hydrologic models. Despite recent progress in hydrological modelling, FFS have remained challenging because they involve consideration of many aspects of hydrological modelling including model structure, modelling approach (continuous vs. event-based, lumped vs. distributed), model calibration data and method, etc.

Extreme flow simulations are very sensitive to model structure (Onyutha 2016). Finding an optimal structure that properly reflects processes and enables fast simulations with minimum data requirements remains one of the greatest challenges in FFS (Cloke & Pappenberger 2009). Structure of a model aimed at FFS should enable simulations of different flood generation mechanisms triggered by various precipitation types, such as short convective rainfall events, prolonged heavy rainfall or snowmelt, although Kauffeldt et al. (2016) in their review noted an absence of a snow routine in many models for flood forecasting. Flow routing in rivers is essential for FFS, but in hydrological models it is often either omitted (especially for smaller catchments; Schaefli et al. 2014) or simulated by simplified methods such as the Muskingum or kinematic-wave methods (for a review see Bierkens (2015) and Kauffeldt et al. (2016)). Some models allow a user to select among several routing methods (e.g. HYPE (Lindström et al. 2010) or HEC-HMS (Feldman 2000)), while application of more complex methods such as full dynamic routing is uncommon because of required data on river cross-sections (Van Der Knijff et al. 2010). Accurate estimation of initial conditions (primarily antecedent soil moisture) is also crucial for effective FFS and continuous simulations are therefore preferred to event-based modelling (Berthet et al. 2009).

Spatial variability of precipitation and catchment properties significantly influences the shape of flood hydrographs (Schuurmans & Bierkens 2007), especially during extreme events (Brigode et al. 2014). This spatial variability can be accommodated by distributed hydrologic models, which simulate hydrologic processes in each computational cell and can therefore capture spatial origin of runoff (Beven 2001). Although distributed models are associated with challenging calibration and high computational demand, finer spatial resolution can improve model performance (Chang & Chao 2014; Cui et al. 2014). However, a balance between finer resolution, i.e. more accurate simulations, and computational time of distributed models is necessary. Fine spatial resolution also requires high-quality input data; Her & Heatwole (2016) showed that errors in topographic data have a significant effect on simulated hydrographs, while the effect of errors in land use type is minor.

Improved results without increasing spatial resolution can be obtained by using irregularly shaped computational cells, which concur with catchment topography better than a raster grid (Paniconi & Putti 2015). Furthermore, flexibility in spatial resolution is advantageous as it allows adjusting resolution for a particular catchment with respect to required computational time.

Numerous models of varying complexity are presently used for FFS, ranging from rather complex physically-based models like MIKE-SHE (Refsgaard & Storm 1995), LISFLOOD (Van Der Knijff et al. 2010) or PIHM (Qu & Duffy 2007) to parsimonious models like GR4 J (Perrin et al. 2003). There are also many models in-between like HYPE, PREVAH (Viviroli et al. 2009), TOPMODEL (Beven & Kirkby 1979; Metcalfe et al. 2015), TOPKAPI (Ciarapica & Todini 2002), VIC (Schulla 2016), SAC-SMA (Zhang et al. 2012), HEC-HMS, HBV (Bergström & Frosman 1973) and its subsequent versions like w_flow HBV (Schellekens 2014). Some models are developed principally for regions with specific dominant runoff generation mechanisms; for example, WaSiM-ETH (Schulla 2017), SEHR-ECHO (Schaefli et al. 2014) or PREVAH are suitable for Alpine catchments with considerable snow cover and glaciers. Brauer et al. (2014) presented a model appropriate for simulations in flat areas influenced by high groundwater levels.

Another challenge related to FFS is to provide simulation of inundation outside the main channel and flooding due to levee overtopping. For this purpose, hydrologic models can be coupled to hydrodynamic ones, which have much finer spatial and temporal resolution (see e.g. an application of MIKE-SHE-MIKE11 chain by Clilverd et al. (2016)). However, simulations with these modelling chains are computationally and data demanding.

Human impact reflected in deforestation, urbanisation or river-engineering measures poses a great problem for hydrological modelling in general (Seibert & van Meerveld 2016). River engineering measures, particularly reservoirs, should be incorporated in FFS since they significantly affect flood wave transformation. However, few hydrological models can accommodate reservoirs: e.g. MIKE-SHE, TOPKAPI, HYPE and HEC-HMS.

In this paper, a novel geographic information system (GIS)-oriented hydroinformatic platform 3DNet with the built-in 3DNet-Catch hydrologic model and their application for FFS is presented. The platform provides a range of geo-spatial data pre-processing, manipulation and visualisation functions with a smooth interface to several simulation models, one of which is 3DNet-Catch. The 3DNet-Catch model is developed (1) to represent key processes in sloped catchments of temperate climate, with a structure that represents a balanced compromise between complex physically-based and parsimonious models, and (2) to enable simulations with minimum data requirements. The majority of model parameters are physically based and can be inferred from topography, land use and soil types data, thereby facilitating model calibration.

The platform provides the model with spatial flexibility so that the model spatial resolution can be adjusted to provide accurate and efficient simulations. Spatial representation of catchment processes and features can range from fully-distributed to lumped, depending on the modelling task and available data. This flexibility is facilitated by irregularly shaped computational cells, which can be grouped into different hydrologic response units (HRUs) of various sizes, contrary to the majority of raster-based distributed models (for a review see Kampf & Burges (2007)).

The platform and the model offer another feature relevant for extreme flood simulations, and that is the ability to include various water management facilities such as reservoirs or diversions. By imposing limitations on the capacity of river reaches, flooding volume can also be estimated. In this way, water retention and flooding can be approximately simulated without employing a demanding hydrodynamic model.

To demonstrate the platform and model abilities, they are applied for simulation of the extreme flood in the Kolubara catchment in Serbia in May 2014, characterised by exceptionally high precipitation and severe inundation and water retention outside the channel network. In addition, hydrologic observations during the event are limited and uncertain. Limited flow observations are tackled by the physical connotation of the parameters, while the flow routing routine enhancements allow information on levee overtopping and local flooding to be incorporated in the model. Satisfactory results across the catchment suggest the importance and applicability of these platform features for extreme flood simulations.

The remainder of this paper is organised as follows. The 3DNet platform and its general characteristics are presented concisely in the next section. The subsequent section describes in more detail the Catch module of the platform, intended for data pre-processing for hydrological simulations, and focuses on the development of catchment computational structure and preparation of meteorological forcing. The 3DNet-Catch hydrologic model structure, inputs and outputs are then briefly described in fourth section, followed by the section on model application for reconstructing the extreme flood in May 2014, and the sections with the results, discussion and conclusions on development of hydrologic models and tools for FFS.

The 3DNet is a comprehensive GIS-oriented platform for hydraulic and hydrological modelling developed at the Faculty of Civil Engineering of the University of Belgrade. Its name stems from its three-dimensional representation of water systems (sewer system, stream network, etc.) and terrain (Pokrajac & Stanić 2010). It has been used mainly for hydrological simulations and for design of water supply and wastewater systems in the Western Balkans region (Vasilić et al. 2012). Its wide application is accompanied by its continuous development.

The platform consists of two groups of modules. The first group includes the Terrain and Geology modules for geo-spatial data pre-processing and manipulation. The Terrain is the core platform module and enables creating/importing the digital terrain model (DTM), mesh triangulation to obtain triangulated irregular network (TIN) and manipulation of raster data such as topographic maps. A set of points that comprise the DTM can be either imported from an external file or digitised in 3DNet from an underlying topographical map. The TIN terrain model is generated by applying the constrained 3D Delaunay triangulation over a given set of control points and structural lines that make a plain straight-line graph. The control points are comprised of the terrain points, while the stream network and optionally the watershed divide represent the structural lines. The size of the triangles is limited to a specified area of maximum circumcircle, which is a parameter that has to be specified prior to triangulation (Cheng et al. 2013). The current version of 3DNet does not contain algorithms for DTM corrections such as sink removal, which are available in other GIS software (e.g. ArcGIS, GRASS GIS, SAGA). The Terrain module also contains functions for visualisation of contour lines and path lines, i.e. the lines of the steepest decent from arbitrarily selected points. The optional Geology module enables creating a finite element mesh over geological layers based on borehole data, needed for groundwater simulations.

The second group of modules performs hydraulic and hydrological modelling and corresponding pre-processing operations. The focus of this paper is the Catch module, which is aimed at hydrological simulations and is elaborated further in the following section. The remaining modules are SewNet, WatNet and UGROW. The SewNet module is used for flow simulations in wastewater and stormwater sewer systems using different simulation engines: SWMM 5.0 (Rossman & Huber 2016), BEMUS (Radojković & Maksimović 1984; Radojković et al. 1989) or SIPSON (Djordjević et al. 2005). Interactions between the sewer system and groundwater can also be simulated. The WatNet module simulates flow through pressurised networks with EPANET 2.0 (Rossman 2000) as the simulation engine. The UGROW module is a groundwater flow model based on a finite element mesh enhanced to capture interactions between groundwater and urban water systems and/or streams (Pokrajac & Stanić 2010). It is a convenient tool for various practical problems of urban water management such as assessing the risk of groundwater pollution by leaking sewers or estimating groundwater recharge from water supply system leakage (Todorović et al. 2011).

To enable fast computations, all algorithms are programmed in C++. Data are stored in an accompanying database with two-way communication with 3DNet through SQL. Database structure is pre-set to accommodate geo-spatial and water-related input data (e.g. elevation-discharge curves, water consumption patterns, meteorological observations, etc.), processed geo-spatial data (e.g. areas, lengths, slopes, etc.) and simulation results.

Communication between 3DNet and end-users is provided through a simple, intuitive and easy-to-use graphical interface (Figure 1). The embedded modules are represented by drop-down menus and by corresponding data groups in the object panel (left panel in Figure 1). The structure of the object panel in 3DNet is pre-set and, differently from other GIS software, cannot be altered by the user. However, some objects can be organised into user-specified layers (e.g. meteorological stations; see bottom of the object panel in Figure 1). Three buttons beside the objects serve to set data and drawing properties (the L button is for the layer properties, while the A and P buttons stand for object attributes and properties). The objects in the workspace can easily be manipulated by tools from the toolbar. The 3DNet platform enables effective 3D visualisation of terrain, water systems and simulation results (e.g. groundwater levels, hydrographs at selected points).

The Catch module comprises algorithms for data pre-processing and for hydrological modelling. This section describes the data preparation algorithms and data structure, while the hydrologic model is described in the next section.

Data preparation for hydrological modelling includes: (1) catchment computational structure (i.e. stream network, hydraulic structures and associated drainage areas), (2) properties of drainage areas (slopes, areas) and stream network (lengths, slopes) and (3) meteorological forcing.

Geo-spatial data for hydrological modelling are stored in the Catch group in the object panel and are organised in a tree-like structure (Figure 1). The Catch group branches into five subgroups: Graph, HRUs, Areas, Polygons and Meteo. The Graph subgroup represents the stream network. The HRUs, Areas and Polygons subgroups represent spatial computational units. The Meteo subgroup represents meteorological stations. Each subgroup (branch) may contain point-, line- and polygon-type objects (leaf data).

The Graph data subgroup representing the stream network contains point objects (nodes, vertices, outlets, hydro-profiles (HPs)) and line objects (links). Nodes and vertices are the basic point objects, while HPs represent computational nodes at which hydrographs are simulated. Typically, HPs are set at stream gauge locations, bifurcations, diversions and at other user-defined points where hydrographs are required. A special type of point objects are sinks and non-linear HPs that represent reservoirs for which the elevation-volume and elevation-discharge curves must be specified. An identifier, coordinates, name, downstream HP and the layer, and the results of hydrologic simulations are stored in the database for each HP. In addition to these data, identifiers of the related elevation-discharge and elevation-volume curves are saved in the database for non-linear HPs. Simulation results at HPs are obtained/stored in regular time intervals, which, together with the HP coordinates, provide adequate data for coupling with a hydraulic routing model (Harpham & Danovaro 2015). The outlets can be user-defined; if not provided, they are automatically inferred from the topologically sorted HPs (similarly to Drainage Point Processing Tool in ArcGIS; Djokic et al. 2011).

The stream network is comprised of links represented in 3DNet by a directed acyclic graph. This means that each link should be pointed towards the downstream node; otherwise, the user will be informed that the stream network is invalid. Geometric properties of a link are automatically calculated from coordinates of nodes and vertices. Based on these values and predefined roughness, initial estimates of mean flow velocity and travel time for the link are computed. Flow downstream of a bifurcation is evenly distributed among the links by default, but the user can specify it differently. Limiting capacity can be imposed on a link in order to estimate the flooding volume for the link. For each link its identifier, identifiers of the upstream and downstream nodes and links, and the layer, length, slope, roughness, travel time, flow fraction (important for links downstream of a bifurcation), flow threshold, and the results of hydrologic simulations are stored in the database.

HRUs represent computational cells of the hydrological model. They are obtained by producing a Voronoi diagram from the triangles of TIN (Tucker et al. 2001) as shown in Figure 2(a). All points within an HRU are assumed to exhibit hydrologically similar behaviour because of geographical proximity. The advantage of the TIN-based HRUs is that smaller units are generated over areas with distinctive topography and/or dense stream network where more topographic features have to be captured, while HRUs over featureless topographies are reasonably larger but are still limited by the size of the triangles. This approach is recommended by Dehotin & Braud (2008) as a balance between accuracy and spatial resolution, but is not common in hydrological modelling (Kampf & Burges 2007). It is implemented in the tRIBS model (Vivoni et al. 2007), while in the PIHM model catchment is discretised into triangular elements of TIN.

Polygons are created by aggregating HRUs to the (sub)catchment level, allowing semi-distributed or lumped hydrologic setup to be created and computational time to be reduced. HRUs are aggregated by means of subcatchment delineation. In this procedure, each HRU is connected to a downstream HRU having the centre of gravity at the steepest descent from the centre of the considered HRU. The HRUs over the stream network are automatically connected according to network topography, thereby creating a weighted oriented graph with slopes acting as weights. Over this graph the Priority First Search (PFS) algorithm is applied (Sedgewick 2002). Initially, HRUs intersected by the stream network upstream of a HP are attributed to the considered HP. Selection continues along the steepest slope, as illustrated in Figure 2(b), until the water divide is reached. As a result, each HP is assigned a drainage area (Figure 2(c)) and each HRU is assigned a drainage HP to which the HRU's runoff is routed. Optionally, a different drainage HP may be assigned for baseflow. This unique feature of this platform is a particularly convenient option for karstic catchments because it enables including soft data on groundwater flow in the model (Jaćimović et al. 2015). Some HRUs may remain unselected if stream network is not in their vicinity in karstic areas, or due to errors in topographic data. In 3DNet, such HRUs are considered karstic sinks and entire precipitation over them is routed to the nearest downstream HP by default, but this can be changed by the user. This distinct feature of the platform makes the catchment spatial representation less sensitive to errors in DTM than in other GIS tools, in which catchment delineation is based on raster DTM and needs terrain pre-processing (Djokic et al. 2011; Fleming & Doan 2013).

The Areas subgroup may contain HRUs or other units (e.g. digitised subcatchments or elevation zones) that are created externally and imported into 3DNet. Both Polygons and Areas subgroups comprise same leaf-data, including points, catch, and cover. The points refer to the constituent points of other objects in the subgroup; the catch refers to HRUs (either aggregated or imported), while cover comprises polygon objects representing land-use, soil, vegetation or similar data.

For each polygon-type object its geometrical properties, model parameters and simulation results are saved in the database. Geometrical properties stored in the database include perimeter, area, slope, aspect, coordinates of the centre of gravity, identifiers for direct runoff and baseflow drainage HPs, length and slope of the longest flow path (overland and channel flow) and calculated time of concentration. 3DNet calculates geometrical properties of the objects automatically, while other GIS applications require running specific tools (Djokic et al. 2011).

Predefined model parameters for each HRU are obtained from cover data. The cover polygons have specified parameter values and priority. If a HRU's centre of gravity is overlaid by several cover polygons (e.g. land-use and soil type) that suggest different parameter values, the HRU is assigned the parameter value defined by the polygon with the highest priority. Unlike for example GeoHMS, the 3DNet platform allows all model parameters to be spatially distributed and to vary across HRUs rather than across subcatchments. Additionally, no pre-processing is required in 3DNet, while transferring parameters from polygon object to a raster cell in GeoHMS requires intersection of these layers (Fleming & Doan 2013).

The Meteo subgroup contains points representing meteorological stations and corresponding data. The stations metadata (object and layer identifiers, name and coordinates) and the observations are stored in the database. Meteorological forcing for each HRU is obtained from the point meteorological observations either by the inverse weighted distance method or the nearest neighbour interpolation method. Interpolated precipitation and temperature series can be adjusted to account for change with elevation with user-defined precipitation and temperature gradients. There is no spatial variation of meteorological variables within an HRU, but according to Tarasova et al. (2016) elevation adjustment provides a reasonable representation of meteorological spatial heterogeneity. Potential evapotranspiration (PET) can be imported as input data, or calculated within 3DNet using the Hargreaves equation (Hargreaves & Samani 1982). Unlike other tools (e.g. GeoHMS), the current 3DNet version does not support raster data (radar observations). Several alternative meteorological series can be defined for one period (e.g. by applying different interpolation methods).

The 3DNet-Catch is a novel distributed hydrologic model developed for continuous hydrologic simulations. For each HRU runoff volume is simulated from the vertical water balance within the interception, snow and soil routines (upper shaded part of Figure 3) and routed to the drainage HP by the response routine (lower hatched part of Figure 3), resulting in direct runoff, fast groundwater discharge and baseflow. All model routines are interconnected in such a way that water conservation is assured.

Vegetation is represented by a canopy reservoir for which rainfall is the input, while throughfall and canopy evaporation E_CAN comprise the output. The canopy reservoir capacity varies over the growing season according to Leaf Area Index (LAI) and is limited by the maximum reservoir capacity CAN_max as a free model parameter.

Snowfall occurs at temperatures below the discrimination one (T_R-S) providing input to the snowpack reservoir with snowmelt and sublimation E_SUB as the outputs. The snowmelt is simulated with the degree-day method (Neitsch et al. 2011). Spatial heterogeneity of the snowpack and intra-annual meltrate variation are accounted for.

Percolation is limited by the residual soil moisture, which is approximated by volumetric soil water content at permanent wilting point (in m³·m⁻³). Evaporation from soil decreases exponentially with a decline in soil moisture. Percolation from the upper layer is the input to any subsurface soil layer, while the output comprises transpiration E_T and percolation to deeper soil layers or groundwater reservoir. Percolation from all soil layers is computed using the same equations, but the parameters may differ across the layers. Transpiration is distributed among subsurface layers according to their thickness and is limited by the available soil moisture (Neitsch et al. 2011).

Percolation from the deepest soil layer W_PERC,GW is routed through a nonlinear groundwater reservoir to a point that does not necessarily coincide with the surface runoff drainage point. Water in this reservoir below the threshold S_max is routed by using a nonlinear outflow equation resulting in baseflow Q_BASE. If the threshold S_max is exceeded, excess water is routed through a liner reservoir resulting in fast groundwater discharge Q_{GW_FAST}. Flow along a river reach is routed using the linear reservoir equation. The reservoir coefficient is equal to the travel time along the section.

Calibration of the model generally involves optimisation of a relatively high number of parameters. However, many parameters can be estimated from the soil, land use and vegetation data (represented by the cover polygon objects). Each HRU may be assigned a unique parameter set, or one set may be attributed to all HRUs with similar properties.

The 3DNet-Catch model is coded in C++, enabling fast computations. The code is implemented as a dynamic-link library and can be used independently of the 3DNet platform (e.g. it can be connected to external programs with optimisation algorithms).

Geo-spatial and hydro-meteorological data should be supplied to 3DNet-Catch for hydrologic simulations. DTM, stream network and water divide are required, along with the stream gauges, meteorological stations, hydraulic structures and other computational points on the stream network. These data can be either imported or created in the 3DNet platform using background maps. Since 3DNet enables subcatchment delineation within a catchment by grouping HRUs, the stream network and the water divide cannot be created in 3DNet from DTM like in other commercial GIS tools (ArcGIS or GRASS GIS). Point elevations of every object can be defined relative to terrain. All objects can be exported to dxf format and used independently of 3DNet. Land-use, soil and vegetation type maps are optional, though they are recommendable to facilitate model calibration.

Meteorological input includes observed precipitation, maximum and mean temperatures, and PET rates at considered meteorological stations. If PET rates are not provided, minimum daily temperatures are required to enable PET calculation with the Hargreaves method. Observed flows are needed for model calibration. The input data series have to be provided in regular intervals that should correspond to the simulation time step.

Hydrologic simulations provide values of state variables for each HRU (e.g. soil moisture, snow water equivalent, actual evapotranspiration, etc.), hydrographs at every HP (direct runoff, fast groundwater discharge and baseflow), and inflow, outflow and elevation/volume for every reservoir. The user can select a simulated variable of interest from a drop-down menu in the results window for visual inspection or export. The results across the catchment can be visualised in 3DNet simultaneously, which can enable for example tracking runoff origin or detecting flood-prone areas. All results related to a specific object (HRU or HP) are stored in the database and can be used externally.

In this paper, the 3DNet-Catch model is applied for reconstruction of the devastating flood in the Kolubara catchment in Serbia in May 2014. The catchment total drainage area is 3,635.5 km² and elevations range from 76 to 1,346 m a.s.l. (Figure 4) with rather steep headwater areas. Floods occur mainly in early spring due to combined rainfall and snowmelt (Todorovic & Plavsic 2016), but can also occur in summer due to intensive convective rainfall events, which can trigger quick response of the torrential headwater tributaries (Plavšić et al. 2014; Petrović et al. 2015). Flood waves from the tributaries reach the main course almost simultaneously due to catchment geometry.

Forests and arable cultivated land prevail in higher altitudes, while the lower part of the catchment is urbanised and heavily populated. The majority of the 335,000 inhabitants are settled in the lowland, where almost the entire industry of this region is located, including major open-pit coalmines and thermal power stations. The coalmine provides 75% of the total coal excavated in Serbia, and it supplies local thermal power stations, which cover over 40% of the total electricity consumption in Serbia (EPS 2015). These facts indicate extreme vulnerability to floods in the lower Kolubara catchment. Retentions and levees constructed to mitigate flood risk are primarily to prevent flooding of the coalmine and power stations.

The devastating flood in May 2014 that hit the west Balkans region was a result of combined heavy prolonged rainfall, nearly saturated soil and high stage in the Sava River (Plavšić et al. 2014). A very rainy period preceded the major flood in May 2014. In April 2014 there were 19 rainy days with total precipitation of 177 mm observed at Valjevo, which is almost triple the average April precipitation of 62.2 mm. On 14th May 2014 a cyclonic storm system reached western Serbia and persisted over the Kolubara catchment producing excessive precipitation amounts. The southwestern and northern parts of the catchment received the largest precipitation amounts (see Figure S1 in the supplementary material, available with the online version of this paper). For example, the town of Valjevo received 189 mm between 14th and 16th May, with 108 mm on 15th May alone. Total precipitation at Valjevo in May 2014 reached 324 mm, which is 52% greater than the previous May maximum. Total precipitation at Valjevo in April and May 2014 was 501 mm, which is almost twice the previous maximum and makes over 60% of the mean annual precipitation (790 mm). Rainfall intensity during this event was modest but relatively steady. The return periods of the largest rain depths over 12, 24, 48 and 72 hours in this event at Valjevo are estimated at 3, 60, 400 and 650 years, respectively (Prohaska & Zlatanović 2015). Peak flow of the Kolubara at Slovac in May 2014 reached 1,100 m³/s, which is substantially greater than the previous maximum of 320 m³/s.

Extremely high flows led to levee breach and overtopping, major landslides and extensive flooding of many towns and villages in the catchment. It is estimated that the open pits of the coalmine received over 2 million m³ of flooding water from Kolubara and its tributaries. The flood caused enormous damage to infrastructure and property, estimated at more than 1.5 billion Euros (Serbia Floods 2014). Subsequent landslides caused further devastation. Twenty-five casualties were reported.

To reduce flood risk in this catchment, the construction of 19 reservoirs is proposed, whose effects are evaluated in ongoing studies. The 3DNet-Catch model of the Kolubara catchment based on the May 2014 flood event is developed as a part of these efforts. The model is required to reproduce flood flows accurately and to integrate the existing and planned reservoirs.

Rain depths and temperatures at 30 meteorological gauges are available between 11th and 18th May 2014 (Figure 5). Observed daily precipitation depths are downscaled to hourly level based on the hourly data from three recording rain gauges within the catchment and another six in its vicinity. Hourly precipitation at non-recording rain gauges is estimated from weighted dimensionless rainfall mass curves at the recording rain gauges, with areas of the Thiessen polygons as weights (Zlatanović & Prohaska 2015).

Hydrologic data are observed at 13 stream gauges within the catchment. Many stream gauges were washed away, while some had gaps in recording (see for example the observed hydrograph at Bogovadja in Figure 6(d)). Peak flows at some sites are underestimated due to upstream flooding (e.g. Beli Brod in Figure 6(c)). Uncertainty in estimated flow rates is significant since the flows are estimated by extrapolating rating curves to the range of stages never observed before. A comprehensive review of data availability and uncertainty during the May 2014 flood is provided by Prohaska & Zlatanović (2015) and Zlatanović & Prohaska (2015).

For reconstruction of the flood in May 2014, HRUs are obtained from the TIN network created from 30 m·30 m DTM and with the specified maximum circumcircle area of 10,000 m². The HRUs are aggregated into 125 polygons, which are then grouped into 14 subcatchments (shaded areas in Figure 5) draining into HPs that represent the stream gauges and open pits. To simulate flood retention immediately upstream of Ćemanov most caused by a bridge on the Tamnava River, this gauge is represented by a non-linear HP and assigned an elevation-volume curve based on the available cross-section data. The elevation-discharge curve is obtained from hydraulic analysis based on an upstream cross-section and the bridge opening geometry with assumed energy losses. Flooding from the stream network to the open pits is routed through transient links indicated by dashed lines in Figure 5. Maximum flows are imposed on the links representing reaches where flooding to the open pits occurred. Other less significant flooding is not included in the model. Non-linear HPs are also set at locations of the reservoirs. Undrained areas in the upper part of the catchment represent karstic sinks (white areas in Figure 5) and rainfall intercepted by them is routed to the nearest downstream HP.

The semi-distributed model setup is used, meaning that each subcatchment is assigned a unique parameter set. Meteorological forcing is obtained for every unit by the inverse-distance interpolation method, without any adjustment with elevation due to the absence of clear correlation between the observations and station altitudes. Evaporation is neglected due to the continuous rainfall and high humidity in the simulation period.

Model structure with one surface and one subsurface soil layer is employed because the effect of more subsurface layers would not be noticeable in simulations over short periods of several days. Both soil layers are assigned the same parameters except the layer thickness. Variation in vegetative cover over short periods is negligible, thus the LAI parameter is constant in time but allowed to vary across subcatchments. Values of the reservoir coefficients obtained from the predefined roughness and slopes of drainage areas and streams are kept during calibration. Two parameters of the interception routine and eight parameters of soil routine for every subcatchment are estimated through calibration, along with the limiting capacities of the links where flooding into the open pits occurred. Dynamics of the inflow to the pits is unknown, so the link capacities are adjusted to reproduce the intercepted flooding volume.

Model calibration against such extreme flows posed a great challenge due to incomplete and uncertain observations. A trial-and-error calibration is therefore exercised in terms of trade-off between hydrological plausibility of the parameters and consistency of simulations with observations. The calibration begins with the head subcatchments and continues towards the catchment outlet, and is summarised as follows:

• Where the observed hydrographs are considered accurate, the model parameters are adjusted to fit the simulated hydrographs to the observed ones in terms of peak magnitude and timing, rising and recession limbs, and total flow volume.

• Where peaks flows were missing or considered inaccurate due to upstream flooding, the parameters are calibrated to reproduce rising and falling limbs as closely as possible, since these hydrograph parts are considered accurately observed.

• In the absence of observations, the parameters are calibrated against flows at the downstream gauge jointly with the parameters of draining areas within the downstream subcatchment. Calibrated parameters from the neighbouring hydrologically similar subcatchments are used as the initial parameter estimates.

The simulation with the hourly time step covers the period between 08:00 on 11th May and 07:00 on 18th May, in total a duration of 168 hours.

The trial-and-error calibration of the 3DNet-Catch model with the semi-distributed setup against the observed May 2014 flood hydrographs, described in the previous section, resulted in two spatially lumped parameters: surface soil layer thickness of 5 cm and the volumetric soil moisture content at a permanent wilting point of 0.005 m³·m⁻³. Frequencies of the remaining spatially varying parameters of the soil and interception routines are shown in Figure S2 in the supplementary material (available with the online version of this paper). The CN takes high values from 70 to 80, which is consistent with high soil moisture due to antecedent precipitation. High soil moisture content prior to the event is reflected in relatively low values of soil conductivity (from 1.5·10⁻⁷ to 6.7·10⁻⁷). The subsurface layer depth exhibits high spatial variation. HRUs in the lower part of the catchment generally have thicker subsurface layer, which confirms the results of Han & Liu (2015).

The simulated vs. observed hydrographs at selected stream gauges are presented in Figure 6. Consistency between the simulated and observed flows at Valjevo and Slovac is apparent (Figure 6(a) and 2(b)), with tolerably underestimated flood peak at Slovac. The rising and falling limbs at Beli Brod (Figure 6(c)) are satisfactorily reproduced, although the observed rising limb near the peak becomes steeper than the simulated one. The observed hydrograph peak cut-off could be reproduced if data on flooding upstream of Beli Brod were available and included in the model. Figure 6(d) presents the results for Bogovadja with an incomplete observed hydrograph where only the peak flow was estimated (represented with the rectangular marker). The simulated hydrograph at Bogovadja has limbs that correspond to the observations. Since the simulated hydrograph at Bogovadja is comprised in the downstream hydrograph at Beli Brod, good performance at Beli Brod also suggests that the Bogovadja hydrograph is well reconstructed.

Flood retention upstream of the Ćemanov most stream gauge is captured owing to the non-linear HP (reservoir) at this gauge (Figure 6(e)). Hydrographs at Degurić, Ub and Koceljeva are also acceptably reproduced (Figure S3 in the supplementary material, available online). Considerable discrepancies at Draževac (also shown in Figure S3) could be attributed to measurement errors since the observed hydrograph at Draževac is significantly delayed compared to others. This lag at Draževac cannot be resolved by calibration because the simulated hydrograph is primarily determined by the upstream hydrographs and flooding to the open pits, while the contribution of its subcatchment is minor.

The performance measures given in Table S1 in the supplementary material (available online) suggest reasonable model performance. Expectedly, the lowest NSE is obtained for Draževac (0.59). The second lowest NSE of 0.78 is at Degurić where the model cannot reproduce flow fluctuations near the hydrograph peak (Figure S3). NSE at other sites varies between 0.86 for Beli Brod and 0.98 for Valjevo, indicating a very good fit. Coefficient of determination R² shows similar results and takes high values at most gauges revealing properly reproduced flood dynamics. High values of VE and low VB indicate accurately reproduced flow volume. Interestingly, the best agreement in terms of volume is achieved at Draževac; this may be attributed to accurately simulated flooding volume in the coalmine. VB at the remaining gauges takes values between −1.6% at Valjevo and 21.7% at Degurić, and is considered tolerable. Mainly positive VB, suggesting overestimated flow volume, is probably due to other floodings that are not included in the model. Bias in flood peak magnitude ranges from −11.4% for Bogovadja to 14.6% for Beli Brod, but its absolute value is below 10% for most gauges. Considering uncertainty in the observations, especially in flood peaks, bias in simulated peak flows is tolerable. There are no differences between timings of simulated and observed peaks at Valjevo and Bogovadja, but at the remaining gauges the simulated peaks are delayed for a few hours (1 h at Ub, 2 h at Degurić, 3 h at Slovac, 5 h at Koceljeva). Better agreement of peak timing with observations is rather difficult to obtain by calibration because it is principally driven by rainfall dynamics and is insensitive to most model parameters.

It should be noted, however, that reaching high performance measures is not the primary objective of calibration in this flood reconstruction due to the short simulation period and data uncertainty. To avoid overfitting of model to noisy data, the calibration is guided by parameter plausibility and by visual inspection of the hydrographs.

Simulated inflows to the open coalmine pits are shown in Figure 6(f). The simulated flood volumes amount to 166, 23 and 17 million m³ for the Tamnava-east, Tamnava-west and Veliki Crljeni pits, respectively, and correspond to the values estimated by Prohaska & Zlatanović (2015).

Simulation of the May 2014 flood event with the semi-distributed 3DNet-Catch model based on 125 polygons and in 168 time steps lasts less than a second, i.e. 3·10⁻⁶ s per polygon and per time step on Intel ^®Core ™ i7 CPU 3.4 GHz with 16 GB RAM. Certainly, impact of spatial aggregation and the fact that all algorithms are implemented in C++ would become evident in longer simulation periods and/or for catchments discretised into more computational units.

Properly reproduced hydrographs across the catchment suggest that simple enhancements to the model routing routine that allow simulation of water retention and flooding are particularly useful for FFS since these processes are common during floods. For example, inclusion of the Rovni reservoir in the model resulted in excellent performance at downstream Valjevo. Similarly, good performance at the Ćemanov most stream gauge is primarily attributed to the non-linear HP that allowed simulation of the water retention. The majority of hydrologic models with a simple flow routing routine do not comprise such enhancements. An exception is HEC-HMS, which enables inclusion of retentions and diversion in the stream network and simulating of retention and flooding (Fleming & Doan 2013). The catchment computational structure in 3DNet is TIN-based, which has proven advantageous in comparison to lumped setup or raster-grid supported by GeoHMS. Comparison of this model setup to alternative ones (e.g. raster-based), as well to other model structures, requires further research.

The 3DNet platform enables TIN-based HRUs aggregation, but the sensitivity of simulated flows to degree of aggregation requires additional analyses. Rainfall was relatively uniform over the catchment during this event, and further research is needed to analyse how the spatial aggregation affects simulated catchment response to, for example, localised convective rainfall. It is therefore recommendable to apply the aggregation carefully, with spatial resolution that is adjusted for each individual catchment with respect to required simulation accuracy.

To improve accuracy of FFS, more robust flow routing methods should be implemented in 3DNet. Additionally, the current model version estimates reservoir outflow according to an elevation-discharge curve, and therefore the platform should be further developed to allow simulations of reservoir operation.

In this paper, the 3DNet platform for hydraulic and hydrological simulations is presented, with the focus on the platform features relevant for FFS. The platform is used to develop a semi-distributed model for reconstruction of the devastating flood in the Kolubara catchment in May 2014. In addition to simulating the hydrologic processes in the catchment, additional features are introduced to simulate flooding. This is achieved by imposing limitations on the channel capacities in selected river sections and by introducing additional point objects that represents reservoirs (i.e. non-linear HPs) at locations with significant flooding, so that water retention could be taken into account. Incomplete and uncertain flow observations necessitated trial-and-error calibration guided by a physical plausibility of model parameters.

The model yields a hydrologically sound catchment response and the simulations proved consistent with reliable observations, such as the rising and falling hydrograph limbs and the estimated volumes of water captured in the coalmine pits. In this way, a comprehensive and effective tool aimed at evaluation of the effect of the proposed reservoirs on flood flows is obtained.

The presented extreme flood reconstruction reveals useful features for FFS:

• Inclusion of hydraulic structures in the model by adding reservoirs (non-linear HPs) for simulating the effect of retention on flood waves enables evaluating the effect of the existing and proposed reservoirs on flood mitigation. Additionally, inundation and water retention at some locations in the catchments are also simulated by means of these reservoirs, thus enabling mass conservation principle to be applied fully.

• Imposing limiting capacity on selected river reaches also facilitates estimation of the flooding volume. This approach cannot replace the hydrodynamic models nor provide inundation maps, but it allows quick approximate estimation of the flooding volume and consequently results in more accurate flow estimates downstream of a considered section.

• Physically meaningful model parameters prevent model overfitting to unreliable noisy data. This feature is important because the models intended for FFS should be mandatorily calibrated against extreme floods, which are inevitably accompanied by limited and uncertain data.

• Model spatial flexibility is important to provide a balance between the accuracy in describing the catchment processes and computational speed in accordance with different modelling aims and different levels of quality of the available data.

These features of the platform and the model might therefore be considered in development/enhancement of other hydrologic models and tools aimed for extreme flood simulations.

The research presented is part of the project ‘Increased Resilience to Respond to Emergency Situations’, funded by the Government of Japan and implemented by the United Nations Development Programme (UNDP). Data were provided by Republic Hydrometeorological Service of Serbia. The authors are also grateful to three anonymous reviewers for providing the comments that helped to better present our work.