Review on Energy Resilience

Resilience of natural gas networks during conflicts, crises and disruptions

Rui Carvalho,1, ∗ Lubos Buzna,2 Flavio Bono,3, 1 Marcelo Masera,4 David K. Arrowsmith,1 and Dirk Helbing5, 6 1School of Mathematical Sciences, Queen Mary University of London, Mile End Road, London E1 4NS, U.K.

2University of Zilina, Univerzitna 8215/1, 01026 Zilina, Slovakia 3European Laboratory for Structural Assessment, Institute for the Protection and Security of the Citizen (IPSC),

Joint Research Centre, Via. E. Fermi, 2749 TP 480, Ispra 21027 (VA), Italy 4Energy Security Unit, Institute for Energy and Transport,

Joint Research Centre, Westerduinweg 3,NL-1755 LE Petten, The Netherlands 5ETH Zurich, Clausiusstrasse 50, 8092 Zurich, Switzerland

6Risk Center, ETH Zurich, Swiss Federal Institute of Technology, Scheuchzerstrasse 7, 8092 Zurich, Switzerland

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.

Keywords: networked risks—congestion control—complex networks

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 con- tracts [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

∗ [email protected]

FIG. 1. Spatial data layers involved in our analysis: population density (Landscan 2012); gas pipeline network and Liquefied Natural Gas (LNG) terminals (Platts 2011); and major urban areas (European Environment Agency and Natural Earth).

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 pro- pose a decentralized model of congestion control that dis- tributes the available network capacity to each route, without sacrificing network throughput [21–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

ar X

iv :1

31 1.

73 48

v1 [

ph ys

ic s.

so c-

ph ]

2 8

N ov

2 01

3

2

has only a local effect and the need for communication is re- duced. 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–17], cannot be solved by technology alone, but do have a relevant behavioural or social component [18–20].

RESULTS

Data set and model

Our data set is organized in four layers (see Supplemen- tary 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 spa- tial 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 Tmn.

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

LN G

Exports

Imports

Norway

Netherlands

Algeria

France

United Kingdom

Germany Austria

Libya Denm ark

G er

m an

y

Ita ly

Un ite

d Ki

ng do

m Fr

an ce

Ukr aine

Spain

Netherland s

Belarus

Belgium

PolandCzech Republic Austria

Slovakia

Ru ss

ia

Italy S

pain

Turkey Belgium

Belarus H

ungary Hungary

Portugal Greece

Ireland

Finland

Lithuania Switzerland

Rom ania

Bulgaria

M oldova

S erbia &

Latvia S

w eden

Luxem bourg

S lovenia

E stonia

D enm

ark

B o

sn ia &

H erzeg

ovin a

M acedonia

M ontenegro

0

20 k

40 k60

k

02 0k

40 k

60 k

80 k

10 0k

12 0k

140 k

160 k

180k 0

20k

40k

60k

80k

0

20k

40k

0

20k

0

20k

40k

60k

80k

0 20 k

40 k

60 k

0

20 k

40 k

0

20 k

40 k

0

20k

40k

0

20k

0

20k

FIG. 2. Natural gas imports by pipeline and via Liquefied Natural Gas (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.

Paths in Methods), we first go through each non-zero entry in the Tmn transport matrix and link each sink node in an importing country n to the Φmn = min(10, sm) 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 sm 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 Tmn of an importing country n from an exporting country m is proportional to the population of country n [24]. We next split the demand Tmn 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 net- work. 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 ci of one single congested link over the bi =

∑ρ j=1 Bi j paths that pass through

the link, where B is the link-path incidence matrix (Bi j = 1 if link i belongs to the path r j and Bi 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

3

that, for each source-sink pair, finds the path with minimum effective path length, where the effective link length is given by l̃i = (〈hi〉/hi)

α li, li is the length of link i, hi = ci/(1 + bi), 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 con- struction 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 hypo- thetically 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 Tmn that are movements of gas originating in Russia. Similarly, the fu- ture 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 sugges- tions 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 consump- tion 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 de- signed. 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 Supplemen- tary 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 alloca- tion 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 net- work as an optimizer and the proportional fairness policy as a distributed solution to a global optimization problem [29, 30].

transit

Countries

exporting

G lo

b a

lMN e

tw o

rk MT

h ro

u g

h p

u t

L N

G

Future

Present

B a

se lin

e

B e

la ru

s

L ib

ya

N o

rw a

y

P o

la n

d

R u

ss ia

U kr

a in

e

N e

th e

rl a

n d

s

A lg

e ri

a

[M ill

io n

Mm 3 /y

e a

r]

2.5

3.0

3.5

4.0

4.5

5.0 A

B

xM10 5

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.

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 Tmn transport matrix. Since the network topology and the flow network Tmn 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 conges- tion control algorithm to the resulting network and paths. We assume that all countries are willing to cooperate, that is, ad- here to the rules of the congestion control policy. To assess the effect of the range of scenarios, we then analyse the throughput

4

R us

si a

U kr

ai ne

N et

he rl

an ds

B el

ar us

LN G

N or

w ay

LN G

N or

w ay

P ol

an d

A lg

er ia

Li by

a

B a

se lin

e

P ol

an d

A lg

er ia

Li by

a

B a

se lin

e

Belgium Ukraine Latvia Greece Finland Switzerland Ireland Belarus Macedonia Bosnia & Herzegovina Estonia Lithuania Denmark Portugal Poland Serbia & Montenegro Austria Czech Republic Bulgaria Slovakia Netherlands Germany Moldova Romania Sweden Hungary Luxembourg Slovenia Italy France United Kingdom Spain

0 Country throughput per capita

[Million m3/year]

0 20

40 60

80 C

ou nt

Present

Future

Coefficient of variation

A

B0.001 0.003

1.0 0

N e th e rl a n d s

R u ss ia

U kr a in e

L N G

N o rw a y

P o la n d

L ib ya

B e la ru s

A lg e ri a

B a se lin e

B e la ru s

P o la n d

A lg e ri a

L ib ya

B a se lin e

London Ruhr Area Dusseldorf Dublin Paris Hamburg Madrid Berlin Brussels Athens Milan Rome Katowice Warsaw Stuttgart Bucharest Copenhagen Kharkiv Kiev Minsk Munich Cologne Naples FrankfurtamMain Vienna Prague Budapest Torino Belgrade Barcelona Manchester Glasgow Birmingham Leeds

0

0 1 0

3 0

5 0

C o u n t

Present

Future

Coefficientof variation

City throughput [Million m3/year]

C

D2000 4000

0.60

FIG. 4. Heat-map [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.

5

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 ex- pected 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 high- lighted 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 lo- cation); and a relatively small population (see discussion in Supplementary Information “Results”). The coefficient of vari- ation, 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 in- dicate 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 ma- jor 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 sce- narios are clustered together when either Russia, Ukraine, the Netherlands, or Belarus are removed from the network, demon- strating 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

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

1000

2000

3000

4000

5000

Present Future

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

0

10000

20000

30000

40000

50000

Present Future

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

0

500

1000

1500

2000

Present Future

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

0

100

200

300

400

500

Present Future

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

0

10000

20000

30000

40000

Present Future

0 0.2 0.4 0.6 0.8 1 Demand Reduction, β

0

200

400

600

800

1000

Present Future

Slovakia

0 10 20 30 40 50 60 70

N or

m al

is ed

Th ro

ug hp

ut [%

]Germany

0 10 20 30 40 50 60 70

N or

m al

is ed

Th ro

ug hp

ut [%

]Czech Republic

0

20

40

60

80

100

120

N or

m al

is ed

Th ro

ug hp

ut [%

]Italy

0

1

2

3

4

5

N or

m al

is ed

Th ro

ug hp

ut [%

]Ukraine

0

4

8

12

16

20

N or

m al

is ed

Th ro

ug hp

ut [%

]Austria

0

10

20

30

40

N or

m al

is ed

Th ro

ug hp

ut [%

]

N et

w or

k Th

ro ug

hp ut

[M ill

io n

m 3 /

ye ar

] N

et w

or k

Th ro

ug hp

ut [M

ill io

n m

3 / ye

ar ]

N et

w or

k Th

ro ug

hp ut

[M ill

io n

m 3 /

ye ar

]

N et

w or

k Th

ro ug

hp ut

[M ill

io n

m 3 /

ye ar

] N

et w

or k

Th ro

ug hp

ut [M

ill io

n m

3 / ye

ar ]

N et

w or

k Th

ro ug

hp ut

[M ill

io n

m 3 /

ye ar

]

A

B

C

D

E

F

FIG. 5. Network throughput of selected countries in a hypothetical crisis with Russia. The right axis shows the country throughput rela- tive 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, Fin- land, 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 de- mand 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 Nether- lands to group I countries, and these routes get a higher allocation of available capacity as Austrian demand decreases (i.e., as β decreases).

of cities into groups, highlighted in gray. We observe in Fig- ure 4D that the coefficient of variation is larger for cities in Eastern Europe than for cities elsewhere (except Berlin, Vienna and Dusseldorf). Note that Dublin is resilient to all scenarios because it is supplied from the UK, which we never removed from the network. Observe also that Austria gets most of its gas from Russia, and only a little from Norway, so Vienna is in a similar situation to Eastern European cities.

Taken together, Figures 4A–D illustrate the resilience of countries, urban areas and the network to the scenarios, by

6

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 coun- tries 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 popu- lous 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 coun- tries. 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 con- gestion 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 scenar- ios 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 sustain- ability of the gas industry makes governments and the indus- try 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 im- plement 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 pan- demics [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 ci and a length li. 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 Bi j = 1 if the link i belongs to the path r j, and set Bi 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 Bi 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 trans- port 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.

7

Proportional fairness congestion control: a primal algorithm

A decentralized algorithm for congestion control (see Sup- plementary Information “Congestion Control”) solves the sys- tem of coupled ODEs:

d dt

f j(t) = 1 − f j(t) η∑

i=1

Bi jµi(t), (1)

where the price on link i is

µi(t) = pi

 ρ∑ j=1

Bi j f j(t)

 , (2) and the price function is given by

pi(y) = max(0, y − ci + �)

�2 , (3)

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”):

d dt µi(t) =

ρ∑ j=1

Bi j f j(t) − qi(µi(t)), (4)

where

f j(t) = 1∑η

i=1 Bi jµi(t) , (5)

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 Tmn 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:

T′(NO)n = aNOT(RU)n T′(N L)n = aN LT(RU)n,

(6)

where aNO = ∑

j T(NO) j∑ j T(NO) j +T(N L) j

and aN L = ∑

j T(N L) j∑ j T(NO) j +T(N L) j

are the normalised proportions of supply from Norway and the Nether- lands, respectively.

ACKNOWLEDGMENTS

We gratefully acknowledge the support of the RAVEN EP- SRC project (EP/H04812X/1). L.B. was partially supported by VEGA (project 1/0339/13) and APVV (project APVV-0760- 11). We thank ETHZ for granting access to the Brutus high- performance cluster.

[1] Yergin D (2012) The Quest: Energy, Security, and the Remaking of the Modern World (Penguin Books, New York).

[2] Levi M (2013) The Power Surge: Energy, Opportunity, and the Battle for America’s Future (Oxford University Press, New York).

[3] The customary gas stand-off. (2009, Jan 5). The Economist. [4] Mouawad J (2005, September 29) Natural gas prices set record,

pointing to costly winter. The New York Times. [5] Chrisafis A, Borger J, McCurry J, Macalister T (2013, January

25) Algeria hostage crisis: the full story of the kidnapping in the desert. The Guardian.

[6] Plimmer G, Chazan G (2013, May 23) UK gas supply six hours from running out in march. Financial Times.

[7] Schiermeier Q (2013) Germany’s energy gamble. Nature 496:156–158.

[8] Agency IE (2012) Natural Gas Information (OECD/International Energy Agency).

[9] Wright S (2012, July 14) An unconventional bonanza. The Economist.

[10] Moniz EJ, et al. (2011) The future of natural gas, Technical report, (MIT).

[11] Helbing D (2013) Globally networked risks and how to respond. Nature 497:51–59.

[12] Hurricane damage to natural gas infrastructure and its effect on the U.S. natural gas market, (2005), Technical report, (Energy and Environmental Analysis).

[13] Voropai NI, Senderov SM, Edelev AV (2012) Detection of ”bottlenecks” and ways to overcome emergency situations in gas transportation networks on the example of the European gas pipeline network. Energy 42:3–9.

[14] Lochner S (2011) Identification of congestion and valuation of transport infrastructures in the European natural gas market. Energy 36:2483–2492.

[15] Brummitt CD, D’Souza RM, Leicht EA (2012) Suppressing cascades of load in interdependent networks. Proceedings of the National Academy of Sciences of the United States of America 109:E680–E689.

[16] Buldyrev SV, Parshani R, Paul G, Stanley HE, Havlin S (2010) Catastrophic cascade of failures in interdependent networks. Nature 464:1025–1028.

[17] Havlin S, et al. (2012) Challenges in network science: Applica- tions to infrastructures, climate, social systems and economics. European Physical Journal-Special Topics 214:273–293.

[18] Clark WC, Dickson NM (2003) Sustainability science: The emerging research program. Proceedings of the National Academy of Sciences of the United States of America 100:8059–

8

8061. [19] Ajmone-Marsan M, et al. (2012) The emerging energy web.

European Physical Journal-Special Topics 214:547–569. [20] Levi S (2013) The mathematics of sustainability. Notices of the

AMS 60:392–393. [21] Kelly FP, Maulloo AK, Tan DKH (1998) Rate control for com-

munication networks: shadow prices, proportional fairness and stability. Journal of the Operational Research Society 49:237– 252.

[22] Bertsimas D, Farias VF, Trichakis N (2011) The price of fairness. Operations Research 59:17–31.

[23] Carvalho R, Buzna L, Just W, Helbing D, Arrowsmith DK (2012) Fair sharing of resources in a supply network with constraints. Physical Review E 85:046101.

[24] Bettencourt LMA, Lobo J, Helbing D, Kuehnert C, West GB (2007) Growth, innovation, scaling, and the pace of life in cities. Proceedings of the National Academy of Sciences of the United States of America 104:7301–7306.

[25] Frischmann BM (2012) Infrastructure: The Social Value of Shared Resources (Oxford University Press, New York).

[26] Dietz T, Ostrom E, Stern PC (2003) The struggle to govern the commons. Science 302:1907–1912.

[27] Kelly F (1997) Charging and rate control for elastic traffic. European Transactions on Telecommunications 8:33–37.

[28] Srikant R (2004) The Mathematics of Internet Congestion Con- trol (Birkhäuser, Boston).

[29] Chiang M, Low SH, Calderbank AR, Doyle JC (2007) Layer- ing as optimization decomposition: A mathematical theory of network architectures. Proceedings of the IEEE 95:255–312.

[30] Kelly F, Raina G (2011) Explicit congestion control: charging, fairness and admission management, in Next-generation Inter- net: architectures and protocols, eds Ramamurthy B, Rouskas GN, Sivalingam KM (Cambridge University Press, New York).

[31] Eisen MB, Spellman PT, Brown PO, Botstein D (1998) Clus- ter analysis and display of genome-wide expression patterns. Proceedings of the National Academy of Sciences of the United States of America 95:14863–14868.

[32] Kerr RA (2012) Technology is turning U.S. oil around but not the world’s. Science 335:522–523.

[33] Levi M (2013) America’s energy opportunity. Foreign Affairs 92:92–104.

[34] Lempert RJ (2002) A new decision sciences for complex sys- tems. Proceedings of the National Academy of Sciences of the United States of America 99:7309–7313.

[35] Ostrom E (2007) A diagnostic approach for going beyond panaceas. Proceedings of the National Academy of Sciences of the United States of America 104:15181–15187.

[36] Colizza V, Barrat A, Barthelemy M, Valleron AJ, Vespignani A (2007) Modeling the worldwide spread of pandemic influenza: Baseline case and containment interventions. PLOS Medicine 4:95–110.

[37] Kirschen DS, Strbac G (2004) Fundamentals of Power System Economics (Wiley-Blackwell).

9

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 popu- lation 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 net- work 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), Be- larus, Moldova, Western Russia, Ukraine (all part of the for- mer 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 trans- mission 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) ter- minals, 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 coun- tries 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 km2.

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 Euro- pean 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 Tmnof gas transported from m to n, where m stands either for a gas exporting country or for Liquefied Natural Gas (LNG) termi- nals 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,

10

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.

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 rm,k,n,l (k = 1, . . . , Φmn and

l = 1, . . . , tn) the tn sinks in an importing country n to the Φmn closest source nodes in the exporting country m, where

Φmn =

{ min(10, sm) if m is a gas exporting country gn if m is LNG

.

(S1) 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).

11

C 5CCC RCCCC R5CCC p7CCC

Number of cities with > 100,000 inhabitants

Number

terminals of LNG

RR57p

Length of transmission network (km) Present Future

Switzerland R5Rp

Lithuania RCR9

Netherlands p85R

Latvia 7RC

Ireland 88p

United6Kingdom 78666A 8pP7

Macedonia 8G

Bulgaria R8576A pp5R

France 7pCP6A 76RP

Hungary GGpR6A G65C

Denmark 9R7 Estonia PG5 Finland GpC

Italy RpP866A RPp5P

Luxembourg R95

Austria R9596A pCC6 Belgium R57G

68P6pp6A

Czech6Republic Gp8R

Croatia RC85 Bosnia6Herzegovina 8C

Ukraine Rp69C Moldova 9Pp

Belarus pRpG

Germany p5RCC6A p6CCP

Greece 8756A R9R7

Sweden G5G Spain 8Gp56A

Slovenia P5C Slovakia p8RG

Romania GRp56A G66P Portugal 6C96A6 99p

Poland GPG8

Serbia6and6Montenegro

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.

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 Tmn 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 Tmn given by

Dmnl = 1

Φmn

ZnlTmn zn

(S2)

where Znl is the population associated with sink node l of importing country n, zn is the population of importing country n, Tmn 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 mini- mum demand on the network. To do this, we note that Dmnl is independent of k, and we replace path rm,k,n,l having de- mand Dmnl by Dmnl paths identical to rm,k,n,l, each with demand min(Dmnl), where

Dmnl = ⌊

Dmnl min(Dmnl)

⌋ , (S3)

where b·c is the largest integer not greater than ·. We note that the demand of a sink node from a source node is now

12

proportional to the number of paths connecting the source and sink pair.

The path notation rm,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 Dmnl source-sink paths identical to rm,k,n,l by the new index. In other words, for each pair of importing-exporting countries, we go through the DmnlΦmntn shortest paths that connect source to sink nodes (see Table 1), and we index all paths in increas- ing 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

∑tn l=1 DmnlΦmnT̂mn, where

T̂mn = 1 if Tmn 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 ci of a congested link i over the 1 + bi = 1 +

∑ρ j=1 Bi j paths that pass through the

link when we add a new path through i, where B is the link-path incidence matrix (Bi j = 1 if the link i belongs to the path r j and Bi j = 0 otherwise). An equitable way to divide the capacity on the link is to assign a path flow of hi = ci/(1 + bi) to each of the 1 + bi paths. Intuitively, hi is the slice of capacity allocated in a fair way to each of the bi paths that share the capacity ci, and the split is viewed as a fair outcome [S6]. Moreover, 1/hi 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:

l̃i = ( 〈hi〉 hi

)α li, (S4)

A

B

A

B

Urban Areas

A

B

1. Identification of urban nodes

2. Selection of pipelines intersecting urban polygons

3. Definition of urban nodes as urban sinks

gas node (not source)

urban node

urban sink

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.

where 〈hi〉 is the average of hi over all network links, li is the length of link i, and 0 ≤ α < 1. Whereas we weight links by their length li 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.

13

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 compet- ing paths so as to manage network congestion? There are two mechanisms at play in such allocation. On one hand, maximiz- ing 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 hi and increases all path flows by hi. Such procedure distributes parsimoniously the capacity of link i among the paths that pass through the link, and increases all unsaturated path flows by hi. The procedure then fixes the path flows of the paths that cross links with capacity ci, 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 con- tribute 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 ci of a single link among N paths in a fair way, so that each path gets a path flow of ci/N, but the two allocations are distinct when operating on a network.

Definition 1 A vector of path flows f ∗ = ( f ∗1 , . . . , f ∗ ρ ) is pro-

portionally 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]:

ρ∑ j=1

f j − f ∗j f ∗j 6 0. (S5)

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 maxi- mizes 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:

∇U( f ) · (g − f ) > U(g) − U( f ). (S6)

Now assume that f is a proportionally fair allocation. Then, ∇U( f )·(g− f ) 6 0 from equation (S5), and thus U( f )−U(g) >

14

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.

0 from equation (S6) for all other feasible g. Hence, f is a global maximum of U. Conversely, assume that the function U has one global maximum at U( f ). Then,

∇U( f ) · (g − f ) = lim t→0+

U( f + t(g − f )) − U( f ) t

≤ 0,

and thus the flow allocation is proportionally fair. �

Theorem 2 If a vector f ∗ = ( f1, . . . , fρ) of path flows is pro- portionally 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 δ = mini∈E(r j ){ci−∑ρ k=1 Bi,k fk} > 0, such that the new vector of path flows is

f ′ = ( f1, . . . , f j + δ, . . . , fρ). Hence, f ′ is not proportionally fair because

∑ρ q=1( f

′ q − fq)/ fq = δ/ 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

15

FIG. S5. Detail of the gas pipeline network in Italy, showing the presence of parallel routes.

1×10 -6

1×10 -4

1×10 -2

1×10 0

α 4.5

4.6

4.7

4.8

4.9

×10 5

G lo

b a

l N e

tw o

rk T

h ro

u g

h p

u t

[M ill

io n

m 3 /y

e a

r]

Future baseline Future (shortest path routing) Present baseline Present (shortest path routing)

α* = 0.03

FIG. S6. Plot of the global network throughput in the present and future baseline scenarios when we apply the heuristic routing algorithm of Equation (S4). Horizontal lines are a guide for the eye and show network throughput in the present and future baseline scenarios with shortest path routing. We choose the value α = 0.03 that maximizes the global throughput.

being feasible, that is:

maximize f

U( f ) = ρ∑

j=1

log( f j)

subject to B f ≤ c f j ≥ 0,

(S7)

where the link-path incidence matrix is defined by Bi j = 1 if the link i belongs to the path r j and Bi j = 0 otherwise, and c = (c1, . . . , 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]:

L( f,µ) = ρ∑

j=1

log( f j) + µ T (c − B f ) (S8)

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

16

as

∂L( f ∗,µ∗) ∂ f ∗j

= 1 f ∗j −

η∑ i=1

Bi jµ ∗ i = 0 ⇔

f ∗j = 1∑η

i=1 Bi jµ ∗ i

, (S9)

and thus

sup f

L( f,µ) = − ρ∑

j=1

log

 η∑ i=1

Bi jµi

 + η∑ i=1

µici −ρ (S10)

After removing the constant term in equation (S10) and con- verting to a maximization problem, we obtain the dual prob- lem [S10, S11]

maximize L

V (µ) = ρ∑

j=1

log

 η∑ i=1

Bi jµi

 − η∑ i=1

µici

subject to µi ≥ 0.

(S11)

The primal problem (S7) is convex and the inequality con- straints 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 may be easier to solve the primal or the dual. In our case, the primal objective function depends on ρ vari- ables (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:

µ∗i

ci − ρ∑ j=1

Bi j f ∗ j

 = 0. (S12) Equation (S12), often referred to as complementary slack- ness [S17], states that the vector µ of the Lagrange multipliers and the vector of residual capacity have complementary spar- sity patterns. To be more specific, either link i is utilized to full capacity (i.e., ci =

∑ρ j=1 Bi j f

∗ j ) and µ

∗ i > 0, or µ

∗ i = 0 and

the capacity of link i is underused (i.e., ∑ρ

j=1 Bi j f ∗ j < ci). 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 pi(y) for the use of link i, when the total load on the link is y =

∑ρ j=1 Bi 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 > ci, then we can generalize problem (S7) to replace the capacity constraints by the link cost, such that

maximize f

Û( f, p, y) = ρ∑

j=1

log( f j) − η∑

i=1

∫ yi 0

pi(z)dz

subject to B f = y f j, yi ≥ 0.

(S13) To derive the dual of problem (S13), we first find its Lagrange dual

ĝ(µ) = ρ∑

j=1

log  1∑η

i=1 Bi jµi

 −ρ+ η∑

i=1

µi p−1i (µi) − ∫ p−1i (µi ) 0

pi(z)dz

 (S14)

To simplify equation (S14), we integrate by parts and then by substitution,

ĝ(µ) = ρ∑

j=1

log  1∑η

i=1 Bi jµi

 −ρ + η∑ i=1

∫ µi 0

qi(x)d x, (S15)

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):

maximize L

V̂ (µ, q) = ρ∑

j=1

log

 η∑ i=1

Bi jµi

 − η∑ i=1

∫ µi 0

qi(x)d x

subject to µi ≥ 0. (S16)

The dual problem (S16) is equivalent to the original dual problem (S11) if qi(x) = ci. However, this function is non- invertible and thus we approximate it by the invertible func- tion [S10, S11]

qi(x) = xci

x + � . (S17)

Problems (S16) and (S11) are thus equivalent in the limit � → 0. The primal problems (S7) and (S13) are equivalent in the limit � → 0 if

pi(y) ∼ q −1 i (y) = y�/(ci − y). (S18)

17

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.

18

FIG. S8. Map of the flow allocation in the following scenarios: (F) present network without Liquefied Natural Gas (LNG) terminals, (G) present network without Norway, (H) present network without Poland, (I) present network without Russia, (J) present network without Ukraine. 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 pi(y) for the use

of link i, when the total load on the link is y = ∑ρ

j=1 Bi 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

d dt

f j(t) = 1 − f j(t) η∑

i=1

Bi jµi(t), (S19)

where the price on link i is

µi(t) = pi

 ρ∑ j=1

Bi j f j(t)

 . (S20)

19

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.

20

FIG. S10. Map of the flow allocation in the following scenarios: (F) future network without Liquefied Natural Gas (LNG) terminals, (G) future network without Norway, (H) future network without Poland, (I) future network without Russia, (J) future network without Ukraine. 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.

User i responds to an underused capacity with a steady in- crease 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 mecha- nism 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:

pi(y) = max(0, y − ci + �)

�2 . (S21)

As � → 0, the price pi 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).

S4. RESULTS

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 Eu- rope, 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 dam- ages drastically the supply of populations in Europe.

21

Table 1. Summary of notation Indexes: i Link index j Path index m Exporting country n Importing country k Node index in an exporting country l Node index in an importing country Extracted at network level: η Number of links in the network di Diameter of link i (units: in; vector with dimension η) ci = 0.56d2.5i Capacity of link i. The exponent 2.5 is found from data [S20], and the prefactor

is obtained by calibrating the present baseline scenario to match the Tmn flow matrix (units: Mm3/year; vector with dimension η)

li Length of link i (units: km; vector with dimension η) Bi j Link-path incidence matrix. Bi j = 1 if the link i belongs to the path r j and

Bi j = 0 otherwise (matrix with dimensions η×ρ) bi =

∑ρ j=1 Bi j Number of paths that pass through link i (vector with dimension η)

Extracted at country level: ν Number of countries sm Number of (source) nodes in an exporting country m tn Number of sink nodes in an importing country n gn Number of Liquefied Natural Gas (LNG) terminals in an importing country n Tmn Volume of gas received by an importing country n from an exporting country m

and LNG (matrix with dimensions: (η + 1) ×η)

T̂mn = {

1 if Tmn is positive 0 otherwise

Binary matrix with zero entries if there is no transport between countries m and n

Znl Population associated with sink node l of importing country n (matrix with dimensions ν× [Number of Voronoi and urban sinks])

zn Population of an importing country n Computed for the routing: ρ Number of paths on the network rm,k,n,l Path connecting source node k in an exporting country m with sink node l

in an importing country n, where k = 1, . . . , Φmn if m is an exporting country and k = 1, . . . , gn if m is LNG

r j Paths are also indexed by an integer, to simplify the notation used in the congestion control algorithm (vector with dimension ρ)

hi = ci/(1 + bi) Share of capacity allocated to each path passing through link i at the beginning of the heuristic rerouting

Dmnl Demand of sink l in an importing country n satisfied by an exporting country m Dmnl =

⌊ Dmnl

min(Dmnl )

⌋ Number of identical paths between a source and sink pair, each having demand min(Dmnl)

l̃i = ( 〈hi〉 hi

)α li Effective length of link i

←→ l j Effective length of path j (vector with dimension ρ)

Parameters: α = 0.03 Exponent of 〈h〉/h in equation (S4) (see also Figure S6)

Φmn =

{ min(10, sm) if m is a country gn if m is LNG

If m is a country, we connect each sink node to the min(10, sm)

geographically closest nodes in an exporting country m (distance measured along network paths); if m is LNG, we connect each sink node to the gn LNG terminals in an importing country n.

Congestion control algorithm: f j Path flow on path j (dimension ρ) µi Price on link i, and dual of f j (vector with dimension η) pi Price function on link i (vector with dimension η) qi Inverse of pi (vector with dimension η) L( f,µ) Lagrangian function of the primal problem V (µ) Lagrangian function of the dual problem

22

A. Detailed interpretation of results at country and urban levels

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 re- moved from the network. The apparent contradiction is solved by realizing that Europe’s supply from Rus- sia 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 diversifi- cation 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 remov- ing 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 Ger- man 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.

[S1] Bhaduri B, Bright E, Coleman P, Dobson J (2002) Landscan: Locating people is what matters. Geoinformatics 5:34–37.

[S2] Angel S, Sheppard SC, Civco DL (2005) The dynamics of global urban expansion, Technical report, (World Bank).

[S3] Schneider A, Friedl MA, Potere D (2009) A new map of global urban extent from modis satellite data. Environmental Research Letters 4:044003.

[S4] International Energy Agency (2011) OECD - Natural Gas imports by origin (OECD/IEA).

[S5] Bettencourt LMA, Lobo J, Helbing D, Kuehnert C, West GB (2007) Growth, innovation, scaling, and the pace of life in cities. Proceedings of the National Academy of Sciences of the United States of America 104:7301–7306.

[S6] Fehr E, Fischbacher U (2003) The nature of human altruism. Nature 425:785–791.

[S7] Carvalho R, Buzna L, Just W, Helbing D, Arrowsmith DK (2012) Fair sharing of resources in a supply network with constraints. Physical Review E 85:046101.

[S8] Schneider CM, Moreira AA, Andrade JS, Havlin S, Herrmann HJ (2011) Mitigation of malicious attacks on networks. Pro- ceedings of the National Academy of Sciences of the United States of America 108:3838–3841.

[S9] Bertsekas DP, Gallager R (1992) Data Networks, second edi- tion, (Prentice Hall, New Jersey).

[S10] Kelly FP, Maulloo AK, Tan DKH (1998) Rate control for com- munication networks: shadow prices, proportional fairness and stability. Journal of the Operational Research Society 49:237– 252.

[S11] Tan DKH (1999) Mathematical Models of Rate Control for Communication Networks, Ph.D. thesis (University of Cam- bridge Statistical Laboratory).

[S12] Guler O (2010) Foundations of Optimization, (Springer-Verlag Berlin, Berlin) Vol. 258.

[S13] Srikant R (2004) The Mathematics of Internet Congestion Con- trol (Birkhäuser, Boston).

[S14] Kelly F (1997) Charging and rate control for elastic traffic. European Transactions on Telecommunications 8:33–37.

[S15] Courant R, Hilbert D (1989) Methods of Mathematical Physics (Wiley-VCH, New York) Vol. 1.

[S16] Ball K (2008) Optimization and Lagrange Multipliers, ed Gow- ers T (Princeton University Press, New Jersey), pp 255–257.

[S17] Boyd S, Vandenberghe L (2004) Convex Optimization (Cam- bridge University Press, New York).

[S18] Jacobson V (1988) Congestion avoidance and control Proceed- ings of SIGCOMM 88, pp 314–329.

[S19] Chiu DM, Jain R (1989) Analysis of the increase and decrease algorithms for congestion avoidance in computer-networks. Computer Networks and ISDN Systems 17:1–14.

[S20] Carvalho R, Buzna L, Bono F, Gutierrez E, Just W, Arrowsmith D (2009) Robustness of trans-European gas networks. Physical Review E 80:016106.