# Biased random walkers and extreme events on the edges of complex networks

Govind Gandhi and M. S. Santhanam  Indian Institute of Science Education and
Research, Dr. Homi Bhabha Road, Pune, India - 411008

###### Abstract

Extreme events have low occurrence probabilities and display pronounced
deviation from their average behaviour, such as earthquakes or power
blackouts. Such extreme events occurring on the nodes of a complex network
have been extensively studied earlier through the modelling framework of
unbiased random walks. They reveal that the occurrence probability for extreme
events on nodes of a network has a strong dependence on the nodal properties.
Apart from these, a recent work has shown the independence of extreme events
on edges from those occurring on nodes. Hence, in this work, we propose a more
general formalism to study the properties of extreme events arising from
biased random walkers on the edges of a network. This formalism is applied to
biases based on a variety network centrality measures including PageRank. It
is shown that with biased random walkers as the dynamics on the network,
extreme event probabilities depend on the local properties of the edges. The
probabilities are highly variable for some edges of the network, while they
are approximately a constant for some other edges on the same network. This
feature is robust with respect to different biases applied to the random walk
algorithm. Further, using results from this formalism, it is shown that a
network is far more robust to extreme events occurring on edges when compared
to those occurring on the nodes.

††preprint: IISER/123

##  I Introduction

Extreme events arise in many situations of practical interest [1]. The most
common ones are the extremes witnessed in geophysical phenomenon such as the
earthquakes, cyclones and drought. These represent extreme events recorded in
a univariate time series. However, networks have emerged as a standard
modeling paradigm to describe such complex systems [2, 3]. For instance,
traffic jams, large scale power blackouts, call drops in cellular phone
networks, overloaded web servers are a result of large fluctuations beyond the
capacity of the system to service the requests. These events take place,
respectively, on a network of roads, power distribution network, cellular
network and on a network of routers. In general, any large excursion of a
dynamical variable from its typical behaviour can be termed an extreme event.
Though such events typically have low probability of occurrence, they are
regarded as significant due to the social and financial losses suffered on
account of most of these extreme events.

While the classical extreme value theory [4], formulated a century ago, deals
with extreme events primarily in univariate time series, the study of extreme
events on complex networks has attracted attention only since the last decade
[5, 6, 7]. In the latter case, a typical setting is that of a topology of a
complex network, a system of nodes connected through edges [8]. In order to
discuss events and extreme events, a dynamical process needs to be
additionally defined on the networks. If we then let $x(t)$ represent the flux
of any dynamical variable on a node of a network, then an event is deemed to
be an extreme event if $x(t)>x_{\text{th}}$, where $x_{\text{th}}$ is a
suitable threshold that indicates, for instance, the processing capacity of a
node. An early work along this direction used non-interacting random walks as
a dynamical process [9]. By assuming $x_{\text{th}}$ to be proportional to the
typical random walker flux on a node, it was demonstrated that the probability
that an extreme event occurs on that node depends on the degree of the node in
question. In particular, quite contrary to intuition, it was observed that, on
an average, hubs have low probability for occurrence of extreme events, while
small degree nodes have higher likelihood for extreme events [9]. Over the
years, this proposition was observed in many different settings. For instance,
this feature been theoretically reported in a network of coupled chaotic maps
[10], in stochastic processes such as Brownian motion in an external potential
[11], and in biased random walks on networks [12]. However, it appears that in
the case of interactions among agents, the extreme events created by the
agents can possibly have a different profile for the probability for the
occurrence of extreme events [13, 14]. The formulation introduced in Ref. [9]
has been used to study network failures as well [15, 16]. In general, the
focus of these works has been on extreme events, specifically on the nodes of
a network.

A natural extension of these works is to consider extreme events occurring at
the interface where nodes interact - the edges of the network. Recently, using
random walks as flows on a network, extreme events occurring on edges and
their properties were studied [17]. Let $i$ and $j$ be two nodes connected by
an edge $e_{ij}$. It was shown that extreme events on nodes $i$ and $j$, and
those on the edge $e_{ij}$ connecting the nodes are nearly uncorrelated and
that they should be regarded as events which do not influence one another
[17]. Hence results for the occurrence of extreme events on nodes do not
characterise those taking place on the edges. Indeed, most of our personal
experience in the case of road traffic tells us that traffic jams on junctions
and at any other point on the roads are not always correlated. Furthermore,
with unbiased walkers occurrence probability of extreme events on edges
depends only on the gross network parameters such as total number of edges and
not on the local details in the vicinity of the edge. This feature differs
appreciably from what we know about extreme events on nodes. Motivated by
these considerations we study extreme events on the edges using the more broad
dynamics of biased random walkers on complex networks. The standard random
walk provides a simple and yet non-trivial setting to explore this question.
In the rest of this paper, a general formalism is developed to study biased
random walks on networks, and it is applied to the degree biased random walks
on scale-free networks, and also to other centrality measure based biases such
as eigenvector centrality, betweenness centrality and page-rank centrality. It
is also applied to a real-life planar street network. In all the cases studied
here, extreme events arise from inherent fluctuations in the number of
walkers, and provides insights into how inherent fluctuations, as opposed to
external shocks, give rise to extreme events.

We describe the formalism in sections II, III, IV. In section V we apply it to
a planar network and subsequently to several different types of biased random
walkers. Finally, in section VII we discuss the robustness of edges against
occurrence of extreme events.

##  II Biased Random Walks on networks

A connected, undirected network with $N$ nodes and $E$ edges is considered.
Its connectivity structure is described by an adjacency matrix ${\mathbf{A}}$
whose elements $a_{ij}$ are either $1$ or $0$ depending on whether an edge
connects nodes $i$ and $j$ or not. On this network structure, diffusion of $W$
independent random walkers is considered whose dynamics is biased by any
network property that can be associated with the node, such as the degree.
This is represented by a single-step transition probability $\pi_{ij}$ for a
walker to go from node $i$ to node $j$, while satisfying the normalization
condition $\sum_{j}\pi_{ij}=1$ for all $i$. In this work, we study a class of
biased random walks defined by the one-step transition probability

$\pi_{ij}=\frac{a_{ij}f_{j}}{\underset{m}{\sum}a_{im}f_{m}},\;\;\;\;i,j=1,2,\dots,N..$
(1)

The biasing function $f_{j}$ here, could represent any structural or non-
structural property of the network as long as it is time-independent and is
associated with a node. The structural properties could simply be the degree
of a node, the betweenness centrality, it could be any amalgam of structural
features of the node and its local neighbourhood, it could be the ground-truth
community the node is externally assigned, and so much more. This can be
generalized to the case of multiplex networks as well [18]. In the case of
unbiased random walk $f_{j}=1$, for all $j$, indicating that the probability
for a walker to jump to any of its neighbouring nodes is a constant.

Let $P_{ij}(t)$ denote the probability that a walker at node $i$ starting at
time $t=0$ reaches node $j$ after $t$ discrete jumps. The time evolution of
probability is given by the master equation

$P_{ij}(t)=\sum_{k=1}^{N}P_{ik}(t-1)~{}\pi_{kj}.$ (2)

By allowing only time independent biasing functions we are assured transition
probabilities $\pi_{ij}$ that are stationary and that any node in the network
is accessible from any other node via a finite number of steps. It then
follows from Perron-Frobenius theorem [19] that for any such biasing function
$f_{j}$, the stationary probability distribution $p_{i}$ exists. This time-
independent occupation probability on node $i$ turns out to be [20]:

$p_{i}=\frac{f_{i}\sum_{j}a_{ij}f_{j}}{\underset{k}{\sum}f_{k}\left(\underset{j}{\sum}a_{kj}f_{j}\right)},$
(3)

in which the normalization condition $\sum_{j}p_{j}=1$ holds. Note that the
occupation probability $p_{i}$ does not restrict the form of the biasing
function or the network structure at this stage, thereby making it
sufficiently general.

##  III Extreme Events on nodes

The occupation probability in Eq. 3 can be used to evaluate extreme event
probabilities on nodes. The probability of finding $w>1$ non-interacting
walkers on the node labelled $i$ is $p_{i}^{w}$ while the rest of the $W-w$
walkers are distributed among the other nodes of the network. The distribution
of number of walkers on node $i$ is then a binomial distribution given by [9]

$\mathcal{P}_{i}(w)=\binom{W}{w}~{}p_{i}^{w}~{}\left(1-p_{i}\right)^{W-w}.$
(4)

Intuitively, an extreme event represents a pronounced deviation from the mean
flux. The flux $w$ represents events in our framework and it is deemed to be
an extreme event if $w\geq q_{i}$, where $q_{i}$ represents a suitable
threshold. Then, using Eq. 4, the EE probability on $i$-th node can be defined
as

$\mathcal{P}_{EE}(q_{i})=\sum_{w=q_{i}}^{W}\mathcal{P}_{i}(w).$ (5)

The threshold $q_{i}$ is chosen to be proportional to the natural variability
of the flux passing through the node. As was done earlier in Ref. [12], in
this work too the threshold for identifying EE is taken to be

$q_{i}=\langle w_{i}\rangle+m\sigma_{i},$ (6)

where $m\geq 0$. The average and the variance of the flux passing through
$i$-th node are given as

$\displaystyle\langle w_{i}\rangle$ $\displaystyle=Wp_{i},$ (7a)
$\displaystyle\sigma_{i}^{2}$ $\displaystyle=Wp_{i}\left(1-p_{i}\right).$ (7b)

Since these quantities depend only on the biasing function and the local
structural information about immediate neighbours of nodes, these results will
hold good for any network irrespective of its degree distribution. This notion
of EE is sufficiently general to encompass systems that accommodate
congestion, e.g., airplane/vehicular traffic in a network of airports/streets.

##  IV Extreme Events on Edges

While extreme events on nodes have been studied earlier, there has not been
much work about EE on the edges of a network. As was recently shown in Ref.
[17], the occurrence of EE on nodes and its edges are nearly uncorrelated.
Hence, it is not possible to infer about EE on edges/nodes based on the
information about EE on nodes/edges. Physically, this is intuitive because
once an EE takes place on a node, the walkers get scattered into a large
number of edges and might not necessarily lead to an EE on the edges. These
events on the edges will depend on the flux through edges and the carrying
capacity of the edges. For instance, in the water supply network, the
reservoirs might have ample capacity to hold the water, though its transport
channels might have a smaller load-bearing capacity. A microprocessor has a
network of logic gates, where the information is carried as electrical
signals. In all these cases, the conduits (edges) can only support a finite
range of fluctuations [21], beyond which it should be designated as an extreme
event. In most cases, the occurrence of EE on the edges can lead to disasters.
Therefore it is imperative that EE on edges be studied on its own right.

###  IV.1 Load and Flux Distribution on edges

Firstly, we seek quantifiable definitions of the situations described above,
i.e., for load and flux. Consider an edge $e_{ij}$ and the nodes it connects
labelled $i$ and $j$. Suppose there are $w_{i}$ walkers on node $i$ and
$w_{j}$ walkers at node $j$. At the $t$-th time step, let $l_{i}$ out of the
$w_{i}$ walkers jump into node $j$, and $l_{j}$ out of the $w_{j}$ walkers
jump to node $j$. The load $l_{ij}$ and flux $f_{ij}$ of the edge $e_{ij}$ at
time $t$ is then defined as

$l_{ij}(t)=l_{i}(t)+l_{j}(t)\,\,\,\,\,\,\,\,\,\,\,\,\,f_{ij}(t)=l_{i}(t)-l_{j}(t)$
(8)

Since the load $l_{ij}$ is stochastic, a natural question of interest is its
probability distribution for a particular edge $e_{ij}$. The details of
obtaining this probability distribution of load for a given edge $e_{ij}$ is
given in the appendix A. The normalized probability distribution comes out to
be

$P_{L}(l_{ij})=\binom{W}{l_{ij}}~{}(2p_{i}\pi_{ij})^{l_{ij}}~{}(1-2p_{i}\pi_{ij})^{W-l_{ij}}.$
(9)

Comparing this with Eq. 4, an effective occupation probability on edge
$e_{ij}$ can be identified as $p_{ij}^{\text{eff}}=2p_{i}\pi_{ij}$.
Physically, this corresponds to the product of the probability of a walker to
be on node $i$ and probability for that walker to move into node $j$.
Therefore, the distribution of load on an edge is dependent on the product of
the biasing functions of the two nodes involved in the edge. It is dependent
on the local network structure and not on the global network parameters. The
mean load and the variance, respectively, on the edge $e_{ij}$ are

$\displaystyle\langle l_{ij}\rangle$
$\displaystyle=2Wp_{i}\pi_{ij}=Wp_{ij}^{\text{eff}}=\frac{2Wa_{ij}f_{i}f_{j}}{{\sum_{l}}({\sum_{k}}(a_{lk}f_{k})f_{l}},$
(10a) $\displaystyle\sigma_{ij}^{2}$
$\displaystyle=Wp_{ij}^{\text{eff}}(1-p_{ij}^{\text{eff}})$ (10b)

In the limit that $p_{ij}^{\text{eff}}\ll 1$, then it is easy to see that
$\sigma_{ij}^{2}\approx\sqrt{\langle l_{ij}\rangle}$, a result known in the
case of unbiased random walks on networks.

For the probability that a node $i$ has $l_{i}$ walkers, we make a similar
calculation done in Appendix A to get

$\mathcal{L}(l_{i};\,i)=\binom{W}{l_{i}}~{}(p_{i}\pi_{ij})^{l_{i}}~{}(1-p_{i}\pi_{ij})^{W-l_{i}}$
(11)

The flux distribution i.e., the probability that the flux equals $f_{ij}$ on
an edge $e_{ij}$ is then product of probabilities of getting different values
$l_{i}$ and $l_{j}$ constrained by Eq. 8. We sum over all such possibilities
to get the flux distribution on an edge to be

$P_{F}(f_{ij})=\sum_{m=0}^{W-f_{ij}}\mathcal{L}(l_{i}=m+f_{ij};\,i)~{}\mathcal{L}(l_{j}=m;\,j).$
(12)

Using the load (Eq. 9) and flux distributions (Eq. 12), we can now compute the
extreme events on edges for a variety of networks.

###  IV.2 Probability of EE on edges for load

Similar to the definition of EE on nodes, if the load on an edge at any time
$t$ exceeds a defined threshold $q_{ij}$ such that $l_{ij}(t)\geq q_{ij}$,
then it indicates the occurrence of load EE on edge $e_{ij}$. As in the case
of nodal extreme events, $q_{ij}$ is taken to be proportional to the natural
variability of the load on the edge. Thus,

$q_{ij}=\langle l_{ij}\rangle+m\sigma_{ij},$ (13)

where $m\geq 0$, and the mean load $\langle l_{ij}\rangle$ and standard
deviation $\sigma_{ij}$ are given by Eqs. 10a and 10b respectively which
inform us that $q_{ij}$ is a function of $p_{ij}^{\text{eff}}$. Hence, the
probability for load EE on an edge $e_{ij}$ also depends on
$p_{ij}^{\text{eff}}$. Then by using Eq. 9 and Eq. 10a we obtain for
probability of load EE on an edge $e_{ij}$

$Q_{L}\left(\langle l_{ij}\rangle\right)=\sum_{l=q_{ij}}^{W}P_{L}(l).$ (14)

The probability for the occurrence of flux EE, $Q_{F}\left(\langle
f_{ij}\rangle\right)$, on an edge is quite similar to Eq. 14, with $P_{L}(l)$
replaced by flux distribution $P_{F}(f)$ from Eq. 12. In general, the EE
probabilities on edge $e_{ij}$ depend only on the mean load or flux, which in
turn carries information about local network structure and strength of bias in
random walk dynamics. Furthermore, $\langle l_{ij}\rangle$ is characteristic
of an edge irrespective of the network topology, and allows us to compare EE
probabilities for different biases and even between networks on a common
footing. This accessible quantity is a candidate for an observable when
studying real life systems.

We are currently in a position to take our pick of biased random walk dynamics
and a network structure and obtain information on the EE probabilities on
edges and nodes for congestion-type systems. In the rest of the paper, this
framework is tested and compared against simulation results for different
random walk dynamics and network topologies.

##  V Biased Random Walks

Figure 1: The 1.5 km2 planar street network around IISER Pune with 1010 nodes
and 1306 edges. The nodes in the two-dimensional plane denote street
intersections and the edges represent streets.

###  V.1 Degree Bias

In this section, we shall consider one of the most applicable biases - the
degree biased random walk [22, 23, 24]. Its biasing function is of the form

$f_{i}=k_{i}^{\alpha},$ (15)

where $k_{i}=\sum_{j}a_{ij}$ is the degree of $i$-th node. The exponent
$\alpha$ is the parameter that determines the strength of bias imparted to the
walkers. Unbiased random walk corresponds to $\alpha=0$, whereas for
$\alpha>0$ the walkers are biased towards hubs (nodes with high degrees), and
for $\alpha<0$ they are preferentially biased towards small degree nodes.
$|\alpha|$ determines the strength of the bias towards or away from the hubs.
On networks, biased random walks of this type had been studied earlier to
enumerate the characteristic time scales, namely, the recurrence time to the
starting point, cover time, and time taken to visit a node for the first time
[25, 26]. The parameter $\alpha$ can be tuned so as to minimize these time
scales. In general, this results in smaller time scales with biased random
walks in comparison to the unbiased random walks [27]. The local navigation
rules for the biased walks are determined by the transition probability (see
Eq. 1)

$\pi_{ij}=\frac{a_{ij}~{}k_{j}^{\alpha}}{\sum_{m}a_{im}~{}k_{m}^{\alpha}}.$
(16)

The stationary occupation probability follows from Eq. 3, and is

$p_{i}=\frac{k_{i}^{\alpha}~{}\underset{m}{\sum}a_{im}k_{m}^{\alpha}}{\underset{l}{\sum}~{}[k_{l}^{\alpha}(\underset{m}{\sum}a_{lm}k_{m}^{\alpha})]}.$
(17)

The index in each of the summations runs over all nodes.

The results are first obtained for a (planar) street network [28] shown in
Fig. 1. The planarity of the network was confirmed using Kuratowski’s theorem
[29]. Fig. 2 displays the probabiity distribution for load and flux obtained
from biased random walk simulations for an arbitrary edge $e_{0,11}$
connecting nodes labelled 0 and 11. For the purposes of this simulation,
$W=2500$ non-interacting walkers executed a degree biased random walk (DBRW)
on the network in Fig.1 for 100,000 time steps. The results shown in Fig. 2
represent an average over 15 realisations with a fixed network, but with
varying initial positions of the walker for each realisation. Evidently, the
simulation results are in good agreement with the analytical results in Eqs.
11 and 12. These results differ substantially from the corresponding
distributions for unbiased random walks, in which case every edge on the
network has the identical probability distribution for load and flux
independent of network topology [17]. In the present case of biased random
walks, the distribution of flux or load depends on nodal properties. Futher,
as the load distributions reveal, for $\alpha>0$ since walkers preferentially
move towards hubs, edges with smaller mean loads do not attract walkers. On
the other hand, if strongly biased towards small degree nodes, it is possible
for walkers to avoid edges which connect fairly high degree nodes. Hence, the
higher probability of null mean loads as seen in Fig. 2 for $\alpha<0$. Thus,
by tuning the bias, an edge can be made to either attract or repel walkers.
The results for $\alpha=0$ (unbiased random walk) agree exactly with the
results of Ref. [17]. However, unlike the case of unbiased random walk where
every edge has the same probability distribution, edges with different mean
loads have distributions according to Eq. 9.

Figure 2: The probability distributions of (a) load and (b) flux for a
randomly chosen edge $e_{0,11}$ connecting nodes 0 and 11 with degrees 100 and
44, respectively. The different curves indicate the distributions obtained by
executing DBRW for varying values of bias strength
$\alpha\in\\{-2,-1,0,1,2\\}$. Each execution had 2500 walkers running for
100,000 time steps on the network in Fig. 1 and was repeated 15 times with
varying initial positions of random walkers.

Next, the probability for the occurrence of load EE, $Q_{L}\left(\langle
l_{ij}\rangle\right)$ on edge $e_{ij}$ is considered. The theoretical result
is obtained from Eq. 14 and is juxtaposed with simulation results for
$Q_{L}\left(\langle l_{ij}\rangle\right)$ in Fig. 3. It is displayed for all
the edges for several threshold values indexed by $m$ at a fixed value of bias
($\alpha=1$). An excellent agreement is observed between the analytical
results and the simulations. In an average sense, the EE probability is higher
for edges with lower values of $\langle l_{ij}\rangle$ in comparison to edges
with larger $\langle l_{ij}\rangle$. This effect is more pronounced for higher
EE thresholds. As $m\to 0$, the mean load itself becomes the threshold value
effectively leading to most events designated as extreme. In this case, EE
probability does not vary much with $\langle l_{ij}\rangle$. Clearly, this
quantity varies with $\langle l_{ij}\rangle$ but does not vary much as $m$ is
varied. Hence, for a fixed value of threshold indexed by $m$, EE probability
has been averaged over different number of extreme events. It was also found
that the density of edges for each value of $\langle l_{ij}\rangle$ differed,
but for a given value of $\langle l_{ij}\rangle$, the relative density
remained the same for all values of $m$. This is not a peculiarity of biased
random walks but is reflective of the inhomogeneity of the network structure.

Figure 3: Probability of load EE vs mean load on edge for different thresholds
of extremeness: $m\epsilon[1,5]$. 2500 walkers performed a DBRW with
$\alpha=1$ for 100,000 time steps on the network mentioned in Fig.1

Now that we have an understanding of the effect of the choice of threshold we
proceed to analyse how changing the bias strength affects the system. Fig. 4
shows $Q_{L}\left(\langle l_{ij}\rangle\right)$ plotted against the mean load
$\langle l_{ij}\rangle$ for varying values of $\alpha$ at fixed EE threshold
of $m=2$. We have run it on a scale-free network with 2000 nodes and 7984
edges. On this network, 15968 non-interacting walkers execute degree biased
random walk for 75000 time steps (note that these are the parameters and
network used for the rest of the graphs, unless mentioned otherwise). An
excellent agreement between analytical and numerical simulations can be seen.
Note that for the unbiased random walk dynamics ($\alpha=0$) shown in Fig.
4(c) all the edges have the same mean load and hence identical EE occurrence
probability, in agreement with the earlier result [17]. The effects of biasing
is visible in the figure. When biased towards hubs with $\alpha>0$, more
walkers moving towards the hubs leads to EE on edges which connect two hubs.
At the other end of the spectrum, when $\alpha<0$, walkers preferentially move
towards small degree nodes and consequently fewer or no EE takes place on
edges connecting hubs. The discontinuity in $Q_{L}\left(\langle
l_{ij}\rangle\right)$ arises due to the discreteness of walkers while the
threshold $q_{ij}$ in Eq. 14 is a real number.

In general, considerable variability in the EE probability can be observed. As
Fig. 4(e) shows for $\alpha=-2$, for mean load $\langle l\rangle<1$, EE
probability varies by 4 orders in magnitude. For mean loads of $\langle
l\rangle>1$ (yellow background), the variability is about one order of
magnitude. Similar wide range of variabilities are observed for other values
of biases shown in Fig. 4. The exception is the case of $\alpha=0$
corresponding to unbiased random walk, for which the mean load is same on all
the edges of the network and consequently the EE probability is identical on
all the edges [17]. The variability in EE probability becomes more pronounced
for $\alpha<0$. This is a new feature not encountered in the EE probabilities
on nodes or edges of networks on which unbiased random walks are executed [9,
12]. It would imply the following: on the edges characterised by large mean
loads, the extreme events are approximately equiprobable irrespective of its
local network structure. On edges with small mean loads, this probability
varies by large orders of magnitude. Hence, two edges whose mean loads do not
appreciably differ from one another can have very different EE probabilities.
So as far as EE are concerned, biased walker dynamics on edges leads to the
network displaying an approximately logical (not necessarily physical) phase
separation in to two parts – one in which EE have only small variations about
an average value, and the other in which large excursions about the mean are
observed. This will have implications for strategies that try to harness
extreme events. Finally, it must be emphasised that this variability is not
due to insufficient averaging but an inherent feature of biased dynamics in
networked system.

###  V.2 Other Biases

Figure 4: Log-log plot of the probabilities of edge EE for load vs mean edge
load for different values of bias parameters (a) $\alpha=2$, (b) $\alpha=1$,
(c) $\alpha=0$, (d) $\alpha=-1$, (e) $\alpha=-2$ and (f) $\alpha=-4$. The red
squares represent analytical results and the green circles indicate numerical
simulations. All runs have $m=2$, and were carried out for 75,000 time steps
each. The network used was a scale-free network with 2000 nodes and 7984 edges
with 15968 non-interacting biased random walkers.

In this subsection, we switch to a scale-free network [30] and apply different
kinds of biased walks based on network centrality measures. A scale-free
network with 2000 nodes and 7984 edges is created. On this network, 15968 non-
interacting walkers executed random walk dynamics for $75000$ time steps. The
walkers are randomly distributed on the nodes of the network initially.
Simulation data is recorded after discarding the first $10^{3}$ time steps on
account of transient behaviour.

Since centrality measures provide different ways to rank the importance of
nodes in a network, they can also be used as biasing functions. We chose four
measures - Betweenness, Eigenvector, Closeness and PageRank centralities [31,
32, 33, 34] and repeated the procedure followed in section V. The Betweenness
centrality $B_{v}$ of a node labelled $v$ is given by

$B_{v}=\sum_{a,b}\frac{\sigma_{a,b}(v)}{\sigma_{a,b}},$ (18)

where $\sigma_{a,b}$ is the total number of shortest paths from node $a$ to
$b$, while $\sigma_{a,b}(v)$ is the number of such paths passing through node
$v$. The summation is over all pairs of nodes $(a,b)$ on the network such that
$a\neq b\neq v$. It represents the degree to which nodes stand between each
other. The eigenvector centrality defines the importance of a node in terms of
the neighbourhood of each node. Let $\mathbf{A}$ be the adjacency matrix, then
the Eigenvector centrality is given by the left eigenvector $\mathbf{x}$
corresponding to the largest eigenvalue $\lambda$:

$\lambda\mathbf{x}=\mathbf{x}\,\mathbf{A}.$ (19)

The $v$-th element of $\mathbf{x}$ (denoted by $X_{v}$) gives the centrality
measure for the node $v$ on a network. The closeness centrality is the average
length of the shortest path between the node and all other nodes in the graph.
If $d(v,u)$ represents the length of shortest path between nodes $v$ and $u$,
then closeness centrality is defined as

$C_{v}=\frac{N}{\sum_{u}d(v,u)}.$ (20)

Thus, if $C(v)$ is large, then it implies that most $d(v,u)$ are small, and
hence node $v$ is ”closer” to other nodes on the network. For the present
purposes, the biasing function for the $v$-th node is taken to be

$f_{v}\propto\Phi_{v}^{\alpha},$ (21)

where $\Phi_{v}=B_{v},\,X_{v}$, $C_{v}$ or PageRank centrality [34]. The
transition probability and the occupation probability can be obtained by using
$f_{v}$ in Eqs. 1 and 3 respectively. From these, the load and flux
distributions as well as their extreme event probabilities can be computed as
was done earlier for degree biased walks.

In Fig. 5(b-e), the load EE probabilities on edge $e_{ij}$ is plotted for the
four different biasing functions. For comparison, the degree bias
($k^{\alpha}$) is shown in Fig. 5(a). In all the cases, the same strength
parameter $\alpha=1$ and the EE threshold is set to $q_{ij}=\langle
l_{ij}\rangle+2\sigma_{ij}$. Even though each of these centrality measures
capture a different facet of the notion of importance for a node, the EE
probabilities on edges of the network do not differ significantly irrespective
of which biasing function was applied. This can be attributed to correlations
found between the centrality measures [35, 36].

Figure 5: Log-log plot of probabilities of edge EE for load against mean edge
load using various biasing functions from (a) Degree, (b) Betweenness, (c)
Closeness, (d) PageRank, and (e) Eigenvector centralities. The red squares
represent analytical results and the green circles indicate numerical
simulations. All runs have $m=2$, $\alpha$=1 and were carried out for 75,000
time steps each on a scale-free network with 2000 nodes and 7984 edges with
15968 non-interacting biased random walkers.

As seen in the case of degree biased random walks, in Fig. 5 too, a large
variability in EE probability is observed for edges whose mean load is
approximately in the range $\langle l\rangle<1$, as compared to the edges with
$\langle l\rangle>1$ (yellow background). Thus, irrespective of the actual
biased walk algorithm employed, pronounced variability in EE probability for
some of the edges appears to be a robust feature. For a given biased dynamics
on a network, let $\Delta l=\langle l\rangle_{\text{m}ax}-\langle
l\rangle_{\text{m}in}$ represent difference between the maximum and minimum
value of mean loads on its edges. It can be observed that $\Delta l$ is a
function of applied bias. As seen in Figs. 4 and 5, $\Delta l$ is largest for
betweenness centrality based walk ($=8.542\times 10\textsuperscript{2}$), and
smallest for closeness centrality based walk ($=2.726\times
10\textsuperscript{0}$). For unbiased random walks, $\Delta l=0$. This doesn’t
preclude the possibility that the same edge under different biases can have
different EE probabilities. Further, since the EE thresholds are defined based
on these measures, it is evident that what might be extreme in one setting
would theoretically not be so in another. This implies that the choice of bias
and its strength are crucial factors in determining extreme event occurrence
rates.

##  VI Correlations between Extreme Events

As reported in [17], the unbiased random walk process on networks is
uncorrelated, the extreme events on nodes are almost temporally uncorrelated
with that on the edges. Yet, the extreme events arising from the biased random
walks can display significant time-lagged correlations, i.e., extreme events
occurring on nodes can be correlated with the extreme events on the
neighbouring nodes or edges after a time lag $d$. Let us denote by
$\mathcal{N}(i)$ the set of nodes that are connected to node $i$. In this
case, two possible scenarios exist, (i) node-node correlation: EE occurs at
time $t$ on a node $i$, and occurs at time $t+d$ on node $j$ such that
$j\in\mathcal{N}(i)$, (ii) edge-node correlation: EE occurs at time $t$ on an
edge $e_{ij}$, and occurs at time $t+d$ on $i$-th or $j$-th node connected to
this edge. To quantify these correlations, we coarse-grain the time series of
walkers as follows. Let $w(t),\,t=1,2,3,\dots,75000$ represent the time
varying walker flux through an edge or node. This series is modified as
follows. If an EE occurs at time $t=\tau$, then $w(\tau)=1$, else $w(\tau)=0$.
This coarse-grained time series is used to calculate the (normalized) time
lagged cross correlation (TLCC) function [37] between any two time series
$x(t)$ and $y(t)$. It is defined as

$r_{d}=\frac{\sum_{t=1}^{T-d}(x(t)-\langle x\rangle)~{}(y(t+d)-\langle
y\rangle)}{\sqrt{\sum_{t=1}^{T-d}(x(t)-\langle
x\rangle)^{2}}~{}\sqrt{\sum_{t=1}^{T-d}(y(t+d)-\langle y\rangle)^{2}}}.$ (22)

Here $r_{d}\in[-1,1]$ and indicates if the two time series are correlated or
anti-correlated at time lag $d$. We calculate node-node and edge-node
correlation $r_{d}$ for the EE time series $w(t)$ for all pairs of neighbors
in the network with DBRW dynamics for lags $d\in[-10,10]$. The results are
shown in Fig. 6 are for five node pairs. The five node pairs are
representative of the different behaviour found, they are labelled with
degrees of the nodes the edges connect. Each colour represents the same node
pair across all subplots. As is evident in Fig. 6(c), the node-node
correlations for unbiased walkers ($\alpha=0$) are barely significant for any
time lag $d$. On the other hand, as the bias strength increases, i.e,
$|\alpha|>0$, strong correlations emerge among connected pairs of nodes.

Figure 6: Node-Node correlations of EE time series for five representative
node pairs labelled by their degrees [$k_{i},k_{j}$] in the legend. Subplots
compare this over different strengths of the degree bias (a) $\alpha=2$, (b)
$\alpha=1$, (c) $\alpha=1$, (d) $\alpha=-1$, (e) $\alpha=-2$ and (f)
$\alpha=-4$ against lag $d$.

From Fig. 6, we see that the node pairs display strong correlations for lags
$d=-1$ and $d=1$. This means that if an EE occurred in a node its neighboring
node is likely to have experienced one in the previous time step and/or
experience one in the next step. The likelihood of this happening increases as
$\alpha$ gets more negative. Furthermore, there exist node pairs (such as the
one labelled [$k_{i},k_{j}$]=[5,4] in Fig. 6(d-f)) whose correlations do not
decay rapidly with lag $d$. This means that EE among neighbors can remain
correlated for long times. This might be useful for predicting such events. As
the bias gets more negative the fraction of pairs with a nontrivial amount of
correlation ($\>0.4$ by convention) increases. This fraction depends on degree
distribution of the network and for this scale-free network it is $<10\%$
until $\alpha=-2$. The strength of this lagged correlation depends on the
product of the degrees of the nodes the edge $e_{ij}$ connects, $k_{i}k_{j}$.
If $k_{i}k_{j}$ is large, the lagged correlation will be small since the EE
tend to get ”scattered” away in to the many neighbours at the next time step.
If $k_{i}k_{j}$ is small, the limited local neighbourhood of poorly connected
nodes could result in a fraction of the negatively biased walkers to be
dispersed again to the same node.

Next we investigate the lagged correlation between the EE for load on an edge
$e_{ij}$ and the EE occurring on one of the nodes it connects to, say, node
$i$. In Fig. 7, lagged correlation $r_{d}$ is shown for DBRW dynamics for
varying bias strength $\alpha$. In general, edge-node correlations are far
weaker than the node-node correlations. For a fixed $\alpha$, the strongest
cross correlation is seen at time lags of $d=-1$ and $d=1$. It implies that if
there is a load EE on an edge then one of the nodes connecting it is likely to
encounter an EE either one time step before and/or one time step later. This
effect gets more pronounced as $\alpha$ gets more negative. In this case,
walkers preferentially explore small degree nodes which have fewer edges
connected to them. Hence, once an extreme event takes place in one such node,
the walkers have no option but to transit through one of its edges creating an
extreme event on the edge as well. Due to this, walkers continually circulate
locally among the few small degree nodes leading to extremely slow decay of
correlations as seen in Fig. 7. As $\alpha$ gets more negative, larger
fraction of edges display such slow decay of correlations with lag $d$.

On the other hand, for $\alpha>0$, the walkers preferentially explore hubs on
the network. If an extreme event takes place in one of the hubs, the walkers
are scaterred away through a large number of edges connected to it. Hence, it
might not exceed the threshold required to designate it as an extreme event.
Hence lagged correlations are not pronounced in the case of $\alpha>0$ as
observed in Fig. 7(a-b).

Figure 7: Edge-Node correlations between EE time series for five
representative edges that were taken and labelled by the degrees of the nodes
[$k_{i},k_{j}$] they connected. Subplots compare this over different strengths
of the degree bias with bias strength (a) $\alpha=2$, (b) $\alpha=1$, (c)
$\alpha=1$, (d) $\alpha=-1$, (e) $\alpha=-2$ and (f) $\alpha=-4$ against lag
$d$. The peak values are observed for lag = -1 and 0.

##  VII Network stability under edge extreme events

Using the extremes related results obtained above, we are now in a position to
discuss the robustness of a network in the context of extreme events faced by
edges. To study this, we adopt two broad classes of node or edge deletion
strategies: (i) targeted node/edge deletion based on decreasing order of EE
occurrence probability, (ii) random edge/node deletion. To illustrate the
effects of these deletion strategies, we monitor the size of the relative
giant component $S(t)=\frac{N_{\text{max}}(t)}{N}$ in the network as a
function of time. It is the ratio of nodes in the largest connected component
$N_{\text{max}}$ to the total number of nodes $N$ initially in a network. If
all the nodes are reachable from any other node,
$S=\frac{N_{\text{max}}}{N}=1$. As we delete nodes/edges $S(t)$ decreases with
time.

In Fig. 8(a), the simulation results for a variety of biased random walks are
considered and for each one the targeted and random node deletion strategies
are shown. In most cases, even if 30-40% of the nodes are removed in a
targeted manner, the giant component in the network falls to less than 50%. In
the case of degree biased random walks, even as 50% of nodes are targeted and
removed, giant component size becomes vanishingly small. However, if the nodes
are randomly removed, a reasonably sized giant component in the network
survives even until about 80% of the nodes are deleted. It is well-known that
scale-free network is resilient to random attacks but not to targeted attacks
[38]. The difference here is that the same viewpoint is reinforced even if the
node deletion strategy is based on occurrence of extreme events.

A somewhat surprising behaviour is seen in Fig. 8(b), in which edges deletion
results are shown for both targeted and random removal strategies. Unlike the
node deletion case, a significant fraction of the giant component survives
well until about 80% of the edges are removed. Indeed, until 30% of edges are
deleted, the size of giant component does not decrease appreciably. Thus, as
far as the survival of network structure is concerned, node and edge deletion
processes are not equivalent. If $q$% of nodes are deleted, the corresponding
giant component size is $S^{n}_{q}$. Similarly, if $q$% of edges are removed,
let the giant component size be denoted by $S^{e}_{q}$. Generally, we observe
that $S^{e}_{q}>S^{n}_{q}$ for nearly all $q$. Physically, this effect arises
because if a node is deleted, typically many edges becomes dysfunctional for
transport of randowm walkers and this tends to quickly reduce the size of
giant component. On the other hand, if an edge is deleted, most other parts of
the network continues to remain functional for transport. Thus, deletion of a
node and an edge lead to different outcomes for network resilience. It might
also be pointed out that, unlike node deletion strategies, there is virtually
no difference between random and targeted edge deletion for the same reason
that scale-free networks are more resilient to edge deletions than node
deletions.

Figure 8: Size of the relative giant component tracked as (a) nodes are
removed in decreasing order of their node EE probability ${P}_{EE}(q_{i})$ (b)
edges are removed in decreasing order of $Q_{L}\left(\langle
l_{ij}\rangle\right)$. In both cases an EE threshold of $m=2$ was used. In
(b), curves for $\alpha=-1$ and $-2$ lie on top of each other as do the curves
for $\alpha=1$ and $2$.

##  VIII Conclusions

In summary, we have proposed a general formalism to study extreme events
taking place on the edges, rather than on the nodes of a complex network. In
this work, we considered the dynamics of biased random walkers on the network,
and extreme events on the edges based on the deviation of flux of walkers
above a defined threshold. Unlike the case of nodes, in the edge $e_{ij}$ of
the network, walkers can go from node $i$ to node $j$ and vice versa. Hence,
it was necessary to distinguish between two distinct dynamical quantities –
the load, which is the total number of walkers traversing the edge. Secondly,
the flux denoting the difference between the number of walkers going $i$ to
$j$ or in the reverse direction. Thus, enabling the consideration of extreme
events using both of these dynamical observables. In this work, we have
applied the proposed formalism to study the occurrence of extreme events in
load and flux on the edges. It must be pointed out that extreme events in this
work represent one possible outcome due to intrinsic fluctuations in the
number of biased random walkers, and is not triggered by external sources.

We obtained analytical and numerical estimates for the occurrence
probabilities of extreme events on the edges of a network. The occurrence
probability of extreme events on edges depends on the local network structure
around the edge of interest. Within the framework of biased random walkers, it
is shown that the edges can be characterised by the mean loads (or mean
fluxes). Then, $e_{ij}=e_{ij}(\langle l\rangle)$, where $\langle l\rangle$
denotes mean load. For the edges with low mean loads ($\langle l\rangle<1$),
it is shown that the variability of the extreme event probability is higher
compared to edges with higher mean loads ($\langle l\rangle>1$). Strongly
biased random walks tend to increase the variability for $\langle l\rangle<1$,
and suppress the variability at $\langle l\rangle>1$. We show these results
for varying strengths of the degree biased random walk on scale-free networks.
We also show similar trends are realised for other biases – Betweenness,
Closeness, Eigenvector, and PageRank centrality based biases. In addition,
these are also demonstrated for a real-life planar network.

Further, we studied how the extreme events on nodes and edges are related.
Earlier, it was shown that these two classes of extreme events are
uncorrelated for unbiased random walkers on networks. However, while using
biased random walkers a non-trivial correlation is found to exist. This
implies that an extreme event on a node tends to influence the occurrence or
non-occurrence of an extreme event on an edge connected to it. We quantify
such correlations for various strengths of the degree bias.

Finally, we also studied the vulnerability of the network as a whole to the
occurrence of extreme events on the edges. We removing nodes/edges in
decreasing order of the probability of their extreme event occurrence. For
instance, the edge with the highest probability will be removed first, and
then the node with next highest probability will be removed and so on. After
each removal, the size of the largest connected component of the network is
computed. It is shown that the largest connected component stays intact almost
until 50% of the edges are removed due to extreme events. This is in contrast
to the case of node removal based on similar criteria outlined above. In which
case, the size of the giant component falls almost by 50% when half of the
nodes are removed. Thus, networks can be expected to maintain their resilience
much more in the face of extreme events on the edges, than if it were to
happen on the nodes.

It would be interesting to study extreme events on real-life network with
measured loads and fluxes and compare them to patterns observed in synthetic
networks. An absorbing exploration would be to map community structure [24]
but when supplied with the EE data. Keeping track of how EE cascades through
the network using higher order versions of the correlations observed among the
nodes and edges is another avenue of research. Using measures that quantify
the local neighbourhood of nodes involved in node-node and edge-node
correlations instead of $k_{i}kj$ can provide better insight about strong
correlations observed. Developing a more accommodating definition of EE can
lead to modeling more than congestion-type systems. Using time evolving
properties of nodes, edges (and if possible even the walkers) can lead to
interesting generalisations of the theory. Trying out other biasing functions,
especially non-structural ones whose values are bestowed externally could lead
to interesting trends. One possible extension to this framework is to get
equations analytically for time-dependent biasing functions as well. In
general, dynamics on the edges have long been neglected and these results, we
hope, might lead to more work on extreme events on the edges of networks.

###### Acknowledgements.

GG thanks Aanjaneya Kumar for many discussions during this work. GG
acknowledges support from IISER Pune and INSPIRE Scholarship for Higher
Education (SHE). MSS acknowledges the support from MATRICS grant of SERB, Govt
of India.

##  Appendix A Distribution of load over an edge

We have $W$ total walkers. The total number of walkers on node $i$ and node
$j$ is $n$. The number of walkers that traverse edge $e_{ij}$ is $l_{ij}$
(load). There are three kinds of events that can happen in a network from the
point of view of an edge $e_{ij}$:

  1. 1.

The walker could be neither on node $i$ nor on node $j$. There are $W-n$ such
walkers with probability $=1-p_{i}-p_{j}\equiv E_{1}$.

  2. 2.

The walker could be on node $i$ and not traversing to node $j$; or the walker
could be on node $j$ and does not travel to node $i$. There are $n-l_{ij}$
such walkers with probability $=p_{i}(1-\pi_{ij})+p_{j}(1-\pi_{ji})\equiv
E_{2}$.

  3. 3.

The walker could be on node $i$ or node $j$ and traverses to node $j$ or node
$i$ respectively. There are $l_{ij}$ such walkers with probability
$=p_{i}\pi_{ij}+p_{j}\pi_{ji}=2p_{i}\pi_{ij}\equiv E_{3}.$

The detailed balance condition holds for ergodic markov chains i.e.,
$p_{i}\pi_{ij}=p_{j}\pi_{ji}$. The probability that there are $W-n$ events of
case 1, $n-l_{ij}$ events of case 2, and $l_{ij}$ events of case 3 happening
for an edge $e_{ij}$ is

$\frac{W!}{(W-n)!(n-l_{ij})!(l_{ij})!}[E_{1}]^{W-n}[E_{2}]^{n-l_{ij}}[E_{3}]^{l_{ij}},$
(23)

where the $\frac{W!}{(W-n)!(n-l)!(l)!}$ factor accounts for repeated events
since the order in which these three independent events occur does not matter.
Now, this is just the joint distribution of three events occurring. Since we
are not interested in the distribution of walkers among cases 1 and 2. The
probability that an edge has load $=l_{ij}$ is simply the marginal probability
of Eq. 23 for case 3 to occur i.e., the probability that there are $l_{ij}$
walkers that hop from node $i$ to node $j$ or vice versa. This is given by

$\displaystyle
P_{L}(l_{ij})=\sum_{n=l_{ij}}^{W}\frac{W!}{(W-n)!(n-l_{ij})!(l_{ij})!}$
$\displaystyle\times[E_{1}]^{W-n}[E_{2}]^{n-l_{ij}}[E_{3}]^{l_{ij}}$ (24)

After some simple manipulations, we find ourselves with the load distribution
from Eq. 9

$\displaystyle
P_{L}(l_{ij})=\binom{W}{l_{ij}}(E_{3})^{l_{ij}}(1-E_{3})^{W-l_{ij}},$
$\displaystyle=\binom{W}{l_{ij}}(2p_{i}\pi_{ij})^{l_{ij}}(1-2p_{i}\pi_{ij})^{W-l_{ij}}.$
(25)

## References

  * Albeverio  _et al._ [2006] S. Albeverio, V. Jentsch, and H. Kantz,  _Extreme events in nature and society_ (Springer Science & Business Media, 2006).
  * Vespignani [2012] A. Vespignani, Modelling dynamical processes in complex socio-technical systems, Nature physics 8, 32 (2012).
  * Barabási [2012] A.-L. Barabási, The network takeover, Nature Physics 8, 14 (2012).
  * Gumbel [1958] E. J. Gumbel,  _Statistics of extremes_ (Columbia university press, 1958).
  * Chen  _et al._ [2013] Z. Chen, W. Hendrix, H. Guan, I. K. Tetteh, A. Choudhary, F. Semazzi, and N. F. Samatova, Discovery of extreme events-related communities in contrasting groups of physical system networks, Data Mining and Knowledge Discovery 27, 225 (2013).
  * Chen  _et al._ [2014] Y.-Z. Chen, Z.-G. Huang, and Y.-C. Lai, Controlling extreme events on complex networks, Scientific reports 4, 1 (2014).
  * Mahecha  _et al._ [2017] M. D. Mahecha, F. Gans, S. Sippel, J. F. Donges, T. Kaminski, S. Metzger, M. Migliavacca, D. Papale, A. Rammig, and J. Zscheischler, Detecting impacts of extreme events with ecological in situ monitoring networks, Biogeosciences 14, 4255 (2017).
  * Albert and Barabási [2002] R. Albert and A.-L. Barabási, Statistical mechanics of complex networks, Reviews of modern physics 74, 47 (2002).
  * Kishore  _et al._ [2011] V. Kishore, M. Santhanam, and R. Amritkar, Extreme events on complex networks, Physical review letters 106, 188701 (2011).
  * Moitra and Sinha [2019] P. Moitra and S. Sinha, Emergence of extreme events in networks of parametrically coupled chaotic populations, Chaos: An Interdisciplinary Journal of Nonlinear Science 29, 023131 (2019), https://doi.org/10.1063/1.5063926 .
  * Amritkar  _et al._ [2018] R. E. Amritkar, W.-J. Ma, and C.-K. Hu, Dependence of extreme events on spatial location, Phys. Rev. E 97, 062102 (2018).
  * Kishore  _et al._ [2012] V. Kishore, M. S. Santhanam, and R. E. Amritkar, Extreme events and event size fluctuations in biased random walks on networks, Physical Review E 85, 056120 (2012).
  * Chen  _et al._ [2015] Y.-Z. Chen, Z.-G. Huang, H.-F. Zhang, D. Eisenberg, T. P. Seager, and Y.-C. Lai, Extreme events in multilayer, interdependent complex networks and control, Scientific reports 5, 1 (2015).
  * Gupta and Santhanam [2021] K. Gupta and M. S. Santhanam, Extreme events in nagel–schreckenberg model of traffic flow on complex networks, The European Physical Journal Special Topics , 1 (2021).
  * Mizutaka and Yakubo [2013] S. Mizutaka and K. Yakubo, Structural robustness of scale-free networks against overload failures, Physical Review E 88, 012803 (2013).
  * Mizutaka and Yakubo [2014] S. Mizutaka and K. Yakubo, Network fragility to overload failures: influence of the scale-free property, journal of complex networks 2, 413 (2014).
  * Kumar  _et al._ [2020] A. Kumar, S. Kulkarni, and M. Santhanam, Extreme events in stochastic transport on networks, Chaos: An Interdisciplinary Journal of Nonlinear Science 30, 043111 (2020).
  * Battiston  _et al._ [2016] F. Battiston, V. Nicosia, and V. Latora, Efficient exploration of multiplex networks, New Journal of Physics 18, 043035 (2016).
  * Ninio [1976] F. Ninio, A simple proof of the perron-frobenius theorem for positive symmetric matrices, Journal of Physics A: Mathematical and General 9, 1281 (1976).
  * Noh and Rieger [2004] J. D. Noh and H. Rieger, Random walks on complex networks, Physical review letters 92, 118701 (2004).
  * Barabási  _et al._ [2004] A. Barabási _et al._ , Fluctuations in network dynamics., Physical review letters 92, 028701 (2004).
  * Fronczak and Fronczak [2009] A. Fronczak and P. Fronczak, Biased random walks in complex networks: The role of local navigation rules, Physical Review E 80, 016107 (2009).
  * Wang  _et al._ [2006] W.-X. Wang, B.-H. Wang, C.-Y. Yin, Y.-B. Xie, and T. Zhou, Traffic dynamics based on local routing protocol on a scale-free network, Physical review E 73, 026111 (2006).
  * Zlatić  _et al._ [2010] V. Zlatić, A. Gabrielli, and G. Caldarelli, Topologically biased random walk and community finding in networks, Physical Review E 82, 066109 (2010).
  * Goldhirsch and Gefen [1987] I. Goldhirsch and Y. Gefen, Biased random walk on networks, Physical Review A 35, 1317 (1987).
  * Lee  _et al._ [2014] Z. Q. Lee, W.-J. Hsu, and M. Lin, Estimating mean first passage time of biased random walks with short relaxation time on complex networks, Plos one 9, e93348 (2014).
  * Bonaventura  _et al._ [2014] M. Bonaventura, V. Nicosia, and V. Latora, Characteristic times of biased random walks on complex networks, Physical Review E 89, 012803 (2014).
  * Boeing [2017] G. Boeing, Osmnx: New methods for acquiring, constructing, analyzing, and visualizing complex street networks, Computers, Environment and Urban Systems 65, 126 (2017).
  * Kuratowski [1930] C. Kuratowski, Sur le probleme des courbes gauches en topologie, Fundamenta mathematicae 15, 271 (1930).
  * Barabási and Albert [1999] A.-L. Barabási and R. Albert, Emergence of scaling in random networks, science 286, 509 (1999).
  * Brandes [2001] U. Brandes, A faster algorithm for betweenness centrality, Journal of mathematical sociology 25, 163 (2001).
  * Bonacich [1987] P. Bonacich, Power and centrality: A family of measures, American journal of sociology 92, 1170 (1987).
  * Freeman [1978] L. C. Freeman, Centrality in social networks conceptual clarification, Social networks 1, 215 (1978).
  * Page  _et al._ [1999] L. Page, S. Brin, R. Motwani, and T. Winograd,  _The PageRank citation ranking: Bringing order to the web._ , Tech. Rep. (Stanford InfoLab, 1999).
  * Oldham  _et al._ [2019] S. Oldham, B. Fulcher, L. Parkes, A. Arnatkevic̆iūtė, C. Suo, and A. Fornito, Consistency and differences between centrality measures across distinct classes of networks, PloS one 14, e0220061 (2019).
  * Valente  _et al._ [2008] T. W. Valente, K. Coronges, C. Lakon, and E. Costenbader, How correlated are network centrality measures?, Connections (Toronto, Ont.) 28, 16 (2008).
  * Gubner [2006] J. A. Gubner,  _Probability and random processes for electrical and computer engineers_ (Cambridge University Press, 2006).
  * Cohen  _et al._ [2001] R. Cohen, K. Erez, D. Ben-Avraham, and S. Havlin, Breakdown of the internet under intentional attack, Physical review letters 86, 3682 (2001).


#  Asymptotic dynamics of high dynamic range stratified turbulence

G. D. Portwood1,2 S. M. de Bruyn Kops1 C. P. Caulfield 3,4 1 Department of
Mechanical and Industrial Engineering, University of Massachusetts, Amherst,
USA 01003  2 Methods and Algorithms (XCP-4), Computational Physics Division,
Los Alamos National Laboratory, Los Alamos, USA 87545  3 BP Institute,
University of Cambridge, Cambridge, CB3 0EZ, UK  4 Department of Applied
Mathematics & Theoretical Physics, University of Cambridge, Cambridge, CB3
0WA, UK

###### Abstract

Direct numerical simulations of homogeneous sheared and stably stratified
turbulence are considered to probe the asymptotic high-dynamic range regime
suggested by Gargett  _et al._ (1984) and Shih  _et al._ (2005). We consider
statistically stationary configurations of the flow that span three decades in
dynamic range defined by the separation between the Ozmidov length scale,
$L_{O}=\sqrt{\epsilon/N^{3}}$, and the Kolmogorov length scale,
$L_{K}=(\nu^{3}/\epsilon)^{1/4}$, up to
$\mathrm{Re_{b}}\equiv(L_{O}/L_{K})^{4/3}=\epsilon/(\nu N^{2})\sim O(1000)$,
where $\epsilon$ is the mean turbulent kinetic energy dissipation rate, $\nu$
is the kinematic viscosity, and $N$ is the buoyancy frequency. We isolate the
effects of $\mathrm{Re_{b}}$, particularly on irreversible mixing, from the
effects of other flow parameters of stratified and sheared turbulence (Mater
and Venayagamoorthy, 2014). Specifically, we evaluate the influence of dynamic
range independent of initial conditions. We present evidence that the flow
approaches an asymptotic state for $\mathrm{Re_{b}}\gtrapprox 300$,
characterized both by an asymptotic partitioning between the potential and
kinetic energies and by the approach of components of the dissipation rate to
their expected values under the assumption of isotropy. As $\mathrm{Re_{b}}$
increases above 100, there is a slight decrease in the turbulent flux
coefficient $\Gamma=\chi/\epsilon$, where $\chi$ is the dissipation rate of
buoyancy variance, but, for this flow, there is no evidence of the commonly
suggested $\Gamma\propto\mathrm{Re_{b}}^{-1/2}$ dependence when
$100\leq\mathrm{Re_{b}}\leq 1000$.

##  I Introduction

Sheared, stratified turbulence, energized by vertical shearing of horizontal
motions in the presence of a statically stable density distribution, arises
throughout the world’s oceans and atmosphere. Dynamical models of such flows,
in particular capturing the vertical transport of heat due to irreversible
mixing, are essential for modeling the global climate system because mixing
occurs on relatively small scales, several orders of magnitude below those
currently resolved in basin-scale models (Ivey  _et al._ , 2008; Ferrari and
Wunsch, 2009; Gregg  _et al._ , 2018).

A key challenge for theoreticians is to develop a robust parameterization
usable in such models for the vertical eddy diffusivity of heat $\kappa_{T}$
defined as

$\kappa_{T}\equiv\frac{B}{N^{2}}\equiv\frac{\langle\frac{g}{\rho_{0}}w\rho\rangle}{N^{2}},\
N\equiv\sqrt{-\frac{g}{\rho_{0}}\frac{d\bar{\rho}}{dz}},$ (1)

where angled brackets denote ensemble averaging, $B$ is the buoyancy flux, $N$
is the buoyancy frequency, $g$ is the constant gravitational acceleration,
$\rho_{0}$ is a constant reference density, $(u,v,w)$ and $\rho$ are the
fluctuating velocity vector and fluctuating density, respectively, in the
$(x,y,z)$ coordinate system, and $\bar{\rho}$ is the mean ambient density with
linear functional dependence in $z$. Here the Boussinesq approximation has
been made such that density variations are sufficiently small so that a linear
equation of state is appropriate and density variations are only significant
in the buoyancy force.

Developing dynamical models of these complex flows has proved exceptionally
difficult and controversial due not least to their potential dependence on a
wide range of parameters (Mater and Venayagamoorthy, 2014; Ivey  _et al._ ,
2018). For simplicity, we fix $\mathrm{Pr}=\nu/\kappa=1$, where $\nu$ is the
kinematic viscosity and $\kappa$ is the molecular thermal diffusivity. Even
with this assumption, stratified sheared turbulence is influenced by at least
four independent length scales, associated with characteristic values of
overall stratification, shear and both the intensity and decay rate of
turbulence (Mater and Venayagamoorthy, 2014). The Ozmidov scale $L_{O}$, the
Corrsin (or shear) scale $L_{C}$, the ‘large-eddy’ turbulent length scale
$L_{L}$, and the Kolmogorov microscale $L_{K}$ are defined here as

$\displaystyle\begin{aligned}
L_{O}&\equiv\left({\frac{\epsilon}{N^{3}}}\right)^{1/2};\
&L_{C}&\equiv\left({\frac{\epsilon}{S^{3}}}\right)^{1/2};\\\
L_{L}&\equiv\frac{E_{k}^{3/2}}{\epsilon};\
&L_{K}&\equiv\left(\frac{\nu^{3}}{\epsilon}\right)^{1/4},\end{aligned}$ (2)

where $E_{k}$ is the averaged turbulent kinetic energy, $\epsilon$ is its
dissipation rate and $S$ is the mean constant vertical shear, i.e. $S\equiv
d\bar{u}/dz$ where $\bar{u}$ is the mean streamwise velocity with linear
functional dependence in $z$. These four scales are determined by the
properties of the fluctuation velocity field (i.e. the turbulence) relative to
the background buoyancy frequency and shear. As noted by Ivey  _et al._
(2018), it is entirely plausible that length scales comparing the perturbation
scalar field to these background quantities may be more physically relevant,
although particularly in such stationary flows such as those considered here,
it is entirely plausible that the various length scales may become closely
coupled (as we investigate further below).

Although in general the scales defined in (2) vary in space and time (see e.g.
Mashayek  _et al._ , 2013), even in statistically stationary flows where
characteristic values of these four length scales can be identified,
parameterization of key properties of the flow could depend on at least three
independent nondimensional parameters determined from these scales. Choices
for these parameters are a characteristic Richardson number, $\mathrm{Ri}$, a
turbulent Froude number $\mathrm{Fr}$ and the activity parameter
$\mathrm{Re_{b}}$ (sometimes called the buoyancy Reynolds number and formally
distinct from a related integral-scale quantity; see Portwood  _et al._ (2016)
for further discussion) defined as

$\begin{gathered}\mathrm{Ri}\equiv\frac{N^{2}}{S^{2}}\equiv\left(\frac{L_{C}}{L_{O}}\right)^{4/3},\\\
\mathrm{Fr}\equiv\frac{\epsilon}{NE_{k}}\equiv\left(\frac{L_{O}}{L_{L}}\right)^{2/3},\mathrm{Re_{b}}\equiv\frac{\epsilon}{\nu
N^{2}}\equiv\left(\frac{L_{O}}{L_{K}}\right)^{4/3},\end{gathered}$ (3)

although alternative parameters can be defined, such as a ‘shear number’
$\mathrm{S^{\ast}}\equiv(L_{L}/L_{C})^{2/3}\equiv SE_{k}/\epsilon\equiv
1/(\mathrm{Fr}\mathrm{Ri}^{1/2})$ and others (Mater and Venayagamoorthy,
2014).

In terms of these parameters, the central challenge regarding $\kappa_{T}$
becomes that of determining the functional dependence of
$\kappa_{T}(\mathrm{Ri},\mathrm{Fr},\mathrm{Re_{b}})$. In a profoundly
influential paper, Osborn (Osborn, 1980) argued from consideration of the
turbulent kinetic energy equation that in stationary flows it is reasonable to
suppose that the buoyancy flux $B$, or equivalently, for the stationary flows
considered here, the dissipation rate of buoyancy variance $\chi$ defined as

$\chi\equiv\kappa\left\langle\frac{g^{2}}{\rho_{0}^{2}N^{2}}{\left|\nabla\rho\right|^{2}}\right\rangle,$
(4)

can be linearly related to the dissipation rate $\epsilon$ through a
‘turbulent flux coefficient’ or ‘mixing efficiency’ $\Gamma$, i.e.
$\chi=\Gamma\epsilon$. Osborn hypothesized an upper bound for $\Gamma\leq 0.2$
on semi-empirical grounds, although in practice, $\Gamma$ is implemented as a
constant saturating the upper bound(Gregg  _et al._ , 2018). The introduction
of $\Gamma$ transforms the fundamental issue of modeling $\kappa_{T}$ into
identifying the functional dependence of
$\Gamma(\mathrm{Ri},\mathrm{Fr},\mathrm{Re_{b}})$, since
$\kappa_{T}=\nu\Gamma\mathrm{Re_{b}}$ provided the flow is statistically
stationary.

Many models have been presented for the functional dependence of $\Gamma$
based on observations, experiments and numerical simulations in a variety of
different flows, with much recent activity (e.g. Ivey and Imberger, 1991;
Barry  _et al._ , 2001; Shih  _et al._ , 2005; Ivey  _et al._ , 2008; Maffioli
_et al._ , 2016; Venaille  _et al._ , 2017; Salehipour  _et al._ , 2016a;
Mashayek  _et al._ , 2017a; Zhou  _et al._ , 2017; Gregg  _et al._ , 2018;
Monismith  _et al._ , 2018) though with little sign of consensus, due
principally to three fundamental issues. Issue I is that disentangling
transient and/or reversible processes from the irreversible processes crucial
for quantifying and parameterizing mixing is extremely challenging. Issue II
is that, even for flows where the assumption of stationarity is appropriate,
stratification in a gravitational field and vertical shear both break
isotropy, and so the extent to which anisotropy is important is difficult to
determine, particularly when there is a relatively small dynamic range between
$L_{K}$ and $\min(L_{C},L_{O})$ (as encountered by Gargett  _et al._ (1984)).
The condition $L_{C}<L_{O}$ is associated with growing and stationary flow
configurations, as necessarily $\mathrm{Ri}<1$, suggesting that shear
instabilities have to be sufficiently vigorous to sustain turbulence. This has
the consequence that $L_{C}\gg L_{K}$ is a typical condition for high dynamic
range, i.e. $\mathrm{Ri}\mathrm{Re_{b}}\gg 1$ by the relations defined in (3)
and as suggested by Itsweire  _et al._ (1993). Finally, issue III is that
although in principle $\mathrm{Re_{b}}$, $\mathrm{Ri}$ and $\mathrm{Fr}$ are
all independent parameters, there is emerging evidence (e.g. Lucas and
Caulfield, 2017) that in many flows the parameters become correlated, and so
an apparent dependency of $\kappa_{T}$ or $\Gamma$ on one parameter is
actually associated with variation in another parameter.

To address all three issues, we consider the model flow of stationary
homogeneous sheared and stratified turbulence (S-HSST)(e.g. Holt  _et al._ ,
1992; Itsweire  _et al._ , 1993; Jacobitz  _et al._ , 1997; Shih  _et al._ ,
2005). In S-HSST, the production of turbulent kinetic energy by uniform mean
vertical shear exactly balances the dissipation rates of kinetic and potential
energy by molecular motion, addressing issue I by ensuring energetic
stationarity by design. When numerically simulated, the turbulence fills the
flow domain provided that $\mathrm{Re_{b}}$ is sufficiently large (Stillinger
_et al._ , 1983). More generally, stratified turbulence has been shown to
exhibit different dynamics depending on the value of $\mathrm{Re_{b}}$ (Shih
_et al._ , 2005; Hebert and de Bruyn Kops, 2006; Bartello and Tobias, 2013; de
Bruyn Kops, 2015; de Bruyn Kops and Riley, 2019). Gargett  _et al._ (1984)
observed in limited two-point statistics from ocean data that Kolmogorov-
Oboukhov-Corrsin (KOC) scaling may be observed in stratified turbulence when
$\mathrm{Re_{b}}\sim O(10^{3})$ and higher. de Bruyn Kops (2015) observes that
KOC scaling is generally not observed in DNS up to $\mathrm{Re_{b}}=220$,
though some statistics may be consistent with such scalings. An important open
question is thus how stratified turbulence behaves when
$\mathrm{Re_{b}}>O(100)$, a parameter space accessible by modern simulations
of S-HSST with sufficiently small $\mathrm{Fr}$ for the stratification to have
a first-order effect on the turbulent dynamics. However, to address issue II
there is the further constraint that $\mathrm{Ri}\mathrm{Re_{b}}\gg 1$.

In fact, in S-HSST $\mathrm{Ri}$ and $\mathrm{Fr}$ cease to be free, but
rather must adjust to ensure statistical stationarity. It has been empirically
observed in simulations that the various dynamic length scales ‘tune’ so that
$\mathrm{Fr}\approx 0.5$ and $\mathrm{Ri}\approx 0.16$, and equivalently
$\mathrm{S^{\ast}}\approx 5$. An immediate consequence of the emergence of an
apparently $\mathrm{Re_{b}}-$invariant stationary value of $\mathrm{Ri}$ in
S-HSST is that issue II can indeed be addressed at sufficiently large
$\mathrm{Re_{b}}$.

Perhaps more importantly, issue III is also addressed, as the dependence of
flow properties on $\mathrm{Re_{b}}$ can be probed independently of the other
parameters. Specifically, although the adjustment of $\mathrm{Fr}\approx 0.5$
and $\mathrm{Ri}\approx 0.16$ occurs at $\mathrm{Re_{b}}\sim O(10)$, other
dynamical scalings have been observed to change as $\mathrm{Re_{b}}$ increases
beyond $\mathrm{Re_{b}}\sim O(100)$ (Shih  _et al._ , 2005; de Bruyn Kops,
2015). The objective of this study is thus to investigate how the key
properties of the turbulence, in particular the energy partitioning, mixing
and small-scale anisotropy depend on variations of $\mathrm{Re_{b}}$. A
specific question to investigate is whether the postulated power law
$\Gamma\propto\mathrm{Re_{b}}^{-1/2}$ (Shih  _et al._ , 2005; Ivey  _et al._ ,
2008; Monismith  _et al._ , 2018) for sufficiently large $\mathrm{Re_{b}}$
occurs in this flow in which dependence on other flow parameters can be
completely eliminated.

##  II Simulations

The incompressible Navier-Stokes equations are considered, subject to the
nonhydrostatic Boussinesq approximation and coupled with equations for
continuity and buoyancy transport. A turbulent decomposition is performed
relative to a mean buoyancy gradient and a mean streamwise velocity gradient,
both in the vertical direction (see also Rogallo, 1981; Holt  _et al._ , 1992;
Jacobitz  _et al._ , 1997; Shih  _et al._ , 2000). We implement the system
with the same Fourier pseudo-spectral scheme described by (Almalkie and de
Bruyn Kops, 2012; de Bruyn Kops, 2015), except that an additional shear term
is handled by an integration factor (Brucker  _et al._ , 2007; Chung and
Matheou, 2012; Sekimoto  _et al._ , 2016).

Stationarity is induced by fixing a value of $\nu$, choosing a target
turbulent kinetic energy $E_{t}$ then adjusting the Richardson number via $g$
(c.f. Taylor  _et al._ , 2016) using a mass-spring-damper control system:

$c_{0}S\mathrm{Ri}^{\prime}(t)+2\alpha\omega\tilde{E}_{k}^{\prime}(t)+\omega^{2}(\tilde{E}_{k}(t)-1)=0$
(5)

where the prime notation denotes a temporal derivative,
$\tilde{E}_{k}(t)\equiv E_{k}/E_{t}$ is the normalized turbulent kinetic
energy, $\omega$ is the characteristic frequency of oscillation and $\alpha$
is a dimensionless damping factor. The control system has been derived by
assuming that the kinetic energy follows a second order linear system (e.g.
Rao and de Bruyn Kops, 2011), and then by applying the first-order
approximation $\tilde{E}_{k}^{\prime}(t)\approx
c_{0}S(\mathrm{Ri}(t)-\mathrm{Ri}_{c})$ such that
$\tilde{E}_{k}^{\prime\prime}(t)\approx c_{0}S\mathrm{Ri}^{\prime}(t)$. The
parameter $c_{0}\approx-1$ is supported by Jacobitz  _et al._ (1997), the
characteristic frequency $\omega$ is determined by the mean shear, and a
damping coefficient $\alpha=1.5$ was found to work well.

Crucially, following this procedure, $\mathrm{Ri}\approx 0.16$ emerges without
presupposition for all our cases, as shown in the table 1 along with other
emergent parameters; flow statistics are averaged over a period of $St\approx
100$ unless noted otherwise.

Case | $\mathrm{Re_{b}}$ | $\mathrm{Ri}$ | $\mathrm{Fr}$ | $N_{x}$
---|---|---|---|---
SHSST-R1 | 36 | 0.163 | 0.46 | 1024
R2 | 48 | 0.159 | 0.47 | 1280
R3 | 59 | 0.162 | 0.48 | 1536
R4 | 81 | 0.154 | 0.50 | 1792
R5 | 110 | 0.155 | 0.52 | 2048
R6 | 160 | 0.157 | 0.48 | 3072
R7 | 240 | 0.156 | 0.48 | 4096
R8 | 390 | 0.146 | 0.46 | 6144
R9 | 550 | 0.163 | 0.45 | 8192
R10 | 900 | 0.152 | 0.42 | 9600
Table 1:  Simulation parameters. $N_{x}$ is the number of grid points in the
$x$-direction and the grid spacing is isotropic.

Furthermore, the dissipation rate $\epsilon$ also ‘tunes’ such that
$\mathrm{Fr}\approx 0.5$. As the dissipation rate is an emergent quantity, the
smallest length scales are resolved by adjusting the resolution such that
$k_{max}L_{K}\approx 2$, where $k_{max}$ is the largest Fourier domain
wavenumber. We also found it necessary to use a relatively large domain with
$L_{x}/L_{y}=2$, $L_{x}/L_{z}=4$ and $L_{x}/L_{L}\approx 40$, where $L_{x}$,
$L_{y}$, and $L_{z}$ are the dimensions of the domain, in order to support the
anisotropic large scales of the flow.

##  III Results and discussion

###  III.1 Energetics

We stress that, although $E_{k}\approx E_{t}$ is enforced, the parameters
$\mathrm{Ri}$, $\mathrm{Fr}$ and $\mathrm{Re_{b}}$ emerge from the
simulations, as do the structure of the turbulence and the scalar. We consider
the ratio of the potential energy to kinetic energy, $\mathrm{R_{PK}}\equiv
E_{p}/E_{k}$, and the ratio of the potential energy to the kinetic energy of
vertical motion, $\mathrm{R_{PV}}\equiv E_{p}/E_{v}$ where the energies are
defined as:

$E_{v}\equiv\frac{1}{2}\langle{w^{2}}\rangle\;\text{and}\;E_{p}\equiv\frac{1}{2}\left<\frac{g^{2}/\rho_{0}^{2}}{N^{2}}{\rho^{2}}\right>\text{
.}$

Energy-partitioning is a critical component to mixing models wherein Reynolds
number, or $\mathrm{Re_{b}}$, dependence is often omitted and the mixing is
assumed to be a function of $\mathrm{Ri}$ (Osborn and Cox, 1972; Schumann and
Gerz, 1995; Pouquet  _et al._ , 2018). Furthermore, in ‘strongly’ stratified
turbulent flow, Billant and Chomaz (2001) suggest that there should be
approximate equipartition between potential and kinetic energy, i.e.
$\mathrm{R_{PK}}\approx 1$, an assumption also used by Lindborg (2006). Figure
1 illustrates that this basic assumption of equipartition is not appropriate.
Perhaps more interesting, the energy ratios decrease by a factor of two
relative to the lowest $\mathrm{Re_{b}}$ case until $\mathrm{Re_{b}}\approx
300$, after which the energy ratios appear to remain constant with anisotropy
of velocity variance explaining the different behaviors of $\mathrm{R_{PK}}$
and $\mathrm{R_{PV}}$.

Figure 1:  Energy partitions as a function of $\mathrm{Re_{b}}$. Bars indicate
upper and lower quartile measurements of instantaneous quantities.

###  III.2 Small-scale anisotropy

The observation of anisotropy in the energy partitions of the previous section
is particularly important, due not least to the inevitable requirement to
estimate energetic dissipation rates in circumstances where each component of
the rate-of-deformation tensor is not available. The existence of anisotropy
at dissipative length scales would require dynamical estimates of $\kappa_{T}$
to account for anisotropy. Therefore, it is necessary to evaluate common
dissipation surrogate models to assess their applicability and to characterize
small-scale anisotropy in these flows. Here, we use single component surrogate
models which rely on a single derivative and the isotropy assumption:

$\displaystyle\tilde{\epsilon}_{ij}$
$\displaystyle=\begin{cases}15\nu\left\langle\big{(}\frac{\partial
u_{i}}{\partial x_{j}}\big{)}^{2}\right\rangle&\text{if $i=j$}\\\
15/2\nu\left\langle\big{(}\frac{\partial u_{i}}{\partial
x_{j}}\big{)}^{2}\right\rangle&\text{if $i\neq j$}\end{cases}$ (6a)
$\displaystyle\tilde{\chi}_{j}$
$\displaystyle=3\kappa\left\langle\frac{g^{2}}{\rho_{0}^{2}N^{2}}\left(\frac{\partial\rho}{\partial
x_{j}}\right)^{2}\right\rangle\text{ .}$ (6b)

In flows that are inherently anisotropic due to mean shear or mean flux, it is
widely assumed that small-scale isotropy is a reasonable assumption when the
scale separation is large (e.g. Gargett  _et al._ , 1984). However, evidence
suggests that isotropy assumptions can be very inaccurate in stably stratified
flows at finite $\mathrm{Re_{b}}$ (Itsweire  _et al._ , 1993; Hebert and de
Bruyn Kops, 2006; de Bruyn Kops, 2015; de Bruyn Kops and Riley, 2019).

Figure 2:  Instantaneous measurements of small-scale anisotropy. The dashed
line represents perfect isotropy at small scales.

The validity of (6) is tested with the help of figure 2 in which is plotted
the relative error $(\tilde{\epsilon}_{ij}-\epsilon)/\epsilon$ of the various
single-component dissipation estimates. As with energy partition, there is an
apparently asymptotic regime for $\mathrm{Re_{b}}\gtrapprox 300$ in which the
assumption of isotropy applied to dissipation rate becomes valid within
approximately 15% error. Nevertheless, there is evidence of small-scale
anisotropy as expected from the analysis in Durbin and Speziale (1991), which
shows that dissipation-range isotropy should only exist when
$\mathrm{S^{\ast}}\ll 1$; in these simulated flows $\mathrm{S^{\ast}}\approx
5$. However, we observe that $\tilde{\epsilon}_{22}$, $\tilde{\epsilon}_{33}$,
and $\tilde{\chi}_{2}$ yield good estimates of the dissipation rates even in
cases in which the small scales are strongly anisotropic.

###  III.3 Mixing coefficient

We plot the turbulent mixing coefficient $\Gamma\equiv\chi/\epsilon$ in figure
3a. Whereas there is some modest decrease with increasing $\mathrm{Re_{b}}$
from the peak $\Gamma$ of approximately 0.19, remarkably close to Osborn’s
suggested upper bound, there is no evidence of the commonly-suggested
$\mathrm{Re_{b}}^{-1/2}$ scaling (e.g. Shih  _et al._ , 2005; Ivey  _et al._ ,
2008; Monismith  _et al._ , 2018). A decrease in mixing efficiency with
respect to $\mathrm{Re_{b}}$ is observed until $\mathrm{Re_{b}}\approx 200$,
above which $\Gamma\approx 0.17$. There is no evidence of the ‘energetic’
regime of Shih  _et al._ (2005) for $\mathrm{Re_{b}}>100$, suggesting that
transient multi-parameter effects may well be relevant in the evolution of
their flows such that $\mathrm{Re_{b}}$ is not always an independent parameter
(see issue III as described in §1). There is evidence that the upper bound
proposed by Osborn is a useful estimate, at least in flows where the
underlying assumption of stationarity is well-justified, as such flows
naturally adjust to $\mathrm{Ri}\approx 0.16$ and $\mathrm{Fr}\approx 0.5$.

Alternative descriptions of the mixing behavior are instructive. In the
asymptotically high $\mathrm{Re_{b}}$ regime, the buoyancy variance
dissipation rate adjusts such that $\chi\approx E_{p}N$ as evidenced by
figures 1 and 3. Although this is a natural scaling, it is of interest that
the $O(1)$ constant is actually extremely close to $1$. A second instructive
description of mixing is the turbulent Prandtl number
$\mathrm{Pr_{T}}=\kappa_{M}/\kappa_{T}$, where $\kappa_{M}$ is the eddy
diffusivity of momentum:

$\kappa_{M}\equiv\frac{\langle uw\rangle}{S}=\frac{P}{S^{2}},$ (7)

and $P$ is the production rate of turbulent kinetic energy. Therefore,

$\mathrm{Pr_{T}}=\frac{N^{2}}{S^{2}}\frac{P}{B}\equiv\frac{\mathrm{Ri}}{\mathrm{Rf}}\approx\frac{\mathrm{Ri}(1+\Gamma)}{\Gamma},$
(8)

where $\mathrm{Rf}$ is the ‘flux Richardson number’, and the relationships
$P\approx B+\epsilon=\chi+\epsilon$ in a stationary flow has been used.
$\mathrm{Pr_{T}}$ is plotted in 3b, and proves to be close to one for all
values of $\mathrm{Re_{b}}$. This is perhaps unsurprising, but makes clear
that the turbulent processes that mix heat and momentum in these flows are
highly coupled, and in particular that the stratification is not sufficiently
‘strong’ to modify the turbulent processes greatly, but rather that the
irreversible conversion of kinetic into potential energy occurs in a balanced,
equilibrated way.

Using mixing length arguments, Odier  _et al._ (2009) defined:

$L_{\rho}\equiv\left(\frac{B}{N^{2}S}\right)^{1/2};\
L_{m}\equiv\left(\frac{P}{S^{3}}\right)^{1/2};$ (9)

such that $L^{2}_{m}/L^{2}_{\rho}=Pr_{T}$. Therefore, since for our flows
$Pr_{T}\approx 1$ and $P\approx(1+\Gamma)\epsilon$, it is apparent that the
density length scale at the heart of the model presented by Ivey  _et al._
(2018) is coupled to the Corrsin scale by $L_{\rho}/L_{C}=\sqrt{1+\Gamma}$,
and the key parameter $D\equiv\chi/\epsilon=1/\Gamma$ has nearly fixed value.

##  IV Conclusions

We have used a canonical controlled flow to evaluate the isolated effects of
high dynamic range in sheared stratified turbulent flow. The energy
partitioning varies non-trivially where $\mathrm{Re_{b}}\lessapprox 300$ above
which an apparently asymptotic regime is entered, consistent with (Gargett
_et al._ , 1984; de Bruyn Kops, 2015). The flow retains measurable anisotropy
at dissipation scales, as suggested by the analyses of Durbin and Speziale
(1991), at all values of $\mathrm{Re_{b}}$ we have considered. Nevertheless,
some single-component surrogates exist which accurately estimate dissipation
rates via isotropy assumptions even for smaller $\mathrm{Re_{b}}$.

Figure 3:  a: Variation of turbulent flux coefficient $\Gamma$ with
$\mathrm{Re_{b}}$. The upper bound proposed by Osborn is indicated by the gray
dashed line where $\Gamma\approx 0.2$, and the $\mathrm{Re_{b}}^{-1/2}$-based
parameterization suggested in Shih  _et al._ (2005) is plotted for their
‘energetic’ regime of $\mathrm{Re_{b}}>100$. b: Variation of turbulent Prandtl
number $\mathrm{Pr_{T}}$ with $\mathrm{Re_{b}}$.

The results presented here seem to indicate that the effects of the large
dynamic-range regimes explored by Gargett et al. are strongly influenced by
asymptotic scalar density dynamics, rather than by the velocity field
independently. Nevertheless, the measured mixing is a much weaker function of
$\mathrm{Re_{b}}$ compared to some proposed scalings (Shih  _et al._ , 2005),
with the turbulent Prandtl number $\mathrm{Pr_{T}}\approx 1$ and the turbulent
mixing coefficient $\Gamma$ near the classical bound of $0.2$ as suggested by
Osborn, although decreasing slightly for $\mathrm{Re_{b}}\lessapprox 300$. The
application of these results to higher Prandtl numbers merits further study,
where dynamic range arguments indicate transition at lower $\mathrm{Re_{b}}$
as $\mathrm{Pr}$ is increased.

We stress that the flow we have considered has, by design, controlled
dependence on all other parameters. $\mathrm{Ri}\approx 0.16$ and
$\mathrm{Fr}\approx 0.5$ naturally emerge to ensure stationarity. Therefore,
we conjecture that observed apparent variation of mixing properties with
$\mathrm{Re_{b}}$ (Monismith  _et al._ , 2018) can be explained by breaking
one or the other of these constraints, i.e. the variation is due to either
transient effects, well-known to lead to strong variation in mixing properties
(e.g Salehipour  _et al._ , 2016b; Mashayek  _et al._ , 2017b), or ‘hidden’
and perhaps correlated variation with other parameters (e.g Zhou  _et al._ ,
2017). Additionally, there could also be as yet unquantified strong dependence
on initial or boundary conditions, whereas the flow we have considered is
isolated from such effects. To develop robust mixing parameterizations, it is
necessary to develop appropriate models to capture the effects of such
conditions, informed by and generalizing from such controlled, idealized flows
as considered here.

###### Acknowledgements.

This work was funded by the U.S. Office of Naval Research via grant
N00014-15-1-2248. High performance computing resources were provided through
the U.S. Department of Defense High Performance Computing Modernization
Program by the Army Engineer Research and Development Center, the Army
Research Laboratory and the Navy DSRC under Frontier Project FP-CFD-FY14-007.
The research activity of C.P.C. is supported by EPSRC Programme Grant
EP/K034529/1 entitled ‘Mathematical Underpinnings of Stratified Turbulence’.
Valuable comments from anonymous reviewers have substantially improved the
clarity of the discussion.

## References

  * Gargett  _et al._ (1984) A. Gargett, T. Osborn, and P. Nasmyth, J. Fluid Mech. 144, 231 (1984).
  * Shih  _et al._ (2005) L. H. Shih, J. R. Koseff, G. N. Ivey, and J. H. Ferziger, J. Fluid Mech. 525, 193 (2005).
  * Mater and Venayagamoorthy (2014) B. D. Mater and S. K. Venayagamoorthy, Phys. Fluids 26 (2014), 10.1063/1.4868142.
  * Ivey  _et al._ (2008) G. N. Ivey, K. B. Winters, and J. R. Koseff, Annu. Rev. Fluid Mech. 40, 169 (2008).
  * Ferrari and Wunsch (2009) R. Ferrari and C. Wunsch, Ann. Rev. Fluid Mech. 41, 253 (2009).
  * Gregg  _et al._ (2018) M. C. Gregg, E. A. D’Asaro, J. J. Riley, and E. Kunze, Ann. Rev. Mar. Sci. 10, 443 (2018).
  * Ivey  _et al._ (2018) G. N. Ivey, C. E. Bluteau, and N. L. Jones, J. Geophys. Res.-Oceans 123, 346 (2018).
  * Mashayek  _et al._ (2013) A. Mashayek, C. P. Caulfield, and W. R. Peltier, J. Fluid Mech. 736, 570 (2013).
  * Portwood  _et al._ (2016) G. D. Portwood, S. M. de Bruyn Kops, J. R. Taylor, H. Salehipour, and C. P. Caulfield, J. Fluid Mech. 807, R2 (2016).
  * Osborn (1980) T. R. Osborn, J. Phys. Oceanogr. 10, 83 (1980).
  * Ivey and Imberger (1991) G. N. Ivey and J. Imberger, J. Phys. Oceanogr. 21, 650 (1991).
  * Barry  _et al._ (2001) M. E. Barry, G. N. Ivey, K. B. Winters, and J. Imberger, J. Fluid Mech. 442, 267 (2001).
  * Maffioli  _et al._ (2016) A. Maffioli, G. Brethouwer, and E. Lindborg, J. Fluid Mech. 794, R3 (2016).
  * Venaille  _et al._ (2017) A. Venaille, L. Gostiaux, and J. Sommeria, J. Fluid Mech. 810, 554 (2017).
  * Salehipour  _et al._ (2016a) H. Salehipour, W. R. Peltier, C. B. Whalen, and J. A. MacKinnon, Geophys. Res. Lett. 43, 3370 (2016a).
  * Mashayek  _et al._ (2017a) A. Mashayek, H. Salehipour, D. Bouffard, C. P. Caulfield, R. Ferrari, M. Nikurashin, W. R. Peltier, and W. D. Smyth, Geophys. Res. Lett. 44, 6296 (2017a).
  * Zhou  _et al._ (2017) Q. Zhou, J. R. Taylor, C. P. Caulfield, and P. F. Linden, J. Fluid Mech. 823, 198 (2017).
  * Monismith  _et al._ (2018) S. G. Monismith, J. R. Koseff, and B. L. White, Geophys. Res. Lett. 45, 5627 (2018).
  * Itsweire  _et al._ (1993) E. C. Itsweire, J. R. Koseff, D. A. Briggs, and J. H. Ferziger, J. Phys. Oceanogr. 23, 1508 (1993).
  * Lucas and Caulfield (2017) D. Lucas and C. P. Caulfield, J. Fluid Mech. 832, R1 (2017).
  * Holt  _et al._ (1992) S. E. Holt, J. R. Koseff, and J. H. Ferziger, J. Fluid Mech. 237, 499 (1992).
  * Jacobitz  _et al._ (1997) F. G. Jacobitz, S. Sarkar, and C. W. Van Atta, J. Fluid Mech. 342, 231 (1997).
  * Stillinger  _et al._ (1983) D. C. Stillinger, K. N. Helland, and C. W. Van Atta, J. Fluid Mech. 131, 91 (1983).
  * Hebert and de Bruyn Kops (2006) D. A. Hebert and S. M. de Bruyn Kops, Geophys. Res. Let. 33, L06602 (2006).
  * Bartello and Tobias (2013) P. Bartello and S. M. Tobias, J. Fluid Mech. 725, 1 (2013).
  * de Bruyn Kops (2015) S. M. de Bruyn Kops, J. Fluid Mech. 775, 436 (2015).
  * de Bruyn Kops and Riley (2019) S. M. de Bruyn Kops and J. J. Riley, J. Fluid Mech. 860, 787 (2019).
  * Rogallo (1981) R. S. Rogallo,  _Numerical experiments in homogeneous turbulence_ , Tech. Rep. TM-81315 (NASA, 1981).
  * Shih  _et al._ (2000) L. H. Shih, J. R. Koseff, J. H. Ferziger, and C. R. Rehmann, J. Fluid Mech. 412, 1 (2000).
  * Almalkie and de Bruyn Kops (2012) S. Almalkie and S. M. de Bruyn Kops, J. Turbul. 13, 1 (2012).
  * Brucker  _et al._ (2007) K. A. Brucker, J. C. Isaza, T. Vaithianathan, and L. R. Collins, J. Comp. Phys. 225, 20 (2007).
  * Chung and Matheou (2012) D. Chung and G. Matheou, J. Fluid Mech. 696, 434 (2012).
  * Sekimoto  _et al._ (2016) A. Sekimoto, S. Dong, and J. Jiménez, Phys. Fluids 28, 035101 (2016).
  * Taylor  _et al._ (2016) J. R. Taylor, E. Deusebio, C. P. Caulfield, and R. R. Kerswell, J. Fluid Mech. 808, R1 (2016).
  * Rao and de Bruyn Kops (2011) K. J. Rao and S. M. de Bruyn Kops, Phys. Fluids 23, 065110 (2011).
  * Osborn and Cox (1972) T. R. Osborn and C. S. Cox, Geophysical & Astrophysical Fluid Dynamics 3, 321 (1972).
  * Schumann and Gerz (1995) U. Schumann and T. Gerz, J. Appl. Meteorol. 34, 33 (1995).
  * Pouquet  _et al._ (2018) A. Pouquet, D. Rosenberg, R. Marino, and C. Herbert, J. Fluid Mech. 844, 519 (2018).
  * Billant and Chomaz (2001) P. Billant and J.-M. Chomaz, Phys. Fluids 13, 1645 (2001).
  * Lindborg (2006) E. Lindborg, J. Fluid Mech. 550, 207 (2006).
  * Durbin and Speziale (1991) P. Durbin and C. Speziale, J. Fluids Eng. 113, 707 (1991).
  * Odier  _et al._ (2009) P. Odier, J. Chen, M. K. Rivera, and R. E. Ecke, prl 102, 134504 (2009).
  * Salehipour  _et al._ (2016b) H. Salehipour, C. C. P., and W. Peltier, J. Fluid Mech. 803, 591 (2016b).
  * Mashayek  _et al._ (2017b) A. Mashayek, C. P. Caulfield, and W. R. Peltier, J. Fluid Mech. 826, 522 (2017b).