# MoSES_2PDF: A GIS-Compatible GPU-accelerated High-Performance Simulation
Tool for Grain-Fluid Shallow Flows

Chi-Jyun Ko
Dept. Hydraulic and Ocean Engineering
National Chen Kung University
70101 Tainan City, TAIWAN
&Po-Chih Chen
Dept. of Computer Science and Information Engineering
National Chen Kung University
70101 Tainan City, TAIWAN
&Hock-Kiet Wong
Dept. Hydraulic and Ocean Engineering
National Chen Kung University
70101 Tainan City, TAIWAN
& Yih-Chin Tai*
Dept. Hydraulic and Ocean Engineering
National Chen Kung University
70101 Tainan City, TAIWAN
<EMAIL_ADDRESS>

###### Abstract

We introduce a GPU-accelerated simulation tool, named Modeling on Shallow
Flows with Efficient Simulation for Two-Phase Debris Flows (MoSES_2PDF), of
which the input and output data can be linked to the GIS system for
engineering application. MoSES_2PDF is developed based on the CUDA structure,
so that it can well run with different NVIDIA GPU cards, once the CUDA vers.
9.2 (or higher) is installed. The performance of the MoSES_2PDF is evaluated,
and it is found that the present GPU-CUDA implementation can enhance the
efficiency by up to 230 folds, depending on the PC/workstations, models of GPU
card and the mesh numbers in the computation domain. Two numerical examples
are illustrated with two distinct initial inflow conditions, which are
included in two modes of MoSES_2PDF, respectively. In the numerical example of
a large -scale event, the 2009 Hsiaolin event, the results computed by two
distinct NVIDIA GPU cards (RTX-2080-Ti and Tesla-V100) are found to be
identical but tiny deviation is figured out in comparison with the results
computed by the conventional single-core CPU-code. It is speculated to be
caused by the different structure in the source codes and some float/double
operation. In addition to the illustration in GIS system, the computed results
by MoSES_2PDF can also be shown with animated 3D graphics in the ANSI-
Platform, where the user can interact with 3D scenes. The feasibility,
features and facilities of MoSES_2PDF are demonstrated with respect to the two
numerical examples concerning two real events.

_Keywords_ GPU-accelerated $\cdot$ CUDA $\cdot$ GIS-compatible $\cdot$ Grain-
fluid shallow flows $\cdot$ Simulatiom tool

##  1 Introduction

Many of the geophysical hazardous flows, such as landslides or debris flows,
usually take place in mountain area and flow on non-trivial topography. It is
clear that the geometry of the topography may have significant impacts on the
flow paths. For the sake of hazard assessment, risk management or evaluation
of disaster mitigation countermeasures, numerical simulation is commonly a
powerful tool for scenario investigation. Hence, it is highly requested to
develop an effective tool of high efficiency and high performance for
providing useful information concerning the flow behavior, such as the paths,
velocity and flow thickness, etc. A successful simulation of hazardous flows
on real topography depends on the availability of the digital elevation model
(DEM) essential. With the rapidly growing popularity of the UAV and the
integrated modern remote sensing technology, the expensive high resolution
data, such as the LiDAR-derived DEM or the aerial photos, has become
affordable and enriched the Geographic Information System (GIS) a lot. From
the viewpoint of scenario investigation, one of the urgent tasks would be the
integration of the simulation tools with the GIS as well as with a high
efficient computational performance. Numerical simulation based on Physics-
based dynamical models can help the quantitative investigation on the impacts
caused by the plausible geophysical hazardous flows. However, in addition to
the high-resolution DEMs, it also requests expensive computer resources,
depending on the resolution and computation domain. The launch of Compute
United Device Architecture (CUDA), a parallel computing platform and
programming model, provides an econamic solution for speeding up the
computational efficiency by means of employing the graphical processing units
(GPUs). The computation duration can therefore be shorten to be less than the
scenario time. Hence, either computation with real-time illustration or
training for model parameters has become achieveable in an efficient way even
without expensive computer resources.

In the past decades, many efforts have been devoted to modeling these highly
destructive flows. Since these flows are generally thin (measured normal to
the basal surface) in comparison with the large extension in the direction
tangential to the topographic surface, shallowness assumption and depth-
integration are often introduced for reducing the complexity of the model
equations as well as of the numerical approach. For the ease of deriving the
governing equations, many of the models are given in Cartesian coordinates
(e.g., Pitman and Le, 2005; Li and Duffy, 2011), where the topography is added
on the horizontal or inclined plane, or in the locally inclined Cartesian-type
coordinates (e.g., Pudasaini, 2012; Pudasaini and Mergili, 2019). For an
appropriate description of the depth-averaged velocity as well as for the
employment of the shallowness assumption, the terrain-fitted coordinate system
has been introduced. In addition to the curved-twisted coordinate proposed by
Pudasaini and Hutter (2003), the remarkable works for general topography are
Bouchut and Westdickenberg (2004); Tai and Kuo (2008); Tai et al. (2012) and
Luca et al. (2016). As well, assuming the non-trivial topographic surface to
be composed of finite number of triangles, local coordinate system is
introduced on each triangular element for reducing the complexity of numerical
implementation in a unstructured mesh system (see Rauter and Tuković, 2018;
Tai et al., 2021).

Regarding the composition of the flow body, models for hazardous flows can be
divided into four categories: a) single-phase model; b) quasi-two-phase
models; c) two-phase models; and d) multi-phase approaches. The key issue by
the single-phase approach is the introduction of an appropriate rheological
relation, which is definitely different to the Newtonian one due to the
complex composition. Taking into account the existence of yield stress, models
with the Bingham plastic relation (e.g., Liu and Mei, 1989; Coussot, 1997) or
with the Herschel-Bulkley law (e.g., Huang and Garcia, 1998; Ancey, 2007) have
been suggested. In addition, a generalized visoplastic model is proposed by
Chen (1988), covering the models by Bagnold’s model (Bagnold, 1954) and
Takahashi’s constitutive relation (Takahashi, 1978). For the ease of numerical
implementation, commercial codes are mainly released based on the single-phase
models for engineering purposes. One of the most popular ones is the FLO-2D,
initially following the rheological model proposed by O’Brien and Julien
(1985, 1988) and improved by many sequential works (e.g. O’Brien et al.,
1993).

The quasi-two-phase or the two-phase models are based on the concept of two-
phase mixture. However, in the quasi-two-phase models, the relative velocity
between the two constituents is assumed to be very small, so that its impacts
are assumed to be insignificant and only one velocity is considered for the
conservation of linear momentum (e.g., Iverson, 1997; Iverson and Denlinger,
2001; Pudasaini et al., 2005; Tai and Kuo, 2012). Since only one velocity is
taken into account, the dynamic behavior of the flow body is similar to a
single-phase fluid (if entrainment or deposit at the bottom is not included),
and the effects of the concentration are then parametrically (instead of
dynamically) considered in the rheological relation. In addition to the type
of Coulomb-mixture (Iverson, 1997; Iverson and Denlinger, 2001), another
approach is the quasi-two-phase grain-flow model, widely employed in Japan for
engineering applications, e.g., Egashira (1997, 2007); Nakatani et al. (2008);
Liu et al. (2013); Nakatani et al. (2016). Emphasizing the the role of the
interstitial fluid on the behavior of flows, the Coulomb-mixture theory has
been improved by introducing the non-hydrostatic pore-fluid pressure for the
dilatancy effects by the solid constituent (Iverson, 2009; Iverson and George,
2014).

In the two-phase models, each constituent is described by individual momentum
equation and has its own velocity. Hence, the key feature of a two-phase model
is the capability of describing the phase separation, which is documented in
field surveys and laboratory experiments. To our limited knowledge, the first
inspiring two-phase model for debris flows would be Pitman and Le (2005),
although the fundamental investigation on two-phase mixture theory for grain-
fluid has been made since years (e.g., Anderson and Jackson, 1967). Based on
Pitman and Le (2005), Pelanti et al. (2008) introduced the well-balancing
property in the numerical implementation, where the variable topography is
included. Pailha and Pouliquen (2009) suggested a two-phase description for
the initiation of underwater granular avalanches. Luca et al. (2012) has
theoretically considered over-saturated mixtures in coordinate of general
topography. The quadratic viscous drag force and virtual mass force are
considered in Pudasaini (2012). A two-phase model in two-layer type and taking
into account the dilatancy effects can be found in Bouchut et al. (2016),
which mainly considers the submarine debris flows, i.e. the upper layer is
pure fluid and not thin. With two-layer approach, the partially unsaturated
grain-fluid mixtures are modeled by Meng et al. (2017), where a fluid-
saturated granular layer is overlaid by a pure granular material. With
thermodynamically consistent derivation, models of granular-fluid mixtures
with non-hydrostatic pore pressure can be found in Heß et al. (2017, 2019),
where the solid phase is assumed hypoplastic. The multi-phase approach can be
found in Pudasaini and Mergili (2019), where the flow body is composed of
three phases (the coarse and fine solid fractions, and the viscous
interstitial fluid), and the model equations are given in the locally inclined
Cartesian-type coordinates for describing the flows over a non-trivial
topography.

Regarding the GIS-compatible simulation tool, Mergili et al. (2017) introduced
a GIS-supported open-source computational framework, r.avaflow, which is
developed to link to the open-source GIS package, GRASS GIS. The r.avaflow is
based on a two-phase solid-fluid mixture model (Pudasaini, 2012), and mass
exchange at the basal surface (entrainment or deposition) are taken into
account. In addition, r.avaflow allows multiple-model-runs that multiple
parameter sets can be simultaneously executed with the parallel processing
function. Although it does not shorten the duration for a single run (a
parameter set), this facility may help in looking for the optimal parameter
set or collecting abundant plausible simulation scenarios. The r.avaflow has
been extended for the multi-phase flow model (Pudasaini and Mergili, 2019),
see Mergili (2014-2020), and several engineering applications are available in
Mergili et al. (2020); Baggio et al. (2021).

Although many of the above-mentioned models are with numerical implementation
and able to provide satisfactory results, the computation processes are
usually very time-consuming and cost a lot of computer resources. In order to
improve the efficiency of computation, one may speed up the computing
applications by harnessing the power of graphical processing units (GPUs).
Although GPU was originally designed for graphics rendering on computers,
especial for the illustration of 3D objects, GPU-accelerated high-performance
computation has become more and more popular in various scientific fields,
such as Artificial Intelligence (AI), machine learning (ML), deep learning
(DL) etc. For example, Castro et al. (2011), Brodtkorb et al. (2012), de la
Asunción et al. (2013), and Aureli et al. (2020) have proposed GPU-accelerated
approach based on the shallow water equations (SWEs), while Dazzi et al.
(2019) integrated the GPU-accelerated 2D SWEs code with levee breach
mechanism. A GPU-enhanced 2D SWEs code for curved triangular meshes is
developed by Wu et al. (2021), of which the scheme is entropy stable and
provides high-order accuracy.

However, there is no GIS-compatible and GPU-accelerated simulation tool for
two-phase grain-fluid flows yet. Based on the two-phase model by Tai et al.
(2019), we developed the GPU-accelerated simulation tool, MoSES_2PDF, which
follows the CUDA structure to integrate the GPGPU technology. In addition, the
MoSES_2PDF is designed to be GIS-compatible. That is, the input data is in the
asc-format and can be directly acquired from GIS systems (ArcGIS, QGIS), and
the output results can be transformed into asc-format for illustration in GIS
systems. On the other hand, in MoSES_2PDF we also use the GPU graphics
rendering functions in OpenGL for 3D numerical results presentation, which
combines the computation results with 3D objects to mimic the effect of
reality simulation.

In the second section, we brief the grain-fluid model and employed numerical
implementation. The introduction of the framework of MoSES_2PDF and its
performance with the GPU parallel computation facility are given in Sect. 3,
two modes for different initial inflow conditions are designed and three
output formats available for results illustration. In the fourth section, two
application examples are illustrated, which are two historical events taking
place in 2009 in Taiwan. The application potential and the key features of
MoSES_2PDF are summarized and highlighted in the concluding remarks.

##  2 Modeling and Numerical Implementation

###  2.1 Modeling

In the GIS system, the topography is described by a set of altitudes of
terrain locations over a regular horizontal grid, the digital elevation model
(DEM). Based on the DEM, a terrain-fitted coordinate system
${O}_{\xi\eta\zeta}$ can be introduced in the vertical-horizontal-oriented
Cartesian coordinates ${O}_{XYZ}$, where the $X$\- and $Y$-axis lie on the
horizontal plane and the direction of $Z$ points upwards. The system
${O}_{\xi\eta\zeta}$ can be chosen in the way that, e.g., the projections of
$\xi$-axis and $\eta$-axis on the horizontal plane coincide with $X$-axis and
$Y$-axis, respectively (cf. Fig. 1). The coordinates between $O_{XYZ}$ and
$O_{\xi\eta\zeta}$ is related by the transformation matrix
${\boldsymbol{\Omega}}$, see Tai et al. (2012) or Luca et al. (2016) for
details.

Figure 1: Coordinate system $O_{XYZ}$ and ${O}_{\xi\eta\zeta}$.

In MoSES_2PDF, the evolution of the flow body is computed by a depth-
integrated two-phase grain-fluid model for debris flows over rugged topography
(Tai et al., 2019), where the conservations of mass and momentum are taken
into account and both of the coordinate systems, ${O}_{XYZ}$ and
${O}_{\xi\eta\zeta}$, are used. Hence, there are 6 model equations for the
evolutions of the 6 physical quantities, $(h^{s,f},v_{X}^{s,f},v_{Y}^{s,f})$.
Here and hereafter, the superscripts “$s$” and “$f$ ” indicate the variables
for the solid and fluid constituent, respectively, and the subscript X,Y,Z
indicates component in $O_{XYZ}$, while ξ,η,ζ for the one in
$O_{\xi\eta\zeta}$. For a concise description, the employed model equations
are not shown here, but presented in A for readers.

Since the flow body is treated as a two-phase mixture, in
$(h^{s,f},v_{X}^{s,f},v_{Y}^{s,f})$, $h^{s}=h\phi^{s}$ and $h^{f}=h\phi^{f}$
are defined, where $\phi^{s,f}$ denoting the depth-averaged volume fraction
and $h$ represents the flow thickness measured in the direction normal to the
basal topography. In addition, it should be noted that the variables
$v_{X}^{s,f}$ and $v_{Y}^{s,f}$ are the $X$\- and $Y$-component of the flow
velocity projected on the horizontal plane. With the help of the
transformation matrix ${\boldsymbol{\Omega}}_{b}$ on the basal surface, one
can get the flow velocities tangential to the topographic surface.

###  2.2 Numerical implementation

MoSES_2PDF is developed based on the CUDA structure (CUDA Toolkit vers. 11.0).
The employed numerical implementation is same as in Tai et al. (2019), that
the balance equations with respect to the solid and fluid constituent are
computed separately at each time step, where the anti-diffusive, non-
oscillatory central scheme (adNOC, see Kurganov and Tadmor (2000); Kurganov et
al. (2007)) is used for solving the six differential equations. For details of
the numerical implementation, we refer the readers to Tai et al. (2019), in
which the code is developed for single-processor (single-core) computation.

The high resolution in space is achieved by utilizing the cell reconstruction,
such as the total variation diminishing (TVD), essentially non-oscillatory
(ENO) or weighted essentially non-oscillatory (WENO) scheme. And the high-
order accuracy in time is carried out by the multi-step Runge-Kutta (RK)
method with respect to a semi-discrete formulation given in Kurganov and
Tadmor (2000). In MoSES_2PDF, the Minmod TVD slope limiter is adopted and the
two-step modified Euler (second-order RK) is employed. As indicated in
Kurganov and Tadmor (2000), the time step should be limited by the Courant-
Friedrichs-Lewy (CFL) condition with a number of 0.125 and the CFL = 0.1 is
set in MoSES_2PDF.

##  3 Framework of MoSES_2PDF and its Performance

###  3.1 Framework

Figure 2: The framework of computation

####  3.1.1 Teamwork of CPU and GPU

With respect to the hardware structure, the computational framework of
MoSES_2PDF can be divided into two sections, the CPU- and GPU-section, where
they have to communicate with each other to complete the operation of
computation. Since the data transfer between the sections could be very time-
consuming, it is therefore designed in the way to minimize the data transfer
as possible, that the CPU-section mainly processes file reading and output,
and the rest computational procedures are operated in the GPU-section, cf.
Fig. 2.

According to the characteristics of applied numerical scheme, the execution
hierarchy of single instruction corresponding to multiple threads (SIMT) is
employed as the GPU computing architecture for the purpose of a high
efficiency in computation. Besides, most of the calculations are matrix
operations, yielding the fact that the memory location of a matrix may cover
several blocks, depending on the size of the computational domain (grid
numbers). It leads to a hurdle in looking for the global extreme (maximum or
minimum) value in the matrices, because the function of global synchronization
is not supported in the CUDA parallel-programming model (library). If one uses
the type of loop for searching over the whole domain in the CPU-section, the
data transfer/collection among the blocks will significantly retard the
computational efficiency. Hence, the extreme-value-searching plays a crucial
role in efficiency enhancement, and one has to find an effective way for data
communication among different blocks. In MoSES_2PDF, we adopt the method of
optimizing parallel reduction, proposed by Harris et al. (2007), in the CUDA
code to compute the maximum/minimum values in the entire domain. That is,
taking the maximum value as example, the maximum value on each block is first
figured out, and then they will be transferred to a specified block, in which
the global maximum value is identified. In this way, the computation
efficiency of the searching process is improved with a significantly higher
efficiency than the computation of using the loop type.

####  3.1.2 Operation of MoSES_2PDF

The operation of MoSES_2PDF consist of four main stages (cf. Fig. 2): (0)
Preparation; (1) Input and parameter setting; (2) Computation; and (3) Results
illustration:
Stage 0: Preparation
Based on different initial inflow types at the boundaries of computational
domain, MoSES_2PDF provides two modes: In the first mode (Mode-I), finite
masses with specified concentration(s) are to be released in the computational
domain, where there is no mass inflow from the boundary in the initial
condition. This mode is mainly used for the simulation of landslide type, and
this type of releasing finite mass(es) can be found in many works (e.g., Gray
et al., 1999; Pitman and Le, 2005; Kuo et al., 2009; Mergili et al., 2017;
Pudasaini and Mergili, 2019; de’Michieli Vitturi et al., 2019). The second
mode (Mode-II) is used when there are materials flowing into the computational
domain through the boundaries. That is, the inflow condition (the hydrograph
of the flow thickness, volume fraction (concentration) of solid/fluid phase,
velocities and the inflow locations) at the specified locations of boundary
can be defined in Mode-II. This kind of inflow condition is often adopted in
debris flow simulations, such as the debris flow simulator Kanako 2D (e.g.,
Nakatani et al., 2008; Liu et al., 2013; Nakatani et al., 2016), developed by
the team in Kyoto university, Japan. In both modes, no reflection boundary
condition is set in the whole process of computation. In addition to choosing
the suitable mode, one has to prepare the DEM of the specified computational
domain and its associate initial condition (the flow thickness, concentration
of solid/fluid phase and velocity distribution of the flow body). The DEM is
in the asc-format, which can be directly extracted from GIS system (ArcGIS or
QGIS). According to the initial inflow type, the data file, the initial
configuration for Mode-I or the inflow condition for Mode-II, should be
prepared (in ASCII-format).
Stage 1: Input and parameter setting
Once the code is initiated, a text-file “par_list” is to be read first, on
which all the model parameters, the total simulate time, the time interval for
results output, and the locations of the DEM file, inflow/initial condition
file can be specified. The model parameters consist of: (a) the angle of basic
friction $\delta_{b}$ between the solid constituent and the basal surface; (b)
the drag coefficient $C_{d}$ between the solid and the fluid constituents,
caused by the related velocity between the two constituents; (c) the viscosity
coefficient $N_{R}$, similar to the Reynolds number with value inversely
proportional to the viscosity; and (d) the fluid friction coefficient
$\vartheta_{b}$. Since the employed model is for two-phase solid-fluid flows,
the initial concentration of solid phase $\phi^{s}_{0}$ is also given in the
text-file “par_list”.
Stage 2: Computation
Once the reading process is completed in the CPU-section, all the data will be
transmitted to the GPU-section for computing the evolutions of the variables.
When it arrives at the specified time level for results output, the data is
transferred from the GPU memory to the host memory for the output operation by
CPU, cf. Fig. 2.

As elaborated in Sect. 2.2, the time interval should be determined according
to the assigned CFL condition, for which the global maximum wave speed is
needed. Since the time interval has to be determined at each time step, this
process is the most time-consuming one. With the method of optimizing parallel
reduction mentioned in Sect. 3.1.1, the efficiency of the model computation
can be enhanced by ca. 20-fold.
Stage 3: Results illustration
The output of the computed results includes the total flow thickness (measured
normal to the basal surface), the fraction of solid phase, the $X$\- and
$Y$-components of velocity for the solid and fluid phases, respectively. In
the data processing, the results can be converted in three different types,
depending on the requirement of illustration platform. The first is the asc-
file format for illustration in GIS system (ArcGIS or QGIS) as well as for the
subsequent applications. The second is in the format for drawing contour maps
by Python, usually for scientific investigation. The third type is for a
three-dimensional scenario illustration tool (Advanced Numerical Scenario
Illustration Platform, ANSI-Platform), developed by Y.C. Tai and his coworkers
with OpenGL. The ANSI-platform provides animated 3D graphics with user-
interactive facility, that one can interact with 3D scenes (zoom in/out,
translate or rotate) for examining the simulation results. In case the
corresponding ortho-photo of satellite image is available, the animated 3D
topography can be integrated with the photograph on its surface for the
purpose of or scenario investigation, hazard assessment or disaster
mitigation.

###  3.2 Performance of GPU parallel computation

The performance of the MoSES_2PDF with/without GPU-CUDA computation is
evaluated by a campaign of the accumulated elapsed time (duration of
computation) for running 10,000 steps with various number of meshes. In this
campaign, computations by single-process of PC (i7-8700<EMAIL_ADDRESS>GHz, 16 GB
memory, Linux OS) and by workstations with two different models of NVIDIA GPU
cards (Geforce RTX 2080 Ti with i7-8700<EMAIL_ADDRESS>GHz and Tesla V100-SXM2-32GB
with Xeon Gold 6154<EMAIL_ADDRESS>are performed. Altogether, there are 6
numbers of meshes (i.e., 10,000, 50,000, 100,000, 250,000, 500,000 and
1,000,000) in the computational domain for evaluating the performance.

The accumulated elapsed time is calculated by the averaged duration of three
repeated runs for minimizing the time deviations by each computation. As shown
in Fig. 3, the duration of computation increases for more meshes, where the
solid lines indicate the accumulated elapsed time and the dashed lines stand
for the speedup ratio. The speedup ratio reflects the increasing rate, of
which the value is determined by the duration for the GPU computation divided
by the elapsed time for the CPU single-processor computation. For computation
with 10,000 meshes, the GPU speedup ratio is about 20-fold, while its
magnitude rapidly increases as the number of meshes grows up to 100,000. When
the number of meshes reaches 500,000, the speedup ratio exceeds 150-fold,
while it is about 160-fold for computation with 1,000,000 meshes.

Figure 3: Performance of CPU and GPU and the associated speedup ratio

It is interesting to find that there is one nick of speedup ratio at mesh
number of 250,000 with respect to the computation by GPU-Tesla V100, although
the elapsed time monotonically increases as the mesh number grows. This nick
is suspected to be caused by the arrangement of memory allocation in the
Tesla-GPU card, e.g., the block size in threads as elaborated in Ryoo et al.
(2007). In general, the enhancement of the efficiency is extremely
significant, and the more the meshes there are, the better the efficiency is
improved. The employment of the GPU computation sheds the light on the
achievement of a real-time computation and illustration for scenario
investigation of hazard assessment or disaster management.

##  4 Application examples

Two application examples are introduced: the 2009 Hsiaolin landslide and the
DF004 debris flow event. Both of the events took place during the typhoon
Morakot in 2009 and locate at southern Taiwan. More than 470 people were
killed in the 2009 Hsiaolin event, and 41 victims are identified during the
debris flow event DF004 (cf. Dong et al., 2011; Soil Water Conserv. Bureau,
2009). In the 2009 Hsiaolin landslide event, all the blocks of released mass
lay in the computational domain (see Fig. 4a). In the computation, the blocks
are released simultaneously, sliding downward and attacking the Hsiaolin
village, so that this event is categorized as the landslide type and the
MoSES_2PDF-I (Mode-I) is utilized. In the DF004 debris flow event, a large
amount of soil, sand and water rushed from upstream into the Nansalu village
along the gully, so that the Mode-II is used for the simulation.

In the two application examples, there are totally five parameters to be set,
of which the four, $(\delta_{b},C_{d},N_{R},\vartheta_{b})$ , are the material
parameters and $\phi^{s}_{0}$ is the solid concentration of the initial/inflow
flow body. For showing the feasibility of MoSES_2PDF, they are set identical
to the ones used in previous studies (Tai et al., 2019, 2020), as listed in
Table 1. In addition, for ease of depicting the dynamic evolution of the flow
body as well as identifying the flow paths, only the area covered by flow
thickness more than 10 cm will be shown in the results illustration, if not
additionally specified.

Table 1: Parameter values used in the application examples Parameter | $\delta_{b}$ | $C_{d}$ | $N_{R}$ | $\vartheta_{b}$ | $\phi^{s}_{0}$
---|---|---|---|---|---
Value | $16^{\circ}$ | $6.0$ | $268$ | $5.0$ | $0.5$

###  4.1 The 2009 Hsiaolin Landslide

The 2009 Hsiaolin Landslide is of an abnormal large-scale, that the largest
scar area caused by the released mass covers about 57 hectares and the total
volume exceeds 20 million cubic meters (Dong et al., 2011; Kuo et al., 2011).
Because of the huge size and the severe casualties, many post-event
investigations have been conducted. Based on the difference of the pre- and
post-event DEMs, landslide blocks are identified for reconstructing the flow
paths, (see Kuo et al., 2011; Tai et al., 2019). The computational domain and
mesh size are identical as in Tai et al. (2019), that it covers $3.7\times
2.2$ km2 and with resolution $\Delta x=\Delta y=10$ m. Hence there are totally
85,956 cells, including 3 ghost cells at each boundary, in the computation.
For the sake of validating the GPU implementation in MoSES_2PDF, identical
parameters are adopted (see Table 1), and the total simulation time is set to
181.82 s, when the whole flow body is nearly at the state of rest.

With identical initial condition and parameters, the single-processor (core)
computation and GPU parallel computation with RTX 2080 Ti and Tesla V100 are
performed. The single-core computation took about 2.5 hours, while the GPU
parallel computations with CUDA implementation took only ca. 80 s and 56 s
with RTX-2080-Ti and Tesla-V100, respectively. Although the GPU computations
are identical (12,427 steps), they are slightly deviated from the CPU results
(12,454 steps). After a tedious checking process between the CPU- and GPU-
codes as well as results, it is found that the difference of the first time
step occurs at the eighth decimal place. The discrepancy is suspected to be
caused by the coding structure and cut-off error in the float/double operation
in the GPU-computations. Nevertheless, regarding the final deposition area the
difference is rather minor. At t = 181.82 s, the deposited mass covers 14,872
meshes where only 25 are in difference, if the meshes with thickness less than
1 cm are isolated. In case only flow thickness more than 10 cm is taken into
account, there are 23 in 13,600 meshes to be deviated. In both cases, the
deviation is less than $0.17\%$. Hence, only the facility of GPU computation
is provided in MoSES_2PDF.

a)
b)
c)

Figure 4: Simulation of Hsiaolin event: (a) The initial flow thickness, where
the red line outlines the Hsiaolin village; (b) Distribution of the solid
concentration at $t=181.83$ s; and (c) The computed flow path with the
satellite orthophoto

Figure 4b depicts the computed distribution of solid phase over the contour
topographic map at $t=181.82$ s. It is found that the high concentration of
solid phase approximately coincides with the locations of valley. In Fig. 4c,
the flow paths are indicated by the cyan shade, which is illustrated on the
orthophoto of satellite image taken about 6 months after the event. The bar
areas in the satellite image can be recognized beneath the cyan shade and they
approximately point out the zones rushed by the landslide mass. The sound
agreement reveals the feasibility of the MoSES_2PDF with CUDA implementation
for scenario simulation in high efficiency.

Figure 5: The 3D terrain illustration in the ANSI-Platform: Simulation of the
2009 Hsiaolin event, where the flow thickness is shown by the color map and
the brown zones indicates the flow paths.

Figure 4ab are plotted by the illustration facility written by Python, and
Fig. 4c is constructed by the toolbox embedded in the QGIS. Figure 5 is one
snapshot of the animated 3D scenes in the ANSI-Platform, where the surface of
the topography is integrated with the satellite image. In the ANSI-Platform,
the dynamic distribution of flow body is depicted with color map, the
accumulated flow paths are given by brown shade, where the partly transparent
shade enables the structure or topographic features in the photographs to be
recognized.

###  4.2 The channel-type debris flow: DF004 event

Figure 6: The applied inflow condition for DF004 debris flow event

In the DF004 event, the debris flow rushed through the channelized valley and
attacked the Nansalu village after a long traveling distance around 17:00 on
Aug. 9th, 2009. Large amount of debris were entrained from the river bed
during the movement. For isolating the complex process of collecting the
debris during the traveling stage and improving the simulation efficiency, the
computational domain focuses on the encountered village, where appropriate
inflow conditions are estimated based on the records of the nearest rainfall
station several hours before the occurrence of the debris flow.

In the computation domain there are 27,720 meshes (including 3 ghost cells at
each boundaries), covering an area of $870\times 740$ m2. The mesh size is
$\Delta x=\Delta y=5$ m. The inflow condition consists of the inflow
locations, the associate fluxes and durations at the boundaries of the
computational domain. In general, the inflow location and fluxes (hydrograph)
are estimated by a 1D computation with empirical relation (based on the
rainfall records and the peak discharge with respect to a return period of 200
years), see Nakatani et al. (2008); Liu et al. (2013); Nakatani et al. (2016).
In the present study, the inflow condition is based on the computed results by
a rainfall run-off model and scaling the 3 hours rainfall record before the
event into 6 minutes (see the top panels of Fig. 6 and report by Soil Water
Conserv. Bureau (2020)), where the locations with high fluxes are set to be
the inflow points of debris flows, cf. the yellow section at the eastern
boundary of the computational domain in Fig. 7. Figure 6 shows the rainfall
records (3 hours before the event) and the estimated inflow flux at the 30
inflow meshes, which coincide with the bar area in the post-event picture (see
Fig. 7b). The scenario period is 12 mins, and it takes 202.8 s and 119.43 s
for the GPU-2080-Ti and GPU-Tesla-V100 computations, respectively.

a)
b)
c)

Figure 7: Simulation of DF004 event. (a,b) The flow thickness. (c) The solid
volume fraction. The area enclosed by the cyan lines is the domain in
computation.

Figure 7 illustrates the flow thickness (panels a and b) and the associated
distribution of volume fraction for the solid constituent (panel c) at $t=12$
mins in the computation. Here we recall that the areas with flow depth less
than 10 cm are isolated. Figure 7a is integrated with the satellite image
taken before the event, and panels b and c are integrated with post-event
images. Both of the images are provided by the Big Geospatial Information
System (BigGIS), a web-based database created by the Soil and Water
Conservation Bureau (SWCB), Taiwan. The brown line marks the plausible debris-
flow-endangered area delineated by the SWCB, Taiwan, using the Ikeya’s
empirical formula (Ikeya, 1981). It is found that the flow body approximately
coincides with the brown-line-marked area. In addition, the purple-line-
circulated area depicts the survey area post disaster (Soil Water Conserv.
Bureau, 2009), in which the volume of the deposited debris was estimated up to
70,000 m3 lying in the circulated area, while it is about 66,066 m3 in the
computation by MoSES_2PDF.

As shown in Fig. 7bc, the inflow locations in the computation coincide with
the bare area in the post-event image, which approximately mimics the
plausible flow paths from upstream during the event. In the computation, the
flow thickness may reach more than 10 cm in the area of village (see Fig.
7ab). Nevertheless, the transparency feature reveals that houses and buildings
in the northern part of the village was neither flushed away nor buried. One
of the possible reasons is that the flow thickness is thin, although it is
more than 10 cm. The other reason is the fact that the applied DEM does not
include the buildings, so that the impacts of the building cannot be
reproduced and less mass was collected in the purple-line-circulated area. In
Fig. 7c, it is found that the variation of the solid concentration
approximately indicates the topography, where the locations of low solid
concentration are in coincidence with the positions of high topography
gradient, especially with the contour of the river bank. The village is mainly
covered by the flow body of high solid concentration, which is speculated to
be caused by the relatively flat and gentle slope in the village area.

##  5 Concluding remarks

In the present study, we have presented a GPU-accelerated simulation tool,
MoSES_2PDF, for two-phase grain-fluid shallow flows, where the input and
output data can be integrated with the GIS system for subsequent application.
In comparison with the conventional CPU computation, the present GPU-CUDA
implementation can enhance the efficiency by up to 230 folds, depending on the
PC/workstations, GPU cards and the mesh numbers in the computation domain.
According to the types of initial condition, there are two modes by
MoSES_2PDF: Mode-I is for the landslide type and Mode-II is for debris flow
type. Two application examples are illustrated: the 2009 Hsiaolin landslide
event and the DF004 debris flow event. Each of them is corresponding to the
associated mode, and sound agreements are delievered with respect to the post-
event satellite images as well as field survey.

The MoSES_2PDF is able to provide the identical results either with the
RTX-2080-Ti or the Tesla-V100 GPU cards, although the results are is slightly
deviated from the ones computed by the CPU-code. This discrepancy is suspected
to be caused by the different structure in the source codes and some
float/double operation. In addition, the utilization of MoSES_2PDF is not
limited by the above-listed GPU-cards. In our tests, MoSES_2PDF can well run
with many other NVIDIA GPU cards, such as models GTX-1060, GTX-1080-Ti,
RTX-2070, RTX-3070, RTX-3090. The requirement is that the CUDA vers. 9.2 (or
higher version) should be installed.

As well, the MoSES_2PDF is also equipped with three output facilities for
illustrating the computed results. In addition to the asc-format by GIS and
conventional contour maps, the ANSI-Platform provides with the animated 3D
graphics, where the user can interact with 3D scenes for demonstrating or
investigating the computed results. Besides, taking benefits of the high
computational efficiency, MoSES_2PDF is able to perform abundant plausible
scenarios in short time. It may help us to figure out the best fitted or
optimal parameters by means of back-calculation of historical events for
hazard assessment or prediction before the event taking place. The features
and facilities of MoSES_2PDF make it suited in an efficient way for
engineering applications as well as for the hazard assessment or evaluation of
plausible disaster mitigation countermeasures.

##### Codeavailability

The codes for producing the results, the associated illustration facilities
and the numerical examples are available on request to th ecorresponding
author. It is remarked that, following the government regulation, only the
DEMs with resolution of $20\times 20$ m can be provided.

##### Acknowledgements

The financial support of the Ministry of Science and Technology, Taiwan (MOST
109-2221-E-006-022-) and the Soil and Water Conservation Bureau, Council of
Agriculture, Taiwan (SWCB-109-269) are acknowledged. In addition, the authors
thank the Taiwan Computing Cloud (TWCC) for the utilization of the GPU
computation (NVIDIA Tesla-V-100-SXM2-32GM).

## References

  * Pitman and Le [2005] E Bruce Pitman and Long Le.  A two-fluid model for avalanche and debris flows.  _Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 363(1832):1573–1601, 2005.
  * Li and Duffy [2011] Shuangcai Li and Christopher J Duffy.  Fully coupled approach to modeling shallow water flow, sediment transport, and bed evolution in rivers.  _Water Resources Research_ , 47(3), 2011.
  * Pudasaini [2012] Shiva P Pudasaini.  A general two-phase debris flow model.  _Journal of Geophysical Research: Earth Surface_ , 117(F3), 2012.
  * Pudasaini and Mergili [2019] Shiva P Pudasaini and Martin Mergili.  A multi-phase mass flow model.  _Journal of Geophysical Research: Earth Surface_ , 124(12):2920–2942, 2019.
  * Pudasaini and Hutter [2003] Shiva P Pudasaini and Kolumban Hutter.  Rapid shear flows of dry granular masses down curved and twisted channels.  _Journal of Fluid Mechanics_ , 495:193, 2003.
  * Bouchut and Westdickenberg [2004] F. Bouchut and M. Westdickenberg.  Gravity driven shallow water models for arbitrary topography.  _Comm. Math. Sci_ , 2:359–389, 2004.
  * Tai and Kuo [2008] Yih-Chin Tai and Ch-Y Kuo.  A new model of granular flows over general topography with erosion and deposition.  _Acta Mechanica_ , 199(1-4):71–96, 2008.
  * Tai et al. [2012] Yih-Chin Tai, Chih-Yu Kuo, and Wai-How Hui.  An alternative depth-integrated formulation for granular avalanches over temporally varying topography with small curvature.  _Geophysical & Astrophysical Fluid Dynamics_, 106(6):596–629, 2012.
  * Luca et al. [2016] Ioana Luca, Yih-Chin Tai, Chih-Yu Kuo, et al.  _Shallow geophysical mass flows down arbitrary topography_.  Springer, 2016.
  * Rauter and Tuković [2018] Matthias Rauter and Željko Tuković.  A finite area scheme for shallow granular flows on three-dimensional surfaces.  _Computers & Fluids_, 166:184–199, 2018.
  * Tai et al. [2021] Yih-Chin Tai, Jeaniffer Vides, Boniface Nkonga, and Chih-Yu Kuo.  Multi-mesh-scale approximation of thin geophysical mass flows on complex topographies.  _Communications in Computational Physics_ , 29(1):148–185, 2021.
  * Liu and Mei [1989] Ko Fei Liu and Chiang C Mei.  Slow spreading of a sheet of bingham fluid on an inclined plane.  _Journal of fluid mechanics_ , 207:505–529, 1989.
  * Coussot [1997] Philippe Coussot.  Mudflow rheology and dynamics.  _Balkema editions, Rotterdam, The Netherlands_ , 1997.
  * Huang and Garcia [1998] Xin Huang and Marcelo H Garcia.  A herschel–bulkley model for mud flow down a slope.  _Journal of fluid mechanics_ , 374:305–333, 1998.
  * Ancey [2007] Christophe Ancey.  Plasticity and geophysical flows: a review.  _Journal of Non-Newtonian Fluid Mechanics_ , 142(1-3):4–35, 2007.
  * Bagnold [1954] Ralph Alger Bagnold.  Experiments on a gravity-free dispersion of large solid spheres in a newtonian fluid under shear.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 225(1160):49–63, 1954\.
  * Takahashi [1978] Tamotsu Takahashi.  Mechanical characteristics of debris flow.  _Journal of the Hydraulics Division_ , 104(8):1153–1169, 1978.
  * O’Brien and Julien [1985] J.S. O’Brien and P.Y. Julien.  Physical properties and mechanics of hyperconcentrated sediment flows.  _Proc. ASCE HD Delineation of landslides, flash flood and debris flow Hazards_ , 1985.
  * O’Brien and Julien [1988] Jim S. O’Brien and Pierre Y. Julien.  Laboratory analysis of mudflow properties.  _Journal of hydraulic engineering_ , 114(8):877–887, 1988.
  * O’Brien et al. [1993] James S O’Brien, Pierre Y. Julien, and W.T. Fullerton.  Two-dimensional water flood and mudflow simulation.  _Journal of hydraulic engineering_ , 119(2):244–261, 1993.
  * Iverson [1997] Richard M Iverson.  The physics of debris flows.  _Reviews of geophysics_ , 35(3):245–296, 1997\.
  * Iverson and Denlinger [2001] Richard M Iverson and Roger P Denlinger.  Flow of variably fluidized granular masses across three-dimensional terrain: 1. coulomb mixture theory.  _Journal of Geophysical Research: Solid Earth_ , 106(B1):537–552, 2001.
  * Pudasaini et al. [2005] SP Pudasaini, Y Wang, and K Hutter.  Modelling debris flows down general channels.  _Natural Hazards and Earth System Sciences_ , 5(6):799–819, 2005.
  * Tai and Kuo [2012] Y.-C. Tai and C.-Y. Kuo.  Modelling shallow debris flows of the Coulomb-mixture type over temporally varying topography.  _Natural Hazards and Earth System Sciences_ , 12(2):269–280, 2012.
  * Egashira [1997] Shinji Egashira.  Constitutive equations of debris flow and their applicability.  In _1st. International Conference on Debris-Flow Hazards Mitigation, ASCE, 1997_ , pages 340–349, 1997.
  * Egashira [2007] Shinji Egashira.  Review of research related to sediment disaster mitigation.  _Journal of Disaster Research_ , 2(1):11–18, 2007\.
  * Nakatani et al. [2008] Kana Nakatani, Takashi Wada, Yoshifumi Satofuka, and Takahisa Mizuyama.  Development of “kanako 2d (ver. 2.00),” a user-friendly one-and two-dimensional debris flow simulator equipped with a graphical user interface.  _International Journal of Erosion Control Engineering_ , 1(2):62–72, 2008.
  * Liu et al. [2013] Jinfeng Liu, Kana Nakatani, and Takahisa Mizuyama.  Effect assessment of debris flow mitigation works based on numerical simulation by using Kanako 2d.  _Landslides_ , 10(2):161–173, 2013.
  * Nakatani et al. [2016] Kana Nakatani, Satoshi Hayami, Yoshifumi Satofuka, and Takahisa Mizuyama.  Case study of debris flow disaster scenario caused by torrential rain on kiyomizu-dera, kyoto, japan-using hyper kanako system.  _Journal of Mountain Science_ , 13(2):193–202, 2016.
  * Iverson [2009] Richard M Iverson.  Elements of an improved model of debris-flow motion.  In _AIP Conference Proceedings_ , volume 1145(1), pages 9–16. American Institute of Physics, 2009.
  * Iverson and George [2014] Richard M Iverson and David L George.  A depth-averaged debris-flow model that includes the effects of evolving dilatancy. i. physical basis.  _Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences_ , 470(2170):20130819, 2014.
  * Anderson and Jackson [1967] T Bo Anderson and Roy Jackson.  Fluid mechanical description of fluidized beds. equations of motion.  _Industrial & Engineering Chemistry Fundamentals_, 6(4):527–539, 1967.
  * Pelanti et al. [2008] Marica Pelanti, François Bouchut, and Anne Mangeney.  A roe-type scheme for two-phase shallow granular flows over variable topography.  _ESAIM: Mathematical Modelling and Numerical Analysis_ , 42(5):851–885, 2008.
  * Pailha and Pouliquen [2009] Mickael Pailha and Olivier Pouliquen.  A two-phase flow description of the initiation of underwater granular avalanches.  _Journal of Fluid Mechanics_ , 633:115, 2009.
  * Luca et al. [2012] I Luca, CY Kuo, K Hutter, and YC Tai.  Modeling shallow over-saturated mixtures on arbitrary rigid topography.  _Journal of Mechanics_ , 28(3):523–541, 2012\.
  * Bouchut et al. [2016] François Bouchut, Enrique Domingo Fernández Nieto, Anne Mangeney, and Gladys Narbona Reina.  A two-phase two-layer model for fluidized granular flows with dilatancy effects.  _Journal of Fluid Mechanics, 801, 166-221._ , 2016.
  * Meng et al. [2017] Xiannan Meng, Yongqi Wang, Chun Wang, and Jan-Thomas Fischer.  Modeling of unsaturated granular flows by a two-layer approach.  _Acta Geotechnica_ , 12(3):677–701, 2017.
  * Heß et al. [2017] Julian Heß, Yongqi Wang, and Kolumban Hutter.  Thermodynamically consistent modeling of granular-fluid mixtures incorporating pore pressure evolution and hypoplastic behavior.  _Continuum mechanics and thermodynamics_ , 29(1):311–343, 2017.
  * Heß et al. [2019] Julian Heß, Yih-Chin Tai, and Yongqi Wang.  Debris flows with pore pressure and intergranular friction on rugged topography.  _Computers & Fluids_, 190:139–155, 2019.
  * Mergili et al. [2017] Martin Mergili, Jan-Thomas Fischer, Julia Krenn, and Shiva P Pudasaini.  r.avaflow v1, an advanced open-source computational framework for the propagation and interaction of two-phase mass flows.  _Geoscientific Model Development_ , 10(2):553–569, 2017.
  * Mergili [2014-2020] M. Mergili.  r.avaflow – the mass flow simulation tool. r.avaflow 2.3 user manual.  https://www.avaflow.org/manual.php, 2014-2020.
  * Mergili et al. [2020] Martin Mergili, Michel Jaboyedoff, José Pullarello, and Shiva P Pudasaini.  Back calculation of the 2017 piz cengalo–bondo landslide cascade with r. avaflow: what we can do and what we can learn.  _Natural Hazards and Earth System Sciences_ , 20(2):505–520, 2020.
  * Baggio et al. [2021] Tommaso Baggio, Martin Mergili, and Vincenzo D’Agostino.  Advances in the simulation of debris flow erosion: The case study of the rio gere (italy) event of the 4th august 2017.  _Geomorphology_ , 381:107664, 2021.
  * Castro et al. [2011] Manuel J Castro, Sergio Ortega, Marc De la Asuncion, José M Mantas, and José M Gallardo.  Gpu computing for shallow water flow simulation based on finite volume schemes.  _Comptes Rendus Mécanique_ , 339(2-3):165–184, 2011.
  * Brodtkorb et al. [2012] André R Brodtkorb, Martin L Sætra, and Mustafa Altinakar.  Efficient shallow water simulations on gpus: Implementation, visualization, verification, and validation.  _Computers & Fluids_, 55:1–12, 2012.
  * de la Asunción et al. [2013] Marc de la Asunción, Manuel J Castro, Enrique Domingo Fernández-Nieto, José M Mantas, Sergio Ortega Acosta, and José Manuel González-Vida.  Efficient gpu implementation of a two waves tvd-waf method for the two-dimensional one layer shallow water system on structured meshes.  _Computers & Fluids_, 80:441–452, 2013.
  * Aureli et al. [2020] Francesca Aureli, Federico Prost, Renato Vacondio, Susanna Dazzi, and Alessia Ferrari.  A gpu-accelerated shallow-water scheme for surface runoff simulations.  _Water_ , 12(3):637, 2020.
  * Dazzi et al. [2019] S Dazzi, R Vacondio, and P Mignosa.  Integration of a levee breach erosion model in a gpu-accelerated 2d shallow water equations code.  _Water Resources Research_ , 55(1):682–702, 2019\.
  * Wu et al. [2021] Xinhui Wu, Ethan J Kubatko, and Jesse Chan.  High-order entropy stable discontinuous galerkin methods for the shallow water equations: curved triangular meshes and gpu acceleration.  _Computers & Mathematics with Applications_, 82:179–199, 2021.
  * Tai et al. [2019] Yih-Chin Tai, Julian Heß, and Yongqi Wang.  Modeling two-phase debris flows with grain-fluid separation over rugged topography: Application to the 2009 hsiaolin event, taiwan.  _Journal of Geophysical Research: Earth Surface_ , 124(2):305–333, 2019.
  * Kurganov and Tadmor [2000] Alexander Kurganov and Eitan Tadmor.  New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations.  _Journal of Computational Physics_ , 160(1):241–282, 2000.
  * Kurganov et al. [2007] Alexander Kurganov, Guergana Petrova, et al.  A second-order well-balanced positivity preserving central-upwind scheme for the saint-venant system.  _Communications in Mathematical Sciences_ , 5(1):133–160, 2007.
  * Harris et al. [2007] Mark Harris et al.  Optimizing parallel reduction in cuda.  _Nvidia developer technology_ , 2(4):70, 2007\.
  * Gray et al. [1999] J.M.N.T. Gray, M. Wieland, and K. Hutter.  Gravity-driven free surface flow of granular avalanches over complex basal topography.  _Proc. R. Soc. Lond. A_ , 455:1841–1875, 1999.
  * Kuo et al. [2009] C.Y. Kuo, Y.C. Tai, F. Bouchut, A. Mangeney, M. Pelanti, R.F. Chen, and K.J. Chang.  Simulation of Tsaoling landslide, Taiwan, based on Saint Venant equations over general topography.  _Engineering Geology_ , 104(3-4):181 – 189, 2009\.
  * de’Michieli Vitturi et al. [2019] Mattia de’Michieli Vitturi, Tomaso Esposti Ongaro, Giacomo Lari, and Alvaro Aravena.  Imex_sflow2d 1.0: a depth-averaged numerical flow model for pyroclastic avalanches.  _Geoscientific Model Development_ , 12(1):581–595, 2019.
  * Ryoo et al. [2007] Shane Ryoo, Christopher Rodrigues, Sam Stone, Sara Baghsorkhi, Sain-Zee Ueng, and Wen-mei W Hwu.  Program optimization study on a 128-core gpu.  In _The First Workshop on General Purpose Processing on Graphics Processing Units_ , pages 30–39. Citeseer, 2007.
  * Dong et al. [2011] Jia-Jyun Dong, Yun-Shan Li, Chyh-Yu Kuo, Rui-Tang Sung, Ming-Hsu Li, Chyi-Tyi Lee, Chien-Chih Chen, and Wang-Ru Lee.  The formation and breach of a short-lived landslide dam at Hsiaolin village, Taiwan—part I: Post-event reconstruction of dam geometry.  _Engineering geology_ , 123(1):40–59, 2011.
  * Soil Water Conserv. Bureau [2009] Soil Water Conserv. Bureau.  Disaster Fast Report: 2009 typhoon Morakot – Namasia township 003 (in Chinese).  Technical report, Soil Water Conserv. Bureau (SWCB), Taiwan, 2009.
  * Tai et al. [2020] Yih-Chin Tai, Chi-Jyun Ko, Kun-Ding Li, Yu-Chen Wu, Chih-Yu Kuo, Rou-Fei Chen, and Ching-Weei Lin.  An idealized landslide failure surface and its impacts on the traveling paths.  _Frontiers in Earth Science_ , 8:313, 2020.
  * Kuo et al. [2011] CY Kuo, Yih-Chin Tai, CC Chen, KJ Chang, AY Siau, JJ Dong, RH Han, T Shimamoto, and CT Lee.  The landslide stage of the hsiaolin catastrophe: simulation and validation.  _Journal of Geophysical Research: Earth Surface_ , 116(F4), 2011.
  * Soil Water Conserv. Bureau [2020] Soil Water Conserv. Bureau.  GIS integrated web-based system for scenario investigation on slope failure and sediment-related disasters (in Chinese).  Technical report, SWCB-109-269, Taiwan, 2020.
  * Ikeya [1981] Hiroshi Ikeya.  A method of designation for area in danger of debris flow.  _Erosion and sediment transport in Pacific Rim Steeplands_ , 132:578–588, 1981.

##  Appendix A Model equations

In MoSES_2PDF, the employed model equations are adopted form the ones in Tai
et al. [2019]. With the shallowness assumption for the flow body, the depth-
integrated dimensionless equations of leading-order read

$\displaystyle\frac{\partial}{\partial
t}\left(J_{b}h^{s}\right)+\frac{\partial}{\partial\xi}\left(J_{b}h^{s}v_{\xi}^{s}\right)+\frac{\partial}{\partial\eta}\left(J_{b}h^{s}v_{\eta}^{s}\right)=0$
(1) $\displaystyle\frac{\partial}{\partial
t}\left(J_{b}h^{f}\right)+\frac{\partial}{\partial\xi}\left(J_{b}h^{f}v_{\xi}^{f}\right)+\frac{\partial}{\partial\eta}\left(J_{b}h^{f}v_{\eta}^{f}\right)=0$
(2)

for the mass conservation, where $J_{b}=\det{\boldsymbol{\Omega}}_{b}$ is the
Jacobian determinant on the basal surface for coordinate transformation,
$v_{\xi,\eta}^{s,f}$ denote the tangential components of the depth-averaged
velocities in $O_{\xi\eta\zeta}$, respectively. In (1) and (2),
$h^{s}=h\phi^{s}$ and $h^{f}=h\phi^{f}$ with $\phi^{s,f}$ the depth-averaged
volume fraction, where $h$ stands for the flow thickness measured normal to
the basal topographic surface.

Regarding the momentum balance equations for the solid and fluid constituents,
the evolutions of
${\mathbf{q}}^{s}=\left(J_{b}h^{s}v^{s}_{X},J_{b}h^{s}v^{s}_{Y}\right)^{T}$
and
${\mathbf{q}}^{f}=\left(J_{b}h^{f}v^{f}_{X},J_{b}h^{f}v^{f}_{Y}\right)^{T}$
are determined by

$\begin{array}[]{l}\displaystyle\frac{\partial{\mathbf{q}}^{s}}{\partial
t}+\frac{\partial{\mathbf{F}}^{s}}{\partial\xi}+\frac{\partial{\mathbf{G}}^{s}}{\partial\eta}=\displaystyle{\mathbf{s}}_{n}^{s}+{\mathbf{s}}_{f}^{s}+{\mathbf{s}}_{d}^{s}+{\mathbf{s}}_{v}^{s}\end{array}$
(3)

and

$\begin{array}[]{l}\displaystyle\frac{\partial{\mathbf{q}}^{f}}{\partial
t}+\frac{\partial{\mathbf{F}}^{f}}{\partial\xi}+\frac{\partial{\mathbf{G}}^{f}}{\partial\eta}=\displaystyle{\mathbf{s}}_{n}^{f}+{\mathbf{s}}_{f}^{f}+{\mathbf{s}}_{d}^{f}+{\mathbf{s}}_{vis}^{f}\,,\end{array}$
(4)

respectively, where ${\mathbf{F}}^{s,f}$ and ${\mathbf{G}}^{s,f}$ stand for
the fluxes,

${\mathbf{F}}^{s}=\left(\begin{array}[]{l}\displaystyle
J_{b}h^{s}v^{s}_{X}v_{\xi}^{s}+\epsilon
J_{b}h{A}_{11}\overline{N}^{s}\\\\[5.69054pt]
J_{b}h^{s}v^{s}_{Y}v_{\xi}^{s}+\epsilon
J_{b}h{A}_{12}\overline{N}^{s}\end{array}\right),$ (5)
${\mathbf{G}}^{s}=\left(\begin{array}[]{l}\displaystyle
J_{b}h^{s}v^{s}_{X}v_{\eta}^{s}+\epsilon
J_{b}h{A}_{21}\overline{N}^{s}\\\\[5.69054pt]
J_{b}h^{s}v^{s}_{Y}v_{\eta}^{s}+\epsilon
J_{b}h{A}_{22}\overline{N}^{s}\end{array}\right),$ (6)
${\mathbf{F}}^{f}=\left(\begin{array}[]{l}\displaystyle
J_{b}h^{f}v^{f}_{X}v_{\xi}^{f}+\epsilon
J_{b}h{A}_{11}\overline{p}^{f}\\\\[5.69054pt]
J_{b}h^{f}v^{f}_{Y}v_{\xi}^{f}+\epsilon
J_{b}h{A}_{12}\overline{p}^{f}\end{array}\right),$ (7)
${\mathbf{G}}^{f}=\left(\begin{array}[]{l}\displaystyle
J_{b}h^{f}v^{f}_{X}v_{\eta}^{f}+\epsilon
J_{b}h{A}_{21}\overline{p}^{f}\\\\[5.69054pt]
J_{b}h^{f}v^{f}_{Y}v_{\eta}^{f}+\epsilon
J_{b}h{A}_{22}\overline{p}^{f}\end{array}\right).$ (8)

Here, we recall that the unit normal vector ${\boldsymbol{n}}$ at the basal
surface reads ${\boldsymbol{n}}=(n_{X},n_{Y},n_{Z})^{T}$ in $O_{XYZ}$, and
$v_{X,Y}^{s,f}$ are the $X$\- and $Y$-components of the phase velocity
projected on the horizontal plane. In Eqs. (5) to (8),
$({A}_{ij})={\boldsymbol{\Omega}}^{-1}_{b}$ is the inverse of the
transformation matrix. $\epsilon$ is the aspect ratio, the ratio of the flow
characteristic thickness to the characteristic length along the basal surface.
$\overline{N}^{s}=c(1-\alpha_{\rho}){h^{s}}/{2}$ stands for the depth-averaged
solid pressure, where $c=n_{Z}$ is the $Z$-component of the unit normal vector
of the basal surface and $\alpha_{\rho}=\rho^{f}/\rho^{s}$ denotes the density
ratio of the fluid constituent to the solid one. Besides,
$\overline{p}^{f}=c\,{h}/{2}$ stands for the mean fluid pressure.

The right hand side of (3) consists of four source terms and they are

${\mathbf{s}}_{n}^{s}=J_{b}p^{s}_{b}\left(\begin{array}[]{l}\\!n_{X}\\!\\!\\\\[5.69054pt]
\\!n_{Y}\\!\\!\end{array}\\!\\!\right),\;\;{\mathbf{s}}_{d}^{s}=-J_{b}p^{s}_{b}\tan\delta_{b}\left(\begin{array}[]{l}\displaystyle\frac{v_{X}^{s}}{||{\boldsymbol{v}}^{s}||}\\\\[5.69054pt]
\displaystyle\frac{v_{Y}^{s}}{||{\boldsymbol{v}}^{s}||}\end{array}\\!\\!\right),$
(9)
${\mathbf{s}}_{f}^{s}=-\epsilon\alpha_{\rho}\phi^{s}\left(\\!\\!\\!\begin{array}[]{l}\displaystyle{A}_{11}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\xi}+{A}_{21}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\eta}\\\\[5.69054pt]
\displaystyle{A}_{12}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\xi}+{A}_{22}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\eta}\end{array}\\!\\!\\!\right),$
(10)
${\mathbf{s}}_{v}^{s}=J_{b}\alpha_{\rho}c_{D}\phi^{s}\phi^{f}h\left(\begin{array}[]{l}\displaystyle
v_{X}^{f}-v_{X}^{s}\\\\[5.69054pt] v_{Y}^{f}-v_{Y}^{s}\end{array}\right),$
(11)

where $p^{s}_{b}$ means the solid basal pressure, $\delta_{b}$ represents the
angle of basal friction and $c_{D}$ the drag coefficient for the drag force
induced by the relative velocity between the constituents. In (9), the solid
basal pressure reads
$p^{s}_{b}=h^{s}\left(c\left(1-\alpha_{\rho}\right)-\epsilon^{\chi}\kappa^{s}\right)$,
where $\kappa^{s}$ means the centripetal acceleration due to the local
velocity of the solid constituent and curvature of the basal surface, and it
is given by [e.g. Tai et al., 2012, 2019]

$\begin{array}[]{l}\displaystyle\kappa^{s}=\phantom{+}v^{s}_{X}\frac{\partial
n_{X}}{\partial\xi}v^{s}_{\xi}+v^{s}_{Y}\frac{\partial
n_{Y}}{\partial\xi}v^{s}_{\xi}+v^{s}_{Z}\frac{\partial
n_{Z}}{\partial\xi}v^{s}_{\xi}\\\
\displaystyle\phantom{\kappa^{s}=}+v^{s}_{X}\frac{\partial
n_{X}}{\partial\eta}v^{s}_{\eta}+v^{s}_{Y}\frac{\partial
n_{Y}}{\partial\eta}v^{s}_{\eta}+v^{s}_{Z}\frac{\partial
n_{Z}}{\partial\eta}v^{s}_{\eta}\,.\end{array}$

The source terms of (4) read

${\mathbf{s}}_{n}^{f}=J_{b}p^{f}_{b}\left(\begin{array}[]{l}\\!n_{X}\\!\\!\\\\[5.69054pt]
\\!n_{Y}\\!\\!\end{array}\\!\\!\right),\;{\mathbf{s}}_{d}^{f}=-\frac{J_{b}\phi^{f}h\vartheta_{b}^{f}}{\epsilon
N_{R}}\left(\begin{array}[]{l}\displaystyle{v_{X}^{f}}\\\\[5.69054pt]
\displaystyle{v_{Y}^{f}}\end{array}\\!\\!\right)$ (12)
${\mathbf{s}}_{f}^{f}=\epsilon\phi^{f}\left(\begin{array}[]{l}\displaystyle{A}_{11}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\xi}+{A}_{21}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\eta}\\\\[5.69054pt]
\displaystyle{A}_{12}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\xi}+{A}_{22}\frac{\partial\left(J_{b}h\overline{p}^{f}\right)}{\partial\eta}\end{array}\right)$
(13)
${\mathbf{s}}_{v}^{f}=-J_{b}c_{D}\phi^{s}\phi^{f}h\left(\begin{array}[]{l}\displaystyle
v_{X}^{f}-v_{X}^{s}\\\\[5.69054pt] v_{Y}^{f}-v_{Y}^{s}\end{array}\right)$ (14)

${\mathbf{s}}_{vis}^{f}=\frac{\epsilon\phi^{f}}{N_{R}}\left(\begin{array}[]{l}s_{vis,X}^{f}\\\\[5.69054pt]
s_{vis,Y}^{f}\end{array}\right)$ (15)

with

$\begin{array}[]{l}s_{vis,X}^{f}=\displaystyle
2\frac{\partial}{\partial\xi}\Bigl{[}J_{b}h\left(A_{11}\partial_{\xi}{v_{\xi}^{f}}+A_{21}\partial_{\eta}{v_{\xi}^{f}}\right)\Bigr{]}\\\\[5.69054pt]
\displaystyle+\frac{\partial}{\partial\eta}\Bigl{[}J_{b}h\left(A_{12}\partial_{\xi}{v_{\xi}^{f}}+A_{22}\partial_{\eta}{v_{\xi}^{f}}+A_{11}\partial_{\xi}{v_{\eta}^{f}}+A_{21}\partial_{\eta}{v_{\eta}^{f}}\right)\Bigr{]}\\\\[11.38109pt]
s_{vis,Y}^{f}=\displaystyle
2\frac{\partial}{\partial\eta}\Big{[}J_{b}h\left(A_{12}\partial_{\xi}{v_{\eta}^{f}}+A_{22}\partial_{\eta}{v_{\eta}^{f}}\right)\Bigr{]}\\\\[5.69054pt]
\displaystyle+\frac{\partial}{\partial\xi}\Bigl{[}J_{b}h\left(A_{12}\partial_{\xi}{v_{\xi}^{f}}+A_{22}\partial_{\eta}{v_{\xi}^{f}}+A_{11}\partial_{\xi}{v_{\eta}^{f}}+A_{21}\partial_{\eta}{v_{\eta}^{f}}\right)\Bigr{]}\end{array}$

for both directions. In (12), $p^{f}_{b}$ stands for the fluid basal pressure
and $\vartheta^{f}_{b}$ represents the fluid friction coefficient. The fluid
basal pressure is determined by $p^{f}_{b}=h^{f}(c-\epsilon^{\chi}\kappa^{f})$
with

$\begin{array}[]{l}\displaystyle\kappa^{f}=\phantom{+}v^{f}_{X}\frac{\partial
n_{X}}{\partial\xi}v^{f}_{\xi}+v^{f}_{Y}\frac{\partial
n_{Y}}{\partial\xi}v^{f}_{\xi}+v^{f}_{Z}\frac{\partial
n_{Z}}{\partial\xi}v^{f}_{\xi}\\\\[5.69054pt]
\displaystyle\phantom{\kappa^{f}=}+v^{f}_{X}\frac{\partial
n_{X}}{\partial\eta}v^{f}_{\eta}+v^{f}_{Y}\frac{\partial
n_{Y}}{\partial\eta}v^{f}_{\eta}+v^{f}_{Z}\frac{\partial
n_{Z}}{\partial\eta}v^{f}_{\eta},\end{array}$

the centripetal acceleration with respect to the fluid velocity and the local
curvature of the topographic surface.


# Confinement and collective escape of active particles

Igor S Aranson  Departments of Biomedical Engineering, Chemistry, and
Mathematics, Penn State University, University Park, Pennsylvania 16802, USA
Arkady Pikovsky  Institute for Physics and Astronomy, University of Potsdam,
Karl-Liebknecht-Strasse 24/25, 14476 Potsdam-Golm, Germany  Department of
Control Theory, Nizhny Novgorod State University, Gagarin Avenue 23, 606950
Nizhny Novgorod, Russia

###### Abstract

Active matter broadly covers the dynamics of self-propelled particles. While
the onset of collective behavior in homogenous active systems is relatively
well understood, the effect of inhomogeneities such as obstacles and traps
lacks overall clarity. Here, we study how interacting self-propelled particles
become trapped and released from a trap. We have found that captured particles
aggregate into an orbiting condensate with a crystalline structure. As more
particles are added, the trapped condensate escape as a whole. Our results
shed light on the effects of confinement and quenched disorder in active
matter.

Assemblies of interacting self-propelled particles, broadly defined as active
matter, continue attracting significant attention Bechinger  _et al._ (2016);
Aranson (2013); Gompper  _et al._ (2020). In the last ten years, notable
progress was achieved in the understanding of the onset of collective behavior
Chaté (2020); Peshkov  _et al._ (2012); Patelli  _et al._ (2019) and
characterization of some collective states Solon  _et al._ (2015); Alert  _et
al._ (2020). Spatial inhomogeneities, surface roughness, or quenched disorder
play a significant role in active systems Peruani and Aranson (2018); Duan
_et al._ (2021); Ro  _et al._ (2021); Olsen  _et al._ (2021). Quenched
disorder, for example, may lead to the onset of trapped states, anomalous
diffusion, and breakdown of ergodicity.

The motion of self-propelled particles on a disordered substrate is relevant
in the context of the “active conductivity”. This situation is realized, for
example, when motile bacteria propagate through a porous environment Waisbord
_et al._ (2021). It is an analog of the equilibrium problem of electrons
migration in random media. If the substrate is approximated as an array of
well-separated traps (impurities), then the process can be viewed as a
sequence of escapes and re-trappings of particles. In the context of electrons
trapped by impurities, a non-zero conductance is due to the overlap of wave-
functions of electrons at individual impurities (Dykhne theorem Bychkov and
Dykhne (1971)). However, no such general result is known for active particles.

In this Letter, we study how self-propelled particles are captured and
released by an isolated trap modeled as a potential well. We start with non-
interacting particles and show that an individual particle typically exhibits
chaotic scattering by a trap. Then we introduce interactions. A Lennard-Jones
potential organizes particles into an orbiting condensate with the hexatic
crystalline order. An alignment coupling brings dissipation and
synchronization in the dynamics, and results in a perpetual capture of active
particles: a trap becomes an analog of a “black hole”. “Bombardment” by active
particles results in particle absorption, condensate melting, and
recrystallization. Above a certain threshold number of captured particles, the
trap storing capacity is exceeded, and the condensate escapes as a whole.

Active particles in a harmonic trap have been studied both experimentally
Dauchot and Démery (2019); Schmidt  _et al._ (2021); Takatori  _et al._ (2016)
and theoretically Pototsky and Stark (2012); Hennes  _et al._ (2014); Wexler
_et al._ (2020); Das  _et al._ (2018); Jahanshahi  _et al._ (2017). The main
focus was on the steady-state distributions or escape of individual particles
due to the combined effect of self-propulsion and thermal fluctuations. In
this Letter, we consider purely deterministic effects of propulsion and
interactions; this aspect of our model manifests a crucial difference not
explored in other publications.

We consider a self-propelled particle moving in two dimensions with a constant
velocity $V$, while a force $\mathbf{f}$ acting on the particle only rotates
its direction of motion:

$\frac{d{\bf r}}{dt}=\tilde{V}{\bf n}\;,\qquad\frac{d{\bf
n}}{dt}=\mathbf{f}-({\bf n}\cdot\mathbf{f}){\bf n}\;.$ (1)

Here ${\bf n}=(\cos\theta,\sin\theta)$ is the unit vector in the direction of
motion. The equation for $\bf n$ ensures that $|{\bf n}|^{2}=1$. Thus, only
the velocity direction changes and obeys $\dot{\theta}=f_{y}\cos\theta-
f_{x}\sin\theta$. Physically, velocity aligns with the potential gradient.
This situation can be realized, for example, for magnetic particles in a
magnetic trap. Models of this type were also discussed in the context of
artificial chemotaxis Liebchen  _et al._ (2016); Liebchen and Löwen (2018);
Stark (2018).

Below we consider different types of forces, but we start with a motion in an
external potential field ${U}(x,y)$, so that $\mathbf{f}_{\rm pot}=-\nabla U$.
If we renormalize variables so that the width and the depth of the potential
are one, the only remaining relevant parameter is the dimensionless velocity
$V$. The limit of strong potential is that of $V\to 0$. The corresponding
equations can be written as

$\displaystyle\frac{d{\bf r}}{dt}$ $\displaystyle=$ $\displaystyle V{\bf
n}\;,$ (2) $\displaystyle\frac{d\theta}{dt}$ $\displaystyle=$
$\displaystyle-\partial_{y}U\cos\theta+\partial_{x}U\sin\theta\;.$ (3)

Equations (2), (3) possess the Hamiltonian

$H({\bf r},{\bf p})=V|{\bf p}|-\exp\left[-\frac{U({\bf r})}{V}\right]=0\;.$
(4)

Because of the energy conservation, the relation $|{\bf p}|=V^{-1}\exp[-U({\bf
r})V^{-1}]$ holds. A substitution ${\bf p}=V^{-1}\exp[-U({\bf r})V^{-1}]{\bf
n}$ reduces the Hamiltonian dynamics to Eqs. (2,3). The Hamiltonian (4)
coincides with that describing ray propagation in geometrical optics, with the
refraction index $\sim\exp[-U({\bf r})V^{-1}]$ Kravtsov and Orlov (1990).

In the small velocity limit, $V\ll 1$, one can separate slow migration of the
particle (Eq. (2)) and its fast alignment with the gradient of the potential
(Eq. (3)). The orientation angle $\theta$ in Eq. (3) adjusts to the gradient
angle $\alpha=\arctan(\partial_{y}U/\partial_{x}U)$ as $\dot{\theta}=|\nabla
U|\sin(\alpha-\theta)$. This fast adjustment of the direction of motion toward
the minimum of the potential is followed by a slow drift (2). Once the
vicinity of the potential minimum is reached, i.e., $|\nabla U|\approx 0$, the
scale separation breaks down, and one has to consider the full Eqs. (2),(3).

In the vicinity of a minimum, a generic potential is approximated by an
asymmetric parabolic well, $U_{\rm pb}=(x^{2}+b^{2}y^{2})/2$. Since the depth
of this harmonic potential is not defined, one can set by virtue of a
renormalization, parameter $V$ in Eqs. (2),(3) to one. Then the only parameter
left is the potential asymmetry $b$. Since the dynamics is Hamiltonian, the
type of motion depends on initial conditions.

Figure 1:  Poincaré map for a parabolic potential $b^{2}=\sqrt{5}-1$ (color
background). Black dots: a chaotic trajectory starting at large $|x|,|y|$. The
chaotic dynamics is not ergodic as it contains voids with closed curves -
images of quasiperiodic trapped motion. Inset: Trapping time $T_{\rm
tr}\sim\exp(1/V)$ in a Gaussian potential well vs velocity $V$,
$b=(\sqrt{5}-1)/2$.

Figure 1 illustrates the dynamics for $b^{2}=\sqrt{5}-1$ (other values of $b$
give qualitatively the same picture). The three-dimensional phase space is
reduced to a two-dimensional Poncaré map (we chose the section $\cos\theta=0$,
$\frac{d}{dt}(\cos\theta)<0$). One sees a large chaotic domain (black dots)
(video 1 in sup ) and regular islands filled with quasiperiodic trajectories.
This Poincaré map is unusual for the Hamiltonian dynamics: the distribution is
very inhomogeneous – density of points at small $\mathbf{r}$ values is much
larger than at large ones. That happens because the “natural” coordinates
$({\bf r},\theta)$ are not the canonical ones. One estimates the density on
the plane $(x,y)$ by assuming a fully developed chaos where the angle $\theta$
is random and uniformly distributed. Then integration of a microcanonical
distribution density for the Hamiltonian (4) $w({\bf r},{\bf
p})\sim\delta\left(V|{\bf p}|-\exp\left[-\frac{U({\bf r})}{V}\right]\right)$
over the angle $\theta$ yields $W({\bf r})\sim\exp\left[-U({\bf
r})/{V}\right]$. This resembles the Gibbs-Boltzmann distribution, with the
velocity $V$ playing a role of the temperature. This result implies that the
chaotic motion in Fig. 1 spreads to arbitrarily large values of the potential,
although it is rather improbable to reach these heights. We conclude that
although a slow particle arrives at the minimum of the potential and moves
chaotically there, after a long time, it returns to the high values of the
potential where it started. Such a return must happen according to the
recurrence of the Hamiltonian trajectories. In addition to a chaotic region in
Fig. 1, there are domains of quasi-periodic dynamics concentrated close to the
potential minimum. This motion can occur for particles starting close to the
minimum.

The Hamiltonian structure of Eqs. (2),(3) implies that capture of particles
falling in a finite-depth potential well is impossible. Only temporary
trapping occurs that can be interpreted as a chaotic scattering: a particle
falls into the well, goes to its minimum, and moves there chaotically like in
Fig. 1, but eventually rises high and escapes (video 2 in sup ). Furthermore,
because of the exponential density $W\sim\exp[-U/V]$, the characteristic
trapping time obeys a Kramers-like law $T_{\rm tr}\sim\exp[-1/V]$. This is
confirmed in Fig. 1(inset) where the mean trapping time on a Gaussian
potential well $U(x,y)=-\exp(-x^{2}-b^{2}y^{2})$ is shown.

We consider two types of interactions between particles below: (i) interaction
by a potential force; (ii) an alignment to an average (over a neighborhood)
orientation of neighboring particles; this latter interaction is of the Vicsek
(or Kuramoto) type Vicsek  _et al._ (1995); Chaté (2020).

Figure 2: Trajectories of three coupled particles in $(x,y)$ plane, for the LJ
potential, for $V=0.1$ (a) and $V=0.2$ (b). Circles show positions of
particles at time $t=100$, squares at time $t=500$. In panel (b) one sees the
breakdown of the bounded state into a couple of particles moving to the left,
and a particle moving to the right. (c) Average lifetime of a crystal of $N$
particles vs inverse velocity $1/V$. The dissociation time is defined as an
event when one particle leaves.

We start with the interaction via a conservative pair-wise force that depends
on the distance between the particles. This case, like the motion in the
external potential above, is Hamiltonian. We explore below the Lennard-Jones
(LJ) potential $U_{\rm LJ}({\bf r}_{1},{\bf r}_{2})=4(\rho^{-12}-\rho^{-6})$,
where $\rho^{2}=|{\bf r}_{1}-{\bf r}_{2}|^{2}$. Several particles having a
small velocity $V$ placed close to each other, form a bounded state due to the
LJ-coupling. For two particles, the bounded state is quasiperiodic, while for
$N\geq 3$ it is typically chaotic. Several particles form a chaotically
vibrating “crystalline molecule” with a hexatic order. The center of mass
(c.o.m.) performs a diffusive motion in the plane $(x,y)$, and the particles
from time to time rearrange their positions, see Fig. 2(a,b) and Video 4 in
sup . Because the LJ potential is finite, there is a non-zero probability for
a particle to escape. Such an event is shown in Fig. 2(b). Figure 2(c) shows
that the lifetime $T_{\rm LJ}$ of a “shaking crystal” without confining
potential has the same Kramers-type dependence on the velocity $T_{\rm
LJ}\propto\exp[V^{-1}]$, as the lifetime of a particle in a potential well
(cf. Fig. 1(inset)). A crystalline molecule has a significant lifetime only
for particles with $V\lesssim 0.5$.

The aligning force $\mathbf{F}_{k}$ acts in the direction of the weighted
average of the velocities of other particles in a neighborhood of $k$-th
particle:

$\mathbf{F}_{k}=\epsilon\sum_{j}S(|\mathbf{r}_{j}-\mathbf{r}_{k}|)\mathbf{n}_{j}\;.$
(5)

The distance-dependent factor (we assume it to be Gaussian
$S(r)=\exp(-r^{2}/r_{0}^{2})$) defines the range $r_{0}$ of the force.
Parameter $\epsilon$ determines the strengths of the alignment. The alignment
force is velocity-dependent and dissipative. In terms of the velocity
direction $\theta$ it has the form of Kuramoto-type coupling
$\dot{\theta}_{k}\sim\sin(\theta_{j}-\theta_{k})$.

We examine next a combination of the alignment and the conservative forces due
to a confining potential or an LJ interaction (for a pure alignment see Kruk
_et al._ (2018)).

We start with a set of chaotic particles in a harmonic trap (Fig. 1),
described by Eqs. (2),(3) with additional alignment (5). The main observation
is that for large values of $\epsilon$ and large ranges of coupling $r_{0}$,
particles always synchronize: after some transient time, all the coordinates
and angles coincide, and the particles form a synchronous point cluster ${\bf
r}_{1}=\ldots={\bf r}_{N}$, $\theta_{1}=\ldots=\theta_{N}$ (video 3 in sup ).
This synchronization is possible because the aligning force is dissipative
Peruani and Aranson (2018). In the final synchronous state dissipation
disappears (the force $\sim\sin(\theta_{j}-\theta_{k})$ vanishes), and the
trajectory of the cluster is described again by the Hamiltonian dynamics
(2),(3). However, in the course of alignment, particles leave the chaotic
domain, and the final dynamics is quasiperiodic (cf. Fig. 1). Thus, strong
alignment synchronizes particles and regularizes their motion. For lower
alignment rates $\epsilon$ and especially for small ranges $r_{0}$, multiple
states with several clusters are observed up to large times. If the alignment
coupling radius $r_{0}$ is small enough, several regular clusters may
effectively stop to interact; then they constitute an “attractor”. Noteworthy,
standard methods of the synchronization theory, like the master stability
function method Pecora and Carroll (1998), are not applicable here because the
type of motion (from chaos to quasiperiodicity) changes over time. In Fig.
3(a) we illustrate the rate of synchronization in dependence on the parameters
$\epsilon$ and $r_{0}$. At large $r_{0}$ there is an “optimal” coupling
strength. We attribute this to the fact that, for large $r_{0}$, the alignment
coupling is a global one. Because the velocities of distant particles are
effectively “de-aligned” by different potential forces they experience at
different positions in the harmonic trap, the alignment slows down. At small
values of $r_{0}$, only neighboring particles interact, and they are much more
“synchronizable” because their trajectories in the potential may easily adjust
as well.

Figure 3:  (a) The average time at which 20 initially randomly distributed
particles in a harmonic potential form 10 clusters due to alignment. For small
$r_{0}$ this time monotonously decreases with coupling strength $\epsilon$,
while large ranges of $r_{0}\gtrsim 0.5$ there is an optimal coupling
strength. (b) Average time to synchronization for chaotic crystals composed of
$N$ particles, for $V=0.1$ and different strengths of alignment forces
$\epsilon$.

In the case of alignment of particles trapped by a finite-depth potential
well, two time scales compete: the trapping time $T_{\rm tr}$ determined by
the velocity $V$ (cf. Fig. 1(inset)); the synchronization time $T_{\rm s}$
(e.g., one can take the clustering time of Fig. 3(a)). If $T_{\rm s}\gg T_{\rm
tr}$, almost all the particles escape from the well and spread. In the
opposite limit, $T_{\rm s}\ll T_{\rm e}$, a complete self-trapping due to
alignment occurs: a cluster of synchronous trapped quasiperiodically moving
particles appears. Thus, the potential well becomes an effective “black hole”:
particles are trapped perpetually due to dissipative alignment. In an
intermediate case $T_{\rm s}\sim T_{\rm tr}$, some particles escape while
others form a perpetually trapped cluster.

A combination of the LJ and the alignment interactions allows for
synchronization of the crystal. Point clusters cannot be formed because
particles cannot come close to each other due to the LJ repulsion at small
distances. Thus, only orientations $\theta_{k}$ can synchronize. Here again, a
relation between the synchronization time and the lifetime due to potential
forces is crucial.

We describe the combined action of the LJ and alignment interactions for
particles in the finite-well potential, because, for small enough velocities,
the particles are practically trapped forever. We chose the width of the well
to be approximately ten times larger than the characteristic spatial scale of
the LJ potential so that only a few particles fit into the well.

Figure 4: Seven particles with a LJ potential and alignment in a potential
well (color coding). Panel (a): the order parameter $R$ vs time. Panels (b,c):
trajectories of the c.o.m on the $(x,y)$ plane in the disordered $t<2300$ (b)
and in the ordered $t>3000$ (c) domains.

Figure 4 illustrates the dynamics of seven particles. Initially, they form a
shaking crystal with a random motion (reminiscent of a confined random walk)
of the c.o.m, Fig. 4(b). The orientational order is characterized by the
Kuramoto order parameter $R=|\langle e^{i\theta}\rangle|$, Fig. 4(a). In the
disordered state, the order parameter fluctuates around $R\approx 0.4$. At
$t\approx 2500$, an abrupt transition to synchrony in the directions of motion
is observed, beyond this transition $R\approx 1$. The crystal becomes ordered,
and the c.o.m. performs a quasiperiodic motion in the well, Fig. 4(c), and all
the particles become perpetually trapped (video 5 in sup ).

The abrupt transition to synchrony from a chaotic crystal should be contrasted
to the clustering transition without the LJ coupling. In the latter case, the
order parameter grows gradually as the particles continuously come closer and
merge. In contradistinction, the process depicted in Fig. 4(a) is
characterized as transient chaos that abruptly ends in an absorbing
synchronized state (see a general exposition of transient chaos Lai and Tél
(2011), and a case of chiral active particles in Pikovsky (2021)).

To examine the dependence of the synchronization time on the size of the
crystal, we consider a set of $N$ active particles. We take the LJ interaction
and additional alignment coupling with $r_{0}=1$ (i.e., the same range as the
LJ potential) but without a confining external potential. To ensure that the
lifetime of the crystal is larger than a characteristic synchronization time
(cf. Fig. 2(c)), we considered particles with small velocity $V=0.1$. Figure
3(b) shows that the dependence of the synchronization time on the size of the
crystal is exponential $T_{\rm syn}\sim g(\epsilon)\exp[h(\epsilon)N]$, with
$\epsilon$-dependent factors $g,h$. That implies that here a super transient
behavior Lai and Tél (2011); Pikovsky (2021) is observed, for which a
characteristic time exponentially grows with the system size.

Figure 5: An example of trapping dynamics for $\epsilon=0.04$ and $V=0.15$ (a
particle is added every 2000 time units). Blue line: evolution of the order
parameter (left axis) showing domains of disorder and synchrony ($R\approx
1$). Red curve: number of particles (right axis). Inset: the average size of
the escaped cluster vs $\epsilon$. The time interval is proportional to the
size of the cluster because a particle is added every 2000 time units.

Finally, we studied how many particles can be trapped in a potential well in
the process illustrated in Fig. 4. We assume for simplicity that
$r_{0}=\infty$, i.e., all particles in the well equally contribute to the
alignment. The particles from outside “bombard” the trap one-by-one at regular
time intervals. They form a crystal, the dynamics of which is followed using
the alignment order parameter $R$, see Fig. 5. One sees that in some cases
adding a particle melts the crystal. In other cases, a particle is absorbed
into the existing order. Finally, the crystal as a whole escapes the potential
well because an incoming particle “kicks” it from a trapped quasiperiodic
regime into an escaping trajectory. We observed that only a synchronized
crystal can escape (video 6 in sup ). Fig. 5(inset) depicts the average size
of the escaping crystal $\langle S\rangle$ vs $\epsilon$. It shows a
significant increase of a typical escaping crystal size for low alignment
rates $\epsilon$. That occurs because for global coupling, the alignment rate
is proportional to the number of particles. Therefore, for small $\epsilon$, a
sufficient number of interacting particles is required to form a synchronous
crystal.

In conclusion, we investigated the trapping of individual and interacting
active particles. The problem is nontrivial for slow particles (deep
potentials), that spend a long time close to the bottom. We have found that
non-interacting particles can be trapped only for a finite time due to the
Hamiltonian structure of the equations of motion. In turn, an alignment of
particles brings dissipation and establishes a time arrow. So, the potential
well becomes a “black hole” with perpetually captured particles on
quasiperiodic trajectories. The particles form either a cluster or a coherent
crystal if additional LJ coupling is taken into account. With the particle
bombardment, we observe a highly nontrivial trapping behavior: melting and re-
synchronization of crystals, and eventually a coherent escape of the entire
assembly.

The problem we have studied is relevant for the understanding of active matter
subject to quenched environmental disorder. For the media, which can be
interpreted as an array of randomly positioned traps, the process can be
represented as a sequence of the described above trapping events. Our study
also shows that the disorder may have a finite trapping capacity. Once the
traps are filled up, the further bombardment may lead to the spontaneous
avalanche-like release of many particles. While this phenomenon is reminiscent
of self-organized criticality (SOC) Bak  _et al._ (1988), there are many
fundamental differences: the trapped states are highly dynamic and form
chaotically or coherently moving crystals. Likely, this will result in a
different distribution for the avalanche sizes and other statistical
characteristics.

###### Acknowledgements.

We thank S. Klapp and H. Stark for useful discussions. A. P. was supported by
the Russian Science Foundation, grant 17-12-01534, and by German Science
Foundation, grant PI 220/22-1. Research of I.S.A. was supported by the U.S.
Department of Energy, Office of Science, Basic Energy Sciences, under Award
DE-SC0020964.

## References

  * Bechinger  _et al._ (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Reviews of Modern Physics 88, 045006 (2016).
  * Aranson (2013) I. S. Aranson, Physics-Uspekhi 56, 79 (2013).
  * Gompper  _et al._ (2020) G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, _et al._ , Journal of Physics: Condensed Matter 32, 193001 (2020).
  * Chaté (2020) H. Chaté, Annual Review of Condensed Matter Physics 11, 189 (2020).
  * Peshkov  _et al._ (2012) A. Peshkov, I. S. Aranson, E. Bertin, H. Chaté, and F. Ginelli, Physical Review Letters 109, 268701 (2012).
  * Patelli  _et al._ (2019) A. Patelli, I. Djafer-Cherif, I. S. Aranson, E. Bertin, and H. Chaté, Physical Review Letters 123, 258001 (2019).
  * Solon  _et al._ (2015) A. P. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri, M. Kardar, and J. Tailleur, Nature Physics 11, 673 (2015).
  * Alert  _et al._ (2020) R. Alert, J.-F. Joanny, and J. Casademunt, Nature Physics 16, 682 (2020).
  * Peruani and Aranson (2018) F. Peruani and I. S. Aranson, Physical Review Letters 120, 238101 (2018).
  * Duan  _et al._ (2021) Y. Duan, B. Mahault, Y.-q. Ma, X.-q. Shi, and H. Chaté, Physical Review Letters 126, 178001 (2021).
  * Ro  _et al._ (2021) S. Ro, Y. Kafri, M. Kardar, and J. Tailleur, Physical Review Letters 126, 048003 (2021).
  * Olsen  _et al._ (2021) K. S. Olsen, L. Angheluta, and E. G. Flekkøy, Soft Matter 17, 2151 (2021).
  * Waisbord  _et al._ (2021) N. Waisbord, A. Dehkharghani, and J. S. Guasto, Nature Communications 12, 1 (2021).
  * Bychkov and Dykhne (1971) Y. A. Bychkov and A. M. Dykhne, Theoretical and Mathematical Physics 6, 307 (1971).
  * Dauchot and Démery (2019) O. Dauchot and V. Démery, Physical Review Letters 122, 068002 (2019).
  * Schmidt  _et al._ (2021) F. Schmidt, H. Šípová-Jungová, M. Käll, A. Würger, and G. Volpe, Nature Communications 12, 1 (2021).
  * Takatori  _et al._ (2016) S. C. Takatori, R. De Dier, J. Vermant, and J. F. Brady, Nature Communications 7, 1 (2016).
  * Pototsky and Stark (2012) A. Pototsky and H. Stark, EPL 98, 50004 (2012).
  * Hennes  _et al._ (2014) M. Hennes, K. Wolff, and H. Stark, Physical Review Letters 112, 238104 (2014).
  * Wexler  _et al._ (2020) D. Wexler, N. Gov, K. Ø. Rasmussen, and G. Bel, Physical Review Research 2, 013003 (2020).
  * Das  _et al._ (2018) S. Das, G. Gompper, and R. G. Winkler, New Journal of Physics 20, 015001 (2018).
  * Jahanshahi  _et al._ (2017) S. Jahanshahi, H. Löwen, and B. Ten Hagen, Physical Review E 95, 022606 (2017).
  * Liebchen  _et al._ (2016) B. Liebchen, M. E. Cates, and D. Marenduzzo, Soft Matter 12, 7259 (2016).
  * Liebchen and Löwen (2018) B. Liebchen and H. Löwen, Accounts of Chemical Research 51, 2982 (2018).
  * Stark (2018) H. Stark, Accounts of Chemical Research 51, 2681 (2018).
  * Kravtsov and Orlov (1990) Y. A. Kravtsov and Y. I. Orlov,  _Geometrical Optics of Inhomogeneous Media_ (Springer, Berlin, Heidelberg, 1990).
  * (27) See Supplemental Material at $http://link.aps.org/supplemental/$ for computational videos .
  * Vicsek  _et al._ (1995) T. Vicsek, A. Czirók, E. Ben-Jacob, I. Cohen, and O. Shochet, Phys. Rev. Lett. 75, 1226 (1995).
  * Kruk  _et al._ (2018) N. Kruk, Y. Maistrenko, and H. Koeppl, Phys. Rev. E 98, 032219 (2018).
  * Pecora and Carroll (1998) L. M. Pecora and T. L. Carroll, Physical Review Letters 80, 2109 (1998).
  * Lai and Tél (2011) Y.-C. Lai and T. Tél,  _Transient Chaos_ (Springer, New York, 2011).
  * Pikovsky (2021) A. Pikovsky, J. Phys. Complexity 2, 025009 (2021).
  * Bak  _et al._ (1988) P. Bak, C. Tang, and K. Wiesenfeld, Physical Review A 38, 364 (1988).