\section{Introduction}
\label{sec:intro}

\noindent \textbf{Motivation.}
Inferring network states from partial observations is an important problem in many applications, such as in epidemics, social networks and Internet/AS graphs. Consider critical infrastructure networks like water distribution system or electric grids --- their rapid functional restoration is a priority right after a disaster such as an earthquake~\cite{gomez_integrating_2019}. The sparse installation of real time monitoring systems installed in these networks (e.g. SCADA systems) imply that only a partial observation of the network is available post disaster i.e., we only get a sample of the functional states of the components (`failed' or not). However, for decision makers to devise an effective recovery plan, the information on functional states of the entire network is required. Hence, it is crucial to infer the states of the components in such networks after a disaster.

Despite its importance, there has been limited research in this space. Most past work has only looked into the so-called \emph{point queries}, i.e. one is given some direct sample of failed components. In epidemics, Shah et al.~\cite{shah2010detecting} proposed inferring missing infections and source detection given a sample of infections.
Indeed, in network tomography, all prior work has focused on independent failures~\cite{jin2006network}.
Recently, there has been some progress on reconstructing network states assuming that the failures are geographically correlated given a sample of the failed components
~\cite{adhikari2018near}. As mentioned earlier, a common theme among these works is that the observed failures are a direct sample of failed components.

A different class of queries, namely, \emph{connectivity queries} also occur naturally. For example,
consider an electric power grid, where transmission line connecting two component may fail~\cite{chen2017hotspots}. Such failures cause disconnections in the network.
SCADA systems (which oversee sensors in a power grid) provide information on
whether or not there is a power supply at a specific demand node (a consumer node). These demand nodes are supplied by identifiable supply nodes (power generators),
thus, a stop of supply to a demand node can be mapped to disconnectivity
between it and the supply nodes. Similarly, the continuous supply on a demand node indicates that there is at least one connected path to the supply nodes.
A motivating example is depicted in Figure~\ref{fig:teaser}. Here, node $1$ is a source node (power plant), the set $\{3,4,6,7,8,9\}$ is the set of demand nodes (consumers), and the remaining nodes $\{ 2,4,5\}$ are transmission nodes. If due to a disaster edges in the set $\{(2,3), (4,7), (5,9) \}$ fail, this would disconnect nodes $7$ and $9$.

\linespread{1.25}
\begin{figure}[!t]
\small

\centering
\begin{tabular}{cc}
      \includegraphics[width=0.4\textwidth]{images/conn_paper_img2.pdf} &  \\
      \vspace{5pt} (a) & \\
     \begin{tabular}[b]{|p{2.8cm}|p{4.0cm}|}
    \hline
    \textbf{Notation} & \textbf{Physical set in Fig. 1(a)} \\
    \hline
    Supply nodes $V_S$ & Electric power plant \\ & $= \{ 1\}$\\\hline
    Transshipment nodes $V_T$ &Transmission towers   $= \{ 2,4,5\}$\\\hline
    Demand nodes $V_D$& Houses/customers $= \{ 3,6,7,8,9\}$\\\hline
    Serviced nodes $\mathcal{S}$ & Houses with electricity \\ & $ = \{ 3,6,8\} \subseteq V_D $\\\hline
    True failure set $I$ & Transmission lines failed \\ & $= \{ (2,3), (4,7), (5,9)\}$\\\hline
    Disaster scenarios $\mathcal{O}$ & Possible epicenters $ = \{ A, B, C\}$ \\\hline
\end{tabular}
\end{tabular}
(b)
\caption{
(a) An example of our problem using an electric power network disrupted by an earthquake. We refer to the demand nodes (houses) still connected to a supply node (power plant) as serviced nodes, i.e. houses with electricity. Thus, given a sample of serviced nodes $\{3,8\}$ (connectivity probes), a smaller sample of the set of failed edges $I$, and knowing that the epicenter can be A, B, or C, we want to infer the rest of true failure set (edges with crosses).
Note that failed edge $\{2,3\}$ does not affect the serviceability of any node; therefore, using only connectivity probes will not allow us to discover this failure.
(b) Description of the physical sets in (a), and its association with our notation.
}
\label{fig:teaser}
\end{figure}

\noindent \textbf{Why study connectivity queries?}
After a natural disaster, connectivity queries are \emph{more practical} to obtain than point queries.
Information about supplied demand and customer service (a natural source of connectivity queries) is often relatively accessible---customers nowadays have multiple ways of reporting issues, including even social media \cite{mcclendon_leveraging_2013}---and has been proven to be useful to localize failures \cite{ji_resilience_2017}. In our example, a sample of the nodes connected to the source are observed i.e, we observe a subset of $\{3,6,8\}$ thanks to signals from sensors and/or consumer reports.

Motivated by the reasons above, we investigate the usefulness of adding connectivity queries to the point queries (collectively referred to as `joint queries') in inferring network failures. We find that using \emph{only} connectivity queries for inference is very challenging. For instance, in our example, no demand node is affected by the failure of edge $\{2,3\}$, which makes it undiscoverable to connectivity queries.
Hence, we propose using `joint queries' where we rely primarily on connectivity queries but also use a small sample of point queries to boost performance. To this end, we propose to the leverage the powerful Minimum Description Length (MDL)~\cite{grunwald04tutorial} concept, to find out both the number of failures and their identity. MDL is especially suited for this task, as it gives us a principled way to formulate our problem.

\noindent \textbf{Contributions.}
Our main contributions in this paper are as follows:
\begin{enumerate}
\item \textit{Formulation and algorithm.} We formulate \prob{} to leverage joint probes for network state inference. We propose \textsc{JointPathMap}\xspace{} algorithm to solve the problem.
\item \textit{Extensive experiments on real datasets}. We use domain-expert simulations by civil engineers of three real network topologies in the aftermath of an earthquake. Our extensive experiments demonstrate that \textsc{JointPathMap}\xspace{} significantly outperforms the baselines.
\end{enumerate}

\noindent To the best of our knowledge, our work is the first to propose connectivity queries and joint queries for network state inference.

\section{Preliminaries}  \label{sec:prelim}

\subsection{Context and Setup}
We are given an undirected graph $G(V,E)$ that represents an infrastructure network. Links in $E$ represent transmission systems and
set $V$ contains three disjoint sets of nodes: supply nodes $V_S$, demand nodes $V_D$, and transshipment nodes $V_T$. Thus, $V=V_S \cup V_D \cup V_T$. Supply nodes are where flow (e.g. electricity, water) is originated, demand nodes are where flow is received and distributed to final users, and transshipment nodes are where flow is transferred to other link(s).
We assume supply and demand nodes are related by many-to-many relationships (i.e. one supply node can service many demand nodes and one demand node can be serviced by many supply nodes) as it is observed in electric and water networks. This kind of structure appears in many different critical infrastructure networks \cite{chen2017hotspots}.

In addition, we are given a set $\mathcal{O}$ containing disaster scenarios that characterize the most likely disaster hazards. For example, in case of earthquakes, they can be described by a magnitude and a geographical epicenter; even though there is a myriad of possibilities for this, typically there is a finite set of earthquake scenarios (our set $\mathcal{O}$) that can characterize most of the seismic hazard~\cite{adachi_comparative_2010}.  Once a scenario is characterized, the effects on the infrastructure network are driven by soil and structural properties in a stochastic manner.

\subsection{Failure Model} \label{sec:failure_model}
The failure model leverages `geographically correlated' failures, which have been extensively used in prior work both in the computer science~\cite{agarwal:infocom11} and civil engineering communities \cite{gomez_integrating_2019}.
At an abstract level, given a disaster scenario $o \in \mathcal{O}$, some set of edges $I \subseteq E$ fail.
Let $p(o)$ be the probability of a disaster scenario $o$ to had happen; hence, by previous definition, $\sum_{o\in\mathcal{O}} p(o) = 1$.
Let $F(e|o)$ be the probability of edge $e$ to fail given $o$, e.g.
probability of an electric transmission line to fail given an earthquake with epicenter $o$.
These probability distributions are calculated considering not only the distance to the epicenter of $o$, but also the type of structure, construction materials, and condition.
In addition, we assume all edge failures depend only on the disaster's effect, i.e. failures are independent from each other.
Both $\mathcal{O}$ and $F(e|o), \forall e\in E, \forall o\in \mathcal{O}$ are given by civil engineering guidelines, as in our experiments.

\section{Formulation and Approach}
\label{sec:problem_formulation}
\noindent Here we systematically formulate our problem based on network connectivity using the MDL principle. See Table~\ref{tab:notation} for notation.

\subsection{Connectivity Probes}
Once a failure has taken place in the network, an operator performs the so-called ``connectivity queries'' on the network, as described next, to gather some information on the status.
Connectivity queries (coming from sensors and/or customer reports) serve as indirect observers for the network states.  As exemplified in Figure~\ref{fig:teaser}, when failures happen, they directly effect the \emph{serviceability} of some demand nodes in the network. Demand nodes require an active path to a nearby supply node to be serviced (e.g. receive electricity).
Thus, due to failures because of the disaster, some demand nodes become non-serviced.
We define a `sensor' set $\mathcal{S} \subseteq V_D$ as the set of demand nodes which are still serviced after removing failure edge set $I$ from $G$.  The operator will finally receive a \emph{connectivity probe} set $\mathcal{Q}_c \subseteq \mathcal{S}$. We can assume that it is chosen uniformly at random (with some given probability, say $\gamma_c$), from $\mathcal{S}$.
Therefore, we have the following dependency relationships:
\setlength\abovedisplayskip{5pt}
\setlength\belowdisplayskip{1pt}
\begin{equation}
o \rightarrow I \rightarrow \mathcal{S} \rightarrow \mathcal{Q}_c
\label{eq:causal}
\end{equation}

\subsection{Joint Probes}
While connectivity queries are more practical to obtain, they are often not enough. If the demand nodes have multiple paths to the supply nodes (i.e., the network is redundant) the connectivity probes $Q_c$ fail to provide enough information to infer $I$ and $o$ correctly. We saw this in Figure~\ref{fig:teaser}, where failure of edge $\{2,3\}$ was not disconnecting any node, thus, it does not affect serviceability.
It turns out that such phenomenon is quite common in real critical infrastructure networks as we empirically demonstrate it in Section~\ref{sec:expts}. As the failures and serviceability are not interdependent in such networks, the connectivity probes do not relay any information regarding the failures; making it near impossible to infer the failures given the connectivity probes.
The problem is compounded by the fact that the number of failure sets which do not effect the serviceability is exponential in terms of number of such edges. As we show in our experiments, using only connectivity probes returns failure set $I$ where recall is low as it fails to identify the edges which have no effect on the serviceability.

Hence, a small sample of point probes of the failures may be helpful in boosting the power of connectivity probes as adding a sample of $I$ informs how many nodes to fail and may help to idenfity the true disaster scenario $o$.
Thus, additionally to $\mathcal{Q}_c$ we assume that a small set of (costlier) point queries $\mathcal{Q}_I \subseteq I$ (a sample of actual components that have failed) is also available.
For the ease of modeling, we assume that edges $ e \in \mathcal{Q}_I$ are sampled uniformly at random with probability $\gamma_I$~\cite{adhikari2018near}. In the rest of the paper, we refer to $\mathcal{Q}_I$ as the failure probe set and refer to both $\mathcal{Q}_I$ and $\mathcal{Q}_c$ together as the joint probes.

\subsection{MDL}
Note that we do not know apriori how many failures are present in the system. Hence, intuitively, to formulate the problem above, we need a model selection framework which finds this automatically in a principled fashion.
To this end, we propose to use the information-theoretic Minimum Description Length (MDL) principle~\cite{grunwald04tutorial}.
We use two-part MDL, also known as the sender-receiver framework as we are specifically interested in the model. Our goal is to transmit the given set of probes $\mathcal{Q} = \mathcal{Q}_c \cup \mathcal{Q}_I$ (the `data') from sender to receiver in the least cost (in bits), by assuming that both of them know the layout of the network $G$, the sampling rates $\gamma_c$ and $\gamma_I$, the possible failure sources $\mathcal{O}$ and the conditional failure parobabilities $F's$.
The MDL cost function consists of two parts: (a) Model cost includes the complexity of the selected model that explains the state of the network; and (b) Data cost that represents the cost of sending the given probe data $\mathcal{Q}$ given the model chosen. Formally, given a set of models $\mathcal{M}$, MDL claims the best model $M^*$ is the one that minimizes $\mathcal{L}(M)+ \mathcal{L}(\mathcal{D} \given M)$, in which $\mathcal{L}(M)$ is the model-cost (length in bits to describe model $M$), and $\mathcal{L}(\mathcal{D} \given M)$ is the data-cost (the length in bits to describe $\mathcal{D}$ using $M$).

\subsection{Model Space and Cost}\label{subsec:sizeformmodelcost}
The first step in our formulation is to identify the right model \textit{space}.
Our failure process states that the failure set $I$ is induced due to a particular disaster scenario $o \in \mathcal{O}$. Hence, we include both $o$ and $I$ in the model.
In real networks different disasters induce failures in totally different regions of the network. Hence, marginalizing over different scenarios and picking the most likely failure set $I$ over all disasters\cite{adhikari2018near} does not work well.

\begin{table}[t]
\centering

\caption{Table of symbols}
\begin{tabular}{|p{1.5cm}|p{5.5cm}|}
    \hline
    \textbf{Symbol} & \textbf{Definition} \\
    \hline
    $V_S, V_T, V_D$ & Supply, transshipment, and demand node set, respectively \\\hline
    $\mathcal{O}$ & Set of disaster scenarios \\\hline
    $I$& Set of failed edges (failure set) \\\hline
    $F(e|o)$& Failure probability of edge $e$ given disaster $o$ \\\hline
    $\mathcal{S}$ & Set of serviced nodes (serviced set) \\\hline
    $\gamma_c, \gamma_I$ & Probability for (uniformly) sampling probes from sets $\mathcal{S}$ and $I$, respectively \\\hline
    $\mathcal{Q}_c$ & Connectivity probe set (connectivity queries) sampled from set $\mathcal{S}$ \\\hline
    $\mathcal{Q}_I$ & Point probe set (point queries) sampled from set $I$ \\\hline
    $h(G, \mathcal{S}, I)$ & Indicates if the concurrence of $G$, $\mathcal{S}$, and $I$ is possible (returns 1) or not (returns 0) \\\hline
\end{tabular}
\label{tab:notation}
\end{table}

We want to primarily infer only $I$, but this set by itself is not enough to be the model because there is no direct relationship between $I$ and $Q_c$; thus, it would not be useful to explain the data. Therefore, we use all three of $o$, $I$ and $\mathcal{S}$ as our model: $\mathcal{M} = (o, \mathcal{S}, I)$. In other words, we first send the disaster scenario $o$, then send the true serviced set $\mathcal{S}$, followed by the true failed set $I$. We then identify the actual joint probes set $\mathcal{Q}_c$ and $\mathcal{Q}_I$ as the data.

\setlength\abovedisplayskip{1pt}
\setlength\belowdisplayskip{1pt}
The MDL model cost
$\mathcal{L}(o, \mathcal{S}, I)$ has the following three components:
\begin{align*}
\mathcal{L}(o, \mathcal{S}, I) = \mathcal{L}(o) +\mathcal{L}(\mathcal{S} \given o) + \mathcal{L} ( I \given \mathcal{S}, o )
\end{align*}

Using the Shannon-Fano code we have the following.
{
\begin{align*}
    \mathcal{L}(o) &= - \log (Pr(o) ) \\
    \mathcal{L}(\mathcal{S} \given o) &= -\log ( Pr(\mathcal{S} \given o) ) \\
    \mathcal{L}(I \given \mathcal{S}, o) &= -\log ( Pr(I \given \mathcal{S}, o) )
\end{align*}

}
Combining and expanding these, we obtain the following.
\vspace{0.05in}
{
\begin{align*}
&\mathcal{L}(o, \mathcal{S}, I) \\
&=- \log (p(o)) -\log \big(Pr(\mathcal{S} \given o)\big) - \log \left( \frac{Pr(\mathcal{S} \given I, o) Pr (I \given o)}{Pr (\mathcal{S} \given o)} \right) \\
&=   - \log (p(o)) - \log \big(  Pr(\mathcal{S} \given I, o) Pr (I \given o) \big)
\end{align*}
}

\noindent
where $Pr(o) = p(o)$.
Note from \eqref{eq:causal} that given $I$, $\mathcal{S}$  does not depend longer on $o$. Hence we have, $Pr( \mathcal{S} \given I, o) =  Pr(\mathcal{S} \given I)$.
Clearly, $Pr(\mathcal{S} \given I)$ is 1 if after removing $I$ from $G$, the serviced set is exactly $\mathcal{S}$, and 0 otherwise.
For $Pr(I\given o)$, recall that we assume that edge failures depend only the disaster's effect and are independent from each other, therefore, we can represent $Pr(I\given o)$ as a product of failure probabilities of the individual failure. Note that $Pr(\mathcal{S} \given I)$  is the probability of set of demand nodes $S$  being  connected to the source the given failure set is $I$. As this probability can only take value of $0$ or $1$ (indicating whether it is feasible for $S$ to be connected or not), we represent it with a indicator function $h(G, \mathcal{S}, I)$, which we refer to as a feasibility function.
Combining all of the above, the final model cost is:

{
\begin{align}
&\mathcal{L}(o, \mathcal{S}, I) = \nonumber \\
&- \log \big(p(o)\big)  -\log \big( h(G, \mathcal{S}, I)  \prod_{e\in I} F( e | o) \prod_{e\in E\setminus I} (1 - F( e | o) ) \big).
\label{eq:final_model_cost_cq}
\end{align}
}

\subsection{Data Cost}
Note that MDL needs a data cost as well (otherwise we will choose a model without considering $\mathcal{Q}_I$ and $\mathcal{Q}_c$).
In the MDL framework, given the model is already transmitted, we now send the data.
For the data cost, we send first the size of $\mathcal{Q}_c$ and $\mathcal{Q}_I$, which are $\norm{\mathcal{Q}_c}$ and $\norm{\mathcal{Q}_I}$, respectively. They followed by the joint queries, i.e. our two probes sets $\mathcal{Q}_c$ and $\mathcal{Q}_I$.
Hence, our data cost consists of four terms.

As $\mathcal{Q}_c$ is sampled uniformly from $\mathcal{S}$ with probability $\gamma_c$, the size $\norm{\mathcal{Q}_c}$ is send as follows:
\vspace{0.08in}
\setlength\abovedisplayskip{5pt}
\setlength\belowdisplayshortskip{5pt}
{
\begin{align*}
\mathcal{L} ( \norm{\mathcal{Q}_c} | \mathcal{S}, I, o )  &= - \log ( Pr(\norm{\mathcal{Q}_c} | \mathcal{S}, I, o)) \nonumber \\
&=  -\log ( \binom{\norm{\mathcal{S}}}{\norm{\mathcal{Q}_c}} \gamma_c^{\norm{\mathcal{Q}_c}} (1 - \gamma_c)^{\norm{\mathcal{S} \setminus \mathcal{Q}_c}})
\end{align*}
}
Similarly, for $\norm{\mathcal{Q}_I}$ we have:
\setlength\abovedisplayskip{5pt}
\setlength\belowdisplayshortskip{5pt}
{
\begin{align*}
\mathcal{L} ( \norm{\mathcal{Q}_I} | \mathcal{S}, I, o )  &= - \log ( Pr(\norm{\mathcal{Q}_I} | \mathcal{S}, I, o)) \nonumber \\
&=  -\log ( \binom{\norm{I}}{\norm{\mathcal{Q}_I}} \gamma_I^{\norm{\mathcal{Q}_I}} (1 - \gamma_I)^{\norm{I \setminus \mathcal{Q}_I}})
\end{align*}
}

Then the  cost of  sending $\mathcal{Q}_c$ and $\mathcal{Q}_I$ in terms of previously sent information is:
\vspace{0.08in}
\setlength\abovedisplayskip{5pt}
\setlength\belowdisplayshortskip{5pt}
{
\begin{align*}
&\mathcal{L}(\mathcal{Q}_c \given[\big] \mathcal{S}, I, \norm{\mathcal{Q}_c}) +
\mathcal{L}(\mathcal{Q}_I \given[\big] \mathcal{S}, I, \norm{\mathcal{Q}_I}) \nonumber \\
= &- \log (Pr(\mathcal{Q}_c \given[\big] \mathcal{S}, I, \norm{\mathcal{Q}_c}) )
- \log (Pr(\mathcal{Q}_I \given[\big] \mathcal{S}, I, \norm{\mathcal{Q}_I}) ) \nonumber \\
= & - \log ( \gamma_c^{\norm{\mathcal{Q}_c}} (1-\gamma_c)^{\norm{\mathcal{S} \setminus \mathcal{Q}_c}} )
- \log ( \gamma_I^{\norm{\mathcal{Q}_I}} (1-\gamma_I)^{\norm{I \setminus \mathcal{Q}_I}} )
\end{align*}
}

Combining the four data cost terms, the total data cost is:
{
\begin{align}
\label{eq:datacost_joint}
= &- \log \binom{\norm{I}}{\norm{\mathcal{Q}_I}} - 2\norm{\mathcal{Q}_I} \log (\gamma_I) - 2(\norm{I} - \norm{\mathcal{Q}_I}) \log(1 - \gamma_I) \nonumber \\
& -\log \binom{\norm{\mathcal{S}}}{\norm{\mathcal{Q}_c}} - 2\norm{\mathcal{Q}_c} \log (\gamma_c) - 2(\norm{\mathcal{S}} - \norm{\mathcal{Q}_c}) \log(1 - \gamma_c).
\end{align}
}

\noindent\textbf{Formal Problem Statement.}
We can now state our problem formally as:
\begin{problem}
\label{prob:problem}
\prob{}: \textit{Given an undirected graph $G(V ,E)$ where $V=V_S \cup V_D \cup V_T$,
a many-to-many relation mapping supply nodes $V_S$ to demand nodes $V_D$,
a set of observed serviced nodes $\mathcal{Q}_c$ sampled uniformly at random with probability $\gamma_c$ from $\mathcal{S}$,
and a set of observed failed edges $\mathcal{Q}_I$ sampled uniformly at random with probability $\gamma_I$ from $I$,
find the complete set of serviced nodes $\mathcal{S}^* \subseteq V_D$ and failed edges $I^* \subseteq E$ (that failed as per the failure model described in Section~\ref{sec:prelim}),
by minimizing the MDL cost function,
i.e.}
{\small
\begin{align*}
&\langle \mathcal{S}^*, I^* \rangle \nonumber\\
&= \underset{\mathcal{S}, I}{\arg\min}  \Big\{ \nonumber\\
&- \log \big(p(o)\big)  -\log \big( h(G, \mathcal{S}, I)  \prod_{e\in I} F( e | o) \prod_{e\in E\setminus I} (1 - F( e | o) ) \big) \nonumber\\
&- \log \binom{\norm{I}}{\norm{\mathcal{Q}_I}} - 2\norm{\mathcal{Q}_I} \log (\gamma_I) - 2(\norm{I} - \norm{\mathcal{Q}_I}) \log(1 - \gamma_I) \nonumber \\
& -\log \binom{\norm{\mathcal{S}}}{\norm{\mathcal{Q}_c}} - 2\norm{\mathcal{Q}_c} \log (\gamma_c) - 2(\norm{\mathcal{S}} - \norm{\mathcal{Q}_c}) \log(1 - \gamma_c)
\Big\},
 \nonumber
\end{align*}
}
i.e.,
{\small
\begin{align}
\label{eqn:objective}
&\langle \mathcal{S}^*, I^* \rangle =  \underset{\mathcal{S}, I}{\arg\min}  \Big\{ \mathrm{Eq.~(\ref{eq:final_model_cost_cq})} + \mathrm{Eq.~(\ref{eq:datacost_joint})} \Big\}
\end{align}
}
\end{problem}

\vspace{6pt}

\noindent\textbf{Algorithm for \prob{} Problem.}
Solving Problem \ref{prob:problem} is very challenging. First of all, the search space of the solution is exponentially large as it involves searching over all possible sets of serviced and failed nodes, and handle the combinatorial explosion of paths.
This makes the problem naturally more challenging than past work in point queries \cite{adhikari2018near}.
Moreover, since the objective is a function of two interdependent sets, designing any algorithm with performance guarantee is of a major challenge.
Formally, we state the following lemma.

\begin{lemma}
\label{lemma:hardness}
\prob{} is NP-hard to approximate within an $\Omega(\sqrt{\log\log{|V|}})$ factor.
\end{lemma}
\vspace{5pt}
\textit{Proof.} (Sketch)
Our proof is by a reduction from the multicut problem, which is defined in the following manner: we are given an
instance $H=(V_H, E_H)$, and a set of source-sink pairs $(s_1, t_1),\ldots, (s_k, t_k)$. The goal is to pick
the smallest subset $E'\subseteq E$ of edges, such that all the source-sink pairs are disconnected in $G[E\setminus E']$,
which was shown by \cite{Chawla2006}
to be NP-hard to approximate within an $\Omega(\sqrt{\log\log{|V_H|}})$ factor,
under standard complexity theoretic assumptions.
We reduce this to an instance of \prob{} in the following manner.
We add a new supply node $s$, and add $m'$ parallel edges $(s, s_i)$ for each $i$ (this can be changed to non-parallel
edges, by adding a new node on each such edge). Finally, we add an additional demand node $t_{k+1}$, and add
the path from $s$ to $t_{k+1}$, which gives us the graph $G$.
We set $V_D=\{t_1,\ldots,t_{k+1}\}$, and $\mathcal{S}=\mathcal{Q}_c=\{t_{k+1}\}$, with $\gamma_c=1$.
Let $M$ be number of edges  in $G$. We set $F(e|o)=1/2^M$ for all edges $e$.
Observe that any feasible solution must make all the terminals $\{t_1,\ldots,t_k\}$ disconnected from $s$.
Next, the cost of any feasible solution $I\subset E$ is $\leq M|I|+1$, because of the choice of $F(\cdot)$.
Finally, due to the parallel edges $(s, s_i)$, none of those edges will be picked in the solution.
This implies, the solution $I$ picked for this instance is also a multicut in $H$.
Since the cost of the solution remains the same, within a fixed multiplicative factor, the lemma follows.
$\blacksquare$
{
\begin{algorithm}[t!]{}
\caption{\textsc{JointPathMap}\xspace{}}
\begin{algorithmic}[1]
\small
    \STATE \textbf{Input:} Instance $(G, V, \mathcal{Q}_c, F, \gamma_c, \mathcal{Q}_I, \gamma_I)$
    \STATE \textbf{Output:} Solution $\hat{I}$ as per Eq. \eqref{eqn:objective}
    \FOR{$o\in\mathcal{O}$}
    \STATE $\hat{I}(o) \gets \mathcal{Q}_I$ \quad \{initialize solution set for $o$\}
  \WHILE{$\exists e \in E \setminus \hat{I}(o)$ that decreases the MDL cost}
    \STATE \{To calculate MDL, first compute $\mathcal{S}$ with BFS in $G \setminus \hat{I}(o)$\}
    \STATE Select $e \in E \setminus \hat{I}(o)$ that decreases the most the MDL cost
        \STATE $\hat{I}(o) \gets \hat{I}(o) \cup \{ e \}$ \quad \{update solution\}
        \STATE $\alpha(o) \gets $ MDL cost for updated $\hat{I}(o)$
    \ENDWHILE
    \ENDFOR
\STATE Return $\hat{I}(o)$ for $o$ that minimizes $\alpha(o)$
\end{algorithmic}
\label{alg:localsearch}
\end{algorithm}
}
\par \noindent
\textbf{Note:} It is worth observing that the corresponding point query only formulation of \cite{adhikari2018near} can be approximated
within an additive $O(\log{n})$ factor. Additive approximations are, in general, easier than multiplicative
approximation guarantees. Therefore, Lemma \ref{lemma:hardness} implies Problem \ref{prob:problem} is a
harder problem than the point query version, from a complexity standpoint.

\prob{} is quite challenging as it is related to the multicut problem with multiple paths between source-sink pairs. Since a related multicut problem leverages a greedy strategy \cite{kortsarts2005greedy}, we adopt this approach. We describe our algorithm \textsc{JointPathMap}\xspace{} for \prob{} next.
It iterates over the disaster scenarios in $\mathcal{O}$ to find a failure set for each $o$.
To obtain $\hat{I}(o)$, it starts by initializing it to $\mathcal{Q}_I$. Then, at each iteration, an edge $e \in E\setminus \hat{I}(o)$, whose addition decreases the MDL cost the most, is added to the failure set $\hat{I}(o)$.
The process is repeated until no such edge can be found. At the end, the $\hat{I}(o)$ that has the lowest MDL score is returned.
The implementation needs the calculation of $\mathcal{S}$ after $\hat{I}(o) \cup \{e\}$ is removed from $G$, which needs to check which demand nodes has at least one path to a supply node,
which is accomplished via breadth-first search (BFS).
The complete pseudocode is in Algorithm~\ref{alg:localsearch}.
\newpage
\begin{lemma}
Algorithm \textsc{JointPathMap}\xspace{} runs in $O(|V|^2 |E|+|V||E|^2)$ time.
\end{lemma}
\textit{Proof.} Calculating total MDL cost for any solution instance $\hat{I}$ requires a BFS (to calculate the $\hat{\mathcal{S}}$ induced by $\hat{I}$) for each demand node, which, takes $O(|V|^2+|V||E|)$. The inner loop in lines 5-10 takes $O(|E|)$ and the loop in 3-11 takes $O(|\mathcal{O}|)$. In practice, $|\mathcal{O}| \lll |V|$ is a small constant; therefore, our algorithm runs in $O(|V|^2 |E|+|V||E|^2)$.
$\blacksquare$

Note that the feasibility function contained in the MDL cost will ensure that the infeasible solutions are not selected. This is because the term $h(G, \mathcal{S}, I) = 0$ for infeasible solutions. Hence, the MDL cost for such solutions will be infinite and they will not be picked by our heuristic. Hence, the algorithm always ensures that $\mathcal{Q}_c \subseteq \hat{\mathcal{S}}$.

\section{Experiments}
\label{sec:expts}
\noindent
In this section, we empirically test our formulations and algorithms in both single and multiple disaster scenarios.

\subsection{Setup}
\par \noindent
\textbf{Setup.}
We implemented the algorithms in Python\footnote{Code and data: \href{https://github.com/arodriguezca/JointPathMap}{https://github.com/arodriguezca/JointPathMap}}\!.
We ran our algorithms on a Linux machine with 8 cores with 3.40 GHz and 16GB of RAM.
Presented results are the average over 30 different \emph{damage scenarios} simulated as per civil engineering standards (explained in detail in our supplement). Each damage scenario is originated from a single \emph{scenario earthquake} $o\in\mathcal{O}$ (each $o$ is a magnitude and an epicenter) which was randomly selected with probability $p(o)$.
\textsc{JointPathMap}\xspace takes roughly 2 seconds to output a solution (for a single damage scenario) for EPN and PWN networks, and roughly 500 seconds for our large network GridKit (see Section \ref{subsec:dataset} for network specifications). Therefore, our algorithm is fast and scalable.

\subsection{Datasets}
\label{subsec:dataset}
\par \noindent
We evaluate the performance of our algorithms on the topologies of three real infrastructure networks. For each network, we follow procedures as described in civil engineering literature~\cite{adachi_comparative_2010} with structural parameters from~\cite{fema_multi-hazard_2013} to generate multiple geographical correlated failure probabilities for nodes. Edge failure probabilities are calculated from these counterparts.

\begin{figure}[h]
    \centering
    \includegraphics[width=0.44\textwidth]{images/shelby_networks_zoom.pdf}
    \caption{Skeletonized power and water network topologies in Shelby County, Tennessee.}
    \label{fig:shelby}
\end{figure}

\par \noindent
\textbf{EPN \& PWN.}
We use two utility networks for Shelby County, TN, USA: electric power network (EPN) and potable water network (PWN), with $\norm{V_{EPN}}=59$, $\norm{E_{EPN}}=73$, $\norm{V_{PWN}}=49$, $\norm{E_{PWN}}=71$, whose topologies are depicted in Figure \ref{fig:shelby}
and are publicly available in~\cite{zhang_probabilistic_2018-1}.
We adopt the 10 scenario earthquakes that characterize the seismic hazard in the area identified in \cite{adachi_comparative_2010}. This means that once an earthquake happens, its effects on structures can be characterized by one of these 10 scenario earthquakes.
Each of the scenario earthquake induce different failure probability-sets ($F$'s) for the network.

\par \noindent
\textbf{GridKit.}
We use the high-voltage power grid for Northern California, obtained from \cite{wiegmans2016gridkit}, from which we extract the subgraph corresponding to Northern California.
It contains $\norm{V_{GK}}=415$ and $\norm{E_{GK}}=543$. We generate 50 scenario earthquakes, each of them as follows: epicenter and magnitude are selected randomly, the first from all node locations,
and the second from the four values considered in \cite{adachi_comparative_2010}.

Table~\ref{tab:network_stats} contains the specific number of demand, source and transshipment nodes for these three networks. As mentioned before, \textsc{JointPathMap}\xspace is fast in these real networks publicly available; therefore, when larger networks become available, our algorithm will apt to run in a suitable time.

\begin{table}[h]
\centering
\caption{Datasets details}
\begin{tabular}{ m{1cm} | m{1cm} | m{1cm} | m{1cm} | m{1cm} | m{1cm} }
Network & $\norm{V}$ & $\norm{V_S}$ & $\norm{V_D}$ & $\norm{V_T}$ & $\norm{E}$\\ \hline
EPN & 59 & 9 & 37 & 13 &  73\\
PWN & 49 & 15 & 34 & 0 & 71 \\
GridKit & 415 & 114 & 82 & 219 & 543
\end{tabular}
\label{tab:network_stats}
\end{table}

We consider the following factual observation in our experiments.
Even though the network topology allows it, in practice, source nodes cannot service \emph{any} demand node, but only the ones close to them that its generative capacity allows. To model this physical constraint, we restrict serviceability of a demand node to the case when there is at least one path of length $L$ to any source. In our experiments, we use $L_{EPN}=3$, $L_{PWN}=2$, and $L_{GK}=4$.

\subsection{Baselines and Metrics}
We use two baselines, namely \modelcost and \pathmap. \modelcost is  an algorithm that greedily minimizes model cost in Eq.~\eqref{eq:final_model_cost_cq}, to demonstrate the importance of considering the data cost. \pathmap is a baseline that only includes connectivity queries.

To evaluate our algorithms, we use
F1-score, for which we need to compute precision and recall. Following the notation in our paper, let $I$ be the true failure set and $\hat{I}$ be our solution set. Our metrics are defined in literature as follows: $\text{precision} = \norm{\hat{I} \cap I } / \norm{\hat{I}}$,
$\text{recall} = \norm{\hat{I}  \cap I } / \norm{I}$, and  $\text{F1-score} = \text{precision} \times \text{recall} / (\text{precision} + \text{recall})$. We add the following limit cases to handle empty $I$ or $\hat{I}$: precision is 1 when $\hat{I}=I=\emptyset$ (we correctly identified no failures), and recall is 1 when $\norm{I}=0$ (there are no failures missed in our solution $\hat{I}$).

\subsection{Results}
\prob{} problem asks to optimize the MDL cost to infer the complete failure set, and \textsc{JointPathMap}\xspace{} succeeds at this optimization task at low computation cost.
This is reflected in our results in Figure~\ref{fig:mdl_u_edges}. We can see that the ratio between optimal MDL and our methods is close to 1.
The optimal MDL is calculated via integer linear programming where all paths are listed.
\begin{figure}[t]
\centering
\includegraphics[width=0.33\textwidth]{images/joint_AAAI_ss_mdl.pdf}
\linespread{1}
\caption{MDL ratio between optimal and our method \textsc{JointPathMap}\xspace. }
\label{fig:mdl_u_edges}
\vspace{-0.1in}
\end{figure}

We present the performance of all the methods in Figure \ref{fig:results} in terms of F-1 score.
As described before, this score is an average over 30 damage scenarios, each of them coming from a different disaster scenarios $o$ sampled from $\mathcal{O}$ (our scenario earthquakes) with probabilities $p(o)$. As failure probabilities $F$ depend on $o$, this set of damage scenarios leads to a variety of ground truth failure sets $I$.
As we can observe, \textsc{JointPathMap}\xspace outperforms all the baselines significantly in all three real world networks for all values of $\gamma_c$. Poor performance of \modelcost can be attributed to the fact that it infers the most likely failures while ignoring the observed data.
Thus, \modelcost's result highlights the necessity of the MDL framework.

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{images/ms_gamma_i03_F1_horiz.pdf}
    \caption{F1-score for  \textsc{JointPathMap}\xspace ($\gamma_I=0.3$), \textsc{OnlyConnectivity}\xspace and \modelcost for multiple disaster scenarios. Datasets from left to right: EPN, PWN, and GridKit.}
    \label{fig:results}
\end{figure*}

Similarly, \pathmap has low F1-score due to low recall, which indicates that while the inferred failures are mostly correct, they are just a fraction of the total failure set. This is due to the fact that, depending on the redundancy of the network topology, some of the failures are not reflected in the connectivity queries, for which it would be very hard to infer by only using $\mathcal{Q}_c$. To get a sense of how many of 'undiscoverable' edges (u-edges) are present in ground truth failure set $I$, we add as many edges from $I$ to $\hat{I}$ as long serviceability is not affected (i.e. $\hat{\mathcal{S}}$ does not change). Then, we calculate the proportion of u-edges in $I$.
In Figure~\ref{fig:u_edges}, we can see that there is a high proportion of u-edges, i.e. edges not reflected in connectivity queries. This indicates there are many failure sets that can produce the same serviced set $\mathcal{S}$; consequently, solely using connectivity probes $\mathcal{Q}_c$ is clearly not sufficient.
As shown by the results, \textsc{JointPathMap}\xspace does not suffer from this issue and works well in real datasets. The improved performance is because \textsc{JointPathMap}\xspace incorporates a small set of probes of the failures $\mathcal{Q}_I$, which provides insight about the location of the epicenter thanks to the geographical correlation among failures.

\begin{figure}
\centering
\includegraphics[width=0.33\textwidth]{images/conn_AAAI_ss_u-edges.pdf}
\caption{
u-edges proportion is an estimate of the proportion of existing edges in the actual failure set $I$ that are not reflected in the connectivity queries. This proportion is especially high for the PWN, which is at the same time the hardest to solve for \pathmap.
}
\label{fig:u_edges}
\end{figure}

For multiple scenario earthquakes, the algorithm has to infer first the culprit scenario $o$ and then the complete failure set $I$.
In a small network, some of scenario earthquakes are similar and missing the right one may not drastically affect further inference; however, in large networks, scenarios earthquakes are completely different (mainly in epicenter), therefore, without a proper inference of $o$ it is very hard to make a good inference of $I$.
Adding a small sample of failure probes helps to properly identify $o$ and the right expected number of failures, and gives a larger improvement in F1-score.
In addition, we can see in Figure \ref{fig:mdl_u_edges_joint} the effect of different sample sizes $\gamma_I$, where each increment has a different effect, being the most relevant the jump from $\gamma_I=0.3$ to $\gamma_I=0.4$.

\begin{figure}
\centering
\includegraphics[width=0.33\textwidth]{images/joint_AAAI_electric_ss_recall_across_gamma_i.pdf}
\caption{The effect of increasing $\gamma_I$ (i.e. increase the sample of failure probes) for EPN dataset.
We present recall here because we want to emphasize the ability of our algorithm to discover elements of the failure set $I$ that are not in $\mathcal{Q}_I$.}
\label{fig:mdl_u_edges_joint}
\end{figure}

As additional results, we also present results for a single disaster scenario.
In practice, after an earthquake occurs, an estimate of its magnitude and epicenter is released within minutes, but this information is then corrected several times over the next hours. Therefore, for this kind of natural disaster, it is preferrable to handle uncertainty by considering multiple earthquake scenarios as we did in our formulation. However, once the estimation about magnitude/epicenter is reliable, we can add this extra information to our problem by reducing $\mathcal{O}$ to a single $o$ and making $p(o)=1$.
For our single disaster scenario experiments, we pick the $o$ that has the largest expected number of failures for EPN and PWN; and, for GridKit, as most disaster scenarios lead to similar damage, we select one at random.
As expected, when know the disaster scenario $o$, this extra information helps the inference as we can see in Figure \ref{fig:results_single}. As we consider the scenario which has the most failures for EPN and PWN (the hardest one), the performance of \textsc{JointPathMap}\xspace may be slightly lower than when aggregating the performance over all scenarios.

\begin{figure*}
    \centering
    \includegraphics[width=0.9\textwidth]{images/ss_gamma_i03_F1.pdf}
    \caption{F1-score for \textsc{JointPathMap}\xspace ($\gamma_I=0.3$), \textsc{OnlyConnectivity}\xspace and \modelcost for a single disaster scenario, i.e. $o$ is known. Datasets from left to right: EPN, PWN, and GridKit. \modelcost returns empty solutions ($\hat{I}=\emptyset$) for PWN because all F's are below 0.5 (i.e. are more likely to remain functional), but the actual expected failure set size is $\approx 10$. Therefore, as its recall is 0, its F1-score is constantly 0.}
    \label{fig:results_single}
\end{figure*}

\begin{figure}
\centering
\resizebox{0.48\textwidth}{!}{\input{images/net_visualization.pgf}}
\caption{Visualization of $I$, $\mathcal{S}$ and $\hat{I}$ in PWN for $\gamma_c=\gamma_I=0.3$.
Source nodes are depicted as squares and demand nodes as circles (there are no transshipment nodes).
Probes $\mathcal{Q}_c$ and $\mathcal{Q}_I$ are represented as green nodes and edges respectively.
Ground truth non-serviced demand nodes have a red X.
Ground truth failed edges are represented as dotted lines.
Blue edges are predicted failures by \textsc{JointPathMap}\xspace.
Red edges are discoverable failures missed by \textsc{JointPathMap}\xspace.
Purple edges are u-edges.
}
\label{fig:net_vis}
\end{figure}

\subsection{Case Study}
Within the single scenario earthquake framework, we study how the failures discovered relate to the ground truth failure set in PWN.
In Figure \ref{fig:net_vis}, we observe that our algorithm misses four failures, from which three of them are actually undiscoverable (u-edges). In other words, these failures do not impede any demand node (circle) to be serviced by source nodes (squares) as it can be recognized from the figure. Serviceability in this network is restricted to available paths of length 2 within source and demand node. Consequently, in this network, there are five demand nodes non-serviced (marked with a red X), for which our algorithms found the right failures to disconnect for four of them. In fact, the two failed edges found (dashed blue edges) are the reason why three demand nodes are non-serviced. Therefore, repairing them is of utmost importance to recover service in most of the non-serviced demand nodes.
This exemplifies how \textsc{JointPathMap}\xspace finds the truly important informative edges to restore serviceability.

\section{Related Work}
\noindent
\par \noindent \textbf{Critical Infrastructure Networks:} Data mining critical infrastructure has been gaining much research interest in recent times~\cite{tabassum2018data}. These include works in modeling~\cite{ouyang2014review}, dependency inference~\cite{chen2016fascinate}, and finding hotspots~\cite{chen2017hotspots}.
Chen et al.~\cite{chen2016fascinate} proposed collective collaborative filtering based approach for inferring dependency between various interconnected critical infrastructure networks. Similarly Chen et al.~\cite{chen2017hotspots} proposed cascade based failure model in heterogeneous critical infrastructure networks to identify most critical nodes in the system.
A more closely related line of work studies failure inference on infrastructure networks. \cite{anwar2015roadrank} proposed an algorithm to infer congested intersections given the road networks and usage data. Similarly \cite{adhikari2018near} proposed a near-optimal MDL-based approach to infer failures, but given point probes.

\noindent\textbf{Network Tomography:} This involves inferring the network states based on end-to-end network measurements. These measurements, available in the Internet and other communication networks,
measure delay across signals transmitted by the probes, which is information that can be used to infer delays on links~\cite{xia2006inference} and network topology~\cite{jin2006network,anandkumar2013topology}.
Xia and Tse~\cite{xia2006inference} propose an expectation maximization algorithm to infer delays on links from performing a great quantity of end-to-end measurements.
Jin et al.~\cite{jin2006network} proposed a heuristic to infer approximate topology based on end-to-end measurements with the goal of reducing the number of such measurements while maintaining accuracy.
Anandkumar et al.~\cite{anandkumar2013topology} study inferring the topology of sparse networks.
They evaluate their algorithm with two routing models and find that they only require a sub‐linear number of participant nodes selected uniformly at random.

\noindent\textbf{Epidemic Inference:}
Another closely related field is reverse engineering an epidemic outbreak over a network given partially observed infections. Shah et al proposed a MLE based approach to detect the source of computer virus spreading over a computer network~\cite{shah2010detecting}. Similarly, \cite{prakash2012spotting,sundareisan2015hidden} leverage MDL based approach to detect both the source and missing infections in networks. These works assume that the infection/virus propagation over the network takes a stochastic form often modelled using a variant of the popular SIR model. On the other hand, a different set of works \cite{rozenshtein2016reconstructing, xiao2018reconstructing, xiao2018robust} model the interaction among individuals as a temporal or static time expanded network and treat the combination of the observed infections and the time as a node in the network. The inference is typically done by solving temporally constrained Steiner-tree problem where the terminals are the observed infections.

\section{Conclusions}
\label{sec:conc}

\noindent We investigated the problem of inferring network states given connectivity queries, a piece of information more practical to obtain than a direct sample of network states (point queries).
Our empirical observations indicate that redundancies, common in infrastructure networks, cause a significant proportion of states to be undiscoverable by connectivity probes. Therefore, we propose adding to the connectivity queries a small sample of point queries to exploit relationships among network states, which proved to be valuable for boosting the power of connectivity queries. We formulated the novel problem with MDL and developed a scalable greedy algorithm that consistently optimizes MDL cost close to optimal in extensive experiments. These experiments are performed in real networks and principled simulations.
Future work can focus on extending to node failures and using noisy probe data. Other interesting future work is extending our formulation to other applications of connectivity probes such as LAN networks and industrial communication networks.

\section*{Acknowledgements}
This paper is based on work partially supported by the NSF (Expeditions CCF-1918770, CAREER IIS-2028586, RAPID IIS-2027862, Medium IIS-1955883, NRT DGE-1545362), CDC MInD program, ORNL, funds/computing resources from Georgia Tech, and the National Institute of Standards and Technology through the Center for Risk-Based Community Resilience Planning [70NANB15H044, 70NANB20H008].
A.V.'s work was partially supported by NSF grants CRISP 2.0 Grant 1832587, CMMI-1745207, IIS-1633028 and DTRA subcontract/ARA S-D00189-15-TO-01-UVA.

\newpage

\bibliographystyle{IEEEtran}

\providecommand{\noopsort}[1]{}


\section*{Plain Language Summary}
We aim to understand the size effect in granular column collapses. Geophysical hazards, such as landslides, debris flows, pyroclastic flows, and rock avalanches, occur frequently worldwide, and it is difficult to predict their initiation and dynamical behavior. Many geophysical flows can be considered as granular materials, where interaction between particles matters. However, it is difficult to upscale a grain-scale particle interaction to a large-scale predicting tool to capture the behavior of many complex granular flows. This paper starts from the collapse of granular columns, a fundamental problem often used to investigate the mobility of granular materials, and carefully navigates to examine the run-out behaviors of columns with different sizes, so that we can obtain a better picture of how different system sizes lead to different collapsing scenarios. We further analyzed the characteristic length scale associated with the collapse of granular columns, and the strong force network presenting at the beginning of the collapse, to understand better what drives the size effect of granular systems, and what we can gain from this research to predict the behavior of geophysical flows.

\section{Introduction}
Granular materials are omnipresent in natural and engineering systems, and the physics and mechanics of them are crucial for understanding geophysical flows, natural hazards, food processing, chemical engineering, and pharmaceutical engineering \cite{Lube2019GenerationOA,guyon2020built}. Granular materials can behave like a solid, a liquid, or a gas in different circumstances \cite{midi2004dense}, which increases the difficulty in capturing their macroscopic behavior. Besides, collective structures may form inside a granular system, and the existence and the size of such collective structures [such as bridging \cite{mehta2007granular}, granular agglomerates \cite{Vo2020evolution}, and contact networks \cite{Zhang2014force}] will introduce a strong size effect, which further increases the complexity of the problem. In recent decades, breakthroughs have been made to understand the basic governing principles, especially the constitutive relationships, of granular materials \cite{midi2004dense,trulsson2012transition,Vo2020additive,man2021granular}, where the behavior of granular materials or granular-fluid mixtures are considered to be described by dimensionless numbers expressed as the ratio between dominating stresses. Even though geophysical flows are rarely dry, investigating dry granular materials is still considered as an important start for more complex systems.

Due to the similarity and potential links between the collapse of granular columns and gravity-driven geophysical flows, previous research investigated the collapse of granular columns to analyze the post-failure behavior of granular systems \cite{thompson2007granular,lacaze2009axisymmetric}. \citeA{lube2004axisymmetric} and \citeA{lajeunesse2005granular} independently determined relationships for both the normalized run-out distance $\mathcal{R}=(R_{\infty}-R_i)/R_i$ (where $R_{\infty}$ is the final radius of the granular pile, and $R_i$ the initial radius of the granular column), and the halt time of a collapsed granular column, which both scale with the initial aspect ratio, $\alpha = H_i/R_i$ (where $H_i$ is the initial height of the column), a parameter drawn out of dimensional analysis. \citeA{zenit2005computer,staron2005study,staron2007spreading,lacaze2009axisymmetric} studied either 2D or axisymmetric granular column collapses using the discrete element method (DEM), and explored the influence of inter-particle properties, energy dissipation, and inertial effects on the run-out and deposition of granular columns. \citeA{Farin2019RelationsBT} even linked the characteristics of granular column collapses to seismic signals. \citeA{warnett2014scalings} and  \citeA{Cabrera2019Granular} studied the collapse of granular columns with experiments and simulations, and argued that the relative size of a granular column, $R_i/d$, where $d$ is the average particle diameter, has strong influence on the run-out distance of a collapsed granular column, and to avoid significant size effects, $R_i/d$ must be larger than 75 for short columns and larger than 50 for tall columns.

Based on dimensional analysis and simulations, we investigated the collapse of axisymmetric granular columns and their resulting deposition with a wide range of inter-granular and particle/boundary frictional coefficients \cite{man2020universality}, and concluded that the normalized run-out distance, $\mathcal{R}$, scales with an effective aspect ratio,
\begin{linenomath*}
\begin{equation}
\begin{split}
    \alpha_{\textrm{eff}} = \left(\mu_w+\beta\mu_p\right)^{-1/2}\left(H_i/R_i\right)\ ,
\end{split}
\end{equation}
\end{linenomath*}
where $\mu_w$ is the frictional coefficient between particles and the boundary, $\mu_p$ is the inter-particle frictional coefficient, and $\beta$ can be seen as the ratio between contributions from inter-particle frictions and those from particle/boundary frictions. They found that the collapse of granular columns can be classified into three different regimes: quasi-static, inertial, and liquid-like. The effective aspect ratio, which was used to describe the run-out behavior, was obtained from dimensional analysis and was shown to represent the ratio between the inertial stress and the frictional stress or, equivalently, the ratio between the kinetic energy and dissipated energy during the collapse of granular columns. In our study, as we change the frictional coefficient, the transition point, where the slope of $\mathcal{R}(\alpha)$ curve changes on a log-log coordinate system, also changes accordingly. Moreover, after using $\alpha_{\textrm{eff}}$ to re-scale the data, this transition point in $\mathcal{R}-\alpha_{\textrm{eff}}$ space becomes universal.

In this paper, after introducing the simulation set-up and the contact law in the DEM, we further conduct a systematic study of granular column collapses with the sphero-polyhedral DEM to investigate the influence of relative column size on the run-out distance to obtain a scaling solution to describe granular column collapses with different sizes. This can not only further explore the physics of granular column collapses but also open a window for linking the behavior of granular columns to the mobility of geophysical flows. Inspired by the work of \citeA{Roche2002ExperimentsOD,Roche2008ExperimentalOO,bougouin2021}, where they presented physical connections between granular dam breaks and pyroclastic flows, to strengthen the tie between granular column collapses and geophysical flows, we present an application of our results to various rock falls around the Soufriere Hills, Montserrat as a result of volcanic eruptions and the resulting pyroclastic flows and granular surges in Section \ref{sec:application} and Figure \ref{discussion}. Finally, we present our concluding remarks.

\section{Simulations}\label{sec:simu}

\subsection{Discrete element method and simulation setup}

We perform DEM simulations of the granular column collapses with Voronoi-based sphero-polyhedra particles\cite{pournin2005generalization,galindo2010molecular}. The sphero-polyhedra technique allows for an easy and efficient definition of contact laws between particles. This is due to the smoothing of the edges of all geometric features by circles (in 2D) or spheres (in 3D). We note that, in our simulation, three types of contacts (vertex-vertex contact, edge-edge contact, and vertex-face contact) are considered. For these types of contacts, we implement a Hookean contact model with energy dissipation to calculate the interactions between particles. The DEM model used in this work has been thoroughly validated with experiments as shown in the supplemental material, where the shape of particles was accurately captured and the frictional coefficients of both inter-particle contact and particle/boundary contact were independently measured.

\begin{figure}[!h]
  \centering
  \includegraphics[scale = 0.30]{figure1.png}
  \caption{Simulation set-up. (a) the initial state of the granular column; (b) the final deposition of the collapsed column in a 3-D view; (c) shows the method for measuring the run-out distance of a collapsed granular column. The $x-$axis is the radial position $r$, and the $y-$axis is the percentage of number of particles, $P_{\textrm{dm}}$, located within $(r-\Delta r/2, r+\Delta r/2)$ divided by the radial position. We determine the final radius of a collapsed granular column when $P_{\textrm{dm}}(r) \leq 5\%P_{\textrm{dm}}^{\textrm{max}}$; (d) Initial state and (e) final state of granular column with $R_i = 5.7$ cm and $H_i = 13$ cm; (f) Initial state and (g) final state of granular column with $R_i = 2.3$ cm and $H_i = 10$ cm.}
  \label{setup}
\end{figure}

In the simulation, we create particles in a certain cylindrical domain with initial height, $H_i$, and initial radius, $R_i$ [Fig. \ref{setup}(a)]. 20\% of the particles are randomly chosen and removed from the simulation to create a granular packing with an initial solid fraction of $\phi_s = $0.8. The average particle size $d$ is 0.2 cm. The coefficient of restitution of particle collisions is 0.1. The material properties are set to be the same as that of quartz sand (particle density 2.65 g/cm$^3$). Using 3D Voronoi algorithm to create sphero-polyhedral particles, we can make sure that, for each particle, the length scales in different directions are kept approximately the same ( as shown in the supplemental material). This paper focuses entirely on simulations with granular columns of circular cross-sections. More general cross-sections are currently under investigation and will be presented in a subsequent publication. The particle/boundary frictional coefficient $\mu_w$ is 0.4, while the inter-particle frictional coefficients $\mu_p$ are 0.2, 0.4, and 0.6.

We treat the relative column size ($R_i/d =$ 2, 2.5, 3.75, 5, 7.5, 10, 12.5, 15, 17.5, 20, 30) as a key parameter. Within one set of simulations, with the same $R_i/d$, we varied the initial height $H_i$ and the inter-particle frictional coefficient $\mu_p$ to obtain the general collapse behavior with a wide range of initial conditions. After the initial Voronoi-based particle packing is generated, we delete the confining tube in the simulation to release all the particles under gravity. After releasing particles to the horizontal plane, the granular material will form a pile of loosely packed grains [Fig. \ref{setup}(b)], with final packing radius, $R_{\infty}$. Thus, the normalized run-out distance, $\mathcal{R}=(R_\infty - R_i)/R_i$, could be obtained. With Voronoi-based spheropolyhedra particles [Fig. \ref{setup}(a) and (b)], we implement the Hookean contact model with an energy dissipation term to calculate the interactions among contacting particles \cite{galindo2010molecular}. In this paper, the contact parameters are the same as previous papers \cite{galindo2013coupled,man2020universality}.

\subsection{\label{2_2}Determination of the run-out distance}

In simulations, the measurement of the final packing radius is complicated due to the lack of a clear mathematical definition of the run-out distance. In cases with small particle/boundary and inter-particle frictional coefficients but large initial aspect ratios, the spread of particles is far-reaching and leads to sparse (single layer) coverage of the area, especially at the front edge. In these cases, it is difficult to determine the edge/boundary and hence the final run-out distance. Thus, we measured the final radius with a histogram of particle distribution for each simulation. Fig. \ref{setup}(c) gives an example of how we measured the run-out distance.

In Fig. \ref{setup}(c), the $x$-axis is the radial position $r$, and the $y$-axis is the percentage of number of particles located within $(r-\Delta r/2, r+\Delta r/2)$ divided by the radial position, $P_{\textrm{dm}}(r) = (1/r)[N(r-\Delta r/2, r+\Delta r/2)/(\Sigma_rN)]$, where $\Delta r$ is the bin width of the histogram, $N(r-\Delta r/2, r+\Delta r/2)$ is the number of particles located between $r-\Delta r/2$ and $r+\Delta r/2$, and $\Sigma_rN$ is the total number of particles in one simulation. FIG. \ref{setup}(c) shows the normalized particle number distribution, in this case, deposition morphology, of a simulation with $\mu_w=0.4$, $\mu_p=0.4$, $R_i = 4$ cm, and $H_i=32$ cm. It shows that most particles locate within $r \leq 27$ cm. Thus, we take the final deposition radius as $R_{\infty} = 27$ cm. After determining the final run-out distance of the collapsed granular column, we obtain the normalized run-out distance, $\mathcal{R}$, the initial aspect ratio, $\alpha$, and the effective aspect ratio \cite{man2020universality}, $\alpha_{\textrm{eff}}$, according to
\begin{linenomath*}
\begin{subequations} \label{eq:normalized_runout}
\begin{align}
    \mathcal{R} &= \left(R_{\infty} - R_i\right)/{R_i}\ ,\ \ \alpha = H_i/R_i\ ,\\
    \alpha_{\textrm{eff}} &= \alpha\left(\mu_w + \beta\mu_p\right)^{-1/2}\ .
\end{align}
\end{subequations}
\end{linenomath*}

\section{Experiments}\label{sec:exp}

\citeA{warnett2014scalings} and \citeA{Cabrera2019Granular} examined the size effect in the granular column collapses with experiments. In this work, we present two sets of experiments to show that changing sizes could influence both the $\mathcal{R}-\alpha$ relationship and the transition point. We perform column collapses of sand particles with grain sizes ranging from 1 mm to 2 mm (the average particle diameter $\bar{d}\approx 1.5$ mm). The experimental setup is the same as that in the simulations. The column diameters are 4.6 cm and 11.4 cm, respectively, resulting in relative column sizes $R_i/d \approx$ 15.3 and 38 [Fig. \ref{setup}(d)-(g)]. We vary the initial height of the column so that the initial aspect ratio varies from $\approx 0.3$ to $\approx 20$. After randomly dropping sand grains into the cylindrical tube, we measure the initial height, $H_i$, of the column, and then quickly lift the tube upward to let the granular column collapse. We measure the final radius of the sand pile in four different directions and take their averages as the final run-out distance, $R_{\infty}$. Then, the relationship between the initial aspect ratio, $\alpha$, and the normalized run-out distance, $\mathcal{R}$, can be obtained accordingly.

We plot the experimental relationship between the initial aspect ratio and normalized run-out distance in Fig. \ref{results}(a). For granular column collapses with both column sizes, as we increase the initial aspect ratio, the normalized run-out distance, $\mathcal{R}$ increases accordingly. We can then observe the clear transition point at $\alpha \approx 1.2$ for cases with $R_i = 5.7$ cm and at $\alpha \approx 2.0$ for cases with $R_i =2.3$. The experiments show that the transition point becomes smaller as we increase the relative column size, and generally, the normalized run-out distance of large columns is larger than that of small columns. The experiments show clear size effect for granular column collapses. It is an experimental proof that changing relative column sizes can not only influence the run-out distance but also affect the transition point on the $\mathcal{R}-\alpha$ relationship. We note that the friction between sand particles and the confining tube while pulling it up might also influence the experimental results, but the influence is not significant, and it does not affect the transition we observe in the $\mathcal{R}-\alpha$ relationship. In the supplemental material, we also present a set of experiments and corresponding simulations of the collapse of cubic particle assemblies. The agreement between experimental results and simulations shows that our DEM model could well capture the behavior of gravity-driven transient granular flows.

\begin{figure}[!ht]
  \centering
  \includegraphics[scale = 0.38]{figure2.png}
  \caption{(a) The relationship between the normalized run-out distance, $\mathcal{R}=(R_{\infty}-R_i)/R_i$, and the initial aspect ratio, $\alpha = H_i/R_i$, of the experiments. The blue dashed line and the blue dash-dot line represent fitting curves before and after \add{the} transition point, when $R_i/d=15.3$. The green dashed line and the green dash-dot line denote the fitting curves before and after \remove{its}\add{the} transition point, when $R_i/d=38$. The black dashed line marks the transition aspect ratio of $R_i/d=15.3$, while the black dash-dot line represents the transition aspect ration of $R_i/d=38$; (b) $\mathcal{R}$ against the effective aspect ratio, $\alpha_{\textrm{eff}}=\alpha\sqrt{1/(\mu_w+2\mu_p)}$ for all the column sizes in the simulations; (c) the relationship between $\mathcal{R}$ and $\alpha_{\textrm{eff}}$ with relative column size $R_i/d =$ 2.5 and 17.5, respectively, to show that the critical inflection point $(\alpha_c,\ \mathcal{R}_c)$ changes as we change the relative size of the column; (d) \remove{shows} the relationship between transitional aspect ratio, $\alpha_c$, and the relative column size, $R_i/d$. (e) \remove{plots} the transitional normalized run-out distance, $\mathcal{R}_c$ against $R_i/d$; (f) \remove{shows} the relationship between $\mathcal{R}(R_i/d)^{\beta_1/\nu}$ and $\mathcal{A}= (\alpha_{\textrm{eff}}-\alpha_{c\infty})(R_i/d)^{1/\nu}$, with $\nu = 1.39$ and $\beta_1 = 0.28$. The inset of the figure plots the same relationship in a log-log coordinate system.}
  \label{results}
\end{figure}

\section{\label{res_discus}Simulation Results and Discussions}

\subsection{Size effects of granular column collapses}

After simulating granular column collapses with various column sizes and frictional coefficients, we obtain the relationship between $\mathcal{R}$ and $\alpha$. Similar to what we have seen previously \cite{man2020universality}, for columns with the same relative size $R_i/d$, less friction leads to larger run-out distances. As we change the $x-$axis to the effective aspect ratio, $\alpha_{\textrm{eff}}$, normalized run-out distance with the same $R_i/d$ collapse onto one curve in $\mathcal{R}-\alpha_{\textrm{eff}}$ space [Fig. \ref{results}(b) and (c)]. The previous work indicates that the effective aspect ratio, which includes the influence of both initial aspect ratio and particle frictional properties, could describe the run-out behavior of the granular column collapse, and most importantly it shows that, for columns with the same relative column size, there is only one transition point for certain column size in the $\mathcal{R}-\alpha_{\textrm{eff}}$ relationship. The collapsing of the transition point of the $\mathcal{R}-\alpha_{\textrm{eff}}$ relationship indicates a possible phase transition related to the collapsing dynamics of granular columns.

However, with various relative column sizes, we see that, in both Fig. \ref{results}(b) and (c), simulation results, especially the transition point in each $\mathcal{R}-\alpha_{\textrm{eff}}$ relationship, varies as we change $R_i/d$. For granular columns with the same initial aspect ratio, granular columns with larger relative column size $R_i/d$ have longer run-out distances. In Fig. \ref{results}(b), we plot simulation results for all the system sizes, and simulation results with $R_i/d = 2.5$ and $R_i/d = 17.5$ in Fig. \ref{results}(c). The latter gives us a clearer picture of how changing relative system size could result in variations in the relationship between $\mathcal{R}$ and $\alpha_{\textrm{eff}}$ and variations in the transition point $(\alpha_c, \mathcal{R}_c)$. This shows that the collapse of the granular column has a significant size effect, which was also observed by \citeA{warnett2014scalings} and \citeA{Cabrera2019Granular}. However, no quantitative studies has been made to universally include different frictional coefficients and boundary conditions to describe the run-out behavior and few researchers pointed out the physical nature of the initial aspect ratio $\alpha$ other than obtaining it from dimensional analysis \cite{man2020universality}.

Further, changing the relative column size also fundamentally influences the shape of $\mathcal{R}-\alpha_{\textrm{eff}}$. For simulations with the same $R_i/d$, a threshold $\alpha_c$ of $\alpha_{\textrm{eff}}$ exists to divide the $\mathcal{R}-\alpha_{\textrm{eff}}$ relationship into two groups [Fig. \ref{results}(b)]. When $\alpha_{\textrm{eff}} < \alpha_c$, $\mathcal{R}$ approximately scales with $\alpha_{\textrm{eff}}$ proportionally. When $\alpha_{\textrm{eff}} > \alpha_c$, $\mathcal{R}$ approximately scales with $(\alpha_{\textrm{eff}})^{0.5}$ with a rather sharp division between the two. Here, the transitional aspect ratio $\alpha_c$ and the corresponding transitional normalized run-out distance $\mathcal{R}_c$ can be seen as the critical aspect ratio and the critical run-out distance, respectively. Both $\alpha_c$ and $\mathcal{R}_c$ vary as we change the size of the granular column. For instance, in Fig. \ref{results}(c), the transition happens at $\alpha_c \approx 8, \mathcal{R}_c \approx 5$ when $R_i/d = 2.5$, and happens at $\alpha_c \approx 2, \mathcal{R}_c \approx 3.5$ when $R_i/d = 17.5$.

We show, in Fig. \ref{results}(b) and (c), that a critical transition point ($\alpha_c$ and the corresponding $\mathcal{R}_c$) exists and also varies with different relative size $R_i/d$. Thus, the $\mathcal{R}-\alpha_{\textrm{eff}}$ curve is dictated by the position of  $\alpha_c$ and $\mathcal{R}_c$ that
\begin{linenomath*}
\begin{equation} \label{eq_general}
    \begin{split}
        \mathcal{R} = f\left( \alpha_{\textrm{eff}}-\alpha_c, \mathcal{R}_c, R_i/d\right).
    \end{split}
\end{equation}
\end{linenomath*}
Both $\alpha_c$ and $\mathcal{R}_c$ decrease with \remove{an} increasing $R_i/d$. In Fig. \ref{results}(d) and (f), based on the relationship between $\alpha_{\textrm{eff}}$ and $\mathcal{R}$ in Fig. \ref{results}(b) and (c), we plot the relationship between $\alpha_c$ and $R_i/d$ and the relationship between $\mathcal{R}_c$ and $R_i/d$. Thus, we could write $\alpha_c$ and $\mathcal{R}_c$ as functions of the relative column size, where, when $R_i/d$ approaches infinity, both $\alpha_c$ and $\mathcal{R}_c$ converge to $\alpha_{c\infty}$ and $\mathcal{R}_{c\infty}$, as shown by
\begin{linenomath*}
\begin{subequations} \label{eq_transition}
    \begin{align}
        \alpha_c = \alpha_{c\infty} + a_1\left( {R_i}/{d}\right)^{b_1}\ , \\
        \mathcal{R}_c = \mathcal{R}_{c\infty} + a_2\left( {R_i}/{d}\right)^{b_2}\ ,
    \end{align}
\end{subequations}
\end{linenomath*}
where $\alpha_{c\infty} = 0.2$ and $\mathcal{R}_{c\infty} = 0.732$ are the fitted critical aspect ratio and the corresponding critical run-out distance, respectively, when the relative column size $R_i/d$ goes to infinity, and both parameters are fitted values. Also, $a_1 = 16$, $b_1 = -0.75$, $a_2 = 5.4$, and $b_2 = -0.3$ are fitting parameters. The fitted curves of $\alpha_c$ and $\mathcal{R}_c$ are the dashed curve in Fig. \ref{results}(d) and the dash-dot curve in Fig. \ref{results}(e).

The power-law decay of both $\alpha_c$ and $\mathcal{R}_c$ with respect to $R_i/d$ inspires us to perform a finite-size analysis of the run-out distance of the collapse of granular columns based on our working hypothesis of a potential phase transition. Since a general scaling analysis often works under thermodynamics limits, where the system size and its corresponding particle number is assumed to be infinity. Finite-size analysis then becomes a powerful tool for investigating real systems with finite sizes. As shown in Fig. \ref{results}(f), the normalized run-out distance, $\mathcal{R}$ indeed exhibits excellent finite-size scaling, suggesting that the transitional aspect ratio, $\alpha_c$, is critical. All the normalized run-out distance data collapse nicely onto a master curve in the form \cite{Torres2015scaling}
\begin{linenomath*}
\begin{equation} \label{eq_fss}
    \begin{split}
        \mathcal{R} = \left({R_i}/{d}\right)^{-\beta_1/\nu}\mathcal{F}_r\left[(\alpha_{\textrm{eff}}-\alpha_{c\infty})\left({R_i}/{d}\right)^{1/\nu}\right]\ ,
    \end{split}
\end{equation}
\end{linenomath*}
with the limiting scaling of the normalized run-out distance being $\mathcal{R}\sim \alpha_{\textrm{eff}} - \alpha_{c\infty}$, where $\nu = 1.39\pm 0.14$ and $\beta_1 = 0.28\pm 0.04$ are obtained to best collapse all the data. In Fig. \ref{results}(f), to better represent the results, we took $\mathcal{A}= (\alpha_{\textrm{eff}}-\alpha_{c\infty})(R_i/d)^{1/\nu}$ as the $x-$axis. Since the system is axisymmetric, the scaling solution for the collapse of granular materials shows similar phenomena in the scaling solutions for connectivity of two-dimensional continuous random networks \cite{Torres2015scaling}, where the scaling parameters $\nu = 1.33$ and $\beta = 0.1389$. This indicates that the flowing behavior of granular column collapse has strong correlations with the connectivity of grain contact networks. This also coincides with the suggestion made by \citeA{mehta2007granular} that the transport of grains takes place percolatively, especially in the context of avalanches.

We then investigate the strong force network of columns with different relative sizes at the very beginning of collapses [Fig. \ref{discussion}(a)-(d)]. In Fig. \ref{discussion}(a)-(d), we only plot the contact forces $f$ that are larger than the mean contact force $\langle f\rangle$ linking the centroids of contacting particle pairs. The blue dashed lines represent the height that the strong force network can reach in the vertical direction, while the red dashed rectangles denote the region which the strong force network could occupy the whole cross-section. As we increase the relative size of the column, the height of the strong force network remains almost unchanged. However, the horizontal occupation of the strong force network is strengthened with the increase of the relative column size. This indicates that larger columns could form better-connected strong force networks in the horizontal direction at the beginning of the column collapse, and therefore might have the capacity to drive the system to flow further than smaller columns.

\begin{figure}[!ht]
  \centering
  \includegraphics[scale = 0.35]{figure3.png}
  \caption{(a) - (d) show strong force networks of granular columns with the same initial aspect ratio but different relative sizes at the beginning of collapse. Solid lines in the figure represent the link between contacting particles with magnitude of contact forces larger than the mean magnitude of contact forces. The blue dashed line represents the height where the strong force network can propagate in the vertical direction, and red dashed lines denote the region where the strong force network could occupy the whole width. The height of the red dashed region is $h_s$, which represent the height of fully occupied strong force networks; (e) The relationship between $\xi/H_t$ and $\alpha_{\textrm{eff}} - \alpha_{c\infty}$ when $\xi$ and $H_t$ are measured at the beginning of the granular collapse, where $\xi$ is the correlation length scale, and $H_t$ is the height of the system at the time we take the measurement. We only measured $\xi$ for four different system sizes and plot the results using circle markers. Meanwhile, the fitted power-law relationships are denoted using dashed lines. (f) shows the relationship between $\xi$ and $h_s/H_t$, where $h_s$ is the height of the fully occupied strong force region in FIG. (a)-(d); (g) Compares simulation data to field data presented in \citeA{Calder1999MobilityOP}. The red \protect\tikzrect{red}{white} markers denote our simulation results, while other markers represent the data compiled from \citeA{Calder1999MobilityOP}.}
  \label{discussion}
\end{figure}

Equation \ref{eq_fss} and the strong force network imply that a length scale exists during the collapse of granular columns in a form that the relative size of the length scale with respect to the column height has a power-law relationship with $\alpha_{\textrm{eff}} - \alpha_{c\infty}$ of the form
\begin{linenomath*}
\begin{equation} \label{eq_length1}
    \begin{split}
        \xi/H_i \sim (\alpha_{\textrm{eff}} - \alpha_{c\infty})^{-\zeta},
    \end{split}
\end{equation}
\end{linenomath*}
where $\zeta$ is the critical exponent for the length scale $\xi$, and $\xi/H_i$ denotes the degree of occupation of sheared granular media along the height of the initial granular column. To further investigate the length scale $\xi$ associated with the collapse of granular columns, we choose four sets of simulations with different relative system sizes. We consider grain pairs along the vertical direction at the very beginning of the collapse when time $t =$ 0.05 s, and define the corresponding correlation $G(z)$, which describes the correlation between particle pairs with vertical distance $z$, based on the equation proposed by Barker and Mehta \cite{Barker1992,mehta2007granular}, by
\begin{linenomath*}
\begin{equation}
    \begin{split}
        G(z) = \frac{\langle\Delta{z_i}\Delta{z_j}\delta\left(|z_{ij}| - z\right)\Theta\left(|d_{ij} - 1/2|\right)\rangle}{\langle|\Delta{z_i}|\rangle^2}\ ,
    \end{split}
\end{equation}
\end{linenomath*}
where $\Delta{z_i}$ and $\Delta{z_j}$ are the vertical displacements of particle $i$ and $j$ after one time step (e.g. $\Delta{z_i} = v_i\Delta{t}$ where $v_i$ is the velocity of particle $i$ and $\Delta{t}$ is the time step), $z_{ij}$ is the vertical distance between two particles, $z$ is the vertical displacement between particles, and $d_{ij}=\sqrt{(x_i - x_j)^2 + (y_i - y_j)^2}$ is the horizontal distance between two particles. Here $\delta(\ )$ is a rectangular function, where the function is equal to one when $z-0.5d < |z_{ij}| \leq z+0.5d$ and equal to zero elsewhere, and $\Theta(\ )$ is the Heaviside step function. This definition ensures that the averages run over all displacements of sphere pairs in the vertical direction. The correlation function $G(z)$ shows how particles influence each other in the $z-$direction.

The displacement correlation function can be fitted with an exponential equation $G(z) = A\textrm{exp}(-z/\xi)$, where $\xi$ can be seen as a correlation length scale associated with the displacement correlation among particles during the collapse of granular columns, and $A$ is a fitting parameter. The occupation of correlated grains across the height of a granular column, represented by $\xi/H_t$, where $H_t$ \remove{here} is the height of the column at the time when we measure the correlation length scale (since $\xi$ is measured at the very beginning of the collapse, $H_t \approx H_i$), indicates the collective motion of particles during a collapse of granular columns. Thus, we plot the relationship between $\xi/H_t$ and $\alpha_{\textrm{eff}} - \alpha_{c\infty}$ in Fig. \ref{discussion}(e), and determine that the ratio between correlation length scale and the system height shows power-law decay as we increase $\alpha_{\textrm{eff}} - \alpha_{c\infty}$, which is consistent with the finite-size scaling of both run-out behavior of granular column collapses, and $\zeta = 0.89$ in Eq. \ref{eq_length1} best fits the power-law decay. When $\alpha_{\textrm{eff}}$ approaches $\alpha_{c\infty}$, $\xi/H_t$ starts to deviate from the power-law relationship. This is because when the initial aspect ratio is too small, only several layers of particles are present along the height of a granular column.

Additionally, we calculate the correlation length scales, $\xi$, of the four cases shown in Fig. \ref{discussion}(a)-(d), and plot them against the ratio between the height of the fully occupied strong force region $h_s$ and the column height $H_t$ [Fig. \ref{discussion}(f)]. We find that, as we increase the relative column size, not only does the height of the fully occupied strong force region increase, but the correlation length scale also increases accordingly. This implies that a larger system could help establish a better correlated moving mass, and as we increase the initial aspect ratio, columns with larger relative size, where well-correlated moving mass is easier to establish, would experience the phase transition from quasi-static collapse to inertial collapse earlier than smaller systems. To further analyze the phase transition associated with granular column collapses, we would have to conduct rigorous investigations on the mean-field theory, which we plan to do in future studies.

\subsection{Application to geophysical flows}
\label{sec:application}
In order to present an application of granular column collapse results to a real geological situation, we also compare our simulation data to the rock falls around Soufriere Hills Volcano, Montserrat as a result of volcanic eruptions and the resulting pyroclastic flows and granular surges [presented in \citeA{Calder1999MobilityOP}]. Although granular column collapses do not behave exactly the same as landslides, rock avalanches or pyroclastic flows, their underlying physics is similar since the mobility of collapsed materials is mainly governed by the ability of transforming potential energy into kinetic energy, by surpassing the dissipative action of friction, a competition of effects that is summarized in the proposed $\alpha_{\textrm{eff}}$ as explained before. \citeA{Roche2002ExperimentsOD,Roche2008ExperimentalOO} also argued that the dam-break experiments could reveal similar behavior as inertial geophysical flows and the results supported the pyroclastic flow modeling approach by \citeA{Levine1991HydraulicsOT}. In \citeA{Calder1999MobilityOP}, $R_\infty/H_i$, which indicates the ability of relocating materials, is plotted against a parameter, $\Psi = \rho gV/H_i^2$, where $V$ is the volume of the material being transported. In Fig. \ref{discussion}(g), we plot our simulation results along with the results presented in \citeA{Calder1999MobilityOP}, where the \tikzrect{red}{white} markers indicate our simulation results, while other markers represent the field data.

We see that most of simulation data do not fall into the region of field observations. This is mainly due to the length scale of simulations being much smaller than a typical geophysical flow. However, as we increase the system size of simulations, we do obtain granular collapses in the same regime as some of the pyroclastic flows and granular surges. The finite-size analysis early in this section implies that increasing the system size could lead to an asymptotic behavior of granular flows, where the transition point $\alpha_c \rightarrow \alpha_{c\infty}$. We can also see this trend in Fig. \ref{discussion}(g) where, as we increase the system size of simulations, our simulation data merge into the field results. More intriguingly, this also shows that almost all the data acquired from dome-collapse pyroclastic flows [\tikztri{black}{white} in Fig. \ref{discussion}(g)] fit in the trend of simulation results; this type of pyroclastic flows is considered to have a flow behavior more similar to granular flows. The size effect of granular column collapses might be potentially connected to predict the behavior of geophysical flows with different length scales. More importantly, the scaling of run-out distances with respect to $R_i/d$ might help explaining the long-range run-out in some pyroclastic flows \cite{Lube2019GenerationOA}, where the initial radius of the pyroclastic ``column" is much larger than the volcanic particles. The analysis of force networks within granular systems, as mentioned in Section \ref{sec:size}, may be helpful to explain some phenomena of geophysical flows. However, it is still difficult to directly obtain the contact information within real geophysical systems.

\section{\label{Conclusion}Conclusions}

Previous research concluded that we could combine the influence of initial and boundary conditions and the initial column aspect ratio to determine the normalized run-out distance with an effective aspect ratio, $\alpha_{\textrm{eff}}$, which introduced a universal relationship to link the behavior of granular column collapses in three different collapsing regimes (quasi-static, inertial, and liquid-like) \cite{man2020universality}. In this paper, we further investigated the size effect associated with the collapse of granular columns. Our research is performed with DEM simulations of Voronoi-based spheropolyhedron particles. We found that the transition point in $\mathcal{R}-\alpha_{\textrm{eff}}$ space, which distinguishes the inertial collapse regime from the quasi-static collapse regime, varies as we change the relative system size $R_i/d$. Both $\mathcal{R}_c$ and $\alpha_c$ experience power-law decay with respect to $R_i/d$, which implies possible finite-size scaling for the normalized run-out distance $\mathcal{R}$.

Similar to the finite-size analysis of the jamming transition of granular materials \cite{Liu2014finite}, where the stress scales with the solid fraction $\phi - \phi_{c\infty}$, we took the previously determined effective aspect ratio $\alpha_{\textrm{eff}}$ as the key parameter, and discovered that both the run-out distance of granular column collapse and the energy consumption of it follows strong finite-size scaling. Interestingly, the scaling parameters we discovered for the finite-size scaling of the run-out distance of granular column collapses are similar to those presented in percolation problems of two-dimensional random networks \cite{Torres2015scaling}, which may imply that the behavior of granular column collapses is strongly influenced by the contact network presented inside the column during the collapse. Additionally, to better understand the scaling of granular column collapses, we further analyze the correlation length scale at the very beginning of the collapse. Simulation results show that, as we increase $\alpha_{\textrm{eff}} - \alpha_{c\infty}$, the length scale $\xi/H_t$ shows a power-law decay with the critical exponent $\zeta\approx0.89$. This study is based on our previous work where we introduced a physics-based dimensionless number (ratio between inertial effects and frictional effects) to describe the behavior of granular column collapses, and we further expand our analysis to include the size effect, which is crucial to applications in engineering, such as chemical engineering, food processing, civil engineering. We also associate the size effect to the strong force networks at the beginning of the column collapse to show that a larger column size results in a well percolated strong force network, which may contribute to relatively larger run-out distances.

However, the granular system we are studying is still axisymmetric, and our preliminary results show intriguing phenomena when the cross-section of a granular column is no longer axisymmetric. Also, we note that the system we are dealing with has no interstitial fluid, and only Voronoi-based particles are considered. In natural and engineering applications, most granular systems are wet, saturated, or submerged in liquid. A rigorous investigation on the fluid-solid interaction and its influence on the bulk behavior, which is beyond the scope of this study, is need to thoroughly describe the behavior of granular-fluid systems. Further investigations related to granular-fluid systems, the influence of cross-section shapes, and the influence of particle shapes will be conducted and presented in future publications.

\section*{Data availability Statement}
The data associated to this work is archived and published as \citeA{Man2021data} in Mendeley Data, which can be accessed freely after registration. Part of the data in Fig. \ref{discussion}(g) is available through \citeA{Calder1999MobilityOP}. The simulation and the corresponding model were described in \citeA{Mechsys} and \citeA{galindo2013coupled}. The software can be obtained from \url{http://mechsys.nongnu.org/}.

\acknowledgments
The authors acknowledge the financial support from Westlake University and thank the Westlake University Supercomputer Center for computational resources and related assistance. H.E. Huppert acknowledges with gratitude the hospitality of his co-authors while he was at Westlake University.

\section{} command to identify level 1 heads;

Math coded inside display math mode \[ ...\]
 will not be numbered, e.g.,:
 \[ x^2=y^2 + z^2\]

 Math coded inside \begin{equation} and \end{equation} will
 be automatically numbered, e.g.,:
 \begin{equation}
 x^2=y^2 + z^2
 \end{equation}

\begin{eqnarray}
  x_{1} & = & (x - x_{0}) \cos \Theta \nonumber \\
        && + (y - y_{0}) \sin \Theta  \nonumber \\
  y_{1} & = & -(x - x_{0}) \sin \Theta \nonumber \\
        && + (y - y_{0}) \cos \Theta.
\end{eqnarray}