# Dynamical stability of water distribution networks

## Abstract

Water distribution networks are hydraulic infrastructures that aim to meeting water demands at their various nodes. Water flows through pipes in the network create nonlinear dynamics on networks. A desirable feature of water distribution networks is high resistance to failures and other shocks given to the system. Such threats would at least transiently change the flow rate in various pipes, potentially mitigating the functionality of the whole water distribution system. Here we carry out a linear stability analysis for a nonlinear dynamical system representing the flow rate through pipes that are interconnected through an arbitrary pipe network with reservoirs and consumer nodes. We show that the steady state is always locally stable and develop a method to calculate the eigenvalue that corresponds to the mode that decays the most slowly towards the equilibrium, which we use as an index for resilience of the system. We show that the proposed index is positively correlated with the recovery rate of the pipe network, which was derived from a realistic and industrially popular simulator. The present analytical framework is expected to be useful for deploying tools from nonlinear dynamics and network analysis to designing, resilience managements and scenario testings of water distribution networks.

## 1 Introduction

Water distribution systems aim to securely provide drinking water to consumers and for firefighting and hence are key infrastructures in our society. They are subjected to drastic changes in demands on the daily and seasonal timescales [1] and to various threats such as contamination [2] and power outage [3]. For safe and secure water supply, it is important that a water distribution system is well designed and operated against potential failures.

In risk management in water engineering, the risk is estimated to be equal to the product of the probability of failures and the consequence of the failure [4, 5]. In contrast to risk management, resilience management is a relatively new framework to aid us to prepare for potential failures in water distribution systems. Resilience is a broad concept and its definition varies according to studies in engineering [6] and in other fields [7]. The concept of resilience is equally wide for water distribution systems; it is construed as the general capacity of a system to resist, absorb, withstand and recover from stresses [6, 8]. It is a useful complement, as many failures that happen in real-life are unforeseeable and are not targeted in the traditional risk management due to the low estimated probability.

Water distribution systems are composed of pipes combined with other functional structures such as reservoirs, pumps and valves. Network disfunctions such as pipe failure are an obvious threat to the functionality of the water distribution system. Network science has accumulated knowledge on how network structure shapes the network’s capacity to withhold random failures and intentional attacks [9, 10]. Therefore, there has been a surge of interest in water engineering community to deploy network analysis to assess the degree of resilience of water distribution networks. However, to the best of our knowledge, most of the existing proposals of resilience measures for water distribution networks are based on the network structure [11, 12, 13, 14, 15, 16, 17, 18], or steady-state or quasi-steady-state flow rates [19, 20, 21, 22, 23, 24]; as such, the dynamic, transient water flows in the pipe network cannot be represented, which is in fact the key to the resilience performance of a water distribution system. Furthermore, one needs to apply stresses to the system to measure resilience [18], which cannot be exhaustive of what failure might occur in a system. Therefore, it is necessary to examine transient dynamics to reveal much of resilience features of the system, in addition to static network structure or flow. In fact, transient dynamics of water flow in pipe networks are not a new issue in water engineering. Examination of transient flows is also a practical requirement because live water distribution networks are incessantly undergoing changes due to, for example, routine operational adjustments, breakdowns of their elements and human error. In turn, transient water flows are known to cause various problems in water distribution systems such as pump failure, collapse of vapor cavities and compromised water quality. Motivated by these needs, various methods to simulate and understand transient flows have been developed [25, 26, 27, 28, 29]. While useful, these developments are mathematical modelling of the transient water dynamics and their numerical simulation methods. These methods do not themselves provide a measure of resilience.

In the present study, we propose an index of resilience via linear stability analysis of a standard set of dynamical equations modelling the steady-state flow in the given water distribution network. The equations are equivalent to a type of a nonlinear electric circuit system. The proposed resilience index is defined in terms of the dominant eigenvalue of the Jacobian matrix around the steady state. Therefore, the index corresponds to the speed at which the transient dynamics decay towards the steady state in response to a small perturbation in the flow rate. The index is calculated from the structure of the network, the diameter and length of the individual pipes, and other constraints such as the energy level at the reservoir nodes and water demand at individual consumer nodes. We test the proposed index against different resilience measures that we previously proposed based on realistic numerical simulations of water distribution systems. MATLAB code for calculating the so-called local stability index is available at Github (https://github.com/naokimas/LSI_water_distribution_net).

## 2 Model

Consider a connected and undirected network of pipes. By definition, a pipe connects two nodes, where nodes are junctions, reservoirs and households. A dead-end node that is only connected to a single pipe (hence, the node’s degree is equal to one) corresponds to a consumer node such as a household. It should be noted that nodes whose degree is larger than one are junction nodes but can also demand water if they are connected to households. Let us denote by and the number of nodes and edges, respectively. Each edge has its own diameter and length, assuming a circular shape of the pipe. These properties affect how much water flows through the pipe under a given condition, thus effectively specifying the weight (i.e., capacity) of the edge. However, we consider the network as unweighted network and explicitly model the effect of the pipe diameter and length.

We consider a rigid water column model on pipe networks under the assumption that flows in a pipe vary slowly in time [30, 25, 31, 32]. The dynamical equation representing transient flow through the th pipe (), which connects the th and th nodes, is given by

(1) |

Nodes and depend on edge . However, unless we state otherwise, we use and here and in the following text to avoid notational abuse. In Eq. (1), is the time-dependent flow rate through the th pipe, in the direction from the th node to the th node; represents the time; represents the total energy at the th node, which depends on the time; is a constant called the pipe coefficient containing the information on the diameter, length and roughness of pipe ; accounts for the nonlinear head loss across the th pipe; is the input term which comes from a valve or other structure located on the th pipe; is the pipe inertia constant. The dynamical system given by Eq. (1) is equivalent to that of a nonlinear resistor-inductor electrical circuit, where is the current, is the voltage, is the inductance, is a nonlinear resistance and is other nonlinear elements placed on the th edge.

The pipe coefficient is given by

(2) |

where is the dimensionless Darcy-Weisbach friction factor for the th pipe and depends on the flow regime (see below for the formula) and other factors, is the length of the th pipe, is the gravitational acceleration and equal to m/s and is the diameter of the th pipe [26]. More generally, if the nonlinear head loss term is given by , one obtains , where is the cross-sectional area of the th pipe, i.e., [33, 26]. The pipe inertia constant is given by .

There are various formula to approximate the friction factor, , as a function of the pipe diameter, , the roughness coefficient of the pipe (which is determined by the pipe’s material, age and other factors), Reynolds number and so forth. Because the Reynolds number is a function of the flow rate, , the friction factor also depends on . Our index developed in section 3 involves differentiation of the right-hand side of Eq (1) with respect to , which requires the derivative of with respect to . In fact, many formula for can be only applicable to particular flow regimes; the flow regime is specified by the value of the Reynolds number. Because may depend on pipes, different pipes may be in different flow regimes and some pipes may be situated near the boundary of two flow regimes. Therefore, we select the formula proposed in Ref. [34] that is applicable in all flow regimes including laminar and turbulent flows and gives as a continuous function of the Reynolds number (and hence continuous in ). The formula is given by

(3) |

In Eq. (3), the Reynolds number for the th pipe is equal to

(4) |

where ms is the kinematic viscosity of water, assuming the temperature of 20 degree celsius. Furthermore, m is the roughness coefficient of cast iron [35] (also see Ref. [26] for a similar value), and

(5) | ||||

(6) |

## 3 Local stability index

In this section, we carry out the linear stability analysis of the steady state of the nonlinear dynamical system given by Eq. (1). The system turns out to be always linearly stable, and we propose a local stability index as the magnitude of the eigenvalue corresponding to the most slowly decaying mode. This task is non-trivial because some nodes act as reservoir of water, other nodes act as consumer of water and the friction factor, , depends on the dynamical variable, .

### 3.1 Removing redundancy from the dynamical equations

Assume that the total energy at nodes is fixed. The nodes are typically reservoir whose altitude is higher than typical consumer and junction nodes. Without loss of generality, we assume that , , , are fixed. Then, the set of Eq. (1) has unknowns, i.e., () and (), whereas Eq. (1) only provides differential equations. The remaining constraints come from the Kirchhoff’s current law (i.e. conservation of water mass) at each of the nodes whose total energy is not fixed. In other words, we obtain

(7) |

where is the flow rate through pipe from the th to th nodes; is the external demand (i.e., withdrawal) of water at the th node; is the element of the adjacency matrix. In other words, if there is a pipe between the th and th nodes. Otherwise, .

Equations (1) and (7) imply that there are unknowns and the same number of equations for solving the steady state and its linear stability. Now we erase () from Eq. (1) using Eq. (7) to derive a self-contained set of -dimensional dynamical system. It should be noted that we do not have to erase () from Eq. (1) because these ’s are constant. To erase (), we differentiate Eq. (7) to obtain

(8) |

By substituting Eq. (1) in Eq. (8), one obtains

(9) |

where is the inertia constant of pipe . Note that Eq. (2) implies that .

We define the -dimensional diagonal matrix by the diagonal elements (). We also denote the head vector by , where denotes the transposition. Furthermore, we denote by , of which the th element is given by

(10) |

To rewrite Eq. (9), we introduce the incidence matrix, denoted by , which is an matrix. By definition, corresponding to each th edge , the entries of the th column of are given by , and (). Although the edges are undirected, we assign and to and , respectively, by arbitrarily selecting either node forming the th edge to which is assigned, for later convenience. We represent in the block form given by

(11) |

where the matrix is the part of the incidence matrix corresponding to the nodes whose total energy is not fixed, and the matrix is the part of the incidence matrix corresponding to the nodes whose total energy is fixed.

We further denote the combinatorial Laplacian matrix of the network by

(12) |

Using

(13) |

we obtain the following block format of the Laplacian matrix:

(14) |

where is the matrix corresponding to the nodes whose total energy is not fixed, is an matrix and so forth.

The combinatorial Laplacian matrix for the network with the same set of edges but with the edge weight given by is called the conductance matrix in electrical circuit theory [36]. We denote the conductance matrix by . The conductance matrix is given by

(15) |

Then, we can rewrite Eq. (9) as

(16) |

With the notation

(17) |

where and , Eq. (16) is equivalent to

(18) |

Matrix is so-called the loopy Laplacian matrix [36]. Because its all eigenvalues are positive if the network is a connected network, which we have assumed, has the inverse and one obtains

(19) |

This procedure is essentially the same as the network Kron reduction [36].

### 3.2 Local stability analysis

In this section, we derive the eigenequation that determines the conventional local linear stability of the steady state of the nonlinear dynamics given by Eq. (1). We use the combination of Eqs. (1) and (7) for solving the steady state and Eq. (20) to perform the local stability analysis around the obtained steady state.

The steady state is given by setting the left-hand side of Eq. (1) to zero. By combining the resulting constraints with the constraints provided by Eq. (7), one can obtain the steady state, which consists of the unknowns, i.e., () and ().

To set up the Newton-Raphson iteration scheme (e.g. [26]), we rewrite Eq. (7) as

(21) |

We set the left-hand side of Eq. (1) to 0 to obtain

(22) |

where and are the two nodes connected by the th pipe; we use and instead of and to avoid possible notational conflict with Eq. (21). The Jacobian of the right-hand side of Eqs. (21) and (22) with unknowns, , , , , , , is given by

(23) |

where is the diagonal matrix whose diagonal elements are given by

(24) |

The update equation for the Newton-Raphson method is given by

(25) |

Denote the obtained steady state by . We now carry out the local stability analysis around this steady state. The Jacobian matrix of the right-hand side of Eq. (20), denoted by , is given by

(26) |

###### Proposition 3.1.

Matrix has -fold zero eigenvalues and negative eigenvalues.

###### Proof.

A direct substitution verifies

(27) |

Therefore, has -fold zero eigenvalues, and the corresponding right eigenvectors are the columns of .

Because is an matrix, its rank is at most . Therefore, the matrix has at most rank and therefore has at least zero eigenvalues.

We denote the eigenvalues of a symmetric matrix by . Using the fact that is positive definite because all its diagonal elements are positive (i.e., ), the Weyl’s theorem on eigenvalues applied to the sum of and , which is equal to , implies that

(28) |

Note that for an (). Therefore, Eq. (28) yields . Because this relationship must be consistent with the fact that has -fold zero eigenvalues, we obtain and . ∎

Let us assume for any pipe that , which means that the head loss owing to friction along the pipe increases as the steady-state flow rate, , increases. We verified held true in all the numerical simulations carried out in the next section. Then, has negative eigenvalues and zero eigenvalues as matrix does (Appendix A).

Therefore, the dynamics given by Eq. (20) are locally neutrally stable in the subspace spanned by the column vectors of . However, perturbations along these directions are infeasible because such a perturbation is not on the hyperplane defined by Eq. (7) on which the dynamics are constrained. Note that the normal vectors of the hyperplane are given by the row vectors of , i.e., column vectors of . Therefore, in the local stability analysis, we should ignore these zero eigenvalues and inspect the other eigenvalues. With this consideration, we conclude that the steady state constrained by the Kirchhoff’s current law (i.e., Eq. (7)) is always locally stable.

To define a local stability index, we calculate the largest eigenvalue of except the -fold zero eigenvalue. Note that all eigenvalues of are real although is not a symmetric matrix in general (Appendix A). We consider the the largest non-zero eigenvalue of (which is negative) and refer to its absolute value as the local stability index, denoted by . A large value implies that dynamics perturbed by a small amount from the steady state would rapidly return to the steady state. Mathematically, the perturbations are given to the flow rate for some pipes. However, we expect that serves as an index that can characterise dynamical stability of water distribution networks in wider contexts.

## 4 Numerical results

In this section, we calculate the local stability index for various networks and examine its relevance to stability and resilience in simulated dynamics of water flow.

We use the synthesized networks used in our previous study [18]. As supplementary materials of Ref. [18], the full network structure and properties of system components (e.g., pipes, nodal demands, elevations) were made open to the public in the EPANET2 format, except for the one network. For completeness, we have made the complete data for all the networks and their system components used in the present study available in the MATLAB format at Github, alongside the code for calculating the local stability index.

Of the 85 networks used in our previous study [18], five networks include pumps, which would pin the flow rate to a constant value. Because our formalism does not directly cover the case of the pinned flow rate, we exclude these five networks. Of the remaining 80 networks, 10 networks have nodes, reservoir nodes among the nodes and pipes; 11 networks have ; 10 networks have ; 10 networks have ; 9 networks have ; 9 networks have ; 5 networks have ; 5 networks have ; 5 networks have ; 5 networks have ; 1 network has . We calculated for each of the 80 networks.

The diameter of a majority of pipes was equal to 400 mm. There were also a few larger values of the diameter (e.g., 900 mm). Such wider pipes were incident (i.e., directly connected) to a reservoir. The diameter of the pipes was automatically generated by the software [37].

We set the total energy of all the reservoirs in each network to 65 m [37]. In addition to these nodes, some nodes in each network had zero water demand (i.e., ) because they are pure junctions connecting different pipes. We set for each pipe to zero because the data are from the 80 networks that do not have valves or pumps.

In our previous study, we numerically simulated dynamics of water flow in these networks using software EPANET2 [38, 23] and measured six strain indices to characterise resilience of water distribution systems. These indices are as follows (see [18] for fuller definitions). First, the time to strain is defined as the time between the application of stress and the start of service failure, which is when the level of service at a node drops below a predefined threshold. Second, the failure duration is defined as the time needed for the system to recover to normal performance since the service failure triggered by the application of stress. Third, the failure magnitude is the most severe drop in the system service performance at a node as a result of the administered stress. Fourth, the failure rate is defined as the failure magnitude divided by the time between the start of the failure and the occurrence of the worst system performance. Fifth, the recovery rate is defined as the failure magnitude divided by the time between the worst system performance and the return to the performance threshold value used in the definition of the time to strain. Sixth, the severity is defined as the threshold minus the system performance (which is positive during the system failure) integrated with respect to time over the period of system failure.

The relationship between each of the six strain indices and the local stability index, , is shown in Fig. 1. In the figure, each circle represents one of the 80 networks, represents the Pearson correlation coefficient and represents the -value for the Pearson correlation coefficient. Figure 1 indicates that is strongly correlated with the failure magnitude, failure rate and recovery rate. In particular, a large value is intuitively associated with a high recovery rate because the recovery rate quantifies how the system returns to the normal performance after being perturbed by external stress. The positive correlation between and the recovery rate (Fig. 1(e)) confirms this intuition. We do not have an active interpretation of the strong correlation results shown in Figs. 1(c) and 1(d) because the failure magnitude and failure rate quantify the magnitude of failure response given the stress, with which the local stability analysis is not directly concerned. Although it may sound contradictory that a large value is associated with a large failure magnitude and failure rate, this result is not due to the definition of . In fact, within our previous numerical results [18], we found a positive correlation between the failure magnitude and the recovery rate (, ) and between the failure rate and the recovery rate (, ). Our interpretation is that the failure magnitude and the failure rate characterise different aspects of the system’s resilience from what the recovery rate does.

In our previous study, we measured eight structural properties of the water distribution network and examined the association between each of them and each of the six strain indices [18]. Here we measured the correlation between each of these structural properties, whose definitions are in Appendix B, and . We found that none of the eight structural properties of the network was as strongly correlated with the recovery rate as was (link density: , ; algebraic connectivity: , ; diameter: , ; average path length: , ; central point dominance: , ; heterogeneity: , ; spectral gap: , ; clustering coefficient: , ). Therefore, the local stability index is better at capturing the recovery rate than these structural properties. This is probably because these structural properties do not directly quantify the speed at which the water distribution network responds to a perturbation to the steady state.

## 5 Discussion

We carried out a local stability analysis of water flow in pipe networks including reservoir nodes and consumer nodes. The steady state was shown to be always linearly stable. We used the eigenvalue corresponding to the slowest relaxation mode as the local stability index. The proposed index was moderately correlated with the recovery rate in response to external stress, which was numerically obtained from an involved simulator commonly used in water engineering research community.

We used a particular formula for the friction factor [34]. However, the present framework accepts other formulae for the friction factor. Because the differentiation of the friction factor with respect to the Reynolds number is used in calculating the local stability index and the value of the friction factor depends on the pipe in general, one should use a formula that is continuous and covers a wide range of the Reynolds number and hence different flow schemes. Some possible choices of the formula different from the one used in the present paper are found in the literature [39, 40, 41]. It is straightforward to extend our local stability index to the case of different functional forms of the friction factor.

In 1972, May proposed to examine the eigenvalue spectra of interaction networks to characterise dynamical stability of complex systems, in particular ecological systems [42]. This is a landmark example of local stability analysis. While the original analysis and many studies that ensued were for random networks [42], the understanding has been extended to the case of interaction networks with general network structure [43, 44]. This family of method is similar to the one proposed in the present study. Main differences are that we have started from a particular set of hydraulic equations modelling water flow in pipe networks and that we have found that our system is always locally stable and hence used the relaxation mode to quantify the extent of the local stability of the system.

The proposed local stability index was significantly correlated with the recovery rate, which was obtained from numerical simulations in which the applied shock was not necessarily small. However, because the present analysis is a local stability analysis, it does not tell us the behaviour of the system when a shock injected to the steady state is not small. In water engineering applications, the magnitude of perturbation is not necessarily small [3]. A recent study in network science has developed a formalism that reduces a class of dynamical systems on interaction networks into a one-dimensional effective dynamical system and assesses the resilience of the dynamical system against perturbation that is not necessarily small [45]. Therefore, adapting their resilience formalism to the case of transient and steady-state dynamics of water flow in pipe networks is of a practical relevance. In Eq. (1), there are usually only a fraction of nodes whose energy value is fixed (i.e., reservoirs) and different nodes have different water demands. In this sense, the flow dynamics on pipe networks are heterogeneous in addition to being heterogeneous in the network structure and edge weight. Other structures such as valves would add more complexity and heterogeneity to the system. For such a heterogeneous system, it is unclear whether one can transform Eq. (1) into a form that is compatible with the resilience function formalism [45, 46]. This issue is left as future work.

## Acknowledgements

N.M. and F.M. acknowledge the support provided through BRIM (Building Resilience into Risk Management) project, EPSRC.

## Appendix A: Eigenvalues of

Here we provide an elementary proof that has zero eigenvalues and negative eigenvalues when , .

Let us introduce a short-hand notation . Then, one obtains . Because is similar to the symmetric matrix , the min-max theorem for the th smallest eigenvalue implies that

(29) |

Similarly, one obtains

(30) |

By combining , and , one obtains

## Appendix B: Structural properties of networks

The eight structural properties of the network that we measured in our previous study [18] and used in the present article for comparison purposes are as follows.

The link density is given by . The algebraic connectivity is the smallest positive eigenvalue of the Laplacian matrix, , of the undirected and unweighted network. The version of the clustering coefficient used in [18] is given by . The average path length is equal to the shortest path length between two nodes, which is averaged over all the possible node pairs. The central point dominance is given by , where is the betweenness centrality of the th node and . The heterogeneity is that in the degree and given by , where is the degree of the th node and is the average degree of the node. The version of the spectral gap used in [18] is given by the difference between the two largest eigenvalues of . The modularity is a quantity representing the quality of community structure detected in the network and approximate modularity maximisation was carried out using Newman’s algorithm [47].

## References

- [1] Butler D, Memon FA. 2006 Water Demand Management. London, UK: IWA Publishing.
- [2] Rasekh A, Brumbelow K. 2014 Drinking water distribution systems contamination management to reduce public health impacts and system service interruptions. Environ. Model. Software 51, 12–25.
- [3] Gullick RW, LeChevallier MW, Case J, Wood DJ, Funk JE, Friedman MJ. 2005 Application of pressure monitoring and modelling to detect and minimize low pressure events in distribution systems. J. Water Supply Res. Tech. Aqua 54, 65–81.
- [4] Jun H, Loganathan GV, Deb AK, Grayman W, Snyder J. 2007 Valve distribution and impact analysis in water distribution systems. J. Environ. Eng. 133, 790–799.
- [5] Giustolisi O, Savic D. 2010 Identification of segments and optimal isolation valve system design in water distribution networks. Urban Water J. 7, 1–15.
- [6] Hosseini S, Barker K, Ramirez-Marquez JE. 2016 A review of definitions and measures of system resilience. Reliability Eng. Syst. Safety 145, 47–61.
- [7] Xue X, Wang L, Yang RJ. 2018 Exploring the science of resilience: critical review and bibliometric analysis. Nat. Hazards 90, 477–510.
- [8] Butler D, Ward S, Sweetapple C, Astaraie-Imani M, Diao K, Farmani R, Fu G. 2017 Reliable, resilient and sustainable water management: the Safe & SuRe approach. Global Challenges 1, 63–77.
- [9] Newman MEJ. 2010 Networks — An Introduction. Oxford: Oxford University Press.
- [10] Barabási AL. 2016 Network Science. Cambridge: Cambridge University Press.
- [11] Yazdani A, Jeffrey P. 2011 Complex network analysis of water distribution systems. Chaos 21, 016111.
- [12] Yazdani A, Otoo RA, Jeffrey P. 2011 Resilience enhancing expansion strategies for water distribution systems: A network theory approach. Environ. Model. Software 26, 1574–1582.
- [13] Yazdani A, Jeffrey P. 2012 Water distribution system vulnerability analysis using weighted and directed network models. Water Resour. Res. 48, W06517.
- [14] Shuang Q, Zhang M, Yuan Y. 2014 Performance and reliability analysis of water distribution systems under cascading failures and the identification of crucial pipes. PLOS ONE 9, e88445.
- [15] Herrera M, Abraham E, Stoianov I. 2016 A graph-theoretic framework for assessing the resilience of sectorised water distribution networks. Water Resour. Manage. 30, 1685–1699.
- [16] Pandit A, Crittenden JC. 2016 Index of network resilience for urban water distribution systems. Int. J. Crit. Infrastruct. 12, 120–142.
- [17] Hwang H, Lansey K. 2017 Water distribution system classification using system characteristics and graph-theory metrics. J. Water Resour. Plann. Manage. 143, 04017071.
- [18] Meng F, Fu G, Farmani R, Sweetapple C, Butler D. 2018 Topological attributes of network resilience: A study in water distribution systems. Water Res. 143, 376–386.
- [19] Todini E. 2000 Looped water distribution networks design using a resilience index based heuristic approach. Urban Water 2, 115–122.
- [20] Zhuang B, Lansey K, Kang D. 2013 Resilience/availability analysis of municipal water distribution system incorporating adaptive pump operation. J. Hydraul. Eng. 139, 527–537.
- [21] Yannopoulos S, Spiliotis M. 2013 Water distribution system reliability based on minimum cut – Set approach and the hydraulic availability. Water Resour. Manage. 27, 1821–1836.
- [22] Zhao X, Chen Z, Gong H. 2015 Effects comparison of different resilience enhancing strategies for municipal water distribution network: A multidimensional approach. Math. Probl. Eng. 2015, 438063.
- [23] Diao K, Sweetapple C, Farmani R, Fu G, Ward S, Butler D. 2016 Global resilience analysis of water distribution systems. Water Res. 106, 383–393.
- [24] Liu H, Walski T, Fu G, Zhang C. 2017 Failure impact analysis of isolation valves in a water distribution network. J. Water Resour. Plann. Manage. 143, 04017019.
- [25] Islam MR, Chaudhry MH. 1998 Modeling of constituent transport in unsteady flows in pipe networks. J. Hydraul. Eng. 124, 1115–1124.
- [26] Boulos PF, Lansey KE, Karney BW. 2006 Comprehensive Water Distribution System Analysis Handbook for Engineers and Planners. Pasadena, CA: NWH Soft second edition.
- [27] Todini E. 2011 Extending the global gradient algorithm to unsteady flow extended period simulations of water distribution systems. J. Hydroinfo. 13, 167–180.
- [28] Nault JD, Karney BW. 2016 Improved rigid water column formulation for simulating slow transients and controlled operations. J. Hydraul. Eng. 142, 04016025.
- [29] Jung BS, Karney B. 2017 A practical overview of unsteady pipe flow modeling: from physics to numerical solutions. Urban Water J. 14, 502–508.
- [30] Wylie EB. 1983 The mircocomputer and pipeline transients. J. Hydraul. Eng. 109, 1723–1739.
- [31] Ghidaoui MS, Zhao M, McInnis DA, Axworthy DH. 2005 A review of water hammer theory and practice. Appl. Mech. Rev. 58, 49–76.
- [32] Chaudhry MH. 2014 Applied Hydraulic Transients. Cambridge, England: Springer third edition.
- [33] Todini E, Pilati S. 1988 A gradient algorithm for the analysis of pipe networks. In Coulbeck B, Orr CH, editors, Computer Applications in Water Supply, Volume 1 — Systems Analysis and Simulation pp. 1–20. New York, NY: John Wiley & Sons Inc.
- [34] Bellos V, Nalbantis I, Tsakiris G. 2018 Friction modeling of flood flow simulations. J. Hydraul. Eng. 144, 04018073.
- [35] Rossman L. 2000 EPANET2 Users’ Manual. Cincinnati, USA: Drinking Water Research Division, Risk Reduction Engineering Laboratory, Office of Research and Development, U. S. Environmental Protection Agency.
- [36] Dörfler F, Simpson-Porco JW, Bullo F. 2018 Electrical networks and algebraic graph theory: Models, properties, and applications. Proc. IEEE 106, 977–1005.
- [37] De Corte A, Sörensen K. 2014 HydroGen: an artificial water distribution network generator. Water Resour. Manage. 28, 333–350.
- [38] Mugume SN, Gomez DE, Fu G, Farmani R, Butler D. 2015 A global analysis approach for investigating structural resilience in urban drainage systems. Water Res. 81, 15–26.
- [39] Swamee PK. 1993 Design of a submarine oil pipeline. J. Transp. Eng. 119, 159–170.
- [40] Cheng NS. 2008 Formulas for friction factor in transitional regimes. J. Hydraul. Eng. 134, 1357–1362.
- [41] Brkić D, Praks P. 2018 Unified friction formulation from laminar to fully rough turbulent flow. Appl. Sci. 8, 2036.
- [42] May RM. 1972 Will a large complex system be stable?. Nature 238, 413–414.
- [43] Okuyama T, Holland JN. 2008 Network structural properties mediate the stability of mutualistic communities. Ecol. Lett. 11, 208–216.
- [44] Allesina S, Tang S. 2012 Stability criteria for complex ecosystems. Nature 483, 205–208.
- [45] Gao J, Barzel B, Barabási AL. 2016 Universal resilience patterns in complex networks. Nature 530, 307–312.
- [46] Tu C, Grilli J, Schuessler F, Suweis S. 2017 Collapse of resilience patterns in generalized Lotka-Volterra dynamics and beyond. Phys. Rev. E 95, 062307.
- [47] Newman MEJ. 2004 Fast algorithm for detecting community structure in networks. Phys. Rev. E 69, 066133.