Non-Local Impact of Link Failures in Linear Flow Networks

The failure of a single link can degrade the operation of a supply network up to the point of complete collapse. Yet, the interplay between network topology and locality of the response to such damage is poorly understood. Here, we study how topology affects the redistribution of flow after the failure of a single link in linear flow networks with a special focus on power grids. In particular, we analyze the decay of flow changes with distance after a link failure and map it to the field of an electrical dipole for lattice-like networks. The corresponding inverse-square law is shown to hold for all regular tilings. For sparse networks, a long-range response is found instead. In the case of more realistic topologies, we introduce a rerouting distance, which captures the decay of flow changes better than the traditional geodesic distance. Finally, we are able to derive rigorous bounds on the strength of the decay for arbitrary topologies that we verify through extensive numerical simulations. Our results show that it is possible to forecast flow rerouting after link failures to a large extent based on purely topological measures and that these effects generally decay with distance from the failing link. They might be used to predict links prone to failure in supply networks such as power grids and thus help to construct grids providing a more robust and reliable power supply.


I. INTRODUCTION
The robust operation of supply networks is essential for the function of complex systems across scales and disciplines. Almost all of our technical and economical infrastructure depends on the reliable operation of the electric power grid [1,2]. Biological organisms distribute water and nutrients via their vascular networks, for instance in plant leaves [3], the human and animal circulatory system [4], or in protoplasmic veins of certain slime molds [5]. Huge amounts of money and assets are exchanged through a complex financial network [6]. Structural damages to such networks can have catastrophic consequences such as a stroke, a power outage or a major economic crisis.
In power grids, large scale outages are typically triggered by the failure of a single transmission or generation element [7][8][9][10][11]. The outages in the United States in 2003, Italy in 2003 and Western Europe in 2006 are very well documented and provide a detailed insight into the dynamics of a large scale network failure [12][13][14]. Each outage was triggered by the loss of a transmission line during a period of high grid load. Subsequently, the power flows were rerouted, causing secondary overloads and eventually a cascade of failures. In these three examples, the cascades propagated mostly locally -overloads took place in the proximity of previous failures. However, this is not necessarily the case during power outages (see, e.g. [15]), raising the question of how networks flows are rerouted after failures [16][17][18][19][20][21][22][23].
In biological distribution networks, robustness against link failure is a critical prerequisite that guards against * JS and FK contributed equally to this work. potentially life-threatening events such as stroke [24] or embolism [25,26], but also to function efficiently in the presence of fluctuations [3,27,28]. Thus, biological networks are often (but not in all cases, such as in the penetrating arterioles of the cortical vasculature [29]) endowed with highly resilient, redundant topologies that optimize rerouting of flow in case of link failure to the network [28] and are generated through adaptive developmental mechanisms [30]. For the understanding of such life-threatening conditions it is therefore crucial to investigate the behavior of the vascular network in the case of failure.
To understand the vulnerability of networks, we here provide a detailed analysis of the impact of link failures in linear flow networks. We focus on how the network topology determines the overall network response as well as the spatial flow rerouting. We consider linear supply network models, where the flow between two adjacent nodes is proportional to the difference of the nodal potential, pressure or voltage phase angle. Linear models are applied to hydraulic networks [31], vascular networks of plants and animals [28,[32][33][34][35], economic input-output networks [36] as well as electric power grids [37][38][39][40][41][42]. The linearity allows to obtain several rigorous bounds for flow rerouting in general network topologies and to solve special cases fully analytically.

II. LINEAR FLOW NETWORKS
Consider a network consisting of N nodes that are connected to each other via lines denoted by (m, n) for a line going from node m to node n. Assign a potential or phase angle θ m ∈ R to each node m in the network. Then we assume the flow F m→n between nodes m and n connected via line (m, n) to be linear in the potential drop along the line (1) Here, b mn = b nm is the transmission capacity assigned to the line (m, n) that describes its ability to carry flow. This equation may for example be used to describe hydraulic networks [31,43] or vascular networks of plants [28], where the θ n denotes the pressure at some node n and the capacity b mn scales with the diameter of a pipe or vein. Our main focus will be its application to electric power engineering, where this linear approximation of the power flow equations is referred to as the DC approximation [38][39][40]. In this case, F m→n refers to the flow of real power along a transmission line (m, n), θ n is the voltage phase angle at node n and b mn is proportional to the line's susceptance.
In addition to that, we assume that Kirchhoff's current law holds at the nodes of the network which states that the inflows and outflows at any node m balance where the right-hand-side denotes the inflow (P m > 0) or outflow (P m < 0) at node m, commonly called the 'power injection' in power engineering. Equations (1) and (2) fully describe the state and the flow of the network once the line parameters b mn and the injections P m are given. These equations may be conveniently written using a vectorial notation. Define the vector θ = (θ 1 , . . . , θ N ) ∈ R N of the nodal potentials or voltage phase angles and the vector P = (P 1 , . . . , P N ) ∈ R N of nodal injections. Here and in the following sections, the superscript ' ' denotes the transpose of a vector or matrix. We further label all lines in the grid by = 1, . . . , L and summarize all line flows in a vector F = (F 1 , . . . , F L ) ∈ R L . Equation (1) may then be rewritten as where B d ∈ R L×L is a diagonal matrix containing the capacities b of all edges. Furthermore, we defined the node-edge incidence matrix K ∈ R N ×L . To define this matrix in an undirected graph, one typically fixes an arbitrary orientation of the graph's edges such that its components read The node-edge incidence matrix also relates the injections to the flows incident at a node. More specifically, Kirchhoff's current law (2) may be rewritten as follows Here, we defined the matrix B = KB d K ∈ R N ×N commonly referred to as the nodal susceptance matrix in power engineering. Mathematically, B is a weighted Laplacian matrix [44,45] with components Here, Λ m is the set of lines which are incident to m.

III. ALGEBRAIC DESCRIPTION AND ANALYSIS OF LINE OUTAGES
An important question in network security analysis is how the flows in the network change if a line fails. Denoting by F the initial flow of the failing line =(r, s), the flow change ∆F e at a transmission line e=(m, n) is written as Adopting the language of power system security analysis [37,38], we call the factor of proportionality the line outage distribution factor (LODF). In the following, we present two alternative derivations as well as a physical interpretation of the linear flow rerouting problem.

A. Self-consistent derivation of line outage distribution factors
To derive an explicit expression for the LODFs one generally starts with a related problem. Consider an increase of the real power injection at node r and a corresponding decrease at node s by the amount ∆P . The new vector of real power injections is then given bŷ where the components of ν rs ∈ R N are +1 at position r, −1 at position s and zero otherwise. Here and in the following, we use a hat to indicate the state of the network after a line outage or a similar change of the topology. The change of the real power injections causes a change of the nodal voltage angles, where B † denotes the Moore-Penrose pseudo-inverse of the Laplacian matrix B. Hence, the real power flows change by The factor of proportionality is referred to as the power transfer distribution factor (PTDF).
The LODFs can be expressed by PTDFs in the following way [38]. To consistently model the outage of line (r, s), one assumes that the line is disconnected from the grid by circuit breakers and that some fictitious real power ∆P is injected at node s and taken out at node r. The entire flow over the line (r, s) after the opening thus equals the fictitious injectionsF rs = ∆P . Using PTDFs, we also know that F rs = F rs + PTDF (r,s),r,s ∆P .
SubstitutingF rs = ∆P and solving for ∆P yields The change of real power flows of all other lines is given by Equation (5) such that we finally obtain LODF (mn),(rs) = PTDF (m,n),r,s 1 − PTDF (r,s),r,s .
For consistency, one usually defines the LODF for the failing line as follows LODF (rs),(rs) = −1. In addition to that, we exclude cases where the failing line is a bridge, i.e., a line whose removal disconnects the graph, from our analysis in the following sections.

B. Algebraic derivation of line outage distribution factors
The LODFs can also be obtained in a purely algebraic way without any self-consistency argument [46]. As the line =(s, r) fails, the nodal susceptance matrix of the network changes as which causes a change of the nodal potentials or voltage phase angles respectively, Equation (3) for the perturbed grid now reads Subtracting Equation (3) for the unperturbed grid, we see that the change of the voltage angles is given by The change of flows after the outage of line (r, s) and thus the LODFs are calculated from the change of the voltage angles which yields In principle, we could now use these equations to calculate the flow changes and the LODFs. However, this would require to invert the matrixB = B + ∆B separately for every possible line (s, r) in the grid, which is impractical. Nevertheless, we can simplify the expression using the Woodbury matrix identity, such that the flow change (9) reads which is identical to Equation (6) obtained using the standard approach.

C. Electrostatic interpretation
A deeper physical insight into the network flow rerouting problem is obtained by the analogy to discrete electrostatics. Equation (8) can be rearranged into a linear set of equations for the change of the nodal potentialŝ Here,B is the Laplacian of the grid after the failure, i.e. the grid where line (r, s) has been removed. Alternatively, we can formulate the equation in terms of the original network topology, substituting Equation (10) into Equation (8). This yields the linear set of equations with the dipole source As noted before, B andB are Laplacian matrices and the right-hand side of both equations (12) and (11) are nonzero only at positions r and s with opposite sign. Hence, these equations are discrete Poisson equations with a dipole source and ψ is a dipole potential, see Ref. [47, ch. 15] for a detailed analysis of this equation. The main complexity of the line outage problem thus arises from the network topology encoded in the Laplacian B, which can be highly irregular. The two equations (12) and (11) yield the same potential ψ, but are formulated on different topologieseither on the original network topology or the topology after the outage. To compare the impact of different failures it is beneficial to use the original topology, such that only the dipole inhomogeneity differs -not the electrostatic problem itself. Then, the strength of the dipole depends on the network topology via the prefac- Using the analogy to electrostatics we can solve the flow rerouting problem for regular network topologies (section IV) and provide some general rigorous results (section VI A). To understand flow rerouting in networks with complex topologies, we thus have to account for the spatial spreading pattern described by B † (see section VI C) as well as the dipole strength, which quantifies the gross response of the grid (see section V).

IV. FAILURES IN REGULAR NETWORKS AND THE CONTINUUM LIMIT
To obtain a first insight into the spatial aspects of flow rerouting, we consider an elementary example admitting a solution in the continuum limit. Consider a regular square lattice embedded in a plane as depicted in Figure 1 and studied in a slightly different form in Ref. [48]. All nodes are labeled by their positions r = (x, y) in this two-dimensional embedding and the lattice spacing is denoted as h. We introduce continuous functions ψ and b such that ψ(x, y) is the potential of the node at (x, y) and b(x + h/2, y) is the weight of the link connecting the two nodes at (x, y) and (x + h, y). The left-hand side of the Poisson equation (12) evaluated at position (x, y) reads Here, we made use of the fact that the components of the gradient ∇ψ = (∂ x ψ, ∂ y ψ) may be expressed as but did not take the limit yet. The derivative with respect to y may be calculated analogously.
Before we proceed to the right-hand side, we remark that the flow changes ∆F according to equation (9) are given by given by such that we obtain in the continuum limit h → 0, Note that the expression ∆F refers to the change in flow due to the link failure here and should not be confused with the continuous Laplace operator.
The right-hand side of the discrete Poisson equation (12) may be calculated similarly noting that only two nodes contribute with opposite signs. Let us assume that the failing link is parallel to the x-axis connecting nodes r and s located at r r = (x r , y r ) and r s = (x s , y s ) = (x r + h, y r ) . The discrete version of the right-hand side reads We will now derive the continuum version of this equation. First, the flow on the failing link before the outage F rs may be calculated as where F (x, y) = b(x, y)∇θ(x, y) is the continuum flow before the outage. Second, the vector ν rs can be formally interpreted in terms of the two-dimensional delta function δ(x, y) and reads for the given link failure Finally, let us assume that a continuum version of the Green's function B † exists. Then the denominator may be calculated as Thus, in total we obtain after expanding the entire right-hand side to lowest order in the continuum limit Here, F (x r , y r ) is assumed to be parallel to the dipole axis, i.e. the direction of the link failure, which is either the x-or the y-direction for the given setting.
We can now formally divide left-hand side (14) and right-hand side (16) by h 2 and take the limit h → 0 to obtain the final continuum limit of the Poisson equation, where the source term is q(x r , y r ) = F (x r , y r ), the unperturbed current field. We note that we obtain the same continuum limit regardless of whether we use Equation (11) or (12) to do the expansions. Thus, the non-locality that is encoded in Equation (12) is lost in the continuum formulation.
If the link weights are homogeneous, b(x, y) = b, and the failing link is assumed to be located at the origin (x r , y r ) = (0, 0) the solution is given by the well-known two-dimensional dipole field We thus obtain a fully analytic solution in the continuum limit. This solution reveals that the impact of link failures decays algebraically in homogeneous lattices. We consider this decay along two different axes. Assume the dipole to be located at the origin in x-direction, such that q = (ε, 0) where ε 1 is some small real number. First, consider the decay in x-direction where r = (x, 0) . In this case, we obtain for the decay of the potential and the flow changes This decay in the flow changes may also be observed in the discrete version of Equation (16) and is shown in Figure 2, (a), for a line failure in a discrete square grid. Along the same lines, we may quantify the decay in ydirection r = (ε, y) for the same dipole orientation. In this case, we obtain Here, we assumed the position vector to be dominated by its y-component ε y such that ||r|| ≈ |y|. In total, we observe a y −3 -scaling in the flow changes in y-direction perpendicular to the dipole source and a y −2 -scaling in y-direction parallel to the dipole source, see Figure 2, (a-c).

V. RIGOROUS BOUNDS ON THE DIPOLE STRENGTH
We now turn to realistic networks with irregular topologies. The change in the nodal potentials or voltage phase angles ψ n and flows ∆F m→n is determined by the discrete Poisson equation (12). We first consider the right-hand side of this equation, the dipole strength, which describes the gross response of the network flows to the outage. This response is proportional to the initial flow of the failing edge F rs and the factor The factor 1 − η(r, s) describes the nonlocality of the network response to a local perturbation at link (r, s). To see this, consider a grid where the real power ∆P is injected and withdrawn at the terminal nodes of the link (r, s). The direct flow over the link is given by whereas the total flow is just given by ∆P . The factor η(r, s) thus measures the fraction of the flow which is transmitted directly and 1 − η(r, s) is the fraction transmitted non-locally via other pathways. Hence, 1 − η(r, s) can also be seen as a measure of redundancy. A high non-local flow indicates that there are strong alternative routes from r to s in addition to the direct link (r, s). If no alternative path exists, the flow must be routed completely via the direct link such that 1 − η(r, s) = 0. We conclude that the properties of alternative and direct paths is decisive for the understanding of flow rerouting. Before we proceed, we thus review the formal definition of a path in graph theory. Definition 1. A path from vertex r to vertex s is defined as an ordered set of vertices LODFs are calculated for (a) links along the x-direction (between (x, 0) and (x + 1, 0)), (b) links along the y-direction parallel to failing link (between (0, y) and (1, y)) and (c) links along the y-direction perpendicular to the failing link (between (0, y) and (0, y + 1)). The dist −2 (a,b) and dist −3 (c) scaling agrees with the dipole scaling predicted using Equation (19) as indicated by black lines. The levels of sparsity considered here do not show any effect on the scaling when considering directions parallel to the dipole axis (a,b), but the scaling becomes more long-ranged with increasing sparsity in direction perpendicular to this axis (c). (d) The dist −2 scaling is not unique to square grids (purple squares, size 1000 × 1000) but may also be observed for the two other regular tilings, namely the hexagonal grid (orange hexagons, 150 × 150 hexagons) and the triangular grid (green triangles, size 1001 × 500). LODFs were again computed along the shortest path in x-direction for links oriented parallel to the dipole. The branching for the hexagonal grid is due to the fact that the path in x-direction is non-unique and non-straight here, such that one of the shortest paths was chosen arbitrarily. Deviations from the scaling occur for large distances due to finite-size effects.
where two subsequent vertices must be connected by an edge and no vertex is visited twice. Two paths are called independent if they share no common edge. The unweighted length of such a path is defined as the number of steps k, while the weighted path length is given by the sum of the edge weights along the path, k j=1 w vj−1vj . In this work, the edge weights are given by the inverse susceptances w ij = 1/b ij . The geodesic or shortest path distance of two vertices r and s is defined as the length of the shortest path from r to s.
The interpretation as a redundancy measure directly relates the factor 1 − η(r, s) to the topology of the network. A first rough estimate can be obtained from the topological connectivity λ T (r, s), which is defined as the number of independent paths from node r to node s. A comparison for several test grids in Figure 3 shows that η(r, s) decreases with λ T (r, s) on average as expected, but that there is a large heterogeneity between the links.
To obtain a better topological estimate for the locality factor we need to take into account the heterogeneity of the link weights. The topological connectivity λ T (r, s) counts the minimum number of edges which have to be removed to disconnect the nodes r and s. We can define a weighted analog λ F (r, s) as the minimum capacity which has to be removed to disconnect the nodes r and s. This is a classical problem in graph theory, where it is referred to as the minimum cut [50]. We will now elaborate this quantity in a definition. An (r, s)-cut can be defined as follows. Let r ∈ S ⊂ V and s ∈ V \ S be two vertices taken from the two disjoint sets. The (r, s)-cut is defined as the set of edges δ(S) = {(u, v) ∈ E | u ∈ S, v ∈ V \ S or v ∈ S, u ∈ V \S} connecting the two disjoint vertexsets. The set of edges δ + (S) = {(u, v) ∈ E | u ∈ S, v ∈ V \ S} is referred to as the forward edges of the cut. The capacity C of a cut δ(S) and the corresponding minimum capacity λ F (r, s) between r and s are then given by By virtue of the max-flow-min-cut theorem [51], λ F (r, s) is equivalent to the maximum flow which can be transmitted from r to s respecting link capacity limits: Numerous efficient algorithms exist to calculate this maximum flow without performing the optimization explicitly [51]. The ratio b rs /λ F (r, s) then gives the ratio of direct flow to total flow from r to s and thus provides an adequate topology-based estimate for the locality factor η(r, s). Indeed, we can prove that it provides a rigorous lower bound.
Proposition 1. The algebraic locality factor η(r, s) is bounded by A proof is given in appendix A. Numerical simulations for several test grids reported in Figure 4 reveal that the topological estimate not only provides a lower bound, but a high-quality estimate for the algebraic locality factor. The Pearson correlation coefficient ρ between η(r, s) and b rs /λ F (r, s) exceeds 0.92 for the three grids under consideration.
We arrive at the conclusion that the dipole strength given by F rs (1 − η(r, s)) −1 generally decreases with the redundancy measures λ T (r, s) and λ F (r, s).
An upper limit for the locality factor η(r, s) can be obtained from an elementary topological distance measure. We consider the weighted geodesic distance of the two nodes r and s after the failure of the direct link (r, s), which we denote by dist w 1 (r, s). The superscript w stands for weighted distance, the subscript 1 for the distance measured in the graph after removal of the link (r, s). We then have the following upper bound. Proposition 2. The algebraic locality factor η(r, s) is bounded from above by A proof is given in appendix B. Numerical simulations for several test grids reported in Figure 5 reveal that the estimate in terms of the shortest path length not only provides an upper bound, but a high quality estimate for the algebraic locality factor. The Pearson correlation coefficient ρ exceeds 0.94 for the three grids under consideration.
We further note that the factor ν rs B † ν rs can also be interpreted as a distance measure -the resistance distance [52,53]. We come back to the quantification of distances in flow networks later in section VI C.

VI. SPATIAL DISTRIBUTION OF FLOW REROUTING
We now turn to the spatial aspects of flow rerouting in general network topologies. We first discuss some rigorous results, showing how the network topology determines the rerouting flows. Then, we return to the regular tilings and study the effect of increasing sparsity in these topologies on the dipole scaling. Finally, we suggest a new measure of distance for flow rerouting and examine its performance on realistic network topologies taken from power grids.

A. Rigorous results
To start off, we first present a lemma due to Shapiro [54], relating the flow changes after a link failure in an unweighted graph solely to the topology of the underlying network. Lemma 1. Consider an unweighted network with a unit dipole source along the edge (r, s), i.e. a unit inflow at node r and unit outflow at node s. Then the flow along any other edge (m, n) is given by where N (r, m → n, s) is the number of spanning trees that contain a path from r to s of the form r, . . . , m, n, . . . , s and N is the total number of spanning trees of the graph.
This lemma exactly gives the LODFs in terms of purely topological properties -the number of spanning trees containing certain paths. A generalization of this theorem to weighted graphs was recently presented in [55].
However, counting spanning trees is typically a difficult task such that these results are of limited use for practical applications. Nevertheless, they reveal the importance of certain paths through networks which we will analyze numerically in more detail below. Before we turn to this issue, we derive some weaker, but more easily applicable rigorous results.
We expect that the flow changes ∆F mn decay with distance as for the case of the square lattice analyzed in section IV. Can we establish some rigorous results on the decay with distance for arbitrary networks? Consider the outage of a single edge and assume that the network remains connected afterwards. We label the failing link as (r, s) such that F r→s > 0 w.l.o.g. We first consider the change of the nodal potential or voltage phase angles ψ n and its decay with distance to the failing link (r, s). More specifically, we define the maximum and minimum values of ψ n attained at a given distance: Here, dist u 0 (n, r) denotes the geodesic distance between two nodes n and r in the initial unweighted graph (indicated by the superscript u for unweighted and subscript 0 for the initial pre-contingency network). We then find the following result. Proposition 3. Consider the failure of a single link (r, s) with F r→s > 0 in a flow network. Then the maximum (minimum) value of the potential change ψ n decreases (increases) monotonically with the distance to nodes r and s, respectively: A proof is given in appendix C. We thus find that potential changes generally decrease with the distance in magnitude and so do the flow changes.
Furthermore, we can exploit the analogy to electrostatics to gain an insight into the scaling of flow changes with distance. As the flows are determine by a discrete Poisson equation, a discrete version of Gauss' theorem follows immediately. We note that we formulate this result in terms of the original network topology, cf. Equation (12).
That is, for each decomposition the total flow between the two parts V 1 and V 2 equals the dipole strength.
For general network topologies this lemma implies that the average flow will decay with distance from the failing link (r, s): Choose V 1 to include all nodes which are closer to r than to s and have a distance to r smaller than a given value With increasing value of d the number of nodes in V 1 increases and typically the number of edges between V 1 and V 2 increases, too. The total flow over these links remains constant according to lemma 2, such that the average flow will generally decrease. The exact scaling of the number of edge between V 1 and V 2 of course depends on the topology of the network. One can furthermore show that a sufficient connectivity is needed for perturbations to spread. Generally, flow can be rerouted via an edge (m, n) only if it can enter and leave the link via two independent paths. One can thus prove the following statement [55,56]. Now that we derived rigorous results on the scaling of LODFs, we want to study the influence of network connectivity on the scaling in more detail.
To do so, we first compare the scaling obtained for the square grid to the one in the other two regular tilings of two-dimensional space, namely the hexagonal grid and the triangular grid. In perfect realizations of these grids, each node has a degree of deg hex = 3 and deg tri = 6, respectively, whereas the degree for the square grid reads deg sg = 4. In Figure 2(d), the LODFs are evaluated for these three topologies with increasing geodesic distance from the failing edge located again in the center of the networks between the nodes at (x r , y r ) = (0, 0) and (x s , y s ) = (1, 0). The quadratic scaling with the geodesic distance in x-direction x −2 (black, dotted line) is preserved for all three topologies, i.e. the triangular grid (green triangles, bottom), the square grid (red squares, center) and the hexagonal grid (blue hexagons, top). The grids used here were of size 1000 × 1000 and 1001 × 500 nodes for the square grid and the triangular grid, respectively, and 150 × 150 hexagons for the hexagonal grid.
Thus, the quadratic scaling is robust throughout different regular networks. However, real networks are in general not regular. For this reason, we proceed by studying the effect of increasing sparsity in these regular tilings. Define the sparsity ξ ∈ [0, 1] ⊂ R as the fraction of edges removed from the original graph. We make use of two different methods to achieve increasing sparsity. Our first method is a completely random removal of edges in the graph followed by measuring the LODFs along a specified path. If an edge along the path does no longer exist, we simply skip the edge. The results obtained from this method are shown in Figure 2 (a-c). There is no change visible in the scaling of LODFs, except for the direction perpendicular to the dipole in panel (c). In particular, only small values of sparsity ξ can be studied using this method, since a random removal of edges may easily result in disconnected graphs. For this reason, we make use of another method.
For the second method, we first construct an arbitrary spanning tree of the network after removal of the failing edge. Then, we subsequently remove random edges from the graph that are not part of the tree until a fraction ξ of its original edges is removed from the graph. This way, we make sure that the whole graph stays connected at all times. We continue by constructing the shortest path from the failing edge ((0, 0), (1, 0)) to the node located at (x max , 0) and quantify the LODFs along this path. Note that using this method to make a graph sparser, we need to take into account the graph-specific maximal sparsity ξ max,G , i.e. the fraction of edges whose removal would disconnect the graph. Assuming the initial tree to be minimal, this fraction may be calculated as ξ max,hex = 1/3, ξ max,sg = 1/2 and ξ max,tri = 2/3 for the hexagonal grid, square grid and triangular grid, respectively.
Using this procedure, we can quantify the scaling of LODFs in grids with increasing sparsity. The direct assessment of a scaling exponent is difficult for sparser graphs due to the large spread in LODF values, see Figure 6(a). This is why we construct a different measure to quantify this scaling. We consider the effective exponent k(ξ), where ξ is the graph's sparsity, and assume a scaling of the form in some region of the geodesic distance r = r from the link failure. This effective exponent is calculated as follows where w ∈ N is a window specifying the range to average over in order to smooth the LODF values considered. We chose a windows size of w = 2 when calculating the effective exponent in practice which we found to result in a good compromise between smoothing and completely removing the trend. However, we did not observe a strong effect of the window size on the results. For a perfect inverse square law |LODF| ∝ r −2 and a vanishing window w = 0, this parameter yields k = − log 5 (5 −2 ) = 2 as required. In Figure 6(b), it can be observed that this effective exponent stays approximately constant at k ≈ 2 over different values of sparsity and the three different topologies considered, where results for each value of sparsity were obtained using 100 random realizations of edge removals and with the same grid sizes as stated previously.
To further quantify the effect of increasing sparsity in regular networks, we make use of another measure which we refer to as the LODF ratio R w (ξ). It is simply calculated as the logarithmic ratio between the LODFs with and without sparsity, again averaged over a fixed window of distances R w (ξ) = log 10 r∈[10 1 −w,10 1 +w] |LODF(r, ξ)| r∈[10 1 −w,10 1 +w] |LODF(r, ξ = 0)| .
Note that we evaluate this parameter at a distance of 10 1 but we found the parameter to yield similar values for all distances considered. A parameter of R w (ξ) = 1 then represents a tenfold increase in the LODFs as compared to the network without any edges removed. In Figure 6(c), this parametzer is shown for the different topologies and sparsities. Here, a window size of w = 5 was used. An increase with increasing sparsity is clearly visible. In particular, the LODFs increase on average more than tenfold close to the highest possible values of sparsity.
In total, we observe that the scaling exponent derived from the dipole analogy in section V holds for the regular networks even when removing a large fraction of their edges. On the other hand, the LODF values at a certain distance from the failing link show an increase with increasing sparsity, such that the actual effect of a link failure can be up to tenfold stronger than for the corresponding regular grid with no links removed. Thus, the overall effect of a link failure is more long-ranged in a sparser network, although no change in the effective exponent can be observed.

C. Scaling with distance
The impact of a link failure generally decays with distances. While the definition of distance is straightforward in regular lattices, different measures are meaningful in networks with complex topologies. The geodesic distance of two links follows from definition 1 for two vertices Here, w rs is the edge weight assigned to the edge (r, s). When considering the unweighted analog, the edge distance is defined analogously setting all edge weights to one. The additional term wrs+wmn 2 ensures that neighboring edges have non-zero distance, e.g. unity distance edist u ge = 1 in the unweighted case. However, this distance is a bad indicator for flow rerouting in real-world irregular topologies. An example shown in Figure 8 demonstrates that this simple distance is only weakly correlated with the magnitude of the LODFs for a real-world power grid test case.
Instead, we need a distance measure based on flow rerouting. If a link (r, s) fails, the flow must be rerouted through other pathways, as described by the electrical lemma 1. However, it is not feasible to take into account all spanning trees which govern the flow rerouting. In order to still be able to estimate the impact on another link (m, n), we will thus consider a path from r to s that crosses this link. The main difference to the ordinary graph theoretical distance is that we have to take into account a path back and forth. We are thus led to the following definition.
Definition 2. A rerouting path from vertex r to vertex s via the edge (m, n) is a path where no vertex is visited twice. The rerouting distance between two edges (r, s) and (m, n) denoted by edist u/w re [(r, s), (m, n)] is the length of the shortest rerouting path from r to s via (m, n) plus the length of edge (r, s). Equivalently, it is the length of the shortest cycle crossing both edges (r, s) and (m, n). If no such path exists, the rerouting distance is defined to be ∞.
The definition of a rerouting path is illustrated in Figure 7. Again, we consider a weighted and an unweighted version of this distance indicated by the superscript w and u , respectively. We note that the length of the edge (r, s) is included in order to make the distance measure symmetric. In Appendix D, we show explicitly that this definition satisfies the axioms of a metric and discuss how to compute the shortest rerouting path.
An example of rerouting distances in comparison to the LODFs is shown in Figure 8 for a small test grid. We observe a much better correlation in comparison to the ordinary geodesic distance defined above. The limitation of geodesic distances becomes especially clear for situations described by proposition 4. If exactly one independent path exists between two links, the rerouting distance is ∞, while the geodesic distance is finite. Hence, the latter fails to explain why the LODF between the two links vanishes.
To further investigate the importance of distance, we simulate all possible link failures in four test grids of different size. For every failing link (r, s) we evaluate the geodesic distance as well the rerouting distance to all other links in the grid. To quantify to which extend the distance predicts the magnitude of the LODFs, we then calculate the Kendall rank correlation coefficient τ [58]. This coefficient is used on ordinal data and assumes values in the interval [−1, 1]. A value of (minus) one indicates perfect (anti)correlation, whereas a zero value implies no correlation between the data. Table I shows the results, averaging over all trigger links (r, s) in the respective grid discarding bridges. The rank correlation is negative as LODFs generally decay with distance. The magnitude of the rank correlation is significantly higher for the rerouting distance. In particular for the test grid 'case1354pegase' we see that the ordinary geodesic distance has a very limited predictive power for the LODFs (|τ | < 0.25), while the rerouting distance is strongly correlated to the magnitude of the LODFs |τ | > 0.83. Figure 9 illustrates this discrepancy in the distribution of τ values for the different distance measures for the test grids 'case118' and 'case1354pegase'. We are thus led to the conclusion that geodesic distances are of limited interest when considering the impact of link failures and should be replaced by other measures such as rerouting distances. Notably, we observe no major difference when comparing weighted and unweighted distances.

VII. CONCLUSION
Link failures represent major threats to the operation of complex supply networks across disciplines. In this article, we examined the impact of such failures in terms of the induced flow changes, which is commonly referred to as Line Outage Distribution Factors (LODFs). We provide mathematically rigorous results and extensive numerical simulations with a focus on the gross network response (i.e. the dipole strength), the scaling of flow changes with distance and the role of network topology. These quantities are crucial to understand the global robustness of supply networks as each failure can trigger a cascade of secondary failures with potentially catastrophic consequences.
First, we demonstrated rigorously that the flow changes created by a single failure in a square lattice corresponds to the field of an electromagnetic dipole. Hence the effects of a failure decay with the distance following an inverse-square law. The dipole analogy developed here allows for an analytical expression describing the spreading of link failures. Although this treatment is rigorously valid only in the continuum limit, we showed that the observed scaling extends to the other regular tilings of two-dimensional space even after removing a fraction of links. Thus, we conclude that the scaling may be expected to hold also for realistic topologies.
Increasing the sparsity of a network promotes more long-ranged effects up to the point when two links are only related by one independent pathway. Then, a rerouting between the two links becomes impossible and a failure of one link does not affect the other. However, this also implies a lack of redundancy such that a link failure can have catastrophic consequences locally. In real-world irregular networks, the gross response of a failure depends on the loading of the link as well as the local network structure. Rigorous upper and lower bounds were given for the dipole strength relating it to the redundancy of the failing link. Furthermore, the common notion of a geodesic graph distance is of limited use to predict flow rerouting. We thus introduced a rerouting distance which we showed to be much more meaningful to predict the impact of failures.
Whereas the classical analysis of link failures relies heavily on simulation results, our results provide heuristic methods and rigorous bounds which allow for an analytical insight into the relationship between the structure of a network and its robustness towards link failures. In particular for large networks where simulations are difficult, our results allow for an a priori analysis of link failures and might also be used to identify critical links, for instance in terms of the locality factor which quantifies the response of a network to a single failure. This type of analysis is aided by the general results on decay of maximal flow changes with geodesic and rerouting distances.
We expect our results to be applicable far beyond power grids since the linearized treatment extends to other phenomena such as hydraulic or biological networks. The rerouting distance along with the bounds on the locality factor may greatly simplify the study of link failures in all kinds of supply networks and makes them more accessible. We expect our results on the scaling of LODFs for networks with increasing sparsity along with this distance measure to help identifying critical parts and paths and improving the overall robustness of supply networks. the joint initiative "Energy System 2050 -a contribution of the research field energy" and the grant no. VH-NG-1025 to D.W.).

Appendix A: Proof of Proposition 1
Proof. By definition, η(r, s) is given by the flow F rs /∆P when the power ∆P is injected at node r and withdrawn at node s, while there is no injection at any other node, We can now use that the potential drop over all other links in the network is smaller than for the link (r, s) Proof. Consider first a reduced network consisting only of the link (r, s) and the shortest alternative path between the two nodes, which we denote as (j 1 = r, j 2 , j 3 , · · · , j n = s). Fixing the nodal potentials such that ψ r − ψ s = 1, the direct flow over the link (r, s) is given by whereas the indirect flow over shortest alternative path is given by Reintroducing all edges to the grid can only increase the total flow from r to s such that n F r→n ≥ F r→s + F r→j2 = b rs + 1 dist w 1 (r, s) . (B3) Thus we obtain (cf. Equation A2) η(r, s) = F r→s n F r→n Appendix C: Proof of Proposition 3 In this appendix we first give the proof for proposition 3 and then show when the decay becomes strictly monotonous.
Proof. The proof is carried out by induction starting from d = d max . We only give the proof for the maximum, the proof for the minimum proceeds in an analogous way. We assume that the network is large enough such that d max ≥ 2, otherwise the statement is trivial anyway.
(1) Base case d = d max : Consider the node n of the network for which dist(n, r) = d max and ψ n assumes its maximum ψ n = u dmax . By assumption we have dist(n, r) ≥ 2 such that the node n cannot be adjacent to the perturbed edge such that q n = 0. The n-th component of Equation and use some important properties of the matrix B: We can furthermore bound the values of ψ m in Equation (C1) by u dmax or u dmax−1 , respectively, such that we obtain (2) Inductive step d → d − 1: We consider the node n of the network with dist(n, r) = d and ψ n = u d . Starting from Equation (12) and using the same estimates as above, we obtain Note that the inhomogeneity q n ≤ 0 for all nodes except for n = r. With the induction hypothesis u d+1 ≤ u d this yields which completes the proof.