Resilience of Natural Gas Networks during Conflicts, Crises and Disruptions

Human conflict, geopolitical crises, terrorist attacks, and natural disasters can turn large parts of energy distribution networks offline. Europe's current gas supply network is largely dependent on deliveries from Russia and North Africa, creating vulnerabilities to social and political instabilities. During crises, less delivery may mean greater congestion, as the pipeline network is used in ways it has not been designed for. Given the importance of the security of natural gas supply, we develop a model to handle network congestion on various geographical scales. We offer a resilient response strategy to energy shortages and quantify its effectiveness for a variety of relevant scenarios. In essence, Europe's gas supply can be made robust even to major supply disruptions, if a fair distribution strategy is applied.

Almost everything we do in the course of a day involves the use of energy.Yet, history has taught us that the threats to the security of supply come in unexpected ways [1,2].Examples of unforeseen energy crises include the recent disputes between Russia and Ukraine over the price of natural gas (2005-2006, 2007-2008, 2008-2009) [3], the disruption of the oil and gas production industry in the US following Hurricanes Katrina and Rita (2005) [4], the terrorist attack on the Amenas gas plant that affected more than 10% of Algerian production of natural gas (2013) [5], and the supply shortage in March 2013, when the UK had only 6 hours worth of gas left in storage as a buffer [6].New vulnerabilities could come from cyber attacks to the infrastructure [1], particularly in the case of state-driven attacks [2]; be the result of prolonged uncertainty or inaction on energy security in the US or Europe [7]; or derive from an extended period of extremely volatile prices due to intense international conflict [2].
Natural gas, a fossil fuel that accounts for 24% of energy consumption in OECD-Europe [8], has been at the heart of these crises.Gas is expensive to transport, and this is done mainly over a pipeline network.The investments are large and are made with long-term horizons, often of decades, and the costs are covered by locking buyers into long-term contracts [9].Moreover, current infrastructure investments in Europe still derive from a historical dependency on supply from Russia and North Africa [10].This dependency leaves the European continent exposed to both a pipeline network that was not designed to transport large quantities of gas imported via Liquefied Natural Gas (LNG) terminals, and to the effects of political and social instabilities in countries that are heavily dependent either on the export of natural gas (e.g., Algeria, Libya, Qatar or Russia) or its transit (e.g., Ukraine).Hence, it is challenging to build infrastructure that will be resilient to a wide range of possible crisis scenarios [11].
In a crisis, less delivery may mean greater congestion.This is due to the breakdown of major transit routes or production losses in affected areas, which cause the supply network to be used in different ways from what it was designed for.Hence, the available resources cannot be distributed well with the remaining transport capacities [12] [13,14].This is why we need a method to handle congestion.
To manage the gas pipeline network during crises, we propose a decentralized model of congestion control that distributes the available network capacity to each route, without sacrificing network throughput [21][22][23].A central controller makes the system vulnerable both to attacks on the control centre and to delays and failures of the lines of communication through the network [21].In contrast, a decentralized method is more resilient to failures because damage to the network arXiv:1311.7348v1[physics.soc-ph]28 Nov 2013 has only a local effect and the need for communication is reduced.To illustrate our model, we analyse the throughput of the present and planned pipeline networks across a range of different crisis scenarios at European, country and urban levels.The most challenging scenario corresponds to a hypothetical crisis with Russia with a complete cut-off of supply to Europe.We analyse how to alleviate the impact of such scenarios, by the identification of country groups with similar interests, which should cooperate closely to manage congestion on the network.This acknowledges that many of the 21st century challenges, such as the management of energy grids and infrastructure networks [15][16][17], cannot be solved by technology alone, but do have a relevant behavioural or social component [18][19][20].

Data set and model
Our data set is organized in four layers (see Supplementary Information "Databases"), three of which are shown in Figure 1.The first layer is the population density, which we compute from the 2012 Landscan global population data set.The second layer is the European gas pipeline network and Liquefied Natural Gas (LNG) terminals, which we extract from the Platts 2011 geospatial data set.This infrastructure is a spatial network, where nodes and links are geographically located, and links have capacity and length attributes.The third layer is defined by the urban areas in Europe with 100, 000 or more inhabitants, and we compile it from the European Environment Agency and Natural Earth.The fourth layer is the network of annual movements of gas via pipelines and of Liquefied Natural Gas (LNG) via shipping routes (see Figure 2).We represent gas flowing from an exporting country m (including LNG) to an importing country n, by a directed network with weighted adjacency matrix T mn .
Gas enters the network at source nodes, is transported over long distances on the pipeline transmission network, and then passed to the distribution network that delivers it to consumers.Here we model only the transport of gas on the transmission network.To model consumption spatially, we first need a tessellation of each country into disjoint sets of urban and non-urban areas, such that the pipeline network in an area is associated with the population it serves.Urban areas are naturally defined by the boundary of their spatial polygons.We partition non-urban areas by a Voronoi tessellation with the gas pipeline nodes as generators, respecting country borders and excluding all urban areas (see Supplementary Information "The Model").
We assume that the flow of gas on each pipeline intersecting an urban polygon (i.e., the border of the urban area) is directed towards the centre of the urban area.For simplicity, we also assume that such pipelines supply the urban area from the closest node to the urban polygon that is located inside the urban area.Moreover, each non-urban area is defined by a Voronoi cell, and we assume that it is supplied by the cell generator node (see Supplementary Information "The Model").
To connect sink to source nodes with paths (see Contract For each exporting country, we show the breakdown of the volumes of gas exported annually, together with the importing countries served.
For each importing country, we show the volumes of gas imported annually, together with the diversity of supply.
Paths in Methods), we first go through each non-zero entry in the T mn transport matrix and link each sink node in an importing country n to the Φ mn = min(10, s m ) closest nodes in an exporting country m, if m is a country, or to all LNG terminals in country n, if m is LNG, where s m is the number of gas pipeline nodes in an exporting country m.
To allocate demand to individual paths, we start with the observation that the demand T mn of an importing country n from an exporting country m is proportional to the population of country n [24].We next split the demand T mn among all source to sink paths between countries m and n, proportionally to the population served by each sink node.We now have a value of demand associated with each path, and therefore with each sink node.Finally, we replace each path by a set of identical paths, each having the minimum demand on the network.This implies that all paths have the same demand, while doubling the demand on a path is equivalent to creating two identical paths with the original demand (see Supplementary Information "The Model").
To begin integrating routing and congestion control, we first consider how to distribute the capacity c i of one single congested link over the b i = ρ j=1 B i j paths that pass through the link, where B is the link-path incidence matrix (B i j = 1 if link i belongs to the path r j and B i j = 0 otherwise), and where ρ is the number of paths on the network (see Table 1 of the Supplementary Information "The Model").To find the exact routing for these paths, we apply an iterative algorithm that, for each source-sink pair, finds the path with minimum effective path length, where the effective link length is given by l i = ( h i /h i ) α l i , l i is the length of link i, h i = c i /(1 + b i ), and α = 0.03 (see Supplementary Information "The Model").
We consider two baseline scenarios: the present and future networks.The present baseline scenario is the network that has been operational since 2011; the future baseline scenario extends the present network by the planned and under construction pipelines.To determine the network effects of crises, we analyse a range of scenarios that consist in hypothetically removing exporting (e.g., Russia) or transit (e.g., Ukraine) countries from the baseline scenarios.The scenarios are, thus, identified by the baseline (present or future) and the hypothetically removed country.For example, the present Russia scenario is given by the present network excluding Russia, that is removing all entries in the transport matrix T mn that are movements of gas originating in Russia.Similarly, the future Ukraine scenario is determined by removing all Ukrainian nodes and links from the future network.
Broadly, there are three strategies to manage congestion [25].First, expanding the network capacity is the most obvious way to lower congestion.The EU has a plan to build major pipelines crossing the continent, that should lower European dependency on Russia (see planned pipelines in Figure S1 of the Supplementary Information).Here, we include these planned pipelines in the future scenarios, but make no suggestions for extra infrastructure because the costs of expanding network capacity are high, and thus our focus is on how to best manage the existing and planned network capacity.Second, implementing congestion pricing is a way to cap the consumption of heavy users that cause network bottlenecks.Finally, by identifying groups of countries that have similar patterns of demand, we map a vast number of consumers to a relatively small number of communities that may be able to cooperate during crises [26].
We are aiming at controlling congestion in situations where the network has to perform a function for which it was not designed.For congestion control, we are using the proportional fairness algorithm (see Methods), which is inspired by the way capacity is managed on the Internet [21,27,28].The main idea behind proportional fairness is to use pricing on the links in order to control congestion (see Methods and Supplementary Information "Congestion Control").Use of non-congested links is free up to a threshold, above which the cost that a path incurs for using a link increases linearly, but steeply, with the difference between link capacity and link utilization.Hence, paths that traverse many congested links pay a high cost for contributing to congestion, and thus get a smaller flow allocation than paths that avoid congestion.A flow is proportionally fair if, to increase a path flow by a percentage ε, we have to decrease a set of other path flows, such that the sum of the percentage decreases is larger or equal to ε.We view the network as an optimizer and the proportional fairness policy as a distributed solution to a global optimization problem [29,30].3. Global network throughput by scenario.(A) A scenario is named after the country that is hypothetically removed from the network, and coloured in blue (orange) if the country is removed from the present (future) baseline scenario.(B) The country removed per scenario is coloured cyan (red) on the map, if it is an exporting (transit) country.The total network throughput increases by 6.3% from the present baseline to the future baseline scenario (i.e., when the future and planned pipelines are added to the present network).The most challenging scenarios are the hypothetical removal of Russia, followed by Ukraine, the Netherlands and LNG.When Russia is removed from the network, the global network throughput falls by 32.7% relative to the present baseline and by 28.1% in relation to the future baseline.

Simulation Results
For each scenario, we hypothetically remove the scenario country from the network and, if m is an exporting country, remove row m in the T mn transport matrix.Since the network topology and the flow network T mn depend on the scenario, we then re-compute the source-sink pairs, the demand of each pair, and we also replace every source-sink path with a number of identical paths, each having the minimum demand in the network.Finally, we apply the proportional fairness congestion control algorithm to the resulting network and paths.We assume that all countries are willing to cooperate, that is, adhere to the rules of the congestion control policy.To assess the effect of the range of scenarios, we then analyse the throughput  [31], illustrating the variation of throughput across various scenarios and the effect of a scenario on the network.The dendrograms are computed using a hierarchical clustering algorithm with the Euclidean norm and average linkage clustering.(A) Heat-map of throughput at country level across various scenarios, allowing for a comparative analysis of the present versus future baseline scenarios, as well as of crises versus baseline scenarios; (B) Coefficient of variation of throughput per capita of a country; (C) Heat-map of throughput at urban level; (D) Coefficient of variation of throughput at urban scale.The gray areas denote groups of countries and urban areas that share common patterns of throughput across scenarios.
at the scales of the European continent, countries, and of urban areas.
We compute the global network throughput, which is the sum of the throughput at all sinks (urban and non-urban), for all the scenarios.Our model reproduces successfully the expected consequences of removing the major source and transit countries from the network (see Figure 3).
We say that a country is resilient to crises if it combines high throughput per capita across scenarios with a low coefficient of variation of throughput.In addition, the network is considered resilient to a scenario if the vectors of country throughput per capita for the scenario and the baseline scenario are similar.To start addressing the resilience of countries and the network to supply and transit crises, we study the signatures in the scenario space given by the country throughput per capita in each of the 20 scenarios.Similarly, a scenario can be seen as a point in the 32-dimensional space of country throughput.The heat-map in Figure 4A shows the throughput per capita for each pair of countries and scenarios [31].
The country groups, determined by dendrograms and highlighted in gray, reflect a similar level of throughput per capita achieved across the scenarios.Countries belong to the high throughput per capita groups (highlighted in dark gray in the figure) due to a combination of effects: diversity of supply; good access to network capacity (strategic geographical location); and a relatively small population (see discussion in Supplementary Information "Results").The coefficient of variation, shown in Figure 4B for present and future scenarios, measures the normalized dispersion of country throughput per capita using the mean as a measure of scale.Larger values indicate that the throughput accessible to a country varies across scenarios.Figure 4B shows that countries in Eastern Europe have high coefficient of variation of throughput per capita in the scenarios where we hypothetically remove Russia or Ukraine.In other words, countries in Eastern Europe are still very much dependent on one single source country (Russia) and one major transit country (Ukraine).Unexpectedly, we observe a spillover effect from countries, such as Germany, which make large investments in infrastructure.These countries themselves seem to benefit less from such investments than some of their smaller neighbours.The reason behind this spillover is that countries with plentiful access to network capacity provide routes for neighbouring countries to also access such capacity.
Figure 4A can be read from left to right: the scenarios that cause the largest disruption appear on the left, and the most benign scenarios are on the right.The present and future scenarios are clustered together when either Russia, Ukraine, the Netherlands, or Belarus are removed from the network, demonstrating that the new pipelines being built will only improve slightly the consequences of a hypothetical crisis with one of the major exporting countries (Russia or the Netherlands), or with a critical transit country (Ukraine or Belarus).It is thus very hard to change the consequences of such scenarios even by building new pipelines.
We illustrate our model at a fine geographical scale in the heat-map of Figure 4C, where we show the throughput for urban areas in Europe with 1.5 million inhabitants or more, as the scenarios vary.The figure suggests possible classifications  showing how countries and urban areas with similar reactions to different types of crises are grouped together by throughput or by its coefficient of variation, and how different scenarios are clustered by their effect on the countries and urban areas.
The most challenging scenario is a hypothetical crisis that would cut-off supply from Russia to Europe.To investigate how Europe could make use of its internal gas production to minimize the impact of such a crisis, we simulate and quantify the effect of replacing gas supply from Russia with supply from Norway and the Netherlands.To do this, we start by creating two groups of countries.Group I is made of the countries that are heavily dependent on Russian gas, and is defined by Eastern Europe (http://eurovoc.europa.eu/100277)together with Estonia, Finland, Greece, Latvia and Lithuania.Group II is defined by all other countries in our study (see Supplementary Information "Databases").We consider a new scenario where Russia is removed from the network and the demand of countries in group I is rerouted to the Netherlands and Norway.To do this, we first create new paths linking each importing country in group I to Norway and the Netherlands (see Supplementary Information "The Model") and we update the matrix of gas flows to T mn (see Methods).Next, we apply a prefactor 0 ≤ β ≤ 1 to the values of the demand T mn of countries in group II.The effect of β is to lower the utilization of the network by countries of group II that do not depend heavily on Russia.These countries typically have a high value of demand, and hence by curtailing their demand, there will be more capacity available to transport gas from Norway and the Netherlands to group I countries.In Figure 5, we observe that group I countries increase their access to network capacity as β decreases.Group II countries, such as Austria, that are geographically on the main routes that link Norway and the Netherlands to group I countries, decrease their throughput as β decreases.These countries are crucial: their throughput decreases as they share their network to benefit the more populous group I countries.In contrast, access to network capacity in routes supplying group II countries, such as Germany and Italy, is broadly unaffected, even as β is lowered considerably, because routes from Norway and the Netherlands to group I countries use little network capacity from these group II countries.Despite the increase in throughput for countries in group I as β decreases, Figure 5 shows the difficulty in replacing Russia by the Netherlands and Norway.Although we can hope to recover between 40 and 50% of the baseline throughput for the Czech Republic and Slovakia, we will only recover up to 5% of the Russian supply to Ukraine and up to 20% of the Austrian supply.

DISCUSSION
Agreed political management processes are needed for crises scenarios, to guarantee supply to the most affected countries and urban areas and minimize the loss of gas by populations.Here, we propose a decentralized algorithm inspired by congestion control on the Internet, which would eliminate the need of improvisation and complicated, lengthy negotiations every time a crisis occurs.Such mechanism has a stabilizing effect because it lowers the resource deficiency of the most affected countries [11,26].We demonstrate how a wide range of scenarios impacts network throughput at global, country and urban levels, and how countries and urban areas react to scenarios of hypothetical crises.We show and quantify how countries that are heavily dependent on Russian supply can lower the impact of a crisis, if other countries accept to reduce their demand.Finally, our model tries to systematically compare alternative policy options during energy crises, using complex system models [34].
In summary, Europe is not necessarily trapped and helpless during energy crises.The long-term interest in the sustainability of the gas industry makes governments and the industry likely to invest in rules and norms to enhance reciprocity and collective efforts during crises.Because the number of governments and companies ultimately involved in taking the decisions in Europe is relatively high, governments could implement decentralized solutions similar to the one we propose here, perhaps with a centralized control solution as backup.At its heart, energy security, like preparedness for future pandemics [36], is about cooperation among nations [1].To avoid European-wide crises, nations must cooperate to share access to their critical infrastructure networks.

METHODS
Let G = (V, E, c, l) be an undirected and connected weighted graph with no loops, node-set V and link-set E = {1, . . ., η}.
Each link i has a capacity c i and a length l i .The network has a set of ρ paths connecting source to sink nodes.All links of a path transport the same path flow.Different paths can share a link, even to perform transport in different directions (e.g.,, during distinct time intervals).
The relationship between links and paths can be described by the link-path incidence matrix B as follows.Set B i j = 1 if the link i belongs to the path r j , and set B i j = 0 otherwise.Matrix B has dimensions η × ρ, and maps paths to the links contained in these paths.When B is applied to a vector of path flows, the resulting vector with components (B f ) i = ρ j=1 B i j f j is the total flow on the links, or link throughput.We say that a link is a bottleneck if the sum of the path flows of paths that pass through it is equal to the link capacity.We assume that flows are elastic, that is that path flows are determined by the available network capacity.

Contract paths
The pipeline contracts are for physical point-to-point transport on a given system over a contract path [37].The contract path is a route between a pair of source and sink nodes, such that gas flows from source to sink along that path and the transport costs are only incurred on links along that route.
Proportional fairness congestion control: a primal algorithm A decentralized algorithm for congestion control (see Supplementary Information "Congestion Control") solves the system of coupled ODEs: where the price on link i is and the price function is given by Proportional fairness congestion control: a dual algorithm Consider a system where the shadow prices vary gradually as a function of the path flows (see Supplementary Information "Congestion Control"): where and q(•) is the inverse of p(•).As → 0, the dual and primal algorithms become equivalent.
Rerouting the demand from Russia to the Netherlands and Norway When Russia is removed from the network, we reroute paths between group I countries and Russia to paths between group I countries and the Netherlands and Norway.To do this, we pair the new source and sink nodes as described in Supplementary Information "The Model", but we modify the T mn matrix of gas flows.The new T mn matrix is found by reallocating the demand from Russia for group I countries to the Netherlands and Norway, proportionally to the production of gas of these two exporting countries: where a NO = j T (NO) j j T (NO) j +T (NL) j and a NL = j T (NL) j j T (NO) j +T (NL) j are the normalised proportions of supply from Norway and the Netherlands, respectively.Supplementary Information: Resilience of natural gas networks during conflicts, crises and disruptions

S1. DATABASES
The data set is illustrated on Figures S1 and S2 and is the result of compiling GIS and population databases into several layers.

A. First layer: population
We use the 2012 LandScan [S1] high-resolution global population distribution data that estimates the population count with a spatial resolution of approximately 1 km, or 30 × 30 seconds of arc (see http://www.ornl.gov/sci/landscan/).

B. Second layer: the gas pipeline network
We compiled the European gas pipeline transmission network and the Liquefied Natural Gas (LNG) terminals from the 2011 Platts Natural Gas geospatial data (see http://www.platts.com/Products/gisdata),including pipelines that are planned or under construction.The data set covers 25 of the 27 EU member states (except Malta and Cyprus), Belarus, Moldova, Western Russia, Ukraine (all part of the former USSR), Bosnia and Herzegovina, Croatia, Macedonia, Serbia (all part of the former Socialist Federal Republic of Yugoslavia), Algeria, Libya, Morocco, Tunisia (all part of the Maghreb), Norway, Switzerland and Western Turkey.
Similarly to electrical power grids, gas pipeline networks have two layers: transmission and distribution.The transmission network transports natural gas over long distances (typically across countries) and has a non-trivial topology.The distribution network is tree-like and comprises pipelines with smaller diameter that deliver gas to consumers.We extract the gas pipeline transmission network as all the important pipelines with diameter d ≥ 15 inches.To finalize the network, we add pipelines interconnecting major branches, so that the resulting network is connected.Network links are weighted by pipeline diameter and length.To simplify, we assume that gas can flow on both directions of a pipeline, although over different time periods.The compiled network has 2, 649 nodes (compressor stations, city gate stations, Liquefied Natural Gas (LNG) terminals, etc.) connected by 3, 673 pipeline segments spanning 186, 132 km.

C. Third layer: urban areas
To avoid the controversy in the definition of an urban area [S2, Ch IV, p 49], we considered only urban areas with 100, 000 or more inhabitants as defined by the Eurostat urban audit (see http://www.urbanaudit.org).We are interested not just in the administrative boundaries of cities, but intend also to capture the surrounding areas that include a substantial share of the commuters into the city, since the gas pipeline infrastructure also supplies these peripheral urbanized districts.Note that the infrastructure network supplies directly the major urban areas, but may not intersect spatially with the built-up area of cities.
Urban areas in the European Union member countries and candidate countries are defined by Eurostat as Larger Urban Zones (http://www.urbanaudit.org),and the GIS files are provided by the European Environment Agency (http://www.eea.europa.eu/data-and-maps/data/urban-atlas).The city levels in non-EU countries are defined from remotely sensed data (see [S3] and http://www.naturalearthdata.com/downloads/10m-cultural-vectors/10m-urban-area/).These city level areas are too small compared with the EU Larger Urban Zones.Hence, we define an urban area in non-EU countries to be the union of the third-level administrative divisions (http://www.gadm.org/)that intersect the corresponding city level polygon.We have found 376 urban areas with a total area of 723, 957 km 2 .

D. Fourth layer: network of gas movements by pipeline and LNG
The fourth layer is the network of annual movements of gas by pipeline and of Liquefied Natural Gas by ship into European terminals, collected from the International Energy Agency Natural Gas Information Statistics for 2011 [S4] (see http: //www.oecd-ilibrary.org/statistics).This directed network is represented by the weighted adjacency matrix T mn of gas transported from m to n, where m stands either for a gas exporting country or for Liquefied Natural Gas (LNG) terminals that supply an importing country n (see Figure 2 of the main paper).We make use of ISO alpha-2 country codes in m and n to denote individual countries (see http://www.iso.org/iso/home/standards/country_codes.htm),so that, for example, T (RU)(FR) is the amount of gas imported annually by France from Russia.

S2. THE MODEL A. Tessellation of urban and non-urban areas and location of source and sink nodes
For simplicity, we consider that all nodes in a gas exporting country are source nodes.The partition of non-urban areas is such that all points within a given Voronoi cell are closer to their corresponding gas pipeline node than to any other node.If a gas pipeline node is inside an urban polygon, we call it an urban node, otherwise, we say the node is non-urban.To simplify, we assume that an urban area is supplied by the pipeline links that cross its border and have a node inside its polygon.This node turns out to be also the urban node on the pipeline that is the closest to the border of the urban polygon, and is thus the first node that gas will cross along the pipeline when entering the urban area.Hence, we naturally say that the node is an urban sink, and we consider no other sink nodes along the pipeline for the given urban area.If an urban area polygon contains no gas nodes, we associate it to the closest gas node (urban or not) and say this node is a sink.We exclude pipeline links that have both end nodes located inside urban areas.In other words, we only consider pipeline links that have one urban and one non-urban node (see Figures S3 and S4).

B. How we pair sink and source nodes
Source nodes are either located in an exporting country m, or at Liquefied Natural Gas (LNG) terminals.When m stands for a country, we connect by a path r m,k,n,l (k = 1, . . ., Φ mn and l = 1, . . ., t n ) the t n sinks in an importing country n to the Φ mn closest source nodes in the exporting country m, where In other words, when m is a country, we connect each sink node in an importing country n to a maximum of ten source nodes in an exporting country m.When m stands for LNG, we assume that sink nodes in an importing country n are supplied from all the LNG terminals in country n (see summary of the notation in Table 1

C. How we define demand
We define the demand of a country to be the amount of gas imported over the gas pipeline network and Liquefied Natural Gas terminals, as given by the T mn matrix (see Figure 2 of the main paper), and the demand of a given geographical area to be the demand of the country weighted by the ratio between the area and the country populations.Since demand for energy is proportional to population [S5], we locate the sink nodes and associate them with the population they supply.When the area is urban, we split its total population equally among the sink nodes inside the urban polygon.If an urban area contains no gas nodes inside its polygon, we add its population to the population associated with the closest gas node.In non-urban areas, we associate the gas pipeline node at centre of a Voronoi cell with the population of the cell.Because each sink node in an importing country n is connected by Φ mn paths to source nodes in an exporting country m, each of these paths has a share of the demand T mn given by where Z nl is the population associated with sink node l of importing country n, z n is the population of importing country n, T mn is the volume of gas imported by an importing country n from an exporting country m, and the number Φ mn of paths from an exporting country m to each sink node is given by equation (S1) (see summary of the notation in Table 1).We next express the demand of a path in units of the minimum demand on the network.To do this, we note that D mnl is independent of k, and we replace path r m,k,n,l having demand D mnl by D mnl paths identical to r m,k,n,l , each with demand min(D mnl ), where where • is the largest integer not greater than •.We note that the demand of a sink node from a source node is now proportional to the number of paths connecting the source and sink pair.The path notation r m,k,n,l has been useful so far to locate the origin and destination of the paths, but the congestion control algorithm uses matrix multiplication, and it is simpler from now on to index paths in the network by an integer.To do this, we loop through all pairs of exporting and importing countries with a non-zero entry in the T matrix and re-label each of the D mnl source-sink paths identical to r m,k,n,l by the new index.In other words, for each pair of importing-exporting countries, we go through the D mnl Φ mn t n shortest paths that connect source to sink nodes (see Table 1), and we index all paths in increasing order of first m, then n and finally k.Now that we have allocated the source to sink paths, we update the number of paths on the network ρ = ν i=1 ν j=1 t n l=1 D mnl Φ mn T mn , where T mn = 1 if T mn is positive and zero otherwise, and we write r j to denote path j, where j = 1, . . ., ρ.

D. The problem with shortest path routing
The pattern of route intersection determines how much the paths condition each other in their sharing of network links, and the capacity of links limits how much can be transported locally.If the network is not congested, transport over the geographical shortest paths minimizes the costs.In contrast, shortest path routing in congested networks can be inefficient, because it may cause congestion at a few overloaded links, while avoiding alternative routes that are only slightly longer but have higher capacity.Moreover, routing over shortest paths in gas pipeline networks makes the effect of congestion even worst.Indeed, parallel routes with similar capacity are often available, but only one of these routes is the shortest path (see Figure S5), and hence the network capacity is largely underused.

E. How we determine the source to sink paths
To begin integrating routing and congestion control, we consider first how to distribute the capacity c i of a congested link i over the 1 + b i = 1 + ρ j=1 B i j paths that pass through the link when we add a new path through i, where B is the link-path incidence matrix (B i j = 1 if the link i belongs to the path r j and B i j = 0 otherwise).An equitable way to divide the capacity on the link is to assign a path flow of h i = c i /(1 + b i ) to each of the 1 + b i paths.Intuitively, h i is the slice of capacity allocated in a fair way to each of the b i paths that share the capacity c i , and the split is viewed as a fair outcome [S6].Moreover, 1/h i can be interpreted as a simple measure of network congestion, since it has a maximum at the most congested link [S7].Hence, we combine routing and congestion through an effective link length: where h i is the average of h i over all network links, l i is the length of link i, and 0 ≤ α < 1.Whereas we weight links by their length l i in the calculation of geographical shortest paths, we now weight each link by l i in the calculation of weighted shortest paths.Thus, a link becomes less attractive (its effective length is increased) if it is more congested than the average.We find that the global network throughput is maximized for α = 0.03 (see Figure S6), and thus we use this value in the simulations.
We define the effective path length ← → l j of path j as the sum of the effective lengths of each of its links.We interpret the sum of link weights on a path as a penalty, which we then use to reroute paths iteratively via the following heuristic [S8].We i) go through each source and sink node pair and find a new path j connecting the two nodes; ii) if this new path has lower value of ← → l j than the previously found path, then it replaces the existing source to sink path; iii) we recompute the weights l i for all links on the new paths and repeat the procedure for all paths, until it has been executed 20 times (we found that the solution is does not change significantly when the number of iterations is larger than 20).

S3. CONGESTION CONTROL
How should we allocate scarce network resources to competing paths so as to manage network congestion?There are two mechanisms at play in such allocation.On one hand, maximizing the flow transported on the network may lead to some paths being assigned a zero share of network capacity, and hence zero path flow.These paths are effectively blocked from using the network, and hence the flow allocation is unfair.On the other hand, allocations that share network capacity fairly are known to deliver low throughput and are thus inefficient [S7].Hence, a good solution to the problem of congestion control aims at a trade-off between efficiency and fairness.
How should we generalize equation (S4) when paths pass through several congested links on the network?Our intuitive notion of fairness breaks down on networks, because paths typically cross several congested links and hence share the capacity of these links with other paths.Roughly, a solution is to allocate path flows iteratively, such that at each iteration we increase all path flows that do not pass through existing bottleneck links by the slice of capacity that is found by sharing equitably the capacity of the most congested links.The slice of capacity available to each path is the smallest on the most congested links, that is, the ratio h is the smallest on these links.A procedure to do this finds one link i, with the smallest ratio h i and increases all path flows by h i .Such procedure distributes parsimoniously the capacity of link i among the paths that pass through the link, and increases all unsaturated path flows by h i .The procedure then fixes the path flows of the paths that cross links with capacity c i , and decreases the capacity of links crossed by these paths by the amount of flow fixed.This creates a residual network, on which the procedure is then repeated, such that path flows are saturated and the capacity available at the links they cross is updated at each iteration.The procedure is repeated until all paths in the network are saturated.Such allocation is known as max-min fair[S7, S9], a name that comes from the way that path flows with minimum allocation are maximized by splitting equitably the capacity at the bottleneck links in an iterative process.The max-min fair allocation is such that to increase a path flow we have to decrease another path flow that is already smaller.
The major limitation of the max-min fair method is that network throughput is low compared to max-flow.To under-stand the mechanism behind this, we have to look at how both max-min fair and max-flow allocate path flows.The efficient allocation (max-flow) privileges short, over long paths that pass through several bottleneck links.Long paths take up capacity from other paths at each bottleneck, but only contribute to network throughput at the sink node.Hence, network throughput is maximized by minimizing the share of capacity to the long paths that pass through many bottlenecks, so that shorter paths can get a higher allocation of capacity and thus provide a higher contribution to network throughput.On the other hand, the max-min fair allocation shares the capacity of bottlenecks among the paths that pass through them.Thus, unlike max-flow, max-min fair allocations do not restrict the amount of network capacity that long paths can consume and are often inefficient.This limitation prompted the search for a trade-off between max-min fairness and max-flow, which would still distribute network capacity in an equitable way, and thus proportional fairness appeared in the late 1990s.

A. Proportional Fairness
Both proportional fairness and max-min fairness share the capacity c i of a single link among N paths in a fair way, so that each path gets a path flow of c i /N, but the two allocations are distinct when operating on a network.Definition 1 A vector of path flows f * = ( f * 1 , . . ., f * ρ ) is proportionally fair if it is feasible and if for any other feasible vector of path flows f , the sum of proportional changes in the path flows is non-positive [S10, S11]: If there were no capacity constraints, equation (S5) would be verified when f * j = ∞ for all j = 1, . . ., ρ.The capacity constraints imply that a flow allocation f * is proportionally fair if all other feasible vector of path flows f j = (1 + δ j ) f * j , for δ ∈ R ρ where j = 1, . . ., ρ, verify that the aggregate of percent changes ρ j=1 δ j is non-positive.Theorem 1 The unique set of feasible paths flows that maximizes the function U( f ) = ρ j=1 log( f j ) is proportionally fair.Proof.The proof given here is a direct application of the properties of convex functions [S12, S13] and global maxima of a function (a sketch of the proof is given in [S14]).First, observe that the set of feasible path flows is compact (closed and bounded) and convex.The functions log( f j ) are strictly concave, and thus U( f ) is strictly concave, since it is the sum of strictly concave functions.Thus U( f ) has a unique global maximum.Second, note that the tangent plane at any point of a convex (concave) function lies below (above) the graph of the function.Hence, since U( f ) is concave: and thus the flow allocation is proportionally fair.
Theorem 2 If a vector f * = ( f 1 , . . ., f ρ ) of path flows is proportionally fair, then each path will pass through a bottleneck.
Proof.To see this, assume that there is one path r j that does not pass through any bottleneck.Consider link i ∈ E(r j ) on the path r j .The path flow f j can be increased by δ = min i∈E(r j ) {c i − ρ k=1 B i,k f k } > 0, such that the new vector of path flows is f = ( f 1 , . . ., f j + δ, . . ., f ρ ).Hence, f is not proportionally fair because ρ q=1 ( f q − f q )/ f q = δ/ f j > 0, and the path flow f j can be increased.

B. A centralized algorithm for Proportional Fairness
Now in order to find the proportionally fair allocation, we need to maximize U( f ), constrained to the vector of path flows being feasible, that is: where the link-path incidence matrix is defined by B i j = 1 if the link i belongs to the path r j and B i j = 0 otherwise, and c = (c 1 , . . ., c η ) is the vector of link capacities.The aggregate utility U( f ) is concave and the inequality constraints are convex, and hence the optimization problem (S7) is convex.Thus, any locally optimal point is also a global optimum and we can use results from the theory of convex optimization to solve problem (S7) (see [S15] and [S16] for a brief introduction to Lagrange multipliers, and [S17] on convex optimization).The Lagrangian associated with the optimization problem (S7) is [S10, S11]: where µ = (µ 1 , . . ., µ η ) is a vector of Lagrange multipliers.The Lagrange dual function [S17] is then given by sup f L( f, µ), which is easily determined analytically by ∂L( f * , µ * )/∂ f = 0 as and thus After removing the constant term in equation (S10) and converting to a maximization problem, we obtain the dual problem The primal problem (S7) is convex and the inequality constraints are affine.Hence, Slater's condition is verified and thus strong duality holds.This means that the duality gap, i.e., the difference between the optimal of the primal problem (S7) and the optimal of the dual problem (S11), is zero [S17].Strong duality has potentially immense implications as, depending on the problem, it be easier to solve the primal or the dual.In our case, the primal objective function depends on ρ variables (the path flows) and is constrained by an affine system of equations, whereas the dual objective function depends on η variables (the links) and is constrained only by the condition that the Lagrange multipliers are non-negative.Taken together, the methods of convex optimization provide us with powerful tools to gain insights into patterns of congestion in networks where the number ρ of transport routes can be considerably larger than the number η of available transport links.Strong duality then states that the optimal path flows f * are related to the optimal Lagrange multipliers µ * by equation (S9).Since f * maximizes the Lagrangian over f , it follows that its gradient must vanish at f * , and thus the following KKT condition is satisfied: Equation (S12), often referred to as complementary slackness [S17], states that the vector µ of the Lagrange multipliers and the vector of residual capacity have complementary sparsity patterns.To be more specific, either link i is utilized to full capacity (i.e., c i = ρ j=1 B i j f * j ) and µ * i > 0, or µ * i = 0 and the capacity of link i is underused (i.e., ρ j=1 B i j f * j < c i ).This gives us a simple and powerful way to identify bottleneck links numerically, as the links with a positive Lagrange multiplier µ * i .

C. A centralized algorithm for Proportional Fairness with link price
Now suppose that the network operator charges a price per unit flow p i (y) for the use of link i, when the total load on the link is y = ρ j=1 B i j f j .This means that the price at one link depends on all the paths that pass through the link [S11].Hence, the problem (S7) can be generalized by adding a cost or penalty that is a function of the price [S10, S11].If the penalty is infinite when the link capacity is exceeded, y > c i , then we can generalize problem (S7) to replace the capacity constraints by the link cost, such that maximize ) To derive the dual of problem (S13), we first find its Lagrange dual To simplify equation (S14), we integrate by parts and then by substitution, ) where q(•) is the inverse of p(•).Following [S10, S11], we now remove the constant term in equation (S15) and covert to a maximization problem to obtain the dual of problem (S13): The dual problem (S16) is equivalent to the original dual problem (S11) if q i (x) = c i .However, this function is noninvertible and thus we approximate it by the invertible function [S10, S11] ) Problems (S16) and (S11) are thus equivalent in the limit → 0. The primal problems (S7) and (S13) are equivalent in the limit → 0 if Link thickness is proportional to the total flow on the link.Links in dark red are bottlenecks and links in blue are not used to their full capacity.The scenario country, which is hypothetically removed from the network, is highlighted in gray on the map.

D. A decentralized algorithm for Proportional Fairness
The key to a decentralized algorithm for Proportional Fairness is to translate the problem (S13) into an autonomous system of coupled differential equation, with a fixed point equivalent to the optimal solution of the optimization problem.To do this, we use the result that the stable fixed point of a system of differential equations is the maximum of the equations' Lyapunov function For each path r j , the network is offering a certain path flow f j with unit rate of change, d f j /dt = 1.Now suppose that the network operator charges a price per unit flow p i (y) for the use of link i, when the total load on the link is y = ρ j=1 B i j f j .This means that the price at one link depends on all the paths that pass through the link.The price causes a reduction in the path flow f j , such that where the price on link i is User i responds to an underused capacity with a steady increase of its path flow, and to congestion with a multiplicative decrease of its path flow at a rate proportional to the congestion price.This additive-increase/multiplicative-decrease mechanism is best known for its use in communication networks, and is implemented in TCP congestion avoidance [S18, S19].
A possible pricing policy consists in charging only for link flows that are close to capacity, with a sharp increase in the price that each path pays as the link becomes saturated: As → 0, the price p i tends to zero for link flows below capac-ity, and to infinity for saturated links.Hence, problem (S13) approximates arbitrarily closely the primal problem (S7). A. Detailed interpretation of results at country and urban levels

Figures S7, S8
At country level: • Greece receives its gas from diverse sources, and thus is resilient to the scenarios we analyse.It gets most of its gas from Russia (65.2%), from Turkey (16.4%) and LNG (18.4%); • Ireland gets its gas from the UK and is unaffected by our scenarios, since the UK is a secure source in our model; • Switzerland acts like a hub between South and Northern Europe, so it has very good access to network capacity; • Ukraine is a major transit route for gas to Europe; • Latvia and Finland have a relatively small population, good access to network capacity and are very close to Russia; • Surprisingly, Belarus does better when Ukraine is removed from the network.The apparent contradiction is solved by realizing that Europe's supply from Russia has been historically built around Ukraine.Hence, Ukrainian routes have higher capacity and shorter routes to central Europe than routes that pass through Belarus.In contrast, when Ukraine is removed from the network, routes through Belarus become the first choice to supply central Europe; • Belgium draws its high energy security from diversification of supply.It gets its gas from the Netherlands (42.5%), from Norway(48.7%,including LNG), Russia (2.8%), Germany (1.2%), and the UK (4.8%).
At urban level: • Surprisingly, Rome seems to gain slightly from removing Libya.Rome is approximately in the middle of the Italy, and the country is supplied both from the South and from the North.When Libya is removed there is no need to transport gas from the South to the North of Italy and this frees capacity to bring more gas from the North to Rome; • Unexpectedly, Berlin gains from the removal of Poland.Germany is transporting gas to Poland.When Poland is removed, the capacity that is freed can be used for German cities located close to the Polish border, including Berlin; • Finally, Dublin is resilient to all scenarios because it gets all of its supply from the UK, and we do not have any scenario affecting the ability of UK to supply gas.

FIG. 2 .
FIG. 2. Natural gas imports by pipeline and via Liquefied NaturalGas (LNG) terminals in Europe during 2011 (million cubic meters).Gas exporting (importing) countries are on the left (right) of the image.For each exporting country, we show the breakdown of the volumes of gas exported annually, together with the importing countries served.For each importing country, we show the volumes of gas imported annually, together with the diversity of supply.
FIG. 3. Global network throughput by scenario.(A)A scenario is named after the country that is hypothetically removed from the network, and coloured in blue (orange) if the country is removed from the present (future) baseline scenario.(B) The country removed per scenario is coloured cyan (red) on the map, if it is an exporting (transit) country.The total network throughput increases by 6.3% from the present baseline to the future baseline scenario (i.e., when the future and planned pipelines are added to the present network).The most challenging scenarios are the hypothetical removal of Russia, followed by Ukraine, the Netherlands and LNG.When Russia is removed from the network, the global network throughput falls by 32.7% relative to the present baseline and by 28.1% in relation to the future baseline.

FIG. 5 .
FIG.5.Network throughput of selected countries in a hypothetical crisis with Russia.The right axis shows the country throughput relative to the present baseline scenario.To minimize the impact of the loss of Russian supply, we re-allocate paths that originate in Russia to Norway and the Netherlands (see Methods).We then partition countries into two groups: group I is composed of Eastern Europe (http://eurovoc.europa.eu/100277)together with Estonia, Finland, Greece, Latvia and Lithuania, and group II includes all other countries in our study.Group II countries have a demand of βT mn , where 0 ≤ β ≤ 1. Panels (A)-(C) show the throughput for selected group I countries (open squares), whereas panels (D)-(F) illustrate the throughput for group II countries (open circles).Panels (A)-(C) demonstrate that countries in group I benefit from curtailing the demand of countries in group II.In contrast, panels (D)-(E) show that some countries in group II are largely unaffected even when their own demand is curtailed considerably.Finally, panel (F) demonstrates that supply to Austria is dominated by the demand reduction prefactor, β.Indeed, Austria is crossed by routes from Norway and the Netherlands to group I countries, and these routes get a higher allocation of available capacity as Austrian demand decreases (i.e., as β decreases).

FIG. S1 .
FIG. S1.European gas pipeline network including part of North Africa.The present network is shown in dark blue, and the planned pipelines are shown in red.The population density is plot in dark green and Larger Urban Zones are indicated in cyan.

FIG. S2 .
FIG. S2.Number of urban areas and Liquefied Natural Gas terminals, and length of the present and planned gas pipeline networks of the countries analysed.

FIG. S3 .
FIG. S3.Schematic figure illustrating the allocation of sink nodes.Urban sink nodes are shown in red, urban nodes that are not sinks are shown in yellow, and non-urban sink nodes are shown in blue.

)
FIG. S4.Urban sink nodes are highlighted in red for the Larger Urban Zones of (A) Hamburg and North West Germany, (B) Madrid, and (C) Milan.

FIG. S5 .FIG. S6 .
FIG. S5.Detail of the gas pipeline network in Italy, showing the presence of parallel routes.
FIG. S7.Map of the flow allocation in the following scenarios: (A) present network, (B) present network without Algeria, (C) present network without Belarus, (D) present network without Libya, and (E) present network without the Netherlands.Link thickness is proportional to the total flow on the link.Links in dark red are bottlenecks and links in blue are not used to their full capacity.The scenario country, which is hypothetically removed from the network, is highlighted in gray on the map.

FIG. S9 .
FIG. S9.Map of the flow allocation in the following scenarios: (A) future network, (B) future network without Algeria, (C) future network without Belarus, (D) future network without Libya, and (E) future network without the Netherlands.Link thickness is proportional to the total flow on the link.Links in dark red are bottlenecks and links in blue are not used to their full capacity.The scenario country, which is hypothetically removed from the network, is highlighted in gray on the map.
Figures S7, S8, S9 and S10 show the load on network links in the present and future scenarios.These figures demonstrate that our model reproduces the main transport corridors in Europe, and show the spatial pattern of bottleneck links for each scenario.It is apparent how a hypothetical removal of either Russia or Ukraine cuts-off the major transport routes, and damages drastically the supply of populations in Europe. ).
RR57pLength of transmission network (km)