# Simulated slidequakes: Insights from DEM simulations into the high-frequency seismic signal generated by geophysical granular flows

Matthew I Arran, A Mangeney, J de Rosny, Renaud Toussaint

<figure><img src="image_1.png" /></figure>

## To cite this version:

Matthew I Arran, A Mangeney, J de Rosny, Renaud Toussaint. Simulated slidequakes: Insights from DEM simulations into the high-frequency seismic signal generated by geophysical granular flows. 2023. [hal-04311906]

HAL Id: hal-04311906

https://hal.science/hal-04311906v1

Preprint submitted on 28 Nov 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulated slidequakes: Insights from DEM simulations into the high-frequency seismic signal generated by geophysical granular flows

M. I. Arran<fnref n="1" />, A. Mangeney<fnref n="1" />, J. De Rosny<fnref n="2" />, R. Toussaint<fnref n="3" />,<fnref n="4" />

## Key Points:

* We used discrete-element simulations to examine the fluctuating forces exerted by steady granular flows on fixed, rough bases

* Without basal slip, results differ from Arran et al. (2021)'s and a new model most accurately predicts high-frequency forces

* A bulk inertial number determines the relative magnitudes of mean, intermediate-frequency, and high-frequency forces

<footnote n="1">¹Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France</footnote>

<footnote n="2">²Institut Langevin, ESPCI Paris, PSL University, CNRS, 75005 Paris, France</footnote>

<footnote n="3">³Université de Strasbourg, CNRS, Institut Terre et Environnement de Strasbourg, UMR 7063, F-67084</footnote>

<footnote>Strasbourg, France</footnote>

<footnote n="4">⁴SFF PoreLab, The Njord Centre, Department of Physics, University of Oslo, Oslo, Norway</footnote>
**Abstract**

Geophysical granular flows generate seismic signals known as ‘slidequakes’ or ‘landquakes’, with low-frequency components whose generation by mean forces is widely used to infer hazard-relevant flow properties. Many more such properties could be inferred by understanding the fluctuating forces that generate slidequakes’ higher frequency components and, to do so, Arran et al. (2021, https://doi.org/10.1029/2021JF006172) (A21) compared the predictions of pre-existing physical models to the forces exerted by laboratory-scale flows. However, A21 was unable to establish whether the laboratory flows exhibited basal slip, and the conditions for applying its results are therefore unclear. Here, we describe discrete-element simulations that examined the fluctuating forces exerted by steady, downslope-periodic granular flows on fixed, rough bases that prevented basal slip. We show that, in its absence, A21’s results do not hold: simulated basal forces’ power spectra have high-frequency components more accurately predicted using mean shear rates than using depth-averaged flow velocities, and can have intermediate-frequency components which we relate to chains of prolonged interparticle contacts. We develop a ‘minimal model’, which uses a flow’s collisional properties to even more accurately predict the high-frequency components, and empirically parametrize this model in terms of mean flow properties. Finally, we demonstrate that the bulk inertial number determines not only the magnitude ratio of rapidly fluctuating and mean forces on a unit basal area, consistent with A21, but also the relative magnitudes of the high and intermediate-frequency force components.

**Plain Language Summary**

Any geophysical granular flow - such as a landslide, rockfall, or debris flow - exerts fluctuating forces that cause the ground to vibrate, in a 'slidequake' that can provide useful information about the flow. Here, we examine simulated slidequakes: computer models of the individual particles within idealized flows, the collisions between them, and the rapidly fluctuating forces they exert on the flow's base. By recording particle and collision properties throughout the flow, we examine pre-existing models for the fluctuating forces; develop, test, and simplify a new model; and relate ratios between forces to an 'inertial number' that characterizes different flows. Our results differ from those of laboratory experiments that previously investigated slidequakes, but the two sets of results can be combined to provide information about real geophysical flows.

# 1 Introduction

## 1.1 Background

Geophysical granular flows such as landslides generate seismic waves, otherwise referred to as ‘slidequakes’ or ‘landquakes’, and the properties of landquakes’ low-frequency components have been linked to the flows’ properties for over a century. Following Galitzine (Golitsyn) (1915)’s association of the 1911 Sarez earthquake with a landslide in the Pamir mountains, articles initially focussed on linking large landslides with their landquakes, as in Weichert et al. (1994)’s review and analysis. But Kanamori and Given (1982)’s detailed analysis of a rock avalanche’s low-frequency seismic waves, using the theoretical framework of a single time-dependent point force, led to other reconstructions of landslides’ force histories (e.g. Hasegawa & Kanamori, 1987; Eissler & Kanamori, 1987). And this single-force framework, described by e.g. Kawakatsu (1989) and Dahlen (1993), has permitted the detection of geophysical flows from seismic signals both locally (Lin et al., 2010; Chao et al., 2017) and globally (Ekström & Stark, 2013), while force inversion with single-block landslide models has allowed approximate reconstruction of flows’ evolution or of certain flow parameters (Brodsky et al., 2003; La Rocca et al., 2004; Gualtieri & Ekström, 2017; Ekström & Stark, 2013; Allstadt, 2013; Yamada et al., 2013; Coe et al., 2016).
However, landquakes' low-frequency components provide information about only the flow processes with long timescales and large lengthscales by which they are predominantly generated. Applying the single-force framework to continuum simulations of landslides, rather than single-block models, permits more accurate flow and parameter reconstruction, and inference about processes with smaller extents than the landslide's (e.g. sub-event evolution, erosion, or spatial variation in friction, in Favreau et al., 2010; Moretti et al., 2012, 2015, 2020; Yamada et al., 2016, 2018; Zhao et al., 2020). But short-lengthscale, short-timescale processes are not modeled by such simulations and cannot be resolved from low-frequency components of landquakes, which therefore cannot be used to infer particle sizes, vertical flow profiles, or many other flow properties. Furthermore, Allstadt et al. (2018) suggests that low-frequency components of landquakes, at frequencies < 0.1 Hz in the field, can typically only be detected for landslide volumes over 10⁷ m³.

Landquakes' high-frequency components are predominantly generated by flow processes with short timescales and small lengthscales, so research into such components offers a means of inferring flow properties that may otherwise be inaccessible. Arran et al. (2021) provides a brief review of the empirical relations observed between the properties of a landslide and those of its landquake's high-frequency components' envelope (by Norris, 1994; Deparis et al., 2008; Schneider et al., 2010; Hibert et al., 2011; Levy et al., 2015; Hibert et al., 2015, 2017; Allstadt et al., 2020), and a detailed review of five different models, based on the 'stochastic impact framework' of Tsai et al. (2012), that predict a landquakes' high-frequency power spectrum from a landslide's properties (those of Kean et al., 2015; Lai et al., 2018; Farin, Tsai, et al., 2019; Bachelet, 2018; Bachelet et al., 2023). The Supplementary Information of Arran et al. (2021) includes a discussion of the conditions under which the stochastic impact framework will be valid, while the bulk of the paper tests the different models using experimental data for which those conditions were satisfied, from laboratory flows of monodisperse glass spheres.

Of the models it compares, Arran et al. (2021) finds the 'thin-flow' model of Farin, Tsai, et al. (2019) to predict most accurately the high-frequency seismic signals produced by the laboratory-scale flows, but cannot discriminate between two alternative explanations for this accuracy. In the 'thin-flow' model, the particles within a flow that impact its base have velocities distributed around the flow's depth-averaged downslope velocity, and the laboratory flows could have exhibited the plug-like velocity profiles this implies. Alternatively, flow velocities could have decreased towards zero towards the flow's base, consistent with the 'no-slip' condition documented, for rough flow bases such as that used in the laboratory experiments, by e.g. GDR MiDi (2004); Jing et al. (2016). But in the latter case, the accuracy of the 'thin-flow' model's predictions can be explained only by non-basal particles making significant contributions to seismic signals, as in the models of Bachelet (2018); Bachelet et al. (2023) and, more recently, Zhang et al. (2021)'s model of single-body impacts by chains of particles.

To determine the applicability of Arran et al. (2021)'s results and link the properties of flows to the properties of the high-frequency seismic signals they produce, it's necessary to establish which explanation is correct. If the 'thin-flow' model's predictions owe their accuracy to basal slip, then the model will apply only to flow regions over which basal friction is insufficient to restrain basal flow velocities, and a different model will be required elsewhere. If, on the other hand, non-basal particles contribute significantly to seismic signals, then conditions must be established under which such contributions compensate for basal velocities lower than those theorized by the 'thin-flow' model.

We therefore investigate the production of high-frequency seismic signals by flows for which the absence of basal slip can be confirmed, and the contributions of non-basal particles quantified. Because we're concerned with particle velocities in the
interior of a granular flow, which physical experiments on flow seismology have been
unable to interrogate (Huang et al., 2007; Taylor & Brodsky, 2017; Farin et al., 2018;
Farin, Mangeney, et al., 2019; Allstadt et al., 2020; Bachelet et al., 2023), we conduct
numerical simulations using a discrete element method (DEM), in which all variables
are accessible.

We structure our paper as follows. After introducing and justifying the variety of DEM we use, in section 1.2, we describe our DEM implementation, in section 2.1, and simulation protocol, in section 2.2. Section 2.3 describes how we calculate different models’ predictions for the seismic signal-determining power spectra $P_F$, of the forces each flow exerts on its base. Then, sections 3.1 and 3.2 describe the properties of the different flows and of their signals, respectively, via the depth profiles of kinematic and collisional properties and via the measured power spectra $P_F$. We compare, in section 3.3, the measured power spectra with the predictions of section 2.3 and introduce, in section 3.4, a semi-empirical parametrization of the most successful model. The following sections, 4.1 and 4.2, then discuss Arran et al. (2021)’s experiments in the context of our measurements. Section 4.1 addresses the two possible interpretations of Arran et al. (2021)’s model comparison, while section 4.2 investigates another of Arran et al. (2021)’s results: that the inertial number $I$, a non-dimensional quantity that describes a granular flow’s flow state (GDR MiDi, 2004; da Cruz et al., 2005; Jop et al., 2006), determines the squared ratio $\delta F^2$ between the mean and fluctuating forces exerted by a granular flow. Finally, section 4.3 discusses the signals we observe at intermediate frequencies, and section 5 concludes.

## 1.2 DEM Simulations

In a DEM simulation, each particle within a flow has a specified position, velocity, orientation, and angular velocity, which evolve over time as a result of interactions with other particles, boundaries, a medium, and/or external forces. Since being proposed by Cundall (1971) and examined by Cundall and Strack (1979), such methods have been used and validated in contexts both industrial (e.g. Daraio et al., 2019; Gao et al., 2021) and geophysical (e.g. Wu et al., 2016; Sun et al., 2016), with modellers having access to all state variables, both of individual particles and of the interactions between them. Local weighted averages of these particle-scale properties, calculated through the coarse-graining process described by e.g. Goldhirsch (2010); Weinhart et al. (2012), can then be compared to the predictions of continuum theories or to similar averages from physical experiments.

To perform simulations in a given situation, interparticle interactions must be
simplified. The experiments of Arran et al. (2021) involved approximately spherical
particles with mean diameter $d = 2$ mm, composed of soda lime glass with density
$\rho \approx 2500$ kg/m$^3$, Young's modulus $E \approx 70$ GPa, and Poisson's ratio $\nu \approx 0.2$ (Sigmund
Lindner, 2018). In such soda lime glass, neglecting its $O(10^5)$-year structural relaxation
at room temperature, the Rayleigh wave velocity would be around $5000$ m s$^{-1}$ and the
strain relaxation time would be of order $10^{-9}$ s (Hunklinger & Arnold, 2012). In
comparison, typical collision velocities will have been of order $\sqrt{gd} \approx 0.1$ m s$^{-1}$, for
gravitational acceleration $g = 9.81$ m/s$^2$. Therefore, noting that a normal deformation
$\delta_n$ is associated with a contact area of order $\delta_n d$, by geometry, and so an elastic force
of order $E\delta_n^2$, we can equate the order-$E\delta_n^3$ elastic energy at the point of maximum
deformation with the order-$\rho g d^4$ kinetic energy at the point of impact, indicating
that typical particle deformation will have been of order $O(\delta_n) = (\rho g d^4/E)^{1/3} \approx$
$2 \times 10^{-6}$ m, during collisions of duration of order $O(\delta_n)/\sqrt{gd} = 10^{-5}$ s. Consequently,
particles will have deformed approximately quasi-statically during each collision, with
deformation independent of other collisions. With peak contact pressures of order
$(\rho g d^4/E)^{1/3} E/d \approx 1$ GPa and the compressive yield stress of soda lime glass reportedly
between 0.25 GPa and 4.5 GPa (Meredith & Swab, 2020), some deformation may have
been plastic, but we show in section S2 of the Supplementary Information that the plastic yield consistent with experimentally observed coefficients of restitution will have had a negligible effect on particle collisions. Particles were sufficiently large, dense, and dry that air resistance, capillary forces, van der Waals forces, and external forces other than gravity will also have been negligible. We believe that similar arguments will apply to many geophysical flows, and therefore consider only binary, contact-mediated, quasi-static interactions between particles, without adhesion or plastic deformation.

Dividing such interactions into those that exert forces at the idealized point of contact and those that exert only torques, each is separable into normal and tangential components, and elastic and viscous contributions, with previous authors describing and using a wide range of interaction models. Kruggel-Emden et al. (2007, 2008); Thornton et al. (2011, 2013) review and examine normal, tangential, elastic, and viscoelastic contact forces, respectively, from simple models with normal and tangential spring constants to models derived from Hertz (1881) and Mindlin and Deresiewicz (1953)'s elasticity-theoretic analyses, in which effective stiffnesses are a function of normal deformation. Notably, Kruggel-Emden et al. (2007) recommends the use of Kuwabara and Kono (1987)'s Hertzian, viscoelastic model for the normal contact force, as theoretically justified by Brilliantov et al. (2015); Kruggel-Emden et al. (2008) demonstrates the importance of representing local slip by limiting tangential forces with Coulomb friction; and Thornton et al. (2011) shows that, to ensure energy conservation when the effective tangential stiffness varies, any tangential force must be updated incrementally and relaxed whenever the normal force decreases.

Among interactions that exert only torques, the separation into normal and tangential components divides torsional resistance from rolling resistance, with Lubkin (1951) deriving Hertz theory's prediction for the former and Dintwa et al. (2005) deriving and validating a quasi-static, viscoelastic extension and simplification. Brilliantov and Pöschel (1998) derives an equivalent quasi-static, viscous resistance to rolling, but Ai et al. (2011)'s review shows purely frictional or purely viscous rolling resistances to be insufficient for realistic, imperfectly spherical particles. Instead, Ai et al. (2011) recommends a model 'C' that uses Iwashita and Oda (1998)'s general framework of visco-elastic, friction-limited rolling resistance; Jiang et al. (2005)'s correction of the rotation velocity for particles of differing radius; incremental calculation to account for rollback; and Bardet and Huang (1993)'s expression for rolling stiffness, which is consistent with the rocking of smooth, imperfectly spherical particles about local minima of curvature.

In most applications of DEM simulations, authors are concerned with granular flows' macroscopic behavior, so use interactions that are conducive to simulating large numbers of particles over long periods of time. Most commonly, the normal and tangential spring constants for contact forces are taken to be constant and torsion and rolling resistance are neglected (Kuhn et al., 2020). Simulations often use effective particle stiffnesses far lower than are realistic, lengthening collision durations and permitting larger timesteps (Lommen et al., 2014).

However, we are concerned with the rapid processes with short lengthscales that generate the high-frequency components of granular flows' seismic signals, so we want to accurately model the microscopic interactions between particles. Only by modeling particles' resistance to torsion and rolling can we accurately simulate such micromechanics, and only by simulating Hertzian interactions with realistic particle stiffnesses can we correctly model impact-dependent collision durations, recovering a realistic separation between the timescales of particle rearrangement and of Hertzian collisions ($\sqrt{d/g} \approx 10^{-2}$ s and $(\rho^2 d^5/E^2g)^{1/6} \approx 10^{-4}$ s, respectively, for the experiments of Arran et al. (2021)).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** *Schematic of simulations' interparticle interactions.* Between two particles in contact, the simulated normal and tangential forces ($\mathbf{F}_n$ and $\mathbf{F}_t$) and couples ($\mathbf{M}_n$, $\mathbf{M}_t$, and $\mathbf{M}_r$) depend on the particles' normal and tangential deformations ($\delta_n$ and $\delta_t$) and angular deformations ($\theta_n$ and $\theta_r$). Behavior is viscoelastic and, except in the case of the normal reaction force, friction-limited, so interactions can be modeled at each timestep by the systems of springs, dash-pots, and friction elements that we represent below each force diagram, with viscosities $\beta\kappa_n$, $\beta k_t$, $\beta k_n$, and $\beta_r k_r$ proportional to the $\delta_n$-dependent spring stiffnesses $\kappa_n$, $k_t$, $k_n$, and $k_r$, and friction coefficients $\mu_t$, $\mu_n$, and $\mu_r$.</figcaption></figure>

# 2 Methods

## 2.1 Simulation conditions

We conduct simulations using a fork (Arran, 2023) of MercuryDPM (Weinhart et al., 2020), which simulates particle interactions and evolution using the formalism described by Luding (2008) and the velocity Verlet algorithm (Allen & Tildesley, 1989). Particles' interactions, the interactions' parameters, and boundary conditions are chosen for consistency with Arran et al. (2021)'s experiments.

Particles are simulated as spheres, with corresponding quasi-static, contact-mediated interactions: normal forces, tangential forces, torsional resistance, and rolling resistance. Specifically, between any two particles in contact, we apply i) the viscoelastic normal force introduced by Kuwabara and Kono (1987); ii) the incremental Mindlin tangential sliding resistance described as the TM model in Thornton et al. (2013), but with consistent, quasi-static viscous dissipation; iii) an incremental, linearized Lubkin (1951) friction-limited torsional elasticity with viscous dissipation as described by Dintwa et al. (2005); and iv) Ai et al. (2011)'s model 'C' for rolling resistance, except with the equation for damping given, for consistency, by Brilliantov and Pöschel (1998). Figure 1 illustrates the four interactions simulated between particles in contact, while section S3 of the Supplementary Information describes them in detail.

Parameters for particle interactions are standard values for 2-mm soda lime glass beads (Seward & Vascott, 2005), other than the damping parameter $\beta$; the stiffness, damping, and friction parameters for rolling resistance $\kappa$, $\beta_r$, and $\mu_r$; the coefficient of sliding friction $\mu_t$; and derived parameters such as $\mu_n$. Simulated particle diameters are drawn from a normal distribution with mean $d = 2$ mm and standard deviation 0.1 mm, with simulated density $\rho = 2530$ kg/m$^3$, Young's modulus $E = 74$ GPa, and Poisson's ratio $\nu = 0.22$, and with gravitational acceleration $g = 9.81$ m/s$^2$. The values of $\beta = 4.5 \times 10^{-8}$ s; $\kappa = 0.742$, $\beta_r = 0.002$ s, and $\mu_r = 0.0175$; and $\mu_t = 0.237$
are calibrated in order that simulations are consistent with, respectively, Foerster et al. (1994)'s measurements of glass beads' coefficient of restitution; Arran et al. (2021)'s particles' steady and transient resistances to rolling in single-bead experiments; and these particles' angle of repose. To accurately resolve collisions, the timestep is $\delta t = 1.75 \times 10^{-5} \sqrt{d/g} = 2.5 \times 10^{-7}$ s. Section S4 of the Supplementary Information describes the details of parameter selection, with Table S1 listing their values. In all calculations, parameters are non-dimensionalized by the mass scale $\rho d^3 = 0.02$ kg, the length scale $d = 0.002$ m, and the time scale $\sqrt{d/g} = 0.014$ s.

We conduct simulations in 3D, with a domain periodic in the downslope and cross-slope directions, $x$ and $y$, and a fixed, rough bed with mean surface the plane $z = 0$. The domain extends distances $n_x d$ and $n_y d$ in the $x$ and $y$ directions, with different simulations having $n_x = n_y = 10$, $n_x = n_y = 14$, and $n_x = n_y = 20$. For each domain, the bed is roughened using a preparatory simulation, in which a number $2n_x n_y$ of non-intersecting particles are initialized above a rigid, frictional, planar base at $z = -d$, with uniformly distributed random offsets in $x$ and $y$ from a pattern of regular lines. These particles are allowed to settle under a uniform, vertical gravitational force, and fixed in place if and when they make contact with the base. When motion has ceased, we record the positions and radii of all those fixed particles, which are kept fixed in future simulations, and henceforth referred to as basal particles. Interactions with basal particles are identical to those with non-basal particles, while interactions with the base are identical to those with a particle of infinite radius and interactions across the periodic boundary are handled using 'virtual' particles: copies of near-boundary particles translated in $x$ and/or $y$ as applicable, for which all interactions across the opposite boundary are calculated and then applied to the original particles.

For comparison with the experiments of Arran et al. (2021), we also simulate flow in a channel. In this case, we maintain periodicity in the $x$ direction but take $n_x = 10$ and $n_y = 50$, and simulate a PMMA wall at cross-slope position $y = 0$ and a plane of symmetry at $y = n_y d$, by setting interaction parameters as described in section S4.4 of the Supplementary Information. The rough bed is prepared in the same manner as for the doubly periodic simulations.

## 2.2 Simulation protocol

To ensure that analysis is performed on a steady-state flow, we perform a preparatory simulation for each domain size, bed angle, and flow depth. To attain a flow of approximate depth $n_z d$ in a domain of size $n_x d$ by $n_y d$, we initialize and fix in place the corresponding basal particles, before initializing $n_x n_y n_z$ non-overlapping particles above the rough bed, with their $x$ and $y$ offsets from a pattern of regular lines taken from a uniform distribution. We define a uniform gravitational field, directed at an angle $\theta$ to the vertical in the $x-z$ plane, and run the simulation with parameters identical to those that will be used for the main experiment. As the simulation progresses, we record the mean particle velocity every $10^{-3}$ simulated seconds and calculate the mean and standard deviation of this mean velocity over consecutive 2-simulated-second intervals. The simulated flow is considered to have attained a steady state at the end of an interval over which this mean and standard deviation are both within 5% of their values over the previous interval. Simulations that attained a steady state, without the flow stopping entirely, are indicated in Figure 2, while a video of steady flow is available as Movie S1 of the supplementary material.

Starting from the final system state of each preparatory simulation, we simulate $\Delta t = 2$ seconds of steady flow and record the properties of the flow, of the collisions between particles within it, and of the forces it applies to the base. Specifically, every $10^{-3}$ simulated seconds, we record each particle's volume $V_j$, position $\mathbf{x}_j$, and velocity
<figure><img src="image_4.png" /><figcaption>**Figure 2.** *Preparation of steady-state flows.* a) Evolution of the mean particle velocity $\bar{u}$ over the 4 simulated seconds before the time $t_0$ at which the flow is considered to be in a steady state, for $n_x = n_y = 10$, $\theta = 23.5^\circ$ and $n_z = 11$. b) Simulated system state at $t_0$, with the fixed basal particles in black and other colors indicating the magnitudes of particle velocities $\mathbf{u}_i$. c) Positions in $\theta - n_z$ parameter space of simulations that are conducted, and considered to attain a steady state.</figcaption></figure>
$\mathbf{u}_j$, as well as the list of contacts between particles and hence the set $C$ of particles that are connected to the base by a chain of contacts. Whenever a particle collides with another particle or with the base, we record the time of the collision, the position $\mathbf{x}_j^c$ of the point of contact, and the normal and tangential relative velocities $\dot{\delta}_{n,j}$ and $\dot{\delta}_{t,j}$ of the colliding particles. Finally, after each timestep of $\delta t = 2.5 \times 10^{-7}$ simulated seconds, we record the total force $\mathbf{F}(t)$ exerted by the flow on the basal particles and, in the few places where these basal particles do not screen the base from the flow, on the base itself.

From these recorded data, for each simulation, we reconstruct base-normal profiles of the flow's kinematic and collisional properties, and calculate the mean values and power spectra of the base-normal force exerted on the base by the flow. With the rough bed's surface at $z = 0$ as origin, spatially and temporally averaged profiles of volume fraction $\phi$, flow velocity $\mathbf{u}$, and granular temperature $T$ are calculated from the particle position and velocity records using MercuryDPM's inbuilt coarse-graining program, with truncated Gaussian weight functions $w_z$ of width $d$ and volume integral $V$ and with $3\rho\phi T/2$ the density of fluctuating kinetic energy. Using a newly written implementation of an equivalent procedure, which is described in section S5 of the Supplementary Information, we calculate the depth-dependent proportion $P_c(z)$ of particles connected to the base by a chain of contacts. We then use collision records to calculate the $w_z$-weighted mean number $n_I(z)$ of interparticle collisions per unit volume and unit time, as well as the similarly weighted mean squared normal and tangential relative velocities of colliding particles, $\langle \dot{\delta}_n^2 \rangle(z)$ and $\langle \|\dot{\delta}_t\|^2 \rangle(z)$, at the point of contact and the instant of collision. Finally, the power spectrum $P_F$ of each component of the force exerted by the flow is calculated by taking the discrete Fourier transform of each component of $\mathbf{F}$, as

$$P_F(j/\Delta t) = \frac{\delta t^2}{\Delta t} \left| \sum_{k=0}^{\Delta t/\Delta t-1} F(k\delta t) e^{-2\pi i j k \delta t / \Delta t} \right|^2, \quad (1)$$

with each component of the mean force $\langle \mathbf{F} \rangle$ having squared magnitude $|\langle F \rangle|^2 = P_F(0)/\Delta t$.

In simulating flow in a channel, we divide the domain cross-slope, into an 'outer' region $y < n_y d/2$ and an 'inner' region $y > n_y d/2$, so that the base of the inner region corresponds to the location of the instrumented plate described by Arran et al. (2021). Arran et al. (2021) studied the power spectra of the base-normal force exerted on this plate by experimental flows, so we both record and calculate the power spectrum of each component of the corresponding net force $\mathbf{F}$: that exerted on each simulation's fixed base via any contact points within the inner region. To assess models for these power spectra $P_F$, we seek profiles of the kinematic and collisional properties only of the flow region responsible for $\mathbf{F}$, and therefore calculate coarse-graining averages only over particles with centers in the inner region or over collisions with contact points in the inner region. This is illustrated in section S4.4 of the Supplementary Information.

## 2.3 Model predictions

From the profiles of kinematic properties $\phi$, $\mathbf{u}$ and $T$, we are able to calculate the predictions of the models of Kean et al. (2015); Lai et al. (2018); Farin, Tsai, et al. (2019); Bachelet (2018); Bachelet et al. (2023) for the high-frequency power spectrum $P_{F_z}(f)$ of the base-normal force $F_z$ exerted by each flow on its fixed base. To do so, we first define quantities analagous to those used in section 2.3 of Arran et al. (2021):
the surface area over which the force **F** is measured

$$A = \begin{cases} n_x n_y d^2 & \text{for doubly periodic domains} \\ n_x n_y d^2/2 & \text{for the simulated channel;} \end{cases} \qquad (2)$$

the mean basal pressure

$$p = \langle F_z \rangle / A; \qquad (3)$$

the flow depth, defined as the distance from the rough bed's surface to the surface at which the flow ceases to be dense (see Figures 3a,g),

$$h = \min\{z > \operatorname{argmax} \phi | \phi(z) < \max(\phi)/2\}; \qquad (4)$$

the bulk volume fraction

$$\bar{\phi} = \frac{1}{h} \int_{-d}^{\infty} \phi(z) \, dz; \qquad (5)$$

and the mean downslope velocity

$$\bar{u} = \frac{1}{\bar{\phi}h} \int_{-d}^{\infty} \phi(z) u_x(z) \, dz. \qquad (6)$$

We also define a non-dimensional function $\zeta$ that describes the frequency dependence of models' predictions for $P_{F_z}$. Writing $F_I(t)$ for the normal force between two particles of density $\rho$ and diameter $d$, undergoing a Hertzian collision with impact velocity $\dot{\delta}_n$, we note that its spectral density $|\tilde{F}_I(f)|^2$ depends on frequency $f$ only via a collision timescale $\tau$ determined by the Young's modulus $E$ and Poisson's ratio $\nu$ of the particles' material (Hertz, 1881):

$$\tau(\dot{\delta}_n) = \left[ \frac{\pi^2 \rho^2 (1 - \nu^2)^2}{4E^2 \dot{\delta}_n} \right]^{1/5} d. \qquad (7)$$

Consequently, we can introduce a non-dimensional frequency $f$ and define

$$\zeta(f) = |\tilde{F}_I(f/\tau(\dot{\delta}_n))|^2 / \left( \frac{\pi \rho d^3 \dot{\delta}_n}{3} \right)^2. \qquad (8)$$

As described in section S6 of the Supplementary Information, $\zeta(f)$ is approximately equal to 1 for $f \ll 1$ and is characterized by its rapid drop-off beyond a corner frequency

$$f_c = \min\{f > 0 | \zeta(2f) \le \zeta(f)/10\} = 0.200, \qquad (9)$$

with $\zeta(f_c) = 0.529$ and $\zeta(f) \ll 1$ for $f \gg f_c$.

From these quantities, we follow Arran et al. (2021) in calculating various models' predictions for the base-normal force's power spectral density $P_F^0$, at frequencies that are high enough for the stochastic impact framework to apply but below those corresponding to the durations of individual collisions. Kean et al. (2015)'s model's prediction for a flow's seismic signal is

$$\hat{P}_F^0 = K A \bar{u} p^2 / d^3, \qquad (10)$$

for arbitrary fitting parameter $K$; Lai et al. (2018)'s model's is

$$\hat{P}_F^0 = \pi^2 \rho^2 A d^3 \bar{u}^3 / 9; \qquad (11)$$

and Farin, Tsai, et al. (2019)'s models' predictions are

$$\hat{P}_F^0 = \pi^2 \rho^2 \bar{\phi} A d^3 (1 + e)^2 \xi(\nu) u_b^3 / 36, \qquad (12)$$
for $e = 0.9$ a constant coefficient of restitution and $\xi(v) \approx 0.053(1 + 5.6v^2)$ a non-dimensional function of fitting parameter $v$. In the 'thin-flow' model, the mean down slope velocity of basal particles $u_b = \bar{u}$ while, in the 'thick-flow' model, $u_b = \chi \bar{u} d/h$ for velocity profile shape factor $\chi = 1.25$. The parameter $v$ corresponds to the relative standard deviation of basal particle velocities, so realism requires that $v < 1$.

We also calculate the more general predictions that Arran et al. (2021) infers for the power spectral density $P_{F_z}(f)$, applicable at all frequencies for which the stochastic impact framework is valid. Extending the work of Farin, Tsai, et al. (2019) recovers predictions

$$ \hat{P}_{F_z}(f) = \frac{\int_{S^2} d^2\mathbf{e}_v \int_{\pi/6} d^2\mathbf{e}_n (\dot{\delta}_n \mathbf{e}_n \cdot \mathbf{e}_z)^2 \zeta(\tau(\dot{\delta}_n)f) \mathcal{H}(\dot{\delta}_n)}{\int_{S^2} d^2\mathbf{e}_v \int_{\pi/6} d^2\mathbf{e}_n (\dot{\delta}_n \mathbf{e}_n \cdot \mathbf{e}_z)^2 \mathcal{H}(\dot{\delta}_n)} \hat{P}_F^0, \quad (13) $$

for unit sphere $S^2$; unit spherical cap $S_{\pi/6}^2$ with maximum polar angle $\pi/6$; normal impact velocity $\dot{\delta}_n = u_b(\mathbf{e}_x + v\mathbf{e}_v) \cdot \mathbf{e}_n$; and Heaviside step function $\mathcal{H}$. Meanwhile, the comparable prediction from Bachelet (2018) is

$$ \hat{P}_{F_z}(f) = \frac{4\pi}{3} \bar{\phi} A \rho^2 d^3 \int_0^h u_x'(z) T(z) \zeta \left( \tau(\sqrt{3T})f \right) e^{-\gamma z} dz, \quad (14) $$

where $u_x'$ is the derivative of $u_x$ with respect to $z$, and $\gamma$ is a constant fitting parameter representing acoustic wave attenuation.

Finally, we develop a new model. We calculate predictions for $P_{F_z}(f)$ derived from the profiles of collisional properties, with the minimal assumption that the high-frequency signal is generated by the uncorrelated normal components of the forces exerted during elastic, Hertzian collisions of individual, diameter-$d$ particles. In this case, using equation (8)'s definition of $\zeta$, the total power spectral density of the base-normal basal force will be

$$ \hat{P}_{F_z}(f) := \frac{1}{\Delta t} |\tilde{F}_z|^2 = \frac{1}{\Delta t} \sum_j (\pi \rho d^3 \dot{\delta}_{n,j}/3)^2 \zeta(f \tau(\dot{\delta}_{n,j})) \iota_j \eta(f)^{\frac{z_j^c}{d}}, \quad (15) $$

where $\dot{\delta}_{n,j}$ and $z_j^c$ are the normal impact velocity and base-normal position of the $j$th collision during a time interval $\Delta t$, $\iota_j$ is a binary variable indicating whether the force exerted during this collision is transmitted to the base via a network of connected particles, and $\eta$ is the fraction of the spectral power transmitted a distance $d$ towards the base via this network. Assuming that attenuation is purely geometric and applying the arguments of the Supplementary Information's section S7, we expect the typical base-normal distance between the centers of two particles in contact to be $d/2$ and the typical fraction of spectral power transmitted between these particles to be $1/3$, so that we take a constant $\eta = (1/3)^{d/(d/2)} = 1/9$. Consequently, noting that $\zeta(f\tau(\dot{\delta}_{n,j}))$ is sufficiently insensitive to $\dot{\delta}_{n,j}$ to be approximable by $\zeta(f\tau(\sqrt{\langle\dot{\delta}_n^2(z_j^c)\rangle}))$, we replace the sum over collisions in equation (15) with an integral over base-normal length elements $dz$, with associated numbers of collisions $n_I(z)\Delta t A dz$ and probabilities of force transmission to the base $P_c(z - \frac{1}{4}d)$, to recover the prediction

$$ \hat{P}_{F_z}(f) = \frac{\pi^2 \rho^2 d^6}{9} A \int_0^\infty \mathcal{I}(f, z) \, dz \quad (16) $$

for

$$ \mathcal{I}(f, z) = n_I(\dot{\delta}_n^2) \zeta \left( f \tau \left( \sqrt{\langle \dot{\delta}_n^2 \rangle} \right) \right) P_c \left( z - \frac{1}{4} d \right) e^{-z \ln 9/d}. \quad (17) $$
# 3 Results

Full results are available in the interactive Jupyter notebook at https://gitlab.com/M_Arran/simulated-slidequakes (Arran et al., 2023). Here, we provide a summary, describing the depth profiles of kinematic properties (volume fraction $\phi$, downslope velocity $u_x$, and granular temperature $T$) and of collisional properties (the probability $P_c$ that a particle is connected to the base, the mean squared normal and tangential impact velocities $\langle \dot{\delta}_n^2 \rangle(z)$ and $\langle \|\dot{\delta}_t\|^2 \rangle(z)$, and the rate $n_I$ of interparticle collisions per unit volume). We describe the power spectra $P_F$ of the components of the basal forces exerted by simulated flows and, for the base-normal components, we compare the high-frequency components of the power spectra to the predictions listed above.

## 3.1 Kinematic and collisional profiles

Figures 3a-c and 3g-i present the depth profiles of $\phi$, $u_x$, and $T$ for two different combinations of slope angle $\theta$ and overburden $n_z$, and illustrate the range of variation among all simulations. Most importantly for our purposes, the profiles demonstrate that basal slip is consistently negligible, with $u_x(z)$ decreasing towards zero as $z \to 0$, and that the flow depth $h$ and bulk volume fraction $\bar{\phi}$ are appropriately defined, with $\phi(z)$ well approximated by $\bar{\phi}$ for $z \in (0, h)$ and decreasing rapidly towards zero for $z > h$. Furthermore, within the bulk of each flow, granular temperature $T$ is reasonably well approximated using the shear rate $u'_x(z)$, by

$$ \hat{T}(z) = \left( \frac{g_T d}{\tan \theta} u'_x(z) \right)^2 \qquad (18) $$

for some $g_T \approx 0.3$, as suggested by e.g. Woodhouse et al. (2010). As section S8 of the Supplementary Information demonstrates, increasing $n_x$ and $n_y$ has no systematic effect on the profiles for a given $\theta$ and $n_z$, indicating that the domain with $n_x = n_y = 10$ is sufficiently large for our results to be generally applicable.

Meanwhile, with increasing $\theta$ and $n_z$, the function $u_x(z)$ transitions from being mostly convex to mostly concave, with $T(z)$ changing in a corresponding fashion. Monotonically increasing for all $z < h - d$ when slope angles and overburdens are low, both $u'_x(z)$ and $T(z)$ become decreasing functions within the bulk of each flow. This change does not apply close to the rough bed, where both $u'_x(z)$ and $T(z)$ are increasing for all $\theta$ and $n_z$, or in a boundary layer at the flow's surface, in which $\phi$ decreases from $\bar{\phi}$ to zero and $T(z)$ diverges from $\hat{T}(z)$. However, increasing $\theta$ and $n_z$ increases the depth of this surficial boundary layer.

Moving to the depth profiles of $P_c$, $\langle \dot{\delta}_n^2 \rangle(z)$, and $\langle \|\dot{\delta}_t\|^2 \rangle(z)$, as illustrated by Figures 3d, e and 3j, k, we observe similar transitions with increasing $\theta$ and $n_z$. $P_c(z)$ decreases more rapidly from $P_c(-d/2) \approx 1$ and is increasingly well approximated by an exponential with least-squares best-fit non-dimensional decay length $\lambda(\theta, n_z)$:

$$ \hat{P}_c(z) = e^{-(2z+d)/2\lambda d}. \qquad (19) $$

As flows become more energetic, the fraction of particles that at a given time are in contact with an underlying particle decreases, and depends less on particles' distances from the base. Meanwhile, decreasing correlation between the velocities of colliding particles is evident in the increasing validity of the 'molecular chaos' approximation, in which the components of colliding particles' mean-flow-corrected velocities are uncorrelated and each have mean squared value $T$, so that $\langle \dot{\delta}_n^2 \rangle(z) \approx \langle \|\dot{\delta}_t\|^2 \rangle(z)/2 \approx 2T(z)$.

As illustrated in Figures 3f and 3l, profiles of $n_I(z)$ do not resemble those of $u'_x(z)/d^3$ (as would be the case if the rate of collisions were determined by particles’
<figure><img src="image_5.png" /><figcaption>**Figure 3.** *Depth profiles of flow properties.* Profiles with base-normal coordinate *z* of (a, g) volume fraction *φ*, (b, h) downslope velocity *u*<sub>x</sub>, (c, i) granular temperature *T*, (d, j) the probability *P*<sub>c</sub> that a particle is connected to the base, (e, k) mean squared normal and tangential impact velocities ⟨*δ*<sub>*n*</sub><sup>2</sup>⟩(*z*) and ⟨||*δ*<sub>*l*</sub>||<sup>2</sup>⟩(*z*), and (f, l) collision rates per unit volume *n*<sub>*l*</sub>, for (a-f) *θ* = 22.0° and *n*<sub>*z*</sub> = 8 and (g-l) *θ* = 24.5° and *n*<sub>*z*</sub> = 16. Dotted lines indicate flow depth *h* and bulk volume fraction *φ̃*, while theoretical approximants are: (c, i) equation (18)'s *T̃*, (d, j) equation (19)'s *P̃*<sub>c</sub>, (e, k) 2*T*, and (f, l) shear rate *u'<sub>x</sub>/d<sup>3</sup> and equation (21)'s *n̂*<sub>*l*</sub>. Red lines have separate axes.</figcaption></figure>
advection, with the mean flow, across underlying layers), but are instead well approximated by considering the kinetic theory of pressure that we introduce here. In a flow, the mean steady-state overburden pressure at a base-normal position $z_0$ is

$$ \hat{p}(z_0) = \rho g \cos \theta \int_{z_0}^{\infty} \phi(z) \, dz \qquad (20) $$

and, if a proportion $P_c(z_0)$ is borne by force chains of particles connected to the base, the vertical impulse that must be transferred by collisions, per unit area and time, through the surface $z = z_0$, is $(1 - P_c(z_0))\hat{p}(z_0)$. Applying the geometric arguments of the Supplementary Information's section S7, perfectly elastic collisions with normal impact velocity $\sqrt{\langle\dot{\delta}_n^2\rangle}$ and with no mean tangential force would each transfer a mean base-normal impulse $J = \pi \rho d^3 \sqrt{\langle\dot{\delta}_n^2\rangle}/12$ over a mean base-normal distance $d/2$. The number of such collisions transferring an impulse through the surface $z = z_0$, per unit area and time, would therefore be equal to both $(1 - P_c(z_0))\hat{p}(z_0)/J$ and $d/2$ times the local rate of collisions per unit volume, which is therefore

$$ \hat{n}_I(z_0) = \frac{24g \cos \theta (1 - P_c(z_0))}{\pi d^4 \sqrt{\langle\dot{\delta}_n^2\rangle(z_0)}} \int_{z_0}^{\infty} \phi(z) \, dz. \qquad (21) $$

For all simulations, $\hat{n}_I(z)$ well approximates the profile $n_I(z)$.

## 3.2 Power spectra of the basal force

At zero frequency, macroscopic force balances accurately predict the power spectra $P_F$ of the downslope, cross-slope, and base-normal components, $F_x$, $F_y$, and $F_z$, of the basal forces exerted by the steady flows. Mean basal pressure

$$ p = \frac{1}{A} \sqrt{\frac{P_{Fz}(0)}{\Delta t}} \approx \hat{p}(0) = \frac{\pi}{6} \rho g d n_z \cos \theta, \qquad (22) $$

while $P_{F_x}(0) \approx P_{F_z} \tan^2 \theta$ and $P_{F_y}(0) \ll P_{F_z}(0)$, with error in the approximations consistent with the at-most-5% variation in each flow's mean downslope velocity.

Turning to higher frequencies, Figure 4 presents the power spectra of the basal forces exerted by two simulated flows, chosen to represent the range of variation over our simulations. Spectra are noisy but, since signal-to-noise ratios are comparable for all flows, components, and frequencies, we calculate $\Delta f = \sqrt{g/d}$ moving averages

$$ \langle P_F \rangle_{\Delta f} = \frac{1}{\Delta f \Delta t} \sum_{|j/\Delta t - f| < \Delta f/2} P_F(j/\Delta t) \qquad (23) $$

to note that, for any given flow, $P_{F_z}(f) > P_{F_x}(f) > P_{F_y}(f)$ at most frequencies $f$, with the ratios between these power spectral densities approximating unity with increasing $f$.

Using the timescale $\sqrt{d/g}$ of particle rearrangement under gravity to classify frequencies as low ($f\sqrt{d/g} < 1$), intermediate ($1 < f\sqrt{d/g} < 100$), and high ($100 < f\sqrt{d/g}$), each power spectrum's spectral density generally decreases from one frequency range to the next, so we view the power spectral density in each frequency range as composed of signals from all frequency ranges at least as high. This can be motivated by considering the power spectrum of a superposition of boxcar functions with uncorrelated parameters and varying interval lengths and, from this perspective, the observed power spectra are composed of three types of signal: low-frequency signals, which contribute significantly to $P_{F_x}$ and $P_{F_z}$ for all flows; intermediate-frequency signals, whose contributions to $P_{F_x}$ and $P_{F_z}$ decrease in significance with increasing
<figure><img src="image_6.png" /><figcaption>**Figure 4.** *Power spectra of basal forces.* Thin lines indicate the power spectral densities $P_F(f)$ of the downslope, cross-slope, and base-normal forces exerted by simulated flows on their base, at various frequencies $f$, for $n_x = n_y = 10$ and (a) $\theta = 22.0^\circ$ and $n_z = 8$ and (b) $\theta = 24.5^\circ$ and $n_z = 16$. Thick lines represent simple central moving averages $\langle P_F \rangle_{\Delta f}$, dotted lines indicate the corner frequency $f_c$ and high-frequency spectral amplitude $P_F^0$ of the base-normal component, and black lines indicate the corresponding Hertzian power spectrum $P_F^0 \zeta(f_c f/f_c)$ and its high-frequency asymptote $0.33 P_F^0 (f/f_c)^{-5}$ (see Supplementary Information section S6).</figcaption></figure>
slope angle $\theta$; and high-frequency signals, which contribute significantly to the power spectra of all force components, for all flows. It is the high-frequency signals on which we shall concentrate our attention.

For each power spectrum $P_F$, the high-frequency signal is characterized by a rapid drop-off beyond a corner frequency and, if we extract the corner frequency $f_c$ and corresponding low-frequency amplitude $P_F^0$ for each spectrum $P_{F_z}$, is very well approximated by the signal of a representative Hertzian collision with corresponding $f_c$ and $P_F^0$. On the basis of section 2.3, we define the corner frequency of a power spectrum $P_F(f)$ as the frequency at which a frequency doubling first corresponds to a 10 dB decrease in $\langle P_F \rangle_{\Delta f}$, so that

$$f_c = \min\{f | \langle P_{F_z} \rangle_{\Delta f}(2f) = \langle P_{F_z} \rangle_{\Delta f}(f)/10\}, \quad (24)$$

and we use this to define, for $\zeta$ defined in equation (8) and $f_c$ in equation (9),

$$P_F^0 = \langle P_{F_z} \rangle_{\Delta f}(f_c)/\zeta(f_c). \quad (25)$$

The corresponding power spectrum, for a signal generated by identical, Hertzian inter-particle collisions, is then $P_F^0(\zeta(f_c f/f_c))$, which Figure 4 shows to be a reasonable approximation to the power spectrum of any force component wherever the high-frequency signal is dominant, and a very good approximation to $P_{F_z}$ at high frequencies.

Having verified that each simulation's power spectra $P_F$ are consistent with a high-frequency signal generated by interparticle collisions, and extracted each signal's low-frequency amplitude $P_F^0$ and corner frequency $f_c$, we can compare these values to the predictions of 'stochastic impact' models for such signals.

## 3.3 Model comparison

To compare the accuracy of the models described in section 2.3, we extract predictions $\hat{P}_F^0$ for the low-frequency amplitude of the 'stochastic impact' signal. For the models of Kean et al. (2015) and Lai et al. (2018), we calculate these predictions directly while, for the models of Farin, Tsai, et al. (2019) and Bachelet et al. and when using the profiles of collisional properties within our minimal model, we extract the predicted corner frequencies $\hat{f}_c$ and values $\hat{P}_F^0$ from the predicted power spectra $\hat{P}_{F_z}$, using the exact equivalents of equations (24) and (25). If a model has a free parameter, we use the parameter value that minimizes the sum over all simulations of $\ln(P_F^0/\hat{P}_F^0)^2$, with the resulting best-fit values listed in Table 1.

Figures 5a-e,g-i compare observed values of $P_F^0$ and $f_c$ to the predictions of previous models, showing that the 'thin-flow' model of Farin, Tsai, et al. (2019) is poor at predicting power spectral amplitudes for our simulated flows and that, instead, the most accurate predictions are those of the same paper's 'thick-flow' model. The 'thin-flow' model and the model of Lai et al. (2018), which suppose that the impulse transferred by a representative basal impact is proportional to the mean downslope velocity $\bar{u}$, overestimate the basal force's spectral amplitude $P_F^0$ for thick (high-$n_z$) flows relative to that for thin flows, and the same is true for the model of Kean et al. (2015), which supposes impact impulses to be proportional to the mean basal pressure $p \propto n_z$. In contrast, the 'thick-flow' model and the model of Bachelet et al. assume representative impulses proportional to the local shear rate, of order $\bar{u}/dn_z$, and make predictions that accurately capture $P_F^0$'s dependence on $n_z$, differing from observations by typical factors of only 1.7 and 2.0 respectively, over two orders of magnitude of variation of $P_F^0$.

However, the 'thick-flow' and Bachelet et al. models' predictions for $P_F^0$ are accurate only if model parameters take unrealistic values. Figures 3b,c,h,i indicate that
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Comparison between models' predictions (x-axes) and simulation outputs (y-axes) for the low-frequency amplitudes (a-f) and corner frequencies (g-j) of high-frequency signals. The gray lines represent perfect agreement, while marker size indicates measurement area A (with $A/d^2$ = 250 for channel flows) and colors indicate flow overburden $n_z$. For each $n_z$, $P_F^0$ and $f_c$ increase with $\theta$.</figcaption></figure>
<table><caption>**Table 1.** *Comparison of models.* For each model whose predictions are described in section 2.3 and for the semi-empirical model of section 3.4, we list the equation used to calculate those predictions, the free parameter value for which such predictions best fit measurements, the geometric standard error of these predictions, $\epsilon = \exp[\sqrt{\frac{1}{N}\sum \ln(P_F^0/\hat{P}_F^0)^2}]$, and the value of the Akaike information criterion (AIC), calculated as in Arran et al. (2021)'s section S9.3.</caption><thead><tr><th>Model</th><th>Equation</th><th>Best-fit parameter</th><th>ε</th><th>AIC</th></tr></thead><tbody><tr><td>Kean et al. (2015)</td><td>(10)</td><td>K = 0.015d<sup>5</sup>/g</td><td>3.6</td><td>102</td></tr><tr><td>Lai et al. (2018)</td><td>(11)</td><td>n/a</td><td>9.0</td><td>132</td></tr><tr><td>Farin, Tsai, et al. (2019)</td><td>(13)</td><td></td><td></td><td></td></tr><tr><td>'thick-flow'</td><td></td><td>v = 21</td><td>1.7</td><td>47</td></tr><tr><td>'thin-flow'</td><td></td><td>v = 1.1</td><td>4.6</td><td>112</td></tr><tr><td>Bachelet et al.</td><td>(14)</td><td>γ = 0/d</td><td>1.9</td><td>59</td></tr><tr><td>Minimal model</td><td>(16)</td><td>n/a</td><td>1.4</td><td>14</td></tr><tr><td>Semi-empirical model</td><td>(27)</td><td>(C, a, b) = (9.8, 2.08, -2.70)</td><td>1.13</td><td>-33</td></tr></tbody></table>

$\sqrt{3T}/u_x$, the relative standard deviation of particle velocities, is typically much less than one. But the 'thick-flow' model's predictions for $P_F^0$ are of the correct order of magnitude only when the corresponding parameter $v > 10$, resulting in strong overestimates of impact velocities and so of corner frequencies $f_c$, as observed in Figure 5g. Similarly, Bachelet et al.'s model's predictions for $P_F^0$ are only of the same magnitudes as observations if the attenuation constant $\gamma$ is taken to be much less than $1/\max(h)$, implying that all forces at a flow's surface are transmitted unaltered to the flow's base. Such values of $\gamma$ are unrealistic, especially given Bachelet et al. (2018)'s measurements of $\gamma = 0.44/d$ in static bead packs, and result in Figure 5i's inflated predictions $\hat{f}_c$ for low-$f_c$ flows, within which any impact velocities associated with near-surface shear are much higher than those associated with basal shear.

To recover accurate predictions without unrealistic fitting parameters, we turn to the profiles of collisional properties. From these profiles, without any fitting parameter, the minimal model introduced in section 2.3 predicts $P_F^0$ to within a factor of 2.0 (Figure 5f), with geometric standard error 1.4, while underestimating $f_c$ by factors between 1.1 and 1.3 (Figure 5j). With all quantities having been non-dimensionalized, this accuracy holds for any values of particle density $\rho$ and mean diameter $d$, as well as over factor-4 variation in surface area $A$ and overburden $n_z$, and the entire range of slope angles $\theta$ for which steady flow can be sustained. We believe the main cause of corner frequencies' underestimation to be that, at each base-normal position $z$, high-velocity collisions make contributions to the very-high-frequency component of $P_F$ that exceed those of even a suitably weighted representative impact at velocity $\sqrt{\langle\delta_n^2\rangle}(z)$.

### 3.4 Parametrization of the minimal model

For the minimal model's accuracy to be useful for physical flows, within which collisional properties are very rarely accessible, the contributions in equation (16) must be described in terms of bulk flow parameters.

Figures 6a,b demonstrate the difficulty of such parametrization, as each predicted high-frequency signal $\hat{P}_F$ is dominated by contributions from a basal layer in which collisional properties are hard to model. For $0 < z_0 < h$, $\int_{z_0}^{\infty} \phi dz$ and the volumetric impact rate $n_I$ are well approximated by $\bar{\phi}h(1 - z_0/h)$ and equation (21)'s kinetic-
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Attempts at parametrization of the minimal model. a, b) Approximants to the inte-grand $\mathcal{I}$ of equation (17), for $n_x = n_y = 10$ and (a) $\theta = 22.0^{\circ}$ and $n_z = 8$ and (b) $\theta = 24.5^{\circ}$ and $n_z = 16$. $\mathcal{I}_{\hat{n}}$, $\mathcal{I}_{\hat{P}_c}$, and $\mathcal{I}_T$ are defined by equation (26) with, respectively and cumulatively: no changes, $\mathcal{P}_c$ replaced with equation (19)'s $\hat{P}_c$, and $\langle \dot{\delta}_n^2 \rangle$ replaced with $2T$. $\mathcal{I}_{\hat{\psi}\hat{V}}$ is defined by re-placing, in equation (26), $\psi(z)\sqrt{\langle \dot{\delta}_n^2 \rangle}$ with equations (29)'s $\hat{\psi}\hat{V}$. c, d) Equivalents of Figures 5f,j, except for the semi-empirical model of equation (27) rather than the minimal model of equation (16).</figcaption></figure>

theoretic $\hat{n}_I$, so that the integrand $\mathcal{I}$ of equation (17) approximately equals

$$ \mathcal{I}_{\hat{n}}(f, z) = \frac{24g}{\pi d^4} \bar{\phi} h \cos \theta \left(1 - \frac{z}{h}\right) (1 - \mathcal{P}_c(z)) \mathcal{P}_c\left(z - \frac{1}{4}d\right) \sqrt{\langle \dot{\delta}_n^2 \rangle} \zeta \left(f \tau \left(\sqrt{\langle \dot{\delta}_n^2 \rangle}\right)\right) e^{-\frac{z \ln 9}{d}}. \quad (26) $$

Furthermore, the probability $\mathcal{P}_c$ that a particle is connected to the base has profiles that are reasonably well parametrized by decay lengths $\lambda$, using equation (19). But velocity correlations prevent $\langle \dot{\delta}_n^2 \rangle / 2$ from being consistently well approximated by granular temperature $T$ (as would be the case for ‘molecular chaos’) or by the estimator $\hat{T} \propto (du_x')^2$ of equation (18). Since the models of which we’re aware approximate even less well the near-base granular temperature and shear rates $u_x'(z)$ (see Arran et al. (2023) for details), the simple theoretical considerations in this paper do not allow $\mathcal{I}_{\hat{n}}$ to be parametrized.

We therefore use a semi-empirical parametrization. Since the integral in equation (16) is dominated by contributions from a basal layer whose width is limited by geometric attenuation, then we can use equation (26) to derive the prediction

$$ \hat{P}_{F_z}(f) = \frac{4\pi}{3\ln 3} \rho^2 d^3 g \bar{\phi} A h \cos \theta \hat{\psi} \hat{V} \zeta(f \tau(\hat{V})), \quad (27) $$

for some representative impact velocity $\hat{V}$ and representative value $\hat{\psi}$ of the profile $\psi(z) = (1 - \frac{z}{h}) (1 - \mathcal{P}_c(z)) \mathcal{P}_c(z - \frac{1}{4}d)$. Then, supposing that these are linked to bulk flow parameters and standard lognormal error $\epsilon_j$ by

$$ \ln \left( \frac{\hat{\psi} \hat{V}}{\sqrt{gd}} \right) = \ln C + a \ln \left( \bar{u} / \sqrt{gd} \right) + b \ln (h/d) + \ln \epsilon_j, \quad (28) $$

we use a linear regression to calculate the parameters $C = 9.8e^{\pm 0.16}$, $a = 2.08 \pm 0.05$ and $b = -2.70 \pm 0.08$ that best fit low-frequency predictions $\hat{P}_{F_z}(0)$ to observations $\mathcal{P}_c^0$. Finally, approximating $\psi(z)$ by $\psi(0)$ and $\mathcal{P}_c$ by equation (19)'s $\hat{P}_c$, we recover
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Profiles (a, c) and colorplots (b, d) of mean downslope velocities in simulated channels, for either (a, b) $\theta = 22^{\circ}$ and $n_z = 16$ or (c, d) $\theta = 24.5^{\circ}$ and $n_z = 4$. In a and c, $u(z)$ is the average over the inner region (right of the dashed line in b, d), while $u_0(z)$ is the velocity profile next to the simulated sidewall (red line in b, d), and $u_w(z)$ is a matching profile observed at the experimental sidewall of Arran et al. (2021), with (a) $\theta = 23.7^{\circ}$ and release gate height $h_g = 20d$ and (c) $\theta = 24.7^{\circ}$ and release gate height $h_g = 5d$. Horizontal dotted lines indicate flow depths $h$.</figcaption></figure>

predictions

$$ \hat{\psi} = (1 - e^{-1/2\lambda})e^{-1/4\lambda}, \quad \hat{V} = \frac{Ce^{1/4\lambda}}{(1 - e^{-1/2\lambda})} \left(\frac{\bar{u}}{\sqrt{gd}}\right)^a \left(\frac{h}{d}\right)^b, \quad (29) $$

which may be substituted into equation (27). The accuracy of this semi-empirical parametrization is demonstrated by Figures 6c,d, with predictions of $P_F^0$ having geometric standard error 1.13. Such accuracy suggests some degree of overfitting, with equation (27) making more accurate predictions than the equation (16) that it approximates, but the additional fitting parameters $C$, $a$, and $b$ are justified by the Akaike information criterion (see Table 1).

# 4 Discussion

## 4.1 Comparison to Arran et al. (2021)

Arran et al. (2021) describes laboratory experiments that closely match our simulations in several respects. Our simulations accurately reproduce the interparticle interactions responsible for seismic signal generation in experiments, with the properties of our simulated particles corresponding to those of the experimental particles. The experiments' range of slope angles $\theta$ and flow depths $h$ overlaps with this study's, leading to velocity profiles at the experimental channel's sidewall that are close to those observed at the corresponding boundary of the simulated channel (see Figure 7a,c). And Arran et al. (2021) measures basal forces' power spectral densities over a frequency range $30 < f\sqrt{d/g} < 1700$, so considers the same high-frequency signals as this study.

However, the relationship this study observes between a flow's bulk parameters and its signal's properties is substantially different from that observed in Arran et al. (2021). At $\theta \le 23^{\circ}$, we observe distinct intermediate-frequency signals, at $f\sqrt{d/g} < 100$, which have no equivalent in Arran et al. (2021) (as illustrated by comparing our
Figure 4a and Arran et al. (2021)'s Figure 4). Concerning the high-frequency signals, Arran et al. (2021) found power spectral densities to be most accurately predicted by Farin, Tsai, et al. (2019)'s 'thin-flow' model, in which basal particles' impact velocities are proportional to the depth-averaged flow velocity $\bar{u}$, whilst this study finds Farin, Tsai, et al. (2019)'s 'thick-flow' model to be more accurate and infers representative impact velocities that scale as $\bar{u}^{2.08\pm0.05}/h^{2.70\pm0.08}$.

These differences suggest that the nature of simulated flows differed from that of the experimental flows, and that basal slip is the explanation for the accuracy of the 'thin-flow' model in Arran et al. (2021). Velocities in the interior of Arran et al. (2021)'s experimental channel could not be directly measured, and Figure 7 demonstrates that measurements at a sidewall can be unrepresentative, so experimental interior velocity profiles could have differed substantially both from velocity profiles at the experimental sidewall and from those in simulated flow interiors. That experimental flows exhibited basal slip is therefore possible, despite its absence in simulations; feasible, given the basal vibrations that were documented in Arran et al. (2021)'s experiments but absent in simulations; and probable, for near-base impact velocities to be consistent with the associated basal forces in each case. If basal particle velocities scaled with the depth-averaged flow velocity in experiments, then the accuracy of the 'thin-flow' model in experiments is entirely concordant with its inaccuracy for our simulations, in which basal particle velocities did not.

## 4.2 The inertial number and the seismic signal

With the 'thin-flow' model inapplicable to our simulations, we examine the applicability of another of Arran et al. (2021)'s results: a strong relationship between a non-dimensionalized shear rate $\hat{I}$, the bulk inertial number, and a non-dimensionalized mean squared magnitude of fluctuating forces $\delta F^2$, on which the seismic signal depends.

The inertial number describes how energetic a granular flow is, and in an idealized flow determines all other local, non-dimensional quantities (da Cruz et al., 2005; GDR MiDi, 2004). For flows such as our simulations', with particle diameter $d$, slope angle $\theta$, mean flow depth $h$, and mean depth-averaged volume fraction $\bar{\phi}$ and downslope velocity $\bar{u}$, Jop et al. (2005) derives the estimate used by Arran et al. (2021):

$$ \hat{I} = \frac{5\bar{u}d}{2h\sqrt{\bar{\phi}gh\cos\theta}}. \qquad (30) $$

Figure 8 indicates that, in simulations, $\hat{I}$ so-defined is an imperfect but strong determinant of other estimates of local, non-dimensional quantities. The slope gradient $\tan\theta$ is the steady-state ratio between local base-tangential and base-normal forces, and has a relationship with $\hat{I}$ that, though dependent on the number $n_z$ of overlying particles, is both strong and consistent with that observed in Arran et al. (2021)'s experiments (Figure 8a). The contact network's non-dimensional decay length $\lambda$, as defined in section 3.1, reflects the local probability that a particle is in contact with an underlying counterpart, and decreases dramatically as this probability tends from 1 at $\hat{I} = 0$ towards 0 as $\hat{I} \to \infty$ (Figure 8b). That $\hat{I}$ almost exactly determines $\lambda$ suggests that the bulk inertial number is appropriately defined, despite deviations from idealized flow behavior.

The quantity $\delta F^2$, meanwhile, requires redefinition when, unlike in Arran et al. (2021)'s experiments, the area $A$ over which it is measured varies. We first recall that $p$ represents basal pressure, $P_{F_z}(f)$ the base-normal basal force's power spectrum, and $P_F^0$ and $f_c$ the amplitude and corner frequency of $P_{F_z}$'s high-frequency signal. Now, the model performance represented by Figure 5 confirms that $P_F^0 \propto A$ with
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Relations to the bulk inertial number $\hat{I}$ of other non-dimensional quantities: the slope gradient $\tan \theta$ (a), contact network decay length $\lambda$ (b), and squared ratio of high-frequency force fluctuations to mean forces $\delta \mathcal{F}^2$ (c). Marker size indicates measurement area A, colors indicate flow overburden $n_z$, and unfilled points represent experimental results from Arran et al. (2021), with marker size reduced and ‘transitional-regime’ flows excluded for clarity.</figcaption></figure>

$f_c$ independent, as predicted by the stochastic impact framework, while equations (3) and (22) confirm the implication of a static force balance, that basal pressure $p$ is independent of A. Therefore, we expect $P_F^0 f_c/(Ap)^2 \propto 1/A$ and, to recover an A-independent local-variable estimate, exchange Arran et al. (2021)'s ratio between squared forces on area A for a ratio between the inferred squared forces on a particle area $d^2$. Consequently, we define

$$ \delta \mathcal{F}^2 = \frac{2P_F^0 f_c}{Ad^2 p^2} \qquad (31) $$

and examine its relation to $\hat{I}$.

Figure 8c demonstrates that, as for Arran et al. (2021), $\delta \mathcal{F}^2$ is strongly correlated with $\hat{I}$. Over an order of magnitude of change in simulations’ $\hat{I}$, the power-law fit $\delta \mathcal{F}^2 \propto \hat{I}^{2.4}$ has a geometric standard error of 1.3. Furthermore, the experimental and simulation-derived relationships agree over much of the range of $\hat{I}$ where they overlap, despite the differing flow behavior discussed in section 4.1. Observed disagreement, at the low end of $\hat{I}$’s experimental range, may be due to that differing flow behavior or to experimental anomalies.

Given the close relationship between a flow’s bulk inertial number $\hat{I}$ and its high-frequency signal, we investigate its relationship with the intermediate-frequency signals introduced in section 3.2.

## 4.3 Signals at intermediate frequencies

The intermediate-frequency signals we observe in e.g. Figure 4a, contributing to basal-force power spectral densities $P_{F_x}(f)$ and $P_{F_z}(f)$ at frequencies $f < 100\sqrt{g/d}$, can be explained neither by the single-particle, basal impacts considered by Kean et al. (2015), Lai et al. (2018) and Farin, Tsai, et al. (2019), nor by the coherent motion with which the single-force framework explains low-frequency signals. The forces exerted during binary collisions make contributions to the power spectra $P_F$ that are approximately constant over all frequencies less than the collision timescale’s associated corner frequency, equations (7) and (9)’s $f_c/\tau(\delta_n)$, which in our case is less than $100\sqrt{g/d}$ only for the rare and insignificant impacts at velocities $\delta_n < 1.3 \times 10^{-8}\sqrt{g/d}$. But signal frequencies significantly greater than $\sqrt{g/d}$ cannot be due to coherent mo-
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Association of force chains and intermediate-frequency signals. a, b) Taking flows with $n_x = n_y = 10$, $n_z = 8$ and either $\theta = 22^\circ$ (a) or $\theta = 24.5^\circ$ (b), at the times $t_0$ at which each is considered to be in a steady state, each pair of particles in contact is indicated by a line segment joining their centres. Opacity and color correspond to the interparticle normal force $F_n$, with case-b forces frequently exceeding the saturation value of $400\rho g d^3$. c, d, e) Ratios of smoothed power spectral densities $\langle P_{F_j} \rangle_{\Delta f}(f)$ of downslope ($j = x$), cross-slope ($j = y$), and base-normal ($j = z$) basal forces, at intermediate ($f = 10\sqrt{g/d}$) and high ($f = f_c$, $P_F^0 \propto \langle P_{F_z} \rangle_{\Delta f}(f_c)$) frequencies, as functions of force chain decay length $\lambda$, bulk inertial number $\hat{I}$, and slope gradient $\tan\theta$. Marker size indicates measurement area A and colors indicate flow overburden $n_z$.</figcaption></figure>

tion, whose minimum timescale is that of gravitational acceleration over the flow's roughness-d base.

Instead, intermediate-frequency signals appear to be associated with the dynamics within a flow of multi-particle chains, along which forces propagate: 'force chains'. Figures 9a,b illustrate, at a given instant, the interparticle forces within two flows with power spectra corresponding to those in Figures 4a,b, demonstrating that the relative magnitude of intermediate-frequency signals corresponds to the prevalence of such force chains. Figure 9c shows that this relation holds across all flows, with a Spearman rank correlation coefficient of 0.98 between the ratio $\langle P_{F_z} \rangle_{\Delta f}(10\sqrt{g/d})/P_F^0$, a measure of the relative magnitude, and the decay length $\lambda$, which by equation (19) approximates the typical base-normal extent of basal force chains.

This relationship is unlikely to be due to the single-body force-chain impacts considered by Zhang et al. (2021). If $n$ density-$\rho$ particles in a force chain impact the base as a single body, at velocity $\hat{v}_n$, then the duration of the collision will be
that which equation (7) specifies for a single density-$\rho n$ particle's impact. The corner
frequency associated with this impact will therefore be lower than that associated
with an equivalent single-particle impact, but even for a 10-particle chain will only be
less than $100\sqrt{g/d}$ for the still-rare impacts with $\dot{\delta}_n < 1.3 \times 10^{-6}\sqrt{gd}$. Furthermore,
no single-body impacts of force chains are evident in Movie S1 and Movie S2 of the
supplementary material, from which Figures 9a and 9b are frames.

Intermediate-frequency signals are more likely induced by force chains’ transmission of force to each flow’s base. The supplementary videos show that contact networks evolve over timescales between those of individual collisions and of particle rearrangement, altering the transmission of prolonged-contact forces, and this could induce basal-force signals at the corresponding intermediate frequencies. Less probably, if transmission attenuates the high-frequency components of forces more than the intermediate-frequency components (contrary to the assumptions of Bachelet et al. (2023) and our minimal model), then the transmission of non-basal, single-particle impact forces could be responsible for relative magnitudes $\langle P_{F_z} \rangle_{\Delta f}(10\sqrt{g/d}/P_F^0 > 1$. Since Hertz theory predicts power spectra’s high-frequency components well, such frequency-dependent attenuation would have to be constant for $f > 100\sqrt{g/d}$, but this explanation cannot be ruled out without further work.

Even without explaining intermediate-frequency signals, we can describe their
properties. Figure 9d shows the relative magnitude of intermediate-frequency signals
to decrease with the bulk inertial number $\dot{I}$ and increase slightly with overburden
$n_z$, while demonstrating that intermediate-frequency signals have no significant cross-
slope component ($\langle P_{F_y} \rangle_{\Delta f}(10\sqrt{g/d})/P_F^0 \approx 0.3$, regardless of intermediate-frequency
signals' magnitude). Figure 9e, meanwhile, examines the relative magnitudes of inter-
mediate signals' downslope and base-normal components: while $P_{F_x} \approx 0.6P_{F_z}$ at high
frequencies, intermediate-frequency signals appear to contribute to $\langle P_{F_x} \rangle_{\Delta f}(10\sqrt{g/d})$
approximately $\tan^2 \theta$ as much as to $\langle P_{F_z} \rangle_{\Delta f}(10\sqrt{g/d})$. This may reflect the difference
between the orientations of single-particle impacts, which are widely distributed, and
those of force chains, which are preferentially aligned with the mean force (Campbell,
2006; Thomas & Vriend, 2019). Furthermore, all aforementioned relative magnitudes
are independent of A, indicating that spatially separated contributions to intermediate-
frequency signals are as uncorrelated as the individual impact forces that contribute
to high-frequency signals.

**5 Conclusion**

With the properties of intermediate-frequency signals described, we conclude by
reviewing our work. We simulated steady flows of spherical, diameter-d particles under
gravitational acceleration *g* and over rough bases towards which downslope velocities
approached zero. Slope angles spanned the range that supported steady flow and flow
depths *h* were between 4*d* and 15*d* (see Figure 2). Writing *F*<sub>*z*</sub> for the base-normal
force exerted by a flow upon its base, we considered its power spectrum *P*<sub>*F*<sub>*z*</sub></sub>(*f*) as the
sum of three components: a low-frequency component which vanishes for frequencies
*f* > √*g*/*d*, an intermediate-frequency component which is constant for *f* < √*g*/*d*
and vanishes for *f* > 100√*g*/*d*, and a high-frequency component which is constant
for *f* < 100√*g*/*d* and vanishes as *f* → ∞ (see Figure 4). Various models use a
flow's properties to predict the properties of this high-frequency component, resulting
from the forces exerted during individual interparticle collisions, and we compared our
observations to these predictions.

Our results differ from those of Arran et al. (2021)'s laboratory experiments, to
which our simulations closely correspond, and we link these differences to differences
in basal flow velocities. While the most accurate predictions in the laboratory were
those of Farin, Tsai, et al. (2019)'s 'thin-flow' model, in which basal impact velocities
vary around the depth-averaged flow velocity $\bar{u}$, we find the most accurate pre-existing model to be the corresponding ‘thick-flow’ model, in which base-impacting particles’ velocities are approximately $\bar{u}d/h$ (see Table 1 and Figure 5). This supports the argument that, in laboratory flows, vibration-induced friction weakening led to basal slip that wasn’t observed in our fixed-base simulations. Instead, at low slope angles, simulated flows exhibited prolonged interparticle contacts, and force chains of multiple such contacts were associated with intermediate-frequency components of $P_{F_z}$ that were absent in the laboratory.

However, we propose generally applicable relations between the properties of a flow and those of the basal forces it exerts. We derive a ‘minimal model’ that, from depth profiles of particle collisions’ properties but without fitting parameters, predicts the high-frequency component of $P_{F_z}$ more accurately than any pre-existing model. Our empirical parametrization, replacing the minimal model’s inputs with $\bar{u}$ and $h$, works well for our simulated flows (see Figure 6), and other parametrizations will apply to flows with slip or slow creep at their base. Furthermore, we demonstrate a relation between the bulk inertial number $\hat{I}$ and the ratio $\delta\mathcal{F}^2$ of mean squared fluctuating and mean forces on a basal area $d^2$, consistent with the relation described in Arran et al. (2021).

These relations are not applicable to geophysical flows without further work. Firstly, geophysical signals are considered ‘low-frequency’ up to an approximately 0.1 Hz threshold, below which spatially separated forces are strongly correlated, whereas we consider the ‘low-frequency’ component of $P_{F_z}$ to vary up to a $\sqrt{g/d}$ threshold, above which such forces are almost uncorrelated, and the ‘high-frequency’ component to be that with a collision-related corner frequency of order given by section 2.3:

$$f_c \approx f_c/\tau(0.1\sqrt{gd}) \approx 0.11[E/\rho(1-\nu^2)]^{2/5}[g/d^9]^{1/10}, \quad (32)$$

for material density $\rho$, Young’s modulus $E$, and Poisson’s ratio $\nu$. These frequencies correspond to 70 Hz and 40 kHz for the experimental particles of Arran et al. (2021) but 3 Hz and 140 Hz for 1 m-diameter blocks of rock with similar properties. Secondly, changes in material properties will change flow parameters and hence signal frequencies, as well as the relation $\delta\mathcal{F}^2(\hat{I})$. And finally, we simulate dense, dry flows of monodisperse spheres, so our results are not directly applicable to saltating particles at or ahead of a geophysical flow’s front, to collisions significantly affected by a fluid medium or by asphericity, or to flows altered by polydisperse particles’ segregation by size.

We nevertheless suggest areas of application. Given the differences between our results and Arran et al. (2021)’s, seismic signals can indicate the basal velocities within an experimental flow, which are otherwise hard to measure. Alternatively, one can combine our results with Arran et al. (2021)’s to estimate the high-frequency seismic signals produced by flows for which the basal boundary condition is known: use Farin, Tsai, et al. (2019)’s ‘thin-flow’ model if there is basal slip and the empirical parametrization of our ‘minimal model’ if not. Finally, even if the conditions for $\hat{I}$ and $\delta\mathcal{F}^2$’s relation are unknown, we show that the same relation applies to purely granular flows with quite different velocity profiles, while Piantini et al. (2023) suggests that an analogous relation underlies a different observation: an inverse relationship between the solid concentrations of laboratory-scale sand-and-cobble sediment flows and the fluctuating forces they exert. Consequently, $\hat{I}$ may be useful for predicting the intermediate and high-frequency force signals generated by a broad range of geophysical granular flows.
# 6 Open Research

V1.0 of the SimSlidequakesMercuryDPM software (Arran, 2023), forked from MercuryDPM (https://www.mercurydpm.org, Weinhart et al., 2020) and used for DEM simulations, is preserved at https://doi.org/10.5281/zenodo.8125567 and available under an attribution licence.

The Jupyter Notebook to execute the analysis in the paper can be found at https://mybinder.org/v2/g1/M_Arran\%2Fsimulated-slidequakes/HEAD?labpath=\%2Fsimulated slidequakes.ipynb and is hosted with its input data, scripts for their creation, and output figures at https://gitlab.com/M_Arran/simulated-slidequakes. This repository is preserved along with datafiles and scripts for using SimSlidequakesMercuryDPM at https://doi.org/10.5281/zenodo.8219510 (Arran et al., 2023).

# Acknowledgments

This work was primarily funded by project ERC-CG-2013-PE10-617472 SLIDEQUAKES, with supporting funds from La Société des Amis de l'ESPCI and from IPGP. RT acknowledges the support of the France-Norway IRP D-FFRACT programs, and of the Research Council of Norway through its Centres of Excellence funding scheme, Project No. 262644. The authors are aware of no conflict of interest. Numerical computations were performed on the S-CAPAD/DANTE platform, IPGP, France, with code using NumPy (Harris et al., 2020), Matplotlib (Hunter, 2007), and pandas (McKinney, 2010).

# References

Ai, J., Chen, J.-F., Rotter, J. M., & Ooi, J. Y. (2011). Assessment of rolling resistance models in discrete element simulations. *Powder Technology, 206*(3), 269-282. Retrieved from http://www.sciencedirect.com/science/article/pii/S0032591010005164 doi: 10.1016/j.powtec.2010.09.030

Allen, M., & Tildesley, D. (1989). *Computer simulation of liquids*. Clarendon Press. Retrieved from https://books.google.co.uk/books?id=032VXB9e5P4C

Allstadt, K. E. (2013). Extracting source characteristics and dynamics of the august 2010 mount meager landslide from broadband seismograms. *Journal of Geophysical Research: Earth Surface, 118*(3), 1472-1490. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/jgrf.20110 doi: 10.1002/jgrf.20110

Allstadt, K. E., Farin, M., Iverson, R. M., Obryk, M. K., Kean, J. W., Tsai, V. C., ... Logan, M. (2020). Measuring basal force fluctuations of debris flows using seismic recordings and empirical green's functions. *Journal of Geophysical Research: Earth Surface, 125*(9), e2020JF005590. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2020JF005590 (e2020JF005590 2020JF005590) doi: 10.1029/2020JF005590

Allstadt, K. E., Matoza, R. S., Lockhart, A. B., Moran, S. C., Caplan-Auerbach, J., Haney, M. M., ... Malone, S. D. (2018). Seismic and acoustic signatures of surficial mass movements at volcanoes. *Journal of Volcanology and Geothermal Research, 364*, 76-106. Retrieved from http://www.sciencedirect.com/science/article/pii/S0377027317306261 doi: 10.1016/j.jvolgeores.2018.09.007

Arran, M. I. (2023). *Simslidequakesmercurydpm*. Zenodo. doi: 10.5281/zenodo.8125567

Arran, M. I., Mangeney, A., De Rosny, J., Farin, M., Toussaint, R., & Roche, O. (2021). Laboratory landquakes: Insights from experiments into the high-frequency seismic signal generated by geophysical granular flows. *Journal of Geophysical Research: Earth Surface, 126*(5), e2021JF006172. Retrieved
from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/
2021JF006172 (e2021JF006172 2021JF006172) doi: 10.1029/2021JF006172

Arran, M. I., Mangeney, A., De Rosny, J., & Toussaint, R. (2023). *Simulated slide- quakes, v1.0*. Zenodo. doi: 10.5281/zenodo.8219510

Bachelet, V. (2018). *Étude expérimentale des émissions acoustiques générées par les écoulements granulaires* (Doctoral dissertation, Institut de Physique du Globe de Paris). Retrieved from https://hal.archives-ouvertes.fr/ tel-03104582

Bachelet, V., Mangeney, A., de Rosny, J., Toussaint, R., & Farin, M. (2018). Elastic wave generated by granular impact on rough and erodible surfaces. *Journal of Applied Physics*, 123(4), 044901. Retrieved from https://doi.org/10.1063/1.5012979 doi: 10.1063/1.5012979

Bachelet, V., Mangeney, A., Toussaint, R., DeRosny, J., Arran, M. I., Farin, M., & Hibert, C. (2023). Acoustic emissions of nearly steady and uniform granular flows: a proxy for flow dynamics and velocity fluctuations. *Journal of Geophysical Research: Earth Surface, (Accepted)*. Retrieved from https://arxiv.org/pdf/2101.04161.pdf

Bardet, J., & Huang, Q. (1993). *Rotational stiffness of cylindrical particle contacts*.

Brilliantov, N. V., Pimenova, A. V., & Goldobin, D. S. (2015). A dissipative force between colliding viscoelastic bodies: Rigorous approach. *EPL (Europhysics Letters)*, *109*(1), 14005. Retrieved from https://doi.org/10.1209/0295-5075/109/14005 doi: 10.1209/0295-5075/109/14005

Brilliantov, N. V., & Pöschel, T. (1998, June). Rolling friction of a viscous sphere on a hard plane. *Europhysics Letters (EPL)*, *42*(5), 511–516. Retrieved from https://doi.org/10.1209%2Fep1%2Fi1998-00281-7 doi: 10.1209/epl/i1998-00281-7

Brodsky, E. E., Gordeev, E., & Kanamori, H. (2003). Landslide basal friction as measured by seismic waves. *Geophysical Research Letters*, *30*(24). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2003GL018485 doi: 10.1029/2003GL018485

Campbell, C. S. (2006). Granular material flows – an overview. *Powder Technology*, *162*(3), 208–229. Retrieved from https://www.sciencedirect.com/science/article/pii/S0032591005005486 doi: 10.1016/j.powtec.2005.12.008

Chao, W.-A., Wu, Y.-M., Zhao, L., Chen, H., Chen, Y.-G., Chang, J.-M., & Lin, C.-M. (2017). A first near real-time seismology-based landquake monitoring system. *Scientific Reports*, *7*, 43510. doi: 10.1038/srep43510

Coe, J. A., Baum, R. L., Allstadt, K. E., Kochevar, J., Bernard F., Schmitt, R. G., Morgan, M. L., ... Kean, J. W. (2016, 04). Rock-avalanche dynamics revealed by large-scale field mapping and seismic signals at a highly mobile avalanche in the West Salt Creek valley, western Colorado. *Geosphere*, *12*(2), 607–631. Retrieved from https://doi.org/10.1130/GES01265.1 doi: 10.1130/GES01265.1

Cundall, P. A. (1971). A computer model for simulating progressive, large-scale movement in blocky rock system. *Proceedings of the International Symposium on Rock Mechanics*, 1971. Retrieved from https://ci.nii.ac.jp/naid/10018723276/en/

Cundall, P. A., & Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. *Géotechnique*, *29*(1), 47–65. Retrieved from https://doi.org/10.1680/geot.1979.29.1.47 doi: 10.1680/geot.1979.29.1.47

da Cruz, F., Emam, S., Prochnow, M., Roux, J.-N., & Chevoir, F. (2005). Rheophysics of dense granular materials: Discrete simulation of plane shear flows. *Physical Review E*, *72*, 021309. Retrieved from https://link.aps.org/doi/10.1103/PhysRevE.72.021309 doi: 10.1103/PhysRevE.72.021309

Dahlen, F. A. (1993). Single-force representation of shallow landslide sources. *Bulletin of the Seismological Society of America*, *83*(1), 130. Retrieved from
http://dx.doi.org/

Daraio, D., Villoria, J., Ingram, A., Alexiadis, A., Hugh Stitt, E., & Marigo, M. (2019). Validation of a Discrete Element Method (DEM) model of the grinding media dynamics within an attritor mill using Positron Emission Particle Tracking (PEPT) measurements. *Applied Sciences, 9*(22). Retrieved from https://www.mdpi.com/2076-3417/9/22/4816 doi: 10.3390/app9224816

Deparis, J., Jongmans, D., Cotton, F., Baillet, L., Thouvenot, F., & Hantz, D. (2008, 08). Analysis of Rock-Fall and Rock-Fall Avalanche Seismograms in the French Alps. *Bulletin of the Seismological Society of America, 98*(4), 1781-1796. Retrieved from https://doi.org/10.1785/0120070082 doi: 10.1785/0120070082

Dintwa, E., Zeebroeck, M. V., Tijskens, E., & Ramon, H. (2005, July 01). Torsion of viscoelastic spheres in contact. *Granular Matter, 7*(2), 169-179. Retrieved from https://doi.org/10.1007/s10035-005-0207-0 doi: 10.1007/s10035-005-0207-0

Eissler, H., & Kanamori, H. (1987, 05). A single-force model for the 1975 Kalapana, Hawaii, earthquake. *Journal of Geophysical Research, 92*, 4827-4836. doi: 10.1029/1987JB092iB06p04827

Ekström, G., & Stark, C. P. (2013). Simple scaling of catastrophic landslide dynamics. *Science, 339*(6126), 1416-1419. Retrieved from http://science.sciencemag.org/content/339/6126/1416 doi: 10.1126/science.1232887

Farin, M., Mangeney, A., de Rosny, J., Toussaint, R., & Trinh, P.-T. (2018). Link between the dynamics of granular flows and the generated seismic signal: Insights from laboratory experiments. *Journal of Geophysical Research: Earth Surface, 123*(6), 1407-1429. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2017JF004296 doi: 10.1029/2017JF004296

Farin, M., Mangeney, A., de Rosny, J., Toussaint, R., & Trinh, P.-T. (2019). Relations between the characteristics of granular column collapses and resultant high-frequency seismic signals. *Journal of Geophysical Research: Earth Surface, 124*(12), 2987-3021. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2019JF005258 doi: 10.1029/2019JF005258

Farin, M., Tsai, V. C., Lamb, M. P., & Allstadt, K. E. (2019). A physical model of the high-frequency seismic signal generated by debris flows. *Earth Surface Processes and Landforms, 44*(13), 2529-2543. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1002/esp.4677 doi: 10.1002/esp.4677

Favreau, P., Mangeney, A., Lucas, A., Crosta, G., & Bouchut, F. (2010). Numerical modeling of landquakes. *Geophysical Research Letters, 37*(15). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2010GL043512 doi: 10.1029/2010GL043512

Foerster, S. F., Louge, M. Y., Chang, H., & Allia, K. (1994). Measurements of the collision properties of small spheres. *Physics of Fluids, 6*(3), 1108-1115. Retrieved from https://doi.org/10.1063/1.868282 doi: 10.1063/1.868282

Galitzine (Golitsyn), B. (1915). Sur le tremblement de terre du 18 février 1911. *Comptes rendus hebdomadaires des séances de l'Académie des sciences, 160*, 810-813.

Gao, Y., De Simone, G., & Koorapaty, M. (2021). Calibration and verification of DEM parameters for the quantitative simulation of pharmaceutical powder compression process. *Powder Technology, 378*, 160-171. Retrieved from https://www.sciencedirect.com/science/article/pii/S0032591020308834 doi: 10.1016/j.powtec.2020.09.019

GDR MiDi. (2004, August 01). On dense granular flows. *The European Physical Journal E, 14*(4), 341-365. Retrieved from https://doi.org/10.1140/epje/i2003-10153-0 doi: 10.1140/epje/i2003-10153-0
Goldhirsch, I. (2010). Stress, stress asymmetry and couple stress: from discrete particles to continuous fields. *Granular Matter*, **12**, 239–252. Retrieved from https://doi.org/10.1007/s10035-010-0181-z

Gualtieri, L., & Ekström, G. (2017). Seismic reconstruction of the 2012 Palisades rockfall using the analytical solution to Lamb's problem. *Bulletin of the Seismological Society of America*, **107**(1), 63–71.

Harris, C. R., Millman, K. J., van der Walt, S. J., Gommers, R., Virtanen, P., Courapeau, D., ... Oliphant, T. E. (2020, September). Array programming with NumPy. *Nature*, **585**(7825), 357–362. Retrieved from https://doi.org/10.1038/s41586-020-2649-2 doi: 10.1038/s41586-020-2649-2

Hasegawa, H. S., & Kanamori, H. (1987). Source mechanism of the magnitude 7.2 grand banks earthquake of november 1929: Double couple or submarine landslide? *Bulletin of the Seismological Society of America*, **77**(6), 1984. Retrieved from http://dx.doi.org/

Hertz, H. R. (1881). Über die berührung fester elastischer körper. *Journal für die reine und angewandte Mathematik*, **92**, 156–171. Retrieved from https://home.uni-leipzig.de/pwm/web/download/Hertz1881.pdf

Hibert, C., Ekström, G., & Stark, C. P. (2017). The relationship between bulk-mass momentum and short-period seismic radiation in catastrophic landslides. *Journal of Geophysical Research: Earth Surface*, **122**(5), 1201–1215. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016JF004027 doi: 10.1002/2016JF004027

Hibert, C., Mangeney, A., Grandjean, G., & Shapiro, N. M. (2011). Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to rockfall characteristics. *Journal of Geophysical Research: Earth Surface*, **116**(F4). doi: 10.1029/2011JF002038

Hibert, C., Stark, C. P., & Ekström, G. (2015). Dynamics of the oso-steelhead landslide from broadband seismic analysis. *Natural Hazards and Earth System Sciences*, **15**(6), 1265–1273. Retrieved from https://www.nat-hazards-earth-syst-sci.net/15/1265/2015/ doi: 10.5194/nhess-15-1265-2015

Huang, C.-J., Yin, H.-Y., Chen, C.-Y., Yeh, C.-H., & Wang, C.-L. (2007). Ground vibrations produced by rock motions and debris flows. *Journal of Geophysical Research: Earth Surface*, **112**(F2). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2005JF000437 doi: 10.1029/2005JF000437

Hunklinger, S., & Arnold, W. (2012). Ultrasonic properties of glasses at low temperatures. In W. P. Mason & R. N. Thurston (Eds.), *Physical acoustics v12: Principles and methods* (p. 155–216). Elsevier Science. Retrieved from https://books.google.fr/books?id=1H\_ksoDvmiMC

Hunter, J. D. (2007). Matplotlib: A 2d graphics environment. *Computing in Science & Engineering*, **9**(3), 90–95. doi: 10.1109/MCSE.2007.55

Iwashita, K., & Oda, M. (1998, 03). Rolling resistance at contacts in simulation of shear band development by dem. *Journal of Engineering Mechanics*, **124**(3), p285. Retrieved from https://doi.org/10.1061/(ASCE)0733-9399(1998)124:3(285)

Jiang, M., Yu, H.-S., & Harris, D. (2005). A novel discrete model for granular material incorporating rolling resistance. *Computers and Geotechnics*, **32**(5), 340–357. Retrieved from http://www.sciencedirect.com/science/article/pii/S0266352X05000649 doi: 10.1016/j.compgeo.2005.05.001

Jing, L., Kwok, C. Y., Leung, Y. F., & Sobral, Y. D. (2016). Characterization of base roughness for granular chute flows. *Physical Review E*, **94**, 052901. Retrieved from https://link.aps.org/doi/10.1103/PhysRevE.94.052901 doi: 10.1103/PhysRevE.94.052901

Jop, P., Forterre, Y., & Pouliquen, O. (2005). Crucial role of sidewalls in granular surface flows: consequences for the rheology. *Journal of Fluid Mechanics*, **541**,
167–192. doi: 10.1017/S0022112005005987

Jop, P., Forterre, Y., & Pouliquen, O. (2006, June). A constitutive law for dense granular flows. *Nature*, **441**(7094), 727–730. Retrieved from http://www.nature.com/nature/journal/v441/n7094/suppinfo/nature04801\_S1.html (10.1038/nature04801) doi: 10.1038/nature04801

Kanamori, H., & Given, J. W. (1982). Analysis of long-period seismic waves excited by the may 18, 1980, eruption of mount st. helens—a terrestrial monopole? *Journal of Geophysical Research: Solid Earth*, **87**(B7), 5422–5432. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB087iB07p05422 doi: 10.1029/JB087iB07p05422

Kawakatsu, H. (1989). Centroid single force inversion of seismic waves generated by landslides. *Journal of Geophysical Research: Solid Earth*, **94**(B9), 12363–12374. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB094iB09p12363 doi: 10.1029/JB094iB09p12363

Kean, J. W., Coe, J. A., Coviello, V., Smith, J. B., McCoy, S. W., & Arattano, M. (2015). Estimating rates of debris flow entrainment from ground vibrations. *Geophysical Research Letters*, **42**(15), 6365–6372. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015GL064811 doi: 10.1002/2015GL064811

Kruggel-Emden, H., Simsek, E., Rickelt, S., Wirtz, S., & Scherer, V. (2007). Review and extension of normal force models for the discrete element method. *Powder Technology*, **171**(3), 157–173. Retrieved from http://www.sciencedirect.com/science/article/pii/S0032591006004360 doi: 10.1016/j.powtec.2006.10.004

Kruggel-Emden, H., Wirtz, S., & Scherer, V. (2008). A study on tangential force laws applicable to the discrete element method (dem) for materials with viscoelastic or plastic behavior. *Chemical Engineering Science*, **63**(6), 1523–1541. Retrieved from http://www.sciencedirect.com/science/article/pii/S0009250907008767 doi: 10.1016/j.ces.2007.11.025

Kuhn, M. R., Suzuki, K., & Daouadji, A. (2020). Linear-frictional contact model for 3d discrete element simulations of granular systems. *International Journal for Numerical Methods in Engineering*, **121**(3), 560–569. Retrieved from https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.6216 doi: 10.1002/nme.6216

Kuwabara, G., & Kono, K. (1987, August). Restitution coefficient in a collision between two spheres. *Japanese Journal of Applied Physics*, **26**(Part 1, No. 8), 1230–1233. Retrieved from https://doi.org/10.1143/2Fjjap.26.1230 doi: 10.1143/jjap.26.1230

La Rocca, M., Galluzzo, D., Saccorotti, G., Tinti, S., Cimini, G. B., & Del Pezzo, E. (2004, 10). Seismic Signals Associated with Landslides and with a Tsunami at Stromboli Volcano, Italy. *Bulletin of the Seismological Society of America*, **94**(5), 1850–1867. Retrieved from https://doi.org/10.1785/012003238 doi: 10.1785/012003238

Lai, V. H., Tsai, V. C., Lamb, M. P., Ulizio, T. P., & Beer, A. R. (2018, 06). The seismic signature of debris flows: Flow mechanics and early warning at Montecito, California. *Geophysical Research Letters*, **45**(11), p5528–5535. doi: 10.1029/2018GL077683

Levy, C., Mangeney, A., Bonilla, F., Hibert, C., Calder, E. S., & Smith, P. J. (2015). Friction weakening in granular flows deduced from seismic records at the soufrière hills volcano, montserrat. *Journal of Geophysical Research: Solid Earth*, **120**(11), 7536–7557. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2015JB012151 doi: 10.1002/2015JB012151

Lin, C.-H., Kumagai, H., Ando, M., & Shin, T. (2010, 11). Detection of land-slides and submarine slumps using broadband seismic networks. *Geophysical Research Letters - GEOPHYS RES LETT*, **37**. doi: 10.1029/2010GL044685
Lommen, S., Schott, D., & Lodewijks, G. (2014). Dem speedup: Stiffness effects on behavior of bulk material. *Particuology*, 12, 107–112. Retrieved from https://www.sciencedirect.com/science/article/pii/S1674200113001387 (Special issue on conveying and handling of particulate solids – Challenges of discrete element simulation, application and calibration) doi: 10.1016/j.partic.2013.03.006

Lubkin, J. (1951). The torsion of elastic spheres in contact. *Journal of Applied Mechanics*, 18(2), p183-187. Retrieved from https://doi.org/10.1115/1.4010273

Luding, S. (2008). Introduction to discrete element methods. *European Journal of Environmental and Civil Engineering*, 12(7-8), 785–826. Retrieved from https://doi.org/10.1080/19648189.2008.9693050 doi: 10.1080/19648189.2008.9693050

McKinney, W. (2010). Data Structures for Statistical Computing in Python. In S. van der Walt & J. Millman (Eds.), *Proceedings of the 9th Python in Science Conference* (p. 56–61). doi: 10.25080/Majora-92bf1922-00a

Meredith, C. S., & Swab, J. J. (2020, Dec). *Compression strength of borosilicate and soda-lime silicate glasses using a dumbbell-shaped specimen* (Tech. Rep.).

Mindlin, R. D., & Deresiewicz, H. (1953, 09). Elastic Spheres in Contact Under Varying Oblique Forces. *Journal of Applied Mechanics*, 20(3), 327–344. Retrieved from https://doi.org/10.1115/1.4010702 doi: 10.1115/1.4010702

Moretti, L., Allstadt, K. E., Mangeney, A., Capdeville, Y., Stutzmann, E., & Bouchut, F. (2015). Numerical modeling of the mount meager landslide constrained by its force history derived from seismic data. *Journal of Geophysical Research: Solid Earth*, 120(4), 2579–2599. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2014JB011426 doi: 10.1002/2014JB011426

Moretti, L., Mangeney, A., Capdeville, Y., Stutzmann, E., Huggel, C., Schneider, D., & Bouchut, F. (2012). Numerical modeling of the mount steller landslide flow history and of the generated long period seismic waves. *Geophysical Research Letters*, 39(16). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2012GL052511 doi: 10.1029/2012GL052511

Moretti, L., Mangeney, A., Walter, F., Capdeville, Y., Bodin, T., Stutzmann, E., & Le Friant, A. (2020, 02). Constraining landslide characteristics with Bayesian inversion of field and seismic data. *Geophysical Journal International*. Retrieved from https://doi.org/10.1093/gji/ggaa056 (ggaa056) doi: 10.1093/gji/ggaa056

Norris, R. D. (1994). Seismicity of rockfalls and avalanches at three cascade range volcanoes: Implications for seismic detection of hazardous mass movements. *Bulletin of the Seismological Society of America*, 84(6), 1925. Retrieved from http://dx.doi.org/

Piantini, M., Gimbert, F., Korkolis, E., Rousseau, R., Bellot, H., & Recking, A. (2023). Solid concentration as a main proxy for basal force fluctuations generated by highly concentrated sediment flows. *Geophysical Research Letters*, 50(1), e2022GL100345. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2022GL100345 (e2022GL100345 2022GL100345) doi: 10.1029/2022GL100345

Schneider, D., Bartelt, P., Caplan-Auerbach, J., Christen, M., Huggel, C., & McArdell, B. W. (2010). Insights into rock-ice avalanche dynamics by combined analysis of seismic recordings and a numerical avalanche model. *Journal of Geophysical Research: Earth Surface*, 115(F4). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2010JF001734 doi: 10.1029/2010JF001734

Seward, T. P., III, & Vascott, T. (Eds.). (2005). *High temperature glass melt property database for process modeling*. The American Ceramic Society.
Sigmund Lindner. (2018). *Product data sheet*. Retrieved online. Retrieved from https://www.sigmund-lindner.com/en/products/glass-beads/dispersing-beads/

Sun, Z., Espinoza, D. N., & Balhoff, M. T. (2016). Discrete element modeling of indentation tests to investigate mechanisms of co2-related chemomechanical rock alteration. *Journal of Geophysical Research: Solid Earth, 121*(11), 7867-7881. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/2016JB013554 doi: 10.1002/2016JB013554

Taylor, S., & Brodsky, E. E. (2017). Granular temperature measured experimentally in a shear flow by acoustic energy. *Physical Review E, 96*, 032913. Retrieved from https://link.aps.org/doi/10.1103/PhysRevE.96.032913 doi: 10.1103/PhysRevE.96.032913

Thomas, A. L., & Vriend, N. M. (2019, Jul). Photoelastic study of dense granular free-surface flows. *Phys. Rev. E, 100*, 012902. Retrieved from https://link.aps.org/doi/10.1103/PhysRevE.100.012902 doi: 10.1103/PhysRevE.100.012902

Thornton, C., Cummins, S. J., & Cleary, P. W. (2011). An investigation of the comparative behaviour of alternative contact force models during elastic collisions. *Powder Technology, 210*(3), 189-197. Retrieved from http://www.sciencedirect.com/science/article/pii/S0032591011000325 doi: 10.1016/j.powtec.2011.01.013

Thornton, C., Cummins, S. J., & Cleary, P. W. (2013). An investigation of the comparative behaviour of alternative contact force models during inelastic collisions. *Powder Technology, 233*, 30-46. Retrieved from http://www.sciencedirect.com/science/article/pii/S0032591012005670 doi: 10.1016/j.powtec.2012.08.012

Tsai, V. C., Minchew, B., Lamb, M. P., & Ampuero, J.-P. (2012). A physical model for seismic noise generation from sediment transport in rivers. *Geophysical Research Letters, 39*(2). Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2011GL050255 doi: 10.1029/2011GL050255

Weichert, D., Horner, R. B., & Evans, S. G. (1994, 10). Seismic signatures of land-slides: The 1990 Brenda Mine collapse and the 1965 hope rockslides. *Bulletin of the Seismological Society of America, 84*(5), 1523-1532.

Weinhart, T., Orefice, L., Post, M., van Schrojenstein Lantman, M. P., Denissen, I. F., Tunuguntla, D. R., ... Thornton, A. R. (2020). Fast, flexible particle simulations — an introduction to mercurydpm. *Computer Physics Communications, 249*, 107129. Retrieved from https://www.sciencedirect.com/science/article/pii/S0010465519304357 doi: 10.1016/j.cpc.2019.107129

Weinhart, T., Thornton, A. R., Luding, S., & Bokhove, O. (2012). From discrete particles to continuum fields near a boundary. *Granular Matter, 14*, 289-294. Retrieved from https://doi.org/10.1007/s10035-012-0317-4

Woodhouse, M. J., Hogg, A. J., & Sellar, A. A. (2010). Rapid granular flows down inclined planar chutes. part 1. steady flows, multiple solutions and existence domains. *Journal of Fluid Mechanics, 652*, 427-460. doi: 10.1017/S0022112010000303

Wu, F., Fan, Y., Liang, L., & Wang, C. (2016, 08). Numerical simulation of dry granular flow impacting a rigid wall using the discrete element method. *PLOS ONE, 11*(8), 1-17. Retrieved from https://doi.org/10.1371/journal.pone.0160756 doi: 10.1371/journal.pone.0160756

Yamada, M., Kumagai, H., Matsushi, Y., & Matsuzawa, T. (2013). Dynamic landslide processes revealed by broadband seismic records. *Geophysical Research Letters, 40*(12), 2998-3002. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1002/grl.50437 doi: 10.1002/grl.50437

Yamada, M., Mangeney, A., Matsushi, Y., & Matsuzawa, T. (2018). Estimation
of dynamic friction and movement history of large landslides. Landslides, 15, 1963–1974. doi: 10.1007/s10346-018-1002-4

Yamada, M., Mangeney, A., Matsushi, Y., & Moretti, L. (2016, 06). Estimation of dynamic friction of the Akatani landslide from seismic waveform inversion and numerical simulation. Geophysical Journal International, 206(3), 1479–1486. Retrieved from https://doi.org/10.1093/gji/ggw216 doi: 10.1093/gji/ggw216

Zhang, Z., Walter, F., McArdell, B. W., de Haas, T., Wenner, M., Chmiel, M., & He, S. (2021). Analyzing bulk flow characteristics of debris flows using their high frequency seismic signature. Journal of Geophysical Research: Solid Earth, 126(12), e2021JB022755. Retrieved from https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2021JB022755 (e2021JB022755 2021JB022755) doi: 10.1029/2021JB022755

Zhao, J., Ouyang, C. J., Ni, S. D., Chu, R. S., & Mangeney, A. (2020). Analysis of the June 2017 Maoxian Landslide Processes with Force Histories from Seismological Inversion and Terrain Features. Geophysical Journal International. Retrieved from https://doi.org/10.1093/gji/ggaa269 (ggaa269) doi: 10.1093/gji/ggaa269