Table Of ContentDominant transport pathways in an atmospheric blocking event
Enrico Ser-Giacomi,1 Ruggero Vasile,2 Irene Recuerda,1 Emilio Hern´andez-Garc´ıa,1 and Crist´obal L´opez1
1)IFISC (CSIC-UIB), Instituto de F´ısica Interdisciplinar y Sistemas Complejos,
Campus Universitat de les Illes Balears, E-07122 Palma de Mallorca, Spain
2)Ambrosys GmbH, Albert-Einstein-Str. 1-5 , 14473 Potsdam, Germany
(Dated: 30 July 2015)
ALagrangianflownetworkisconstructedfortheatmosphericblockingofeasternEuropeandwesternRussia
in summer 2010. We compute the most probable paths followed by fluid particles which reveal the Omega-
block skeleton of the event. A hierarchy of sets of highly probable paths is introduced to describe transport
pathwayswhenthemostprobablepathaloneisnotrepresentativeenough. Thesesetsofpathshavetheshape
of narrow coherent tubes flowing close to the most probable one. Thus, even when the most probable path is
not very significant in terms of its probability, it still identifies the geometry of the transport pathways.
5
1
PACS numbers: 92.60.-e, 47.27.ed, 89.75.Hc
0
2
l Eastern Europe and Western Russia experienced particularly in geosciences, stochastic components are
u
J astrongheatwavewithdevastatingconsequences addedtothemotionstobetterrepresentunresolvedspa-
in the summer of 2010. This was due to an atmo- tialscales16–18. Also,impreciselystatedinitialconditions
9
2 spheric blocking episode that lasted during sev- will develop into a divergent set of possible trajectories,
eral weeks. Despite these type of events have because of the inherently chaotic character of advection
] been well-investigated over the years, a complete by nearly any nontrivial fluid flow, particularly when it
h
understanding and prediction is still missing. In is time-dependent. In fact in real experiments such as in
p
this work we present a characterization of this the deployment of buoys or balloons the trajectories of
-
o flow pattern based on the study of fluid transport closely released objects diverge soon19–21. The so-called
a asaLagrangianflownetwork,sothatthemethod- spaghetti plots18 provide a visual representation of this
.
s ology of complex networks can be applied. In dispersion. Buttheybecome,whenmanytrajectoriesare
c
particular, the most probable paths linking nodes represented,clutteredandunclear. Sometypeofcluster-
i
s of this atmospheric network reveal the dominant ing or the selection of relevant trajectories is needed to
y
pathways traced by atmospheric fluid particles. highlight which are the dominant routes among a large
h
set of possible trajectories.
p
[ We have recently developed22 a formalism that com-
putes,inunsteadyflows,theoptimalfluidpathsstarting
2
I. INTRODUCTION at given initial conditions and also optimal paths con-
v
necting pairs of points. By optimal we refer to the paths
6
1 Lagrangiananalysisoftransportinfluids,inparticular which are more likely to be followed, in a well-defined
5 in geophysical and time-dependent contexts, has experi- sense made explicit below, by the fluid particles initial-
4 enced intense developments in the last decades. These ized in a finite neighborhood of the initial locations. By
0 can be roughly classified in three classes: Some of the this reason they are called most probable paths. The
1. approaches search for geometric objects –lines, surfaces, methodology builds on the set-oriented techniques8–12
0 usually related to invariant manifolds – which bound which discretize space to provide a coarse-grained de-
5 fluid regions with different properties1–3. In the second scription of transport, and draws analogies with network
1 type of approaches one computes different types of Lya- theory22–26, for which tools to compute optimal paths in
:
v punov exponents and other stretching-like fields in the graphs are well developed. A related formalism address-
i fluid domain4–7. Finally, set-oriented methods8–12 ad- ing optimal paths in time-independent flows in continu-
X
dress directly the motions of finite-size regions. ous time has been developed by Metzner et al.27. The
r
a Most of these techniques focus in identifying proper optimal paths provide the main pathways or skeleton of
Lagrangian Coherent Structures13–15,understoodasbar- the transport process in a given geographical area. Be-
riers to transport or coherent regions with small fluid cause of the implicit stochastic ingredient in the coarse-
exchange with the surrounding medium. Much less is grainingprocedureofset-orientedmethods,thismethod-
known about the actual routes of transport, the domi- ology,atvariancewithotheronesmoretiedtothetheory
nantpathways alongwhichfluidparticlestravelandfluid ofsmoothdynamicalsystems,canbeappliedequallywell
properties are interchanged. to cases of deterministic transport and to strongly diffu-
In principle, the pathways are simply given by trajec- sive situations.
tories starting from the desired initial conditions. This In this paper we compute optimal transport paths for
is true when the advection dynamics is represented by the atmospheric circulation during a blocking event oc-
a deterministic dynamical system and the initial condi- curringinSummer2010(inparticularwefocusourstudy
tion is precisely fixed. In many applications however, for the period 20th July - 30th July) over Eastern Eu-
2
ropeandRussia. Thisatmosphericflowhasverydifferent stance, using time-ordered graphs22,28.
temporal and spatial scales, and is much more diffusive, AfundamentalassumptionwemakeisthatofaMarko-
than the oceanic flow analyzed in Ser-Giacomi at al.22. vian dynamics, i.e. at each time interval the ideal fluid
We give a more detailed description of the methodology particlesareinitializedwithuniformdensityineachbox,
sketched in that reference, and generalize it to extend thus without keeping track of the trajectories at the pre-
the concept of most probable path to a hierarchy of sets vious time step. The effect of such assumption is to in-
of paths characterized by an increasing probability. The troduce diffusive effects in the dynamics even when the
spatial coherence of these sets is also discussed. original equations of motion are fully deterministic29. In
The paper is organized as follows: In Sect. II we the limit of very small boxes and very short time steps,
summarize the definition and construction of the opti- this computational diffusion is suppressed and we ap-
mal pathways as most probable paths in a flow network. proach the perfect Lagrangian motion under the given
In Sect. III we extend this concept to sets of highly velocityfield(whichitselfcancontaindiffusiveorfluctu-
probable paths and give rules to establish their signifi- ating terms).
cance and spatial coherence. Sect. IV describes the at- In our network approach spatio-temporal particle tra-
mosphericblockingevent,thedataandmodelsweuseto jectories are mapped into discretized paths between the
compute the Lagrangian trajectories, and construct the network nodes. We define an M-step path µ between
flow network from them. Sect. V contains our results: nodesI andJ astheorderedsequenceof(M+1)nodes,
optimal pathways for different dates and locations, and µ={I,k ,...,k ,J}, crossed to reach node J at time
1 M−1
also a discussion of the statistical representativeness of t starting from node I at time t . Under the Marko-
M 0
the optimal paths on the sets of highly probable paths. vian hypothesis we can associate a probability to each of
The final Section summarizes our Conclusions. An Ap- these paths as
pendix applies our formalism to a simple model flow, an
analyticdouble-gyresystem,sothatthepropertiesofthe (cid:20)M−1 (cid:21)
optimal and highly probable paths computed for the at- (pMIJ)µ =T(I1k)1 (cid:89) T(kll)−1kl T(kMM−)1J, (1)
mospheric dynamics could be more easily understood in l=2
this simplified framework.
where
A(l)
II. OPTIMAL PATHS FROM LAGRANGIAN FLOW T(l) = kl−1kl (2)
NETWORKS kl−1kl s(l) (k )
out l−1
is the probability of a fluid particle to reach node k at
Ourapproachtofindoptimalpathsintime-dependent l
time t if it was initialized at time t in node k ,
fluid flows first represents the fluid transport dynamics l l−1 l−1
asatime-dependentflownetwork26 andthenusesgraph- estimated as the ratio of the number of particles doing
so to the total number of particles released at the initial
theorytechniquestoextractfromittheseoptimalpaths.
Followingtheset-orientedmethodology8–12,26weproceed node and time. The quantity s(olu)t(k)=(cid:80)jA(klj) is called
first by a discretization of the spatial domain of interest, out-strength of node k during the l-th time step.
dividing it into N non-overlapping boxes. In terms of Among all possible M-step paths between node I and
the network-theory approach to transport22,23,26 each of J the one associated with the highest probability in
these boxes will represent a single network node. A large Eq. (1) is called the most probable path (MPP) and
number of ideal fluid particles is released in each box. is denoted by ηM. Since this path depends explicitly
IJ
Under advection by a given velocity field, links between on the number M of steps considered, it could be also
nodesareestablishedbystudyingtheLagrangiantrajec- named “fixed-time most probable path”. Its probability
tories of the particles exchanged among each pair of net- is denoted by PIMJ = maxµ{(pMIJ)µ}. To find the MPP
work nodes. This is conveniently done with a temporal and its probability we use an adaptation of the Dijk-
discretization, i.e. we consider the dynamics restricted stra algorithm30 which takes into account the layered
to a time interval [t ,t ] and divide it in time steps of and directed structure of our time-ordered flow graph.
0 M
length τ, t = t + lτ, l = 0,1...,M. For each time The simplest implementation of the algorithm would in-
l 0
interval [t ,t ] we integrate the equations of motion of volve finding maxima by searching over the full network,
l−1 l
each ideal fluid particle and keep track of each trajec- which can be a computationally expensive task. This
tory. The transport dynamics will then be described by is greatly facilitated by using the concepts of accessibil-
adjacency matrices A(l), (l = 1...M), in which a matrix ity and accessibility matrices31. Thus, for given I and
element A(l) is given by the number of particles initial- J, our implementation of the algorithm consists of two
izedattimIeJtl−1 innodeI thatendupattimetl innode mainparts. InthefirstpartonebuildsthetablesU(IlJ) of
J. Sincethevelocityfieldwillvaryintimetheadjacency nodes accessible from I and J at time step l, i.e. the set
matriceswilldependonthetimeintervalconsidered. The of nodes which can be crossed at t = tl coming from I
weightednetworkwebuildwillthereforehaveanexplicit and proceeding towards J (see Fig. 1). Technically, this
time-dependent character and can be analyzed, for in- is done by including in U(l) the nodes k for which the
IJ l
3
two following conditions are satisfied:
a) b)
(cid:20)(cid:89)l (cid:21) (cid:20) (cid:89)M (cid:21) A J A1
A(i) (cid:54)=0 and A(i) (cid:54)=0. (3) I I
i=1 Ikl i=l+1 klJ
B
Inthesecondpartofthealgorithmonerecognizesthat C A2
the structure of expression (1) allows to maximize it by
recursivelymaximizingoverk ,k ,...,k . Thisisdone 0 t1 t2 tM-1 tM 0 t1 t2 t3
1 2 M−1
by finding, for each accessible node k ∈ U(l) (and only c) d)
for them, without the need of scannling thIeJ remaining A1 A1 B3
nodes in the full network), the highest probability PIlkl I B1 I B1
of the path connecting I and k and the actual path as-
l
C
sociated. For l = 1, i.e. for the first time step, trivially A A B 1
2 2 2
we have P1 =T(1) . For l=2,3,...,M −1 we apply
Ik1 Ik1
recursively the formula 0 t1 t2 t3 0 t1 t2 t3
Pl+1 =max(cid:0)Pl T(l+1) (cid:1). (4)
Ikl+1 kl Ikl klkl+1
FIG. 1. Schematics of the algorithm to find the MPP of M
until the final point k = J is reached, and the maxi- steps between I and J. a) First part: determination of the
M
mum probability, together with the associated path, are accessible nodes. Point A is reachable from I at t = t2 but
obtained (See Fig. 1). The same procedure can then be it is not possible to reach J from it in the rest of the time
applied to any other pair of nodes (I(cid:48),J(cid:48)). interval. Point C is not reachable from I at t = tM−1 even
if J can be reached from it. Point B satisfies both accessi-
Raising the number M of steps we observe a fast
bility conditions, therefore, in contrast to points A and C, it
increase in the number of paths connecting two given
belongs to the accessibility set and it will be considered in
nodes. It is thus crucial to understand how much the
the calculation of the MPP. Systematic identification of all
MPP is representative of the large set of possible paths accessible nodes is done my applying the criteria in Eq. (3).
joining two nodes. To assess in a quantitative way this The rest of the figure illustrates the recursive maximization
issue we introduce the following quantity procedure given by Eq. (4): b) In the first time step one as-
signs to the links towards the nodes A and A (considered
1 2
λMIJ = (cid:80)P(pIMJM) , (5) taobibliteietsheT(o1n)ly onaensdinT(t1h)e ac,ceresssipbeiclittiyveslye.t Uc)(I1JF)o)rthneodperoBb-
µ IJ µ IA1 IA2 1
one considers the links from A and A , evaluates the path’s
1 2
which determines the fraction of probability carried by probabilities T(1)IA1T(2)A1B1 and T(1)IA2T(2)A2B1, and se-
the MPP with respect to the sum of probabilities of all lects the maximum one (in the figure the corresponding to
paths connecting nodes I and J. Note that the denom- the path I,A2B1, red lines). One repeats this for all nodes
inator can be simply computed as the matrix-product B1,B2,B3 in the accessibility set U(I2J) to obtain the MPPs
entry (cid:16)(cid:81)M T(l)(cid:17) . between I and these nodes, and then the procedure can be
l=1 IJ iterated again for the accessible nodes at time t3.
III. SETS OF HIGHLY PROBABLE PATHS
relevant paths. Given the initial (I) and final (J) points
we fix r nodes at intermediate times and look for paths
For large values of M, the MPP progressively loses
betweenI andJ madeofsegmentswhichareMPP’scon-
dominanceand,onaverage,doesnotcarryasignificantly
nectingtheseintermediatenodes,byusingthealgorithm
highfractionofprobability. Howeverthedynamics,char-
above. Different locations and times for these r inter-
acterized by a high number of paths connecting initial
mediate nodes are scanned and paths with probability
and final points, can be still described by a few of them,
larger than (cid:15)PM are retained and incorporated into the
which together have a non-negligible probability. To see IJ
set KM(r,(cid:15)). For (cid:15) → 1, independently on the rank (or
thiswecanrelaxthedefinitionofMPPanddefineafam- IJ
for r = 0) only the MPP is retained. KM(r = M −1,(cid:15))
ily of subsets of highly probable paths (HPP) holding IJ
contains all the paths with probability larger than (cid:15)PM.
most of the probability. In our formulation each subset IJ
However, evaluation of these sets of HPPs can be com-
KM(r,(cid:15)) is characterized by a rank 0 ≤ r ≤ M −1 and
IJ putationally costly for high values of r, since the algo-
a threshold parameter 0 ≤ (cid:15) ≤ 1. Ideally the sets would
rithm scales exponentially with r. Nevertheless interest-
contain all the paths whose probability is larger than
ingresultscanbeobtainedconsideringalreadylow-order
(cid:15)PM. Butsinceexhaustivesearchingofallsuchpathsbe-
IJ HPPs, i.e. r =1 and r =2.
comes computationally prohibitive except for very small
M, the second parameter r is introduced to determine Once one of the subsets is computed we can estab-
thenumberofconstraintsimposedinthesearchforthese lishitssignificancebydefininganextensionofexpression
4
Eq. (5):
(cid:80) (pM)
λM(r,(cid:15))= ν IJ ν, (6)
IJ (cid:80) (pM)
µ IJ µ
where the sum in the numerator is over all the paths in
the subset KM(r,(cid:15)) and the one in the denominator is
IJ
over all paths connecting I to J.
Another important aspect of the sets of HPPs is to
establish how close, spatially, are they with respect to
the corresponding MPP. This is obtained with an aver-
age distance function. Given two generic paths between
initial and final points I and J, µ = {I,k ...,J} and
1 1
µ ={I,l ...,J} we define their average distance as
2 1
M−1
1 (cid:88)
d(µ1,µ2)= M −1 d(ki,li), (7) FIG. 2. Geopotential height at 500hPa (contours, in m) and
i=1 temperature (color code, in degrees C) over the region of in-
terest, on July 24th, 12:00 UTC.
where d(k ,l ) is a metric determining the distance be-
i i
tween two given nodes of the network. For a geophysi-
cal transport network the geographical distance (on the A. Event description
sphere)betweenthecentersofthenodesisthemostnat-
ural choice. For a given pair of nodes (I,J) the average
Eastern Europe and Western Russia experienced a
distance between the subset KM(r,(cid:15)) and the MPP con-
IJ strong, unpredicted, heat wave during the summer of
necting them in M time steps is defined as
2010. Extreme temperatures resulted in over 50000
deathsandinflictinglargeeconomiclossestoRussia. The
DM = 1 (cid:88)d(µ,ηM), (8) heatwavewasduetoastrongatmosphericblockingthat
IJ NM IJ persistedovertheEuro-RussianregionfromlateJuneto
IJ µ
earlyAugust32. DuringJulythedailytemperatureswere
nearoraboverecordlevelsandtheeventcoveredWestern
where NM is the number of paths µ in the subset
IJ Russia, Belarus, Ukraine, and the Baltic nations. Physi-
KM(r,(cid:15)), and the sum is extended over all paths in the
IJ cally, the origins of this heat wave were in a atmospheric
subset(rememberthatηM denotestheMPP).Thisquan-
IJ block episode that produced anomalously stable anticy-
tity provides an estimation of how much paths in the
clonic conditions, redirecting the trajectories of migrat-
subset deviate spatially from the correspondent MPP. A
ing cyclones. Atmospheric blocks can remain in place
large deviation means that the probability to reach J
for several days (sometimes even weeks) and are of large
from I is spatially spread in a large region and indicates
scale (typically larger than 2000 km). In particular, the
furthermoretheimportanceofconsideringtheHPPsub-
Russian block of summer 2010 was morphologically of
set instead of only the MPP. Small values of DM imply
IJ the type known as Omega block that consists in a combi-
HPPsetswiththeshapeofcoherentnarrowtubesaround
nation of low-high-low pressure fields with geopotential
theMPP,sothattheMPPalreadycharacterizesthespa-
lines resembling the Greek letter Ω (see Fig.2). Omega
tial pathways, even if its probability is not large.
blocks bring warmer and drier conditions to the areas
In the next Sections we apply the above formalism to
thattheyimpactandcolder,wetterconditionsintheup-
the atmospheric flow occurring over Eastern Europe in
streamanddownstream33. Westudytheconcreteperiod
Summer 2010. Computations of optimal paths and their
extended from the July the 20th to July 30th.
sets in an analytic double-gyre system, a much simpler
flow in which path properties could be more easily ap-
preciated, are contained in the Appendix.
B. Data
Atmospheric data were provided by the National Cen-
IV. A NETWORK OF ATMOSPHERIC FLOW OVER tersforEnvironmentalPrediction(NCEP)ClimateFore-
EASTERN EUROPE IN SUMMER 2010
castSystemReanalysis(CFSR)throughtheGlobalFore-
cast System (GFS)34. This reanalysis was initially com-
In this section we describe the physical characteristics pleted over the 31 year period from 1979 to 2009 and
of the atmospheric event, the data used and the model extended to March 2011. Data can be obtained with a
we employ to obtain the air particle trajectories. temporal resolution of 1 hour and a spatial horizontal
5
resolution of 0.5◦ ×0.5◦. The spatial coverage contains
70°N
a range of longitudes of 0◦E to 359.5◦E and latitudes of
90◦S to 90◦N.
The variables needed as input to the Lagrangian dis-
persion model described in the next section include dew
pointtemperature,geopotentialheight,landcover,plan-
etary boundary layer height, pressure and pressure re- 60°N
duced to mean sea level, relative humidity, temperature,
zonalandmeridionalcomponentofthewind,verticalve-
locity and water equivalent to accumulated snow depth.
AllthesefieldsareprovidedbyCFSRdataon26pressure 50°N
levels.
40°N
0° 20°E 40°E 60°E 80°E
C. Lagragian Particle Dispersion Model FLEXPART
Asmentioned, theideaistoobtaintheeffectiveveloc- FIG. 3. The geographical domain considered and the dis-
cretization grid defining the nodes of our flow network.
ity field felt by any fluid particle. Then the Lagrangian
dispersionmodel(seenextsubsection)willintegrateitto
provide as output the three-dimensional positions of the
using a sinusoidal projection. The latitudinal extension
particle at every time step.
of each node-box is 1.5◦, the longitudinal one varies de-
The numerical model used to integrate particle veloci-
pendingonthelatitude(seeFig.3). Theareaofeachbox
tiesandobtaintrajectoriesistheLagrangianparticledis-
persionmodelFLEXPARTversion8.235,36. FLEXPART is 27722 km2, so that the typical horizontal size is of the
order of 166.5 km. This is a moderate coarse-graining
simulates the long-range and mesoscale transport, diffu-
of the resolution (0.5◦ ×0.5◦) of the NCEP data used
sion, dry and wet deposition, and radioactive decay of
for particle integration. We take τ = 12 hours as time
tracers released from point, line, area or volume sources.
discretization, which is enough to follow the dynamics of
It most commonly uses meteorological input fields from
the blocking event. It has been shown in an oceanic flow
thenumericalweatherpredictionmodeloftheEuropean
network22thatthevalueofτ hasaminorinfluenceonop-
CentreforMedium-RangeWeatherForecasts(ECMWF)
timalpaths,beingmoreimportantthetotaltime-interval
aswellastheGlobalForescastSystem(GFS)fromNCEP
considered Mτ. We uniformly fill each node with 800
(the one used in our study). Trajectories are produced
ideal fluid particles releasing them at 5000 m of height,
by integrating the equation (the input velocity data are
a representative level in the middle troposphere. FLEX-
interpolated on the present particle position):
PARTtrajectoriesarefullythreedimensional,butbyini-
dX tializing at each time-step particles in a single layer we
=v(X(t)), (9)
dt are effectively neglecting the vertical dispersion (which
is of the order of 800 m in the τ = 12 h time step) and
with t being time, X the vector position of the air par-
focussingonthepathwaysoflargescalehorizontaltrans-
ticle, and v = v¯+vt +vm is the wind vector. FLEX-
port. Fully three-dimensional flow networks will be the
PART takes the grid scale wind v¯ from the CFSR, but
subject of future work.
complements it with stochastic components vt and vm
tobettersimulatetheunresolvedturbulentprocessesoc-
curring at small scales. The turbulent wind fluctuations
vt areparametrizedbyassumingaMarkovprocessviaa V. RESULTS
Langevin equation, and the mesoscale wind fluctuations
vm are implemented also via an independent Langevin A. Optimal paths
equation by assuming that the variance of the wind at
the grid scale provides information on the subgrid vari- Equipped with the tools developed above we can now
ance. Variables entering the parametrizations are ob- compute pathways of transport during the atmospheric
tained from the meteorological CFSR fields. For addi- eventdescribedinSect. IV. Figure4ashowsalltheopti-
tional details we refer to Stohl et al.35,36. mal paths leaving a node in the Scandinavian Peninsula
at July 25 and arriving to all nodes which are reached
in M = 9 steps (i.e. 4.5 days). The graphical repre-
D. Network construction sentation joins with maximal arcs the center of the grid
boxes identified as pertaining to the MPP. The actual
We focus our analysis on the domain in between 0◦E - particle trajectories between two consecutive boxes are
80◦Eand40◦N-70◦N.Inordertodefinethenodesofthe not necessarily such arcs. The paths are colored accord-
networkwediscretizethisregionin626equal-areaboxes ingtotheirprobabilityvaluePM. TheMPPswithhigh-
IJ
6
FIG. 4. Paths of M = 9 steps of τ = 12 hours in our flow
network with starting date July 25th 2010 (panel a)) and
July 20th 2010 (panel b)), represented as straight segments
(in fact, maximal arcs on the Earth sphere) joining the path
nodes. MPPsoriginatingfromasinglenode(blackcircle)and
endinginallaccessiblenodes. ColorgivesthePM valueofthe
IJ
paths in a normalized log-scale between the minimum value
(deep blue) and the maximum (dark red). Panel a): proba- FIG.5. Optimalpathsof9stepsofτ =12hourswithstarting
bilities ranging from 10−3 to 10−14. Panel b): probabilities date July 20th 2010, entrained in the high- and in the two
ranging from 10−3 to 10−15. low-pressure areas of the blocking. Same coloring scheme as
inFig. 4. Panela): probabilitiesrangingfrom10−3 to10−16.
Panelb): probabilitiesrangingfrom10−2 to10−16. Panelc):
probabilities ranging from 10−3 to 10−13.
estprobability(reddishcolors)followadominantanticy-
clonic (i.e. clockwise) route bordering the high pressure
region (see Fig. 2, but note that this is at a particular presence of second hyperbolic structure is visible at the
time, whereas the trajectory plots span a range of dates end of the westward branch, close to 50◦N 30◦E. Figure
of more than four days) without penetrating it. There is 5 displays additional MPPs starting also at July 20th,
also a branch of MPPs with much smaller probabilities but initialized inside the main anticyclonic region of the
(yellow and bluish colors) that are entrained southward blocking,andintwolow-pressureregionsflankingit. Fig.
by a cyclonic circulation. 5a clearly shows the main anticyclonic circulation, high-
Despite the persistent character of the Eulerian block lighting also the escape routes from the high-pressure
configuration, sets of Lagrangian trajectories become zone, associated with the hyperbolic regions described
highlyvariableintime. SeeforexamplethesetofMPPs above. The other two panels show the cyclonic circula-
starting from the same initial location but five days ear- tions at each side of the high, in a characteristic Omega-
lier(Fig. 4b). Thesouthwardcyclonicbranchisnowab- blockingconfiguration. Itisremarkablethecompactness
sent, all MPPs following initially the anticyclonic gyre. of the trajectories inside the eastern low-pressure area,
Remarkably, the set of trajectories bifurcates into two which form a very localized and coherent set with prac-
branches when approaching what seems to be a strong tically no escape in the 4.5 days time-interval displayed.
hyperbolic structure close to 40◦N 60◦E. A hint of the We stress that the plots in Figs. 4 and 5 are different
7
100 7
a)
10−5 6
5
10−10
MPIJ 4
10−15 M = 6 M = 8 λf()
M = 10 3
10−20 2
0 2 4 6 8 10 12 14 16 1
path ranking x 104
0
0 0.1 0.2 0.3 0.4 0λ.5 0.6 0.7 0.8 0.9 1
FIG.6. RankingplotinwhichthePM valuesofallMPPsob-
IJ
tainedforM =6,8,and10startingonJuly25thinthewhole
areaareplottedindecreasingorder. Therangeofprobability 0.8
b)
values of the MPPs can be read from the vertical axis (from 0.7
a few percent to 10−15 for M = 6 or to less than 10−20 for
0.6
M = 10). The total number of optimal paths can also be
0.5
read-off from the horizontal axis.
λ〉 0.4
〈
0.3
fromspaghettiplots forwhichmanyavailabletrajectories
0.2
are plotted from different or related initial conditions.
0.1
For our set of particles this will give 800 trajectories em-
anating from each box. Here we are plotting just one 05 6 7 8 9 10 11 12 13
M
path, the MPP, for each initial and final box pair, which
strongly limits the number of paths from each box but,
as we will see more thoroughly, it is still representative FIG. 7. a) Normalized probability density f(λ) of the merit
of the trajectories of many released particles. figure λMIJ(r,(cid:15)) of paths started on July 20th 2010 for M =
9 and (cid:15) = 0.1, with r = 0 (only the MPPs, black curve),
r=1(blue)andr=2(green). b)MeanvalueoftheλM(r,(cid:15))
IJ
distributions (paths’ starting date July 25th) as a function
B. Relevance of the MPPs
of the number of time steps M for r = 0 (only MPPs, blue
squares), r=1 (red circles) and r=2 (single black star).
The range of colors in Figs. 4 and 5 indicates that,
given an initial box, not all MPPs leading to different
locations are equally probable. This is quantified by the
different is it from the MPPs? We plot in Fig. 8 exam-
probability PM which gives a weight to each MPP. In-
IJ plesofallHPPswithr =1and(cid:15)=0forparticular(I,J)
dead PM takes a very large range of values. Figure 6
IJ values and dates. In all the cases the sets remain coher-
shows a ranking plot in which the values of all MPPs of
entandnarrowtubesoftrajectoriesdefiningroughlythe
a given M and started at a particular date are plotted
same pathway as the MPP.
in decreasing order. We see a huge spread on the values
of PM. Very low probability values arise because of the A quantification of the width of the tubes can be done
expoInJential explosion of the number of paths between with the distance measure DIMJ in Eq. (8). An average
two nodes with increasing M. Given these low values of of it over pairs of locations is shown in Fig. 9. Although
PM except for the smallest values of M, one should ask the tube width increases with M, it remains always be-
IJ
howrepresentativearetheMPPsforthefullsetofpaths. low the typical linear box size of approximately 166.5
Figure 7a shows distributions of the parameter λM(r,(cid:15)) km (see Sect. IVD) indicating that the tubes remain
IJ
giving the relative importance of the different types of narrow. Thus we conclude that, despite the decreasing
paths. We see that λ-values are small when considering probability of the MPPs for increasing M, they remain
only the MPPs (r = 0), but the distributions shift to- good indicators of the dominant pathways in the trans-
wardshighervaluesforpathssetsofincreasingr. Figure port network.
7b gives mean values of the λ distributions. They de- As a final description of properties of the dominant
crease with M, reflecting the lack of representativeness transport paths, we present in Figure 10 (compare with
of the smallest sets of paths for large M. However, al- Fig. 8d) an example on how the MPP and the HPPs be-
ready for r =1 the set of HPPs has a mean value higher tween a fixed pair of nodes change when considering dif-
than 0.5 for a relevant range of time steps. ferent values of M, defining the temporal interval. Typ-
Thus,forthevaluesofM and(cid:15)discussedhere, theset ically, the probability of the MPP shows a maximum at
ofHPPswithr =1seemstoberichenoughtorepresent some intermediate value of M in between shorter values
the transport pathways. But how different is the geom- ofM forwhichveryfewparticlesconnectthetwonodes,
etry of the different paths in this HPP set? And how and larger values of M for which the increasing number
8
FIG. 8. All paths in KM(r = 1,(cid:15) = 0.1) for different I,J, initial point I marked by a circle and final point J marked by a
IJ
square. The color bar gives in logarithmic scale values ranging from the maximum one PM (dark red), corresponding to the
IJ
MPP, to the minimum of 0.1PM. Panel a): M = 8 steps, with starting date July 25th 2011; PM = 7.8×10−5. Panel b):
IJ IJ
M = 12 steps, with starting date July 25th 2011; PM = 2.7×10−5. Panel c): M = 11 steps, with starting date July 20th
IJ
2011; PM =7.4×10−7. Panel d): M =11 steps, with starting date July 20th 2011; PM =1.5×10−7.
IJ IJ
of factors smaller than one in the product (1) defining
160
(pM) makes this quantity to decrease again until van-
IJ µ
ishing. For the example shown in Figs. 8d and 10, the
150
value of M giving the maximum PM is around M ≈ 9,
IJ
i.e. Mτ = 4.5 days. Note that the HPP trajectories
m)140 change length but keep a similar shape in the range of
k
〉D ( M considered, indicating that in this time interval the
〈 130 blocking atmospheric structures evolve slowly.
120
VI. CONCLUSIONS
110
5 6 7 8 9 10 11 12
M
WehaveintroducedMPPsandsetsofHPPsastoolsto
FIG.9. PlotofthemeandistanceDIMJ (Eq. (8))asafunction visualize and analyze dominant pathways in geophysical
of M for r=1 and (cid:15)=0.1. The quantity is further averaged
flows. Wehavecomputedthemforanatmosphericblock-
overalltheHPPsstartingonJuly25th. Unitsarekilometers.
ing event involving eastern Europe and Western Russia.
The computed optimal paths give a Lagrangian view of
9
APPENDIX A: OPTIMAL PATHS IN A SIMPLE MODEL
SYSTEM
In this Appendix we display optimal paths and sets
of optimal paths for an analytic model flow, the double-
gyre. Seeforexample37,38 forbasicpropertiesofthissys-
tem and computations of its Lagrangian coherent struc-
tures and Lyapunov fields. Because of the simplicity
of this flow as compared with the atmospheric situation
studied in the main text, characteristics of the optimal
paths could be appreciated more easily.
The double-gyre is a two-dimensional time-periodic
flow defined in the rectangular region of the plane x =
(x,y) ∈ [0,2]×[0,1]. It is described by the streamfunc-
tion
ψ(x,y,t)=Asin(πf(x,t))sin(πy) , (A1)
with
f(x,t)=a(t)x2+b(t)x (A2)
a(t)=γsin(ωt) , (A3)
b(t)=1−2γsin(ωt) . (A4)
From these expressions, the velocity field is
∂ψ
FIG. 10. All paths in KIMJ(r = 1,(cid:15) = 0.1) for the same I,J x˙ =− =−πAsin(πf(x,t))cos(πy) (A5)
and starting date as in Fig. 8d. Same coloring scheme as in ∂y
Fig.8. Panel a) M = 7 steps; PM = 2.2×10−5. b) M = 9 ∂ψ ∂f(x,t)
IJ
steps; PM =2.3×10−4. y˙ = =πAcos(πf(x,t))sin(πy) . (A6)
IJ ∂x ∂x
For γ =0, this flow is steady. Ideal fluid particles fol-
the Omega-block configuration, with a central anti- cy- low very simple trajectories: theyrotate following closed
cloniccirculationflankedbytwocyclonicones. Moreover streamlines,clockwiseinthelefthalfoftherectangle,and
they give additional insight on it, such as the variability counterclockwiseintherightone. Thecentralstreamline
of the dominant pathways, and the identification of es- x = 1, a heteroclinic connection between the hyperbolic
caping and trapping regions. The statistical significance point at (1,1) and the one at (1,0), acts as a separatrix
of single MPPs decreases with the time interval consid- between the two regions. When γ > 0, more complex
ered, but we find always that the MPPs remain repre- behavior including chaotic trajectories arises. The pe-
sentative of the spatial geometry of the pathways, in the riodic perturbation breaks the separatrix, so that now
sense that the sets of HPPs are coherent narrow tubes someinterchangeoffluidispossiblebetweentheleftand
providing transport paths always close to the optimal the right part of the rectangle. The geometric structures
path. This spatial coherence of transport between pairs involved in this interchange have been studied with a
of locations was already noticed in an ocean flow22 and variety of techniques37,38 but the framework of optimal
it is also present in the model flow discussed in the Ap- pathsdevelopedinthispaperseemsquitenaturalforthis
pendix. Then, it seems to be a general characteristic of purpose.
flow networks. We take the parameters A = 0.1 and ω = 2π/5,
and compute paths in our network framework for two
qualitatively different situations, namely the steady case
ACKNOWLEDGMENTS γ = 0, and the periodically perturbed case (of period
2π/ω = 5) with γ = 0.3. We discretize the fluid domain
We acknowledge financial support from FEDER and into 100×50=5000 square boxes, defining the nodes in
MINECO (Spain) through the ESCOLA (CGrant no. our flow network, and compute the adjacency matrices
TM2012-39025-C02-01) and INTENSE@COSYP (Grant A(l), l = 1...M, by releasing 400 particles from each of
no. FIS2012-30634) projects, and from the Euro- theboxes. Inallthecasesshownbelowwecomputepaths
pean Commission Marie-Curie ITN program (FP7-320 of M =6 steps of duration τ =1, starting at t =0.
0
PEOPLE-2011-ITN) through the LINC project (Grant FigureA1considersthesteadyflow(γ =0)andshows
no. 289447). all optimal paths emanating from two particular initial
10
1.0 1.0
1.0 1.0
0.9 0.9
0.8 0.8
0.8 0.8
0.7 0.7
0.6 0.6 0.6 0.6
y 0.5 y 0.5
0.4 0.4 0.4 0.4
0.3 0.3
0.2 0.2
0.2 0.2
0.1 0.1
0.0 0.0
0.0 0.0
0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0
x x
FIG. A1. Double-gyre 6-step paths in the steady case γ =0, FIG. A2. 6-step paths for the periodically perturbed double
fromtwodifferentstartingnodes(theblackcircles)toallthe gyre at γ = 0.3, from three different starting nodes (black
accessibledestinations. Thenetworknodespertainingtoeach circles) to all the accessible destinations. Color coding as in
patharejoinedbystraightlinesegmentscoloredaccordingto Fig. A1, with probability ranges which are [0.0327,3.355×
the path probability. The color code is logarithmic in the 10−7](bottom-leftnode,66paths),[0.0245,1.879×10−9](top-
full probability range, which is [0.054,7.73×10−10] for the left node, 108 paths), and [0.0128,1.335×10−8] (right node,
leftnode(thereisatotalof94pathsemanatingfromit)and 106 paths).
[0.0234,10−7] for the right one (87 paths).
1.0
nodes and reaching all nodes accessible from them af-
1.0
ter the 6 steps. We see the general clockwise and anti-
clockwisecirculationsateachsideoftheseparatrix. The 0.9
0.8
two halves of the domain remain isolated. Note that the 0.8
paths are different from the closed streamlines. This is 0.7
so because the discretization of the fluid domain into fi- 0.6 0.6
nite boxes, together with the Markov assumption, intro- y 0.5
duces an stochastic component equivalent to an effective
0.4 0.4
diffusivity29 andleadstodispersionoftheparticlesstart-
0.3
ing from a single node. In our atmospheric velocity flow
0.2
therewereinadditionexplicitstochastictermsmodeling 0.2
0.1
turbulentdiffusionandmesoscalefluctuations. Notealso
that, as in the atmospheric case, a huge range of values 0.0
of PIMJ is present. 0.00.0 0.5 1.0 1.5 2.0
Figure A2 shows optimal paths for the periodically x
perturbed flow (γ = 0.3). The general clockwise and
counterclockwise rotations still remain, but now there FIG. A3. Three sets of HPPs in the double gyre for γ =
are pathways connecting the two halves of the domain. 0.3. We show all HPPs in KI6J(r = 1,(cid:15) = 0.05), with
starting nodes I at the black circles and destinations J at
Note the strong divergence of close pathways when they
the black squares. Color coding as in Figs. A1 and A2,
approach the hyperbolic region at the bottom of the do-
with probability ranges [0.0282,0.00208] (bottom-left node,
main, and how is this geometric structure what allows
12 paths), [0.0321,0.00179] (top-left node, 7 paths), and
transport of fluid between the two regions that were iso-
[0.0128,0.000946] (right node, 10 paths).
lated in the steady case.
InFig. A3wedisplaysetsofHPPsbetweenthreepairs
of nodes at γ = 0.3. More specifically we compute the
the left and right regions, using the “opening” around
paths obtained with r =1 and a probability larger than
the hyperbolic region at the top of the domain.
5% of the PM for these pairs of nodes (i.e. the paths
IJ
in the set K6 (r = 1,(cid:15) = 0.05)). The HPPs arrange in
IJ
very narrow tubes around the MPP, which is the same 1A.M.Mancho,D.Small, andS.Wiggins,“Atutorialondynam-
icalsystemsconceptsappliedtoLagrangiantransportinoceanic
behavior observed in the atmospheric paths and also in
flowsdefinedasfinitetimedatasets: Theoreticalandcomputa-
oceancalculations22. ThecentralpathinFig. A3clearly
tionalissues,”PhysicsReports437,55–124(2006).
identifies the pathway followed by particles to connect 2G. Haller and F. J. Beron-Vera, “Geodesic theory of transport
