Assessment of water resources system resilience under hazardous events using system dynamic approach and artificial neural networks

Water resources systems are designed to withstand demands imposed by their service requirements (Simonovic 2021). However, their flood protection facilities are designed by existing standards and under the present climate. Considering the aging process and rapid changes in the environment (e.g., climate variability, climate change, natural hazards), they do not necessarily guarantee an adequate level of service and safety. Hence, major investment is required to renew and upgrade these aging water systems to adapt to the rapidly growing population, whose future is affected by changing climate and natural disasters (e.g., floods, earthquakes).

Climate change leads to flood events characterized by greater frequency and more severity (Liang et al. 2019), and therefore greatly impacts the probabilistic behavior of extreme annual streamflows (Xu et al. 2020). Given that, one of the main roles of water systems is responding to natural disasters, such as floods, by controlling the outflow from the water system (Men et al. 2019; Yang et al. 2022). To cope with such hazards under present and future climate, water resources systems have to be revitalized for several reasons, such as poor initial design of the water system, rapid changes in the downstream river sections (e.g., urbanization), and safety requirements by increasingly risk-averse societies (Bocchiola & Rosso 2014).

Alongside climate-related risks, earthquakes also present a significant hazard that can reduce the safety of the water resources system (ICOLD 2018). Earthquakes commonly trigger landslides producing consequences ranging from small soil cover failures to devastating rock avalanches (Fan et al. 2019). As a part of a water system, a large dam can fail under severe seismic conditions leading to uncontrolled release from the system, and require optimal operational rules for gated spillways (Haktanir et al. 2013). Moreover, a water resources system located within a high seismic area possesses a substantial risk for the downstream river sections since active faults may cause damaging deformation of the embankment leading to uncontrolled water system behavior (Tosun 2015).

Risk analysis of the water resource systems needs to consider a variety of hazardous events (e.g., floods, earthquakes), or a combination of events that could lead to the failure of the system (Bowles et al. 2013). It can be performed through the qualitative or quantitative description of potential hazardous scenarios including external threats (e.g., earthquakes, floods), internal threats (e.g., failure of the system elements), dam failures, etc. The identification of potentially hazardous events is commonly quantified using expert opinions regarding the occurrence probability of hazards, their severity as well as detectability (Ardeshirtanha & Sharafati 2020). Generally, the four-step procedure shown in Table 1 is usually applied to identify core hazardous scenarios leading to loss of strength and water tightness, or even a breach of the dam (Bowles et al. 2013).

As an essential treatment of the water resources system functionality, a procedure for the identification of potentially hazardous events highlights a combination of hazards that can jeopardize the water system services and safety. In contrast to a single hazard risk assessment, quantification of multiple hazards represents a challenging issue within risk analysis, including hazard assessment, vulnerability characterization, and risk assessment (De Angeli et al. 2022). Multiple hazards arise from the spatial and temporal coincidence of potentially interconnected hazards affecting the water system elements and level of safety (Kappes et al. 2012). Their relationships primarily consider a situation where the first hazard triggers or increases the probability of secondary hazards occurring (Gill & Malamud 2014). This can be achieved by using visualization techniques based on the large amounts of the recorded hazards (Kappes et al. 2012).

Furthermore, a multi-hazard modeling approach can be implemented within the mathematical framework, consisting of risk analysis, analytic hierarchy processing, and fuzzy logic to examine multi-hazard interactions and to provide a comprehensive risk assessment of the water resources systems (Jiang & Zhang 2008; Greco et al. 2012; Hwang et al. 2015; Assumma et al. 2019; Nozari et al. 2021; Datola et al. 2022; De Angeli et al. 2022). A general framework for multi-hazard modeling is based on the use of component operating states database of the water resources system, where the failure scenarios are simulated via the application of the deterministic Monte-Carlo (MC) approach to determine a wide range of water system behaviors for each scenario (King & Simonovic 2020). The generated scenarios vary the timing and severity of hazardous events as well as the inflows in the water resources system capturing the system's response to the full range of potential operating conditions (King et al. 2019).

However, the flood risk assessment frameworks do not research in detail the coincidence of hazardous events with low joint frequencies and their impacts on the water system behavior, mainly because a wide range of potential hazardous events and the system's responses require machine learning techniques to deduce the knowledge from a significant amount of data. Additionally, the system response to hazards, like earthquakes, imposes the use of a numerical dam safety model to precisely extract the relationship between the system element functionality and earthquake magnitude (Rakić et al. 2022). To examine the impacts of natural hazards with a high magnitude and low joint frequencies on the water resources system functionality, the presented research proposes the usage of a system dynamics modeling approach underpinned by artificial neural networks (ANNs) for flood-related risk prediction. It incorporates the developed water resources system and hazardous models, as well as outputs from the dam safety model. The outputs from system dynamics models are used for flood dynamic resilience evaluation since it is selected as a time-dependent flood risk metric (Ignjatović et al. 2021). Implemented within the system dynamics modeling approach, hazardous events include the spatial and temporal coincidence of flood and earthquake under stochastically generated inputs, followed by the variable initial state in the water resources system. An ANN is then trained over the generated inputs and outputs from the system response to quantify the flood-related risks in real-time and to support the decision-making process within the river basin under hazardous events.

Therefore, the goals of the presented research are given as follows: (1) proposition of a generalized methodological framework to quantify flood dynamic resilience of water system using a significant amount of generated data related to multiple hazards (e.g. floods, earthquakes, variable state of the system), (2) implementation of the results derived from the dam safety model for accurate estimation of the earthquakes' impacts on the water resources system functionality, and (3) application of ANN to predict the main characteristics of the flood-risk related matrices (flood dynamic resilience) based on the water system response.

The paper is structured in the following way: the first section envelopes the methodological approaches for developing the system dynamic model, hazardous models, and dynamic resilience assessment of the water resource system. Next, the case study illustrates the selected water resource system located in Serbia (southeast Europe). Afterwards, the results of the proposed methodology are depicted to illustrate model development and simulations performed. The next section discusses the results and links them to relevant research from the literature. In the final section, the implications of the presented research are discussed and directions for future research are proposed.

The main objective of this research is to devise a generalized risk-based methodology for the evaluation of flood dynamic resilience of water resources systems using machine learning. To achieve this objective, a system dynamics (SD) model is developed to mimic the non-linear behavior of the water resource system. It provides information required for flood dynamic resilience assessment as a wide range of possible outcomes depending on the initial reservoir levels, inflows, and timing of events (King & Simonovic 2020). The SD model is supported by hazardous models which introduce the time coincidence of several unexpected hazardous events (earthquakes, floods, and changeable water content in a reservoir). Specifically, the earthquake hazardous model simulates the impact of earthquakes on the element of the system reservoir (e.g., the impact of the decreased capacity of the gated spillway). A flood dynamic over the simulation period is represented by the flood hydrograph model, as a part of hazardous modeling. Moreover, the changeable initial water content in the reservoir at the time of system failure is generated since the gate's operation depends on inflows and the initial state in a system. The hazardous models introduce the randomized model parameters enabling the generation of hazardous events with different magnitudes, especially in terms of low-frequency hazardous events. Alongside the SD model, hazardous models support the assessment of flood dynamic resilience for a variety of possible scenarios. Then, a number of the generated hazardous events and estimated dynamic resilience are used as a solid basis for the extraction of the knowledge from the generated data via ANN. The ANN captures the important characteristics of dynamic resilience and enables its application for real-time decision-making in a water resources system.

A deterministic, distributed-parameter, physical process-based hydrologic model (Precipitation-Runoff Modeling System) is used to generate the long-term flow sequences (Ignjatović et al. 2021). It enables the evaluation of the hydrological response based on various climate combinations including the generated long-term sequences of climate data. As the inputs for the hydrological model, the stochastically generated climate data by weather generator is utilized (Ignjatović et al. 2021). The non-parametric K-nearest neighbor (K-NN) weather generator is employed by reshuffling the recorded value of precipitation as well as minimal and maximal air temperature as a basis for daily streamflow assessment (King et al. 2015). Additionally, the reshuffled climate records are perturbed with the aim of generating climate data beyond the historically observed values. Such generated data are used to force a hydrological model to produce daily streamflows over long-term periods.

In this study, the annual maximum series approach is applied to estimate low flood frequencies from the generated annual streamflows (Karim et al. 2017). The general extreme value (GEV) distribution is fitted to annual maximum streamflows since it provides the best fit for extreme flow as well as precipitation data within the region analyzed (Stojković et al. 2017; Gocić et al. 2021). By application to the hydrological data, the flood magnitudes with different return periods (⁠⁠) are estimated. In addition, the shape of a theoretical flood hydrograph is approximated by the Gumbel distribution (Manfreda et al. 2021). Using the flood magnitude and shape of the flood hydrograph, flood hydrographs are defined for the set of return periods.

The varying earthquake events are generated as a stress test of the water resources system addressing the adverse impacts on the water resources system elements considering the fact that earthquakes have frequently occurred within the analyzed area and wider (Maršanić et al. 2021).

The modeling of the water resources system involves an initial water state (⁠⁠) in the reservoir at the first stage of decision making at (Kistenmacher & Georgakakos 2015). The state in the reservoir at the next simulation stage (⁠⁠) is estimated by applying the first stage decisions (e.g., gates opening and hydropower release at ⁠) and first stage inflow (⁠ at ⁠).

Fluctuations of the water state in the system, expectedly, range from the states at minimal operation level to the corresponding state at maximal operation level. It is more likely that stage decisions keep the water level in the reservoir around the normal water level which enables the water resource system to mitigate and attenuate flood hydrograph.

The SD simulation model with a causal loop is developed to mimic the behavior of the water resources system. It is characterized by the following blocks (Ahmad & Simonovic 2000; Stojkovic & Simonovic 2019; Ignjatović et al. 2021): stocks, flows, and connectors. Stocks represent the level of variable balancing among flows, for instance, the inflows and outflows from the water resources system. Variables are used to analytically define the functions needed for system reservoir operations connectors, while connectors direct the variables, inflows, and levels enabling the SD model for real operational usage.

The SD model incorporates causal loops providing a link between the system characteristics, inputs, and operational rules separately defined for flood protection facilities (spillway and bottom outflow) and hydropower plants (Stojkovic & Simonovic 2019). System characteristic describes the main parameters of the water system such as water storage curve, spillway, bottom outflow hydraulic characteristics, and hydropower plant parameters. The SD model incorporates an additional model structure defined as functional indicators since they implicitly describe the reduced capacity of the water system elements caused by natural or man-made hazardous events (Ignjatović et al. 2021).

The SD model is developed in this research to mimic a multipurpose water resources system highlighting the flood management reservoir roles. It releases outflows from the spillway (⁠⁠) and the bottom outflow facility (⁠⁠). Additional water management roles are also incorporated within the SD model allowing for the control of the discharge to suit the demand for hydropower generation (⁠⁠) and ecological flow improvements during low-flow periods (⁠⁠). In addition, the losses from the water resources system are taken into account as leakage from the dam (⁠⁠).

A simulated earthquake coinciding with high-flow river conditions and variable initial state in the water resources system influences the system elements by reducing the system's capacity to respond and withstand during the hazardous event. The loss of reservoir capacity is defined as the loss of the volume in the water system needed to receive the flood hydrograph. On the contrary, the degree of flood peak attenuation is the reservoir's capacity to respond to the flood event.

Flood dynamic resilience represents a time-dependent risk metric (Equation (7)) controlled by the capacity of the spillway to evacuate an excess water volume from the water system maintaining the system below acceptable water levels, for instance, maximal operational level. However, the transformation of flood hydrographs within the system is demanded since it coincides with a reduced capacity of the spillway gates caused by the earthquake. Moreover, the variable initial water state can increase the flood-related risk as the initial water level exceeds the expected water level in the reservoir.

Taking into account the fact that the flood dynamic reliance depends on hazards with variable magnitude, time of occurrence, and duration, there is a need to generalize the assessment of dynamic resilience using an ANN. Application of the ANN provides information on the dynamic resilience in a timely manner and enables the prediction needed to perform adequate strategies aimed at reducing the adverse effect on the environment caused by external hazardous inputs. The previously developed ANN algorithm (Beale et al. 2020) is adopted to extract the information from the flood dynamic resilience of the system (Equation (7)) expressed as the robustness and rapidity, using randomly generated hazardous events with different impacts on the system elements (Equations (1), (3), and (4)).

The primary purpose of the Pirot water system is the mitigation of floods at the Nišava river, hydropower generation, and downstream water quality control by regulation of the outflows from the reservoir over the low-flow seasons. The management of the Pirot water system depends on the actual volume of water stored in the reservoirs, inflows, and energy demand (Ignjatović et al. 2021). The active storage of the Zavoj reservoir is 140 × 10⁶ m³ (Table 2). Three gated spillways are located at the left part of the dam with a capacity of 1,820 m³/s. The power plant has two turbines (40 MW) for power generation, with an installed capacity of 45 m³/s (Table 2). The Zavoj reservoir is hydraulically connected with diversion-type turbines at the HPP Pirot by the 9 km pressure tunnel with a 4.5 m radius. The pressure tunnel conveys the water from the Visočica river to the Nišava river providing a significant contribution to the total annual flow of the downstream river.

The variability of flood events was simulated by applying a random selection of return periods. Next, the generated earthquake magnitudes were introduced as additional pressure on the system functionalities by reducing the spillway capacity. Depending on the earthquake magnitudes, the decreased capacity of the spillway was simulated by employing the empirical law (Equation (2)). Moreover, a variable initial state (volume) in the Zavoj reservoir was also randomly generated respecting the physical limitations of the reservoir (referent operation levels).

The deterministic hydrological model provided the long-term daily flows for a 1,000-year time frame at the site of the Zavoj reservoir. For the sake of the hydrological simulations, the generated precipitation and air temperature from the K-NN water generator were used. The annual maximum series approach had chosen the most extreme annual flows, and then the GEV distribution was applied to estimate the design annual flows with different return periods (from 100- to 10,000-year return periods). To define the flood hydrographs, the Gumbel distribution was used to shape the flood hydrographs and respect the design of annual flows.

Once the flood hydrographs had been estimated, a generation of flood hydrographs was performed via the MC simulations as the input for the SD model simulations. For the purpose of hydrograph generation, a uniform distribution was applied to form a dataset of 1,000 members (Table 3).

The impact analysis of the earthquake magnitudes on the spillway capacity reduction was carried out considering the range of possible seismic acceleration within the analyzed region (Rakić et al. 2022). Given that, the dependence between the spillway capacity and earthquake magnitude was deduced from the numerical stability model of the Zavoj dam resulting in empirical values of (Equation (2)).

To generate the dataset of 1,000 values, the MC approach (Table 3) was performed using a uniform distribution of the standardized magnitude of the earthquake (⁠⁠). Furthermore, the spatial and temporal coincidence of natural hazards within the upper Nišava river basin was taken into account with respect to flood hydrograph occurrence. Having obtained the values, start time (⁠⁠) and recovery time (⁠⁠), the functionality indicators of the spillway were determined for each generated earthquake following Equation (3).

Initialization of the SD model of the Pirot water resources system was performed using the same number of MC simulations. Namely, each simulation has a statically generated initial water state in the Zavoj reservoir varying between maximal and minimal operational levels (Table 3). Moreover, the most probable initial state in the reservoir corresponds to the normal operating level. For these reasons, the simulation was performed following the triangular distribution whose median matches the normal operational level, while the distribution is limited by two extreme operational levels (maximal and minimal).

System dynamics simulation was used to model interactions between the water resources system inputs, system elements, and outputs. In this research, the SD model of the multipurpose Pirot water resources system represents the following system elements: spillway and bottom outflow facilities, hydropower plant, and leakage from the Zavoj Dam. Spillway (⁠ and bottom outflow releases (⁠ serve to control and reduce the devastating impact of the flood events at the downstream river sections of the Visočica and Nišava rivers. Moreover, bottom outflow releases improve the water quality of the Visočica river during the low-water periods, usually during the summer and autumn months. Depending on the energy market's demands, the operator of the Pirot water system controls the hydropower generation releases (⁠ considering the physical and environmental limitations to supply the energy market over peak demand hours.

It can be seen from Figure 6 and Table 4 that the most severe flood event was generated in simulation 1. The hydrograph peak reached the value of 6,183 m³/s (Table 4), while the generated earthquake magnitude caused a relatively moderate decrease in the spillway functionality indicator (Figure 6(b)). The initial state in the reservoir, at the beginning of simulation 1, corresponds to the lowest water level respecting the result from all simulations (Figure 6(d) and 6(c)). Considering the generated inputs for simulation 1, especially in the terms of the extremely high flood hydrographs, the water levels in the Zavoj reservoir and spillway outflows reached the highest possible values (Figure 6(d) and 6(e)). Within this simulation, the bottom outflow had an important role since it assists in the discharging of the reservoir, alongside the spillway. Note that the total capacity of the bottom outflow (80 m³/s) is substantially lower than the spillway capacity (1,820 m³/s). Over the highest reached water level in the Zavoj reservoir, the bottom outflow facility released outflows up to its total capacity (Figure 6(e)).

For the case of simulation 4, the water levels in the Zavoj reservoir also achieved a maximal possible level. Such behavior of the Pirot corresponded to the highest value of the generated earthquake magnitude equal to 0.948 (Table 4), and consequently a hugely decreased capacity of the spillway (Figure 6(b)). However, a moderate peak of the generated hydrograph (Table 4) did not significantly affect the changes in the water reservoir levels.

The rest of the simulations (simulations 2, 3, and 5) produced a conservative system behavior since the increase of the water levels in the Zavoj reservoir was moderate (Figure 6(d)). This response of the Pirot water resources system is mainly attributed to the high capacity of the spillway (up to 1,820 m³/s) and corresponding values of spillway functionality indicators to discharge the reservoir during the flood events with the peaks in the range from 706 to 1,201 m³/s (Table 4). The most indicative are simulations 2 and 5 where the generated earthquake magnitudes did not cause any reduction of the spillway capacity, considering Equation (2). Therefore, the functionality indicators over the whole simulation time kept maximal values enabling the efficient transformation of flood hydrographs to downstream river sections (Figure 6(b)).

Flood dynamic resilience captures the Pirot water resources system capacity over simulation time to respond to and withstand hazardous events. The main properties of dynamic resilience are robustness and rapidity. Robustness is defined as the minimal value of dynamic resilience reached at the time step (Figure 2). The difference between the recovery time (⁠⁠) and the time when a hazardous event occurs is the system's rapidity.

A system performance level of the Pirot water resources system is defined by using Equation (8). When the water levels in the Zavoj reservoir are below the maximal operation levels, the system performance is equal to 1. By exceeding the predefined level, the performance level decreases to zero value. Following this, the flood dynamic resilience (Equation (7)) is quantified for each simulation time step with respect to the hazardous start time (⁠⁠).

The results from Figure 7 suggest that the Pirot water system was less resilient (⁠0.04) in the case of simulations 1 and 4. The former is attributed to the extraordinary values of the flood hydrograph resulting in the prolonged high water levels in the Zavoj reservoir. Although the latter is not characterized by the extremely high peak of the flood hydrograph, the reduced spillway capacity produced a low flood dynamic resilience. The important difference between these simulations lies in the fact that the dynamic resilience rapidly increased for simulation 4 since the duration of high-water levels was less when compared to simulation 1 (Figure 6(d)). Moreover, the Pirot water resources system implied a low resistance level for simulation 3 (0.143). For this simulation, the low flood resilience level is attributed to a substantial decrease in the spillway capacity (Table 4).

On the contrary, the most resilient system behavior was exhibited in simulation 5, maintaining the resilience value at the maximal level (Figure 7). Despite relatively high water levels at the beginning of the simulation (Figure 6(c)), the Pirot water resources system was not disturbed by significant hazardous events.

In the case of simulation 2, the high water levels in the Zavoj reservoir led to a rapid decrease of the flood dynamic resilience equal to 0.50 (Figure 7). Fortunately, the Pirot water resources system coped with the unfavorable starting conditions and promptly reached the full system performance (Figure 6(c)).

For the most extreme hazard combination, the flood dynamic resilience decreased roughly to zero value. A decrease in robustness (with a minimum value of 0.03 given in Figure 8(a)) is brought about by combinations of earthquakes with significant magnitude and flood hydrographs with a low frequency of occurrence. To retrieve the system functionality under such severe hazards, the rapidity of the Pirot water resources system reached 1,000 h (Figure 8(c)). In the case of the moderate hazards, the system robustness had a median of 0.2 (Figure 8(a)) while rapidity was equal to 162 h (Figure 8(b)).

To examine the flood dynamic resilience in more detail, ANN was employed to model dependencies between hazardous events with different impacts on the system elements and flood-risk metrics (robustness and rapidity). For this purpose, a feedforward neural network was used, with three input nodes, three hidden layers of four nodes each, and two output nodes. The network was fully connected, and the hidden, as well as the output layer activation functions, were all sigmoid. The network was trained using the LM algorithm (with 1,000 epochs). The ANN's inputs were: (1) the set of maximal values of generated flood hydrographs (⁠⁠); (2) the reduction of the functional indicators (⁠⁠), and (3) the initial state (⁠⁠) in the Zavoj reservoir. The dataset was split into three subsets: the training set, validation set, and test set using a 70%–15%–15% split, respectively. Additionally, the data for these subsets were selected randomly to preserve generality. The network was trained 200 times and the best of those networks was selected to avoid the problem of local minimum which can adversely affect the network's further application. The selection of the best network was based on the Root Mean Square Error (RMSE) metric. Note that when selecting the best network, robustness and rapidity were valued with equal importance.

The results imply that the ANN accurately reproduces key dynamic resilience parameters (Figure 8(b) and 8(d)) compared to their simulated values via the SD model (Figure 8(a) and 8(c)). For instance, the first and third quartiles on the boxplot diagrams expose similar values for both modeling approaches (SD and ANN models). However, slight discrepancies are indicated in the case of the second quartile (median) equal to +3.3% and −8.1% for robustness and rapidity, respectively (Figure 8). Maximal values of those parameters indicate that the ANN model slightly underestimates the extreme values derived by the SD approach (Figure 8).

The metric used to quantify the efficacy of the ANN approximation was the relative error; the average relative error was equal to 2.14% for robustness and 1.77% for rapidity. Moreover, the RMSE values for robustness and rapidity were equal to 0.004 (–) and 12.7 h, respectively. Considering the timespan of rapidity, such an error is not substantial.

Robustness and rapidity diverge from the corresponding observed parameters for higher values of cumulative distribution function (e.g., higher than 0.95), where robustness diverges more, the quantiles of the approximation distribution are more aligned with the observed parameters in the case of rapidity (Figure 9).

However, the results (Figures 8 and 9) indicate that ANNs are capable of adequately approximating the wide range of the robustness and rapidity values needed for real-time operations of the Pirot water resources system. In addition, results indicate several accessible improvements: (1) simulations of the models are performed with the generated hydrographs instead of using a limited number of the observed flood events, (2) hazardous models envelop the most important system disturbances (floods and earthquakes) while other hazards and are not considered (e.g. failures of monitoring, electrical, or mechanical equipment), (3) sigmoid activation function may not catch the extremely high values of dynamic resilience-related parameters, and application of the ReLU activation function can improve prediction of the flood-related parameters (Bui et al. 2020).

Natural hazards encompass several physical phenomena such as floods, severe storms, earthquakes, landslides, and tornadoes (Gill & Malamud 2014). The southwest European region is sensitive to natural hazards since floods and earthquakes considerably exceeded the largest foreseen events and caused billions in damage (Stojković et al. 2017; Maršanić et al. 2021).

Therefore, this research is focused on the impacts of floods and earthquakes on the functionalities of the Pirot water system (located in Serbia, Southeast Europe) that prevent and reduce loss and damage due to extreme events. To examine the combination of a wide range of less and more severe hazards, we proposed a hazards modeling approach founded upon MC simulations that enable the simulation of hazardous events with low joint probability. Hazardous modeling reduces the uncertainty caused by short recording periods (Koutsoyiannis & Montanari 2007) using the stochastic weather generator creating long climate data (1,000 years span) needed for hydrological modeling (King et al. 2015). Hazardous modeling also encompasses the outputs from the Zavoj dam safety model to precisely deduce the dependence between the earthquake magnitude and functionality of the spillway (Rakić et al. 2022). In addition, MC simulations are used to derive variable water content in the water system since its initial condition significantly impacts on prevention of uncontrolled and rapid release from the system (Kistenmacher & Georgakakos 2015).

The behavior of the Pirot water resources system is simulated under a system dynamics approach enabling the system to respond to the generated hazardous events. A system approach deals with the difficulties introduced by the complexity of the water resources solutions and environmental issues (Simonovic 2020). It is especially highlighted when the system fails to execute the service requirements because of the reduced capacity of the flood-defense facilities under hazardous events. To measure the flood-related risk of the water system, dynamic resilience is used as a time-dependent parameter since it overperforms the static risks measures such as system reliability or vulnerability (Ignjatović et al. 2021). Based on the output from system dynamics simulations of the Pirot water system, the dynamic resilience captures the most important characteristics (Simonovic & Arunkumar 2016), specifically system robustness and rapidity (Figure 2). Rapidity describes the recovery time of the water resources system forced by an external or internal hazard (e.g., floods, failure of electrical or mechanical system facilities) (Bowles et al. 2013). A maximal reduced capacity of the water system in terms of delivering the service requirements is explained by robustness.

Considering an increase in available hydroclimatic data, alongside more computing power, machine learning algorithms become a useful method for flood risk and impact assessments (Wagenaar et al. 2020). In the presented research, flood-related risk characteristics (robustness and rapidity) were captured by postprocessing a huge amount of data generated via the system dynamics model. Then, the generated data, including the robustness and rapidity, was utilized for the supervised training of the ANN (Beale et al. 2020).

ANN was used to make a model for predicting robustness and rapidity based on characteristic values related to floods and earthquakes, and in the network training procedure, large data obtained from simulations of system dynamics models are used.

As mentioned before, the ANN was used to approximate flood-related risk characteristics, through the system modeling simulations. The results suggest that it is possible to train a network to approximate the rapidity and robustness of a system. Note that rapidity highly matches the target values since it continuously accumulates the flood-related risk characteristics over the simulation time. While the real system is more susceptible to measurement errors (e.g., measured water levels in reservoirs), the rapidity timescale allows for such errors to occur without significantly affecting resilience predictions. In contrast, the network cannot capture robustness with the same efficiency as rapidity because it changes rapidly after a disturbance is introduced.

There is an interest to consider the application of ANNs on a real system using actual data (Kostić et al. 2016; Allawi et al. 2018; Allawi et al. 2019; Jenny et al. 2020; Osman et al. 2020). However, its application is possible for the observed period with available records (inflows and water levels in the reservoir), while the usage of the developed system dynamics model cannot be avoided for predicational purposes where the inputs to the system and its behavior are unknown. Moreover, such an approximation is useful for managing water resources systems, as it can give an idea of when the system will return to a predisturbance level of functionality, helping operators to adequately allocate or reallocate resources to other systems, and improving overall functionality.

This paper proposes a novel framework for the evaluation of flood dynamic resilience for a water resources system under hazardous events, where the link between the hazardous events and the parameters used to evaluate flood dynamic resilience is defined by trained ANNs.

The proposed framework is tested on the Pirot water resources system located in Serbia. The framework is based on an SD modeling approach to mimic the non-linear behavior of the water resources system. The first of the hazardous models is a flood hydrograph model that is used to generate data sets using MC simulations based on uniform distribution. The second hazardous model is the earthquake model, which utilizes a functionality indicator (α) that is simulated using MC simulation in the range 0 to 1; this functionality indicator is derived from a simplified analytical form of ESWC values of the water resources system. The third hazardous model is the initial state model which uses a triangular distribution as input for MC simulations. Using the outputs from the SP modeling approach, flood dynamic resilience is estimated to capture several significant numerical characteristics (rapidity and robustness).

The framework relies on ANN which uses characteristic values of the hazardous models' outputs as inputs into the network to approximate flood dynamic resilience, exactly rapidity, and robustness. Quantifying the efficacy of the ANN approximation using the relative error metric shows an average relative error of 2.14% for robustness and 1.77% for rapidity. This is indicative of the fact that the ANN is properly trained. It is important to note that the robustness and rapidity diverge from the corresponding observed parameters for higher values of the cumulative distribution function; the robustness diverges more, while the quantiles of the approximated cumulative distribution are more aligned with the observed parameters in the case of rapidity. These diverging extreme values are to be expected given the sigmoid activation function, which may not catch the extreme values of flood dynamic resilience-related parameters. These results indicate that the ANNs are capable of adequately approximating the robustness and rapidity values caused by the hazards in the Pirot water resources system.

Future research may include several possible directions. For instance, the sigmoid activation function is used as an important part of the ANNs for flood dynamic resilience prediction. An improvement in the resilience prediction could be achieved by using the Rectified Linear Unit (ReLU) activation function instead of the sigmoid function. It helps to avoid the vanishing gradient phenomenon which can significantly hamper the learning progress (Goodfellow et al. 2016).

The authors express their gratitude to the Science Fund of the Republic of Serbia for the support through the project of the PROMIS call, 6062556, DyRes_System: ‘Dynamics resilience as a measure for risk assessment of the complex water, infrastructure and ecological systems: Making a context’.

Data cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.