Evaluation of the vulnerability in water distribution systems through targeted attacks

Water is a scarce and indispensable resource for maintaining human life, and access to it has been undertaken as a sustainable development goal by the United Nations. To accomplish this goal, a water distribution system (WDS) must be extended and optimized to reliably deliver the required water quantity and quality to a larger percentage of the human population. Furthermore, these systems are fundamental to food security through crop irrigation and health through sanitary services (United Nations 2020). WDSs are interdependent and complex structures formed by water sources, treatment plants, pipes, junctions, tanks, reservoirs, pumps, and valves, among other elements, designed to satisfy the water demand at consumption points (Christie et al. 2011). Nevertheless, these structures might be vulnerable in front of natural disasters and other threats such as droughts and ill-intended attacks.

Furthermore, climate change has generated extreme weather conditions that affect water reservoirs and other elements of the WDS. Therefore, it is important to research how to guarantee the stability and continuity of the water supply by understanding the vulnerability of the WDS when they face ‘attacks’ which in this work are defined as the outages of a set of elements that might be caused by different natural or artificial conditions.

Mathematical models are fundamental in researching the vulnerability of the WDS, and typically a WDS model consists of establishing a relationship between the variables that intervene in the water flow problem (pressures and volume of water flow) and the parameters of each element (Cabrera & García-Serra 1999). The complex network theory provides an invaluable framework for the representation and quantitative study of the WDS, through analysis of subsystem connections, quantifying topological redundancy, and identifying critical components whose failure may impact the overall performance of the system (Barabási & Albert 1999; Yazdani & Jeffrey 2012; Barabási 2013; Albarakati et al. 2017; Albarakati & Bikdash 2021; Albarakati 2021a, 2021b).

In recent years, water distribution networks have been studied as complex networks (Murthy & Murthy 2012; Alexiou & Tsouros 2017). In Giudicianni et al. (2018), the authors analyzed the topology of several water distribution networks, concluding that they tend to be sparse and resemble regular square grids with homogenous node degree values. Some vulnerability studies regarding the WDS have been carried out considering topological features of such networks, and additional metrics have been proposed to identify critical nodes (Yazdani & Jeffrey 2012). The relationship of the betweenness centrality of the network's nodes and the vulnerability of a WDS was demonstrated by Agathokleous et al. employing a study case under continuous and intermittent conditions (Agathokleous et al. 2017). In Gutiérrez-Pérez et al. (2013), a methodology based on spectral measurements was introduced to establish the importance of areas in the WDS, the authors employed two ranking algorithms: PageRank and HITS.

More recently, a methodology for calculating an index of vulnerability that represents the tendency of an injected contaminant to spread over the network was introduced, this index compares the behavior of the network to a scale-free network and allows different networks to be ranked in terms of their resistance to contamination (Nicolini 2020). Mortula et al. conducted an integrated network and hydraulic analysis for identifying critical nodes that can be impacted by water shortages and proposing managerial action to reduce network vulnerability, diversify water sources, and reduce possibilities of contamination (Mortula et al. 2020). The adjacency matrix from the WDS network has been employed to detect unintended isolation segments that can cause vulnerability, helping to determine the best location of valves for avoiding the interruption of the service for pipe failures (Jeong et al. 2021). Also, the nodal vulnerability of water distribution networks has been evaluated under cascading failures through a study case, where there was a high probability of the nodal failures triggering a cascading failure (Shuang et al. 2014).

Wéber et al. explored the behavior and topology of the segments resulting from an accidental pipe burst, especially the criticality of such segments from the point of view of the whole system. They defined the network vulnerability as a measure to quantify the expected value of relative consumption that cannot be served due to a single accidental pipe break (Wéber et al. 2020). In Giudicianni et al. (2021), a novel vulnerability index called the cut-vulnerability is introduced and several networks with different topological structures were analyzed under the proposed framework, while the vulnerability measure is based on the loss of connectivity.

Those research works represent important advances in the application of complex network techniques and centrality measures to the study of the WDS. However, such metrics need to be further investigated to understand their reliability by using benchmarks to establish appropriate ways to measure the vulnerability of water distribution networks. In this sense, this research conducts a vulnerability analysis in different WDS benchmarks, considering a graph model that integrates topological features and hydraulic characteristics to compare different centrality measures in terms of their ability to predict the shortage of water supply using the graph model. This framework is different from that proposed in the literature because it is based on the topology and the operating condition of the system, helping to predict the vulnerability under different operating conditions.

The main components of the WDS are water sources, treatment plants, pipes, junctions, tanks, reservoirs, pumps, and valves. For modeling the network, these elements can be represented by nodes and links: nodes are the points of the system where there can be an input or output of water, and links convey water from one node to another.

Each node of the model is associated with an elevation and there are different types of nodes: for example, reservoir nodes where the head is forced and water is provided to the system, consumption nodes where there is a determined water demand, and the junction nodes that are mainly connections between pipes. Links in the network model represent pipes, which are the most abundant components of a WDS, allowing the water to move from one point to another of the system. In these links, it is important to consider the friction losses causing energy to dissipate across them.

Since these equations are nonlinear, a numeric solver is required to obtain a fast solution. In this research, the water flow is determined with the EPANET-MATLAB Toolkit under quasi-static conditions (Eliades et al. 2016).

Centrality metrics are the most widely employed indicators for identifying important nodes in a network, and depending on the type of network and the study case, the importance might have a different meaning (Agathokleous et al. 2017). In this research, the centrality measures the importance in terms of the water flow that runs through a node or a pipe.

In previous works, different centrality measures have been evaluated for similar applications in electric power systems, and the most consistent were degree, eigenvector, and Katz centralities based on power traffic, which can be considered an equivalent of water flow. Furthermore, betweenness, closeness, and PageRank centralities have been considered relevant in current computational applications. Betweenness, closeness, and degree centralities have been cataloged as suitable metrics for studying centrality in spatial networks (Giustolisi et al. 2019).

The complexity of calculating centrality metrics depends on the network size, connectivity, and the weights associated. In this work, the mentioned centralities are compared in terms of their efficiency for assessing the vulnerability of the WDS. Most of the centralities are calculated using the weighted graph model previously described as detailed in the following subsections. For a better comparison, the centralities are normalized in such a way that the sum of the measures for all the nodes or links is equal to one.

Degree centrality is a popular centrality measure that gives some insight into the connectivity of a node⁠; however, it can miss potentially important aspects of the network's architecture and the position of the node within it (Bloch & Jackson 2021).

The Freeman's betweenness centrality of node i measures the relevance of such a node in terms of its role connecting the rest of the nodes in the network. Betweenness is meant to capture the importance of the node as an intermediary in the association between other nodes, in the case of the WDS, in the transportation of water among the rest of the nodes.

It is important to note that the closeness works for fully connected networks; otherwise, the distance between two nodes belonging to different components is infinite, and then, the metric of Equation (7) is null, and harmonic closeness (Equation (8)) is employed to avoid this problem (Giustolisi et al. 2019).

In this work, the PageRank centrality calculation is based on the MATLAB© algorithm, which is the result of a random walk of the network. At each node in the graph, the next node is chosen with a probability (by default 0.85) from the set of successors of the current node: if the node has no successors, the next node is chosen from all the nodes. The centrality score is the average time spent at each node during the random walk (The Mathworks Inc. 2020).

This measure has been previously applied to the WDS (Giudicianni et al. 2020) and in social networks (Riquelme et al. 2018). In some large-size network cases, it did not converge.

To exemplify the calculation of the centrality measures, Figure 1 shows an example network with five nodes and its weighted adjacency matrix ⁠, and this graph models a WDS example obtained from the EPANET-MATLAB Toolkit (EPANET-MSX Example Network) (Eliades et al. 2016). Table 1 contains all the centrality measures calculated as described in the previous subsections, and the same functions are applied through all the research work to determine the centrality measures.

The previous section was based on node (junction) centralities; however in this work, the link (pipe) centralities are also evaluated. To obtain the pipe centralities, the weighted adjacency matrix of graph G(v, e) is converted to an equivalent line graph adjacency matrix by applying the transformation to the line graph (Yoshida 2013). is a matrix with e equal to the number of pipes in the grid, the elements are the line weights, which represent a water flow from pipe i to pipe j.

After converting the graph to a line graph, the centrality measures are calculated in the same way that they are calculated for a node graph, with as the adjacency matrix.

Figure 2 shows the results of transforming the example network graph from Figure 1 to a line graph where the five pipelines are mapped into nodes, and is the new adjacency matrix of the line graph.

Table 2 shows the results of the centrality measures calculated for the equivalent line graph.

The vulnerability assessment framework is based on previous works from the authors focused on electric power systems' vulnerability implementing individually targeted attacks and sequential attack profiles on nodes and links of the corresponding graph model. An attack profile is a sequence of elements (nodes or links) selected from the network for their consecutive removal.

The framework proposes using the vulnerability curve, which represents a physical damage measure in the horizontal axis and a functional damage measure in the vertical axis. This graph is helpful to determine the amount of functional damage caused by a fault profile.

The vulnerability prediction measure (VPM) describes the amount of damage caused in terms of the physical damage measure, and it is calculated with the area under the vulnerability curve. In this sense, the higher the VPM, it is said that the attack sequence causes more severe electrical damage to the power system; thus, the power system is more vulnerable to such attack.

For calculating the unsatisfied demand in each stage, the WDS modeled in EPANET is converted to a MATLAB graph object (The Mathworks Inc. 2020), and then is calculated and recorded for each stage depending on the type of attack, as it is described in the following subsection.

Figure 3(a) shows the simplified flow diagram for the RMCEF algorithm, and in this type of attack, the elements are removed from the network starting from the highest to the lowest centrality one element per stage.

The IMCEF is an attack profile that also relies on centrality measures, but in this type of attack, the centrality measures are recalculated after the removal of an element. The idea under this attack profile is that the centrality measures change after the removal of the most central element; thus, the second most central element in the initial ranking may not be the most central, once the most central element is removed. Then, the vector of centralities must be recomputed to obtain the new ranking of elements according to centrality (Figure 3(b)).

Furthermore, the IMDEF fault profile is defined in this work as an attack based on removing the element that generates the largest raise in the of the WDS in the current stage. The simplified algorithm for obtaining the vulnerability curve using the IMDEF fault profile is shown in Figure 3(c). This algorithm has an internal cycle to determine the unsatisfied demand of each element of the network under current conditions; afterward, it selects the most damaging element (MDE) to remove it from the network.

To compare the different centrality measures and attack strategies, four benchmarks were selected from previous research and reports. The benchmark models are called: the New York Tunnel (NYT), the Hanoi (HAN), the Modena (MOD), and the Balerma (BAL) networks, the former two are considered medium-size networks and the latter two are large-size networks.

Table 3 shows the characteristics of each benchmark selected for this research, from the smallest to the largest based on the data from Centre for Water Systems University of Exeter (2014). The NYT network was originally presented as an extended design problem (Schaake & Lai 1969); the HAN WDS consists of a structure of three loops and one reservoir with a fixed head of 100 m (Kim et al. 1994); the MOD WDS has both minimum and maximum pressure requirements for demand nodes with pipes made of cast iron (Bragalli et al. 2008); the BAL system is an adaptation of an existing irrigation network in the Sol-Poniente irrigation district, located in Balerma in the province of Almeria (Spain), and it has a fixed water consumption across all the demand nodes (Reca & Martínez 2006).

The selection of benchmarks was based on having diverse characteristics, i.e., medium and large systems with different sizes of pipes and materials. The topologies of the WDS benchmarks are shown in Figure 4, where they are meshed networks, which allow the redirection of the power flow when the links are broken.

This section presents the results of the experiments performed using the vulnerability evaluation framework for the benchmarks selected. The experiments are based on fault profiles for junctions and pipes, using the centrality measures described in Section 3. The results are analyzed in terms of the vulnerability curve, the VPM, and the execution time for each centrality measure to detect the centrality that provides the best prediction for the most damaging attack profile. The calculation of the water flows was made assuming an initial condition, and the variability of the demand is not taken into account. This operating condition can be easily changed to predict the VPM under future circumstances.

For the attack profiles performed in the junctions of the WDS, the histogram of the centrality measures calculated in the initial operative point is shown in Figure 5. It is noted that degree, eigenvector, and Katz centralities present a similar behavior, where most of the junctions have a small centrality but there are a few that have larger centralities, while for the rest of the centralities there is not a large deviation among all the junctions, which have small and similar centrality measures.

Figure 6 presents three examples of the vulnerability curve for the MOD benchmark, where is it noted that the IMDEF and the IMCEF attack strategies have similar behavior with a steep curve that rises to 100% of unsatisfied water demand after removing less than 20% of the junctions. Meanwhile, the RMCEF strategy for this network had a less steep curve, only reaching 100% of unsatisfied demand after removing more than 80% of the junctions.

Results from VPM and execution times are shown in Table 4 for each experiment using the different fault strategies. It is noted that the IMDEF strategy presented the highest values of VPM in most of the benchmarks, which is expected since this algorithm finds the nodes of the network with the highest impact in the unsatisfied demand. In terms of the strategies based on centrality measures, degree, eigenvector, and Katz centralities also presented the highest VPM values. Furthermore, it is evident in Figure 7 that the VPM for the experiments with degree and Katz centralities were always close to the IMDEF VPM.

Regarding the execution times of each strategy, it is noted that the IMDEF had the highest times, which is expected given the quadratic nature of the algorithm. The RMCEF and the IMCEF presented similar execution times for the medium-size problems (NYT and HAN), but for the larger systems, the IMCEF was more effective. Although the IMCEF algorithm recalculated the centrality measures in each iteration, it reaches the maximum unsatisfied water demand much faster. The shortest execution times were obtained with the degree centrality.

The three attack strategies (IMDEF, RMCEF, and IMCEF) were also applied to pipes, by employing the transformation to the line graph detailed in Section 4. The histograms of the pipe centrality measures for the benchmarks is shown in Figure 8, where it can be noted that the behavior is similar to the junction centralities, having a larger deviation between values of degree, eigenvector, and Katz centralities than for the rest.

Examples of vulnerability curves are shown in Figure 9 for the MOD benchmark, where it is noted that the steepest curve was obtained for the IMCEF attack strategy, which reaches 100% of unsatisfied water demand before removing 20% of the pipes. This behavior differs from the attacks on junctions, since this time the IMDEF attack strategy does not guarantee that the network reaches the total unsatisfied demand faster, which is because different pipes cause the same damage, and the algorithm chooses one of them without a particular criterion that allows predicting the overall most damaging path.

Examining the VPM values for each experiment, it is noted that the IMCEF attack strategy overall presented equal or higher VPM values than the rest of the strategies for the same WDS. Furthermore, the RMCEF strategy in all the studied cases presented the second highest VPM. The degree, eigenvector, and Katz centralities presented the highest VPM in most of the cases in the same way as for the junction attacks (Table 5). The IMCEF strategy allows high VPM values to be obtained, so it predicts the most damaging attacks on pipes and has the lowest execution times for the largest benchmarks (MOD and BAL).

Analyzing the results obtained for attack profiles on junctions and pipes of the benchmarks, it is noted that the IMDEF strategy worked better for attacks on junctions by predicting the most damaging attacks; however, it presented high execution times and it did not predict the highest VPM for the attacks on junctions. Therefore, the IMCEF strategy was more efficient for attacks in pipes and junctions having high VPM and lower execution times for the larger WDS.

In terms of the centrality measures, the behavior was similar to that obtained for electric power systems in previous works, where the degree, eigenvector, and Katz centralities are representative of the elements that generate the highest impact in the unsatisfied demand. These centralities allow the identification of important elements in the grid considering the information on the calculated water flow.

These experiments assume that the junctions and pipelines can be disconnected with isolation valves, but in practical scenarios, the disconnection will occur through segment networks. That condition can be taken into account for extending the work to more realistic conditions and evaluating the vulnerability in the segment networks using models with data from isolation valves.

The main contribution of this work is to apply the water flows to the calculation of centrality measures for evaluating the vulnerability in different WDS benchmarks. The vulnerability framework is based on a graph model integrating the topological structure of the network and the operational water flow calculated with EPANET and MATLAB.

The vulnerability framework determines the vulnerability curve and the VPM according to an attack profile based on three different strategies. The RMCEF consists of removing the elements according to the centrality measured calculated in the initial operating point, the IMCEF removes the most central element after recalculating the centrality for each stage, and the IMDEF strategy consists of removing the MDE first. These strategies were applied to junctions and pipes of four different WDS benchmarks: the NYT, the HAN, the MOD, and the BAL networks.

After analyzing the results of all the experiments performed, it is concluded that the IMCEF had the best behavior for evaluating the most damaging attacks on the WDS, since it presented high VPM values both for removing junctions and pipes. Moreover, it presented low execution times for the largest benchmarks. This indicates that the recalculation of the centrality measures after removing the most central elements helps to identify new important elements and the initial centrality measures lose relevance after such removal.

Another important contribution is identifying the centralities that are more predictive in terms of the MDEs, which were the degree, eigenvector, and Katz centralities under the conditions of this research. The degree centrality additionally presents the advantage that its calculation is simpler; thus, it is recommended for the larger WDS.

It was demonstrated that the vulnerability framework can be applied to the WDS in the same way in which it was previously applied to electric power systems. This work can be extended to consider the variability of the demand, allowing us to obtain the vulnerability under different operating conditions and the segment graph. Future work will be oriented to the design of the WDS by optimization by minimizing the vulnerability of the system under attacks that can be generated by droughts and other extreme weather conditions.

The authors extend their appreciation to the deputyship for the Research & Innovation Ministry of Education in Saudi Arabia for funding this research work through the project number (IFP-202O.13).

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