Journal of Geophysical Research: Solid Earth

Friction weakening in granular flows deduced from seismic
records at the Soufrière Hills Volcano, Montserrat

Clara Levy1,2,3, Anne Mangeney1,4, Fabian Bonilla1,2, Clément Hibert1,3, Eliza S. Calder5,
and Patrick J. Smith6,7

1Institut de Physique du Globe de Paris, CNRS UMR 7154, Université Paris Diderot-Paris 7, Paris, France, 2IFSTTAR, Université
Paris-Est, Paris, France, 3Bureau des Recherches Géologiques et Minières, Orleans, France, 4ANGE team, CEREMA, INRIA,
Laboratoire Jacques-Louis Lions, Paris, France, 5School of GeoSciences, University of Edinburgh, Edinburgh, UK,
6Montserrat Volcano Observatory, Flemings, Montserrat, 7Seismic Research Centre, University of the West Indies,
St. Augustine, Trinidad and Tobago

Abstract Accurate modeling of rockfalls and pyroclastic flows is still an open issue, partly due to a lack
of measurements related to their dynamics. Using seismic data from the Soufrière Hills Volcano, Montserrat,
and granular flow modeling, we show that the power laws relating the seismic energy Es to the seismic
duration ts and relating the loss of potential energy ΔEp to the flow duration tf are very similar, like the
power laws observed at Piton de la Fournaise, Reunion Island. Observations showing that tf ≃ ts suggest
a constant ratio Es∕ΔEp ≃ 10−5. This similarity in these two power laws can be obtained only when the
granular flow model uses a friction coefficient that decreases with the volume transported. Furthermore,
with this volume-dependent friction coefficient, the simulated force applied by the flow to the ground
correlates well with the seismic energy, highlighting the signature of this friction weakening effect in
seismic data.

1. Introduction

Volcanic activity often generates rockfalls and pyroclastic flows. The areas reached by these natural hazards
as well as their dynamics remain difficult to predict because of the lack of precise measurements on the prop-
agation phase. Yet as the mass accelerates and decelerates over the topography, it generates seismic waves
due to the fluctuations of the basal stress applied to the ground. The emitted signal covers a large range of fre-
quencies related to the mean friction and loading/unloading of the flowing mass and to the rebounds of the
involved grains. Seismic monitoring has therefore been used increasingly over the years to collect valuable
information on the flow dynamics at different scales [Hibert et al., 2014a; Moretti et al., 2015; Allstadt, 2013].
These data overcome the difficulty of making direct visual observations in the presence of clouds, darkness,
or an incomplete view of the volcano.

However, deriving the characteristics of the gravitational flows from seismic properties is not straightforward
because the generated seismic signal is a complex combination of the effects of volume, flow dynamics, inter-
actions with the topography, wave propagation properties, etc., [Suriñach et al., 2005; Deparis et al., 2008;
Hibert et al., 2011]. Low-frequency seismic signals have been inverted to find the seismic source mechanism
and have been simulated using a combination of landslide and wave propagation models [Brodsky et al., 2003;
Favreau et al., 2010; Moretti et al., 2012, 2015; Yamada et al., 2013]. High-frequency seismic signals are, however,
very difficult to invert and simulate, as the Earth model at these frequencies and the associated Green’s func-
tions are usually unknown. As a result, these signals are mostly analyzed using empirical laws that attempt
to relate the signal properties to the granular flow characteristics [Suriñach et al., 2005; Deparis et al., 2008;
Helmstetter and Garambois, 2010; Lacroix and Helmstetter, 2011; Dammeier et al., 2011; Hibert et al., 2011,
2014b]. The challenge is to extract information on the rheological properties of the natural granular flows from
these laws. Indeed, one of the main issues in understanding landslide dynamics is the apparent low friction
acting during the downslope motion of large landslides [Heim, 1932; Dade and Huppert, 1998; Pouliquen, 1999;
Legros, 2002; Lucas and Mangeney, 2007; Lucas et al., 2011, 2014]. Field observations show that large landslides
travel over unexpectedly long distances when comparing to smaller landslides, suggesting lower dissipation
during their flow than for smaller landslides.

RESEARCH ARTICLE
10.1002/2015JB012151

Key Points:
• The friction weakening of granular

flows was evidenced using
seismic data

• Conversion of potential to seismic
energy is almost constant for
granular flows

• Seismic energy correlates with the
force applied by granular flows to
the ground

Correspondence to:
C. Levy,
c.levy@brgm.fr

Citation:
Levy, C., A. Mangeney, F. Bonilla,
C. Hibert, E. S. Calder, and P. J. Smith
(2015), Friction weakening in
granular flows deduced from seis-
mic records at the Soufrière Hills
Volcano, Montserrat, J. Geophys.
Res. Solid Earth, 120, 7536–7557,
doi:10.1002/2015JB012151.

Received 22 APR 2015

Accepted 5 OCT 2015

Accepted article online 10 OCT 2015

Published online 6 NOV 2015

©2015. American Geophysical Union.
All Rights Reserved.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7536

http://publications.agu.org/journals/
http://onlinelibrary.wiley.com/journal/10.1002/(ISSN)2169-9356
http://dx.doi.org/10.1002/2015JB012151


Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Lucas et al. [2014] have quantified this friction weakening effect on the basis of observations of runout
distances (i.e., static measurements) of a wide range of small to large landslides. We address here the question
as to whether the signature of this friction weakening can be identified in seismic data, i.e., in data related to
the flow dynamics.

Since the 1990s, many studies have attempted to relate the characteristics of seismic signals (amplitude,
energy, and duration) to rockfall properties (volume, duration, and path [Norris, 1994; Brodsky et al., 2003;
Suriñach et al., 2005; Vilajosana et al., 2008; Deparis et al., 2008; Helmstetter and Garambois, 2010; Lacroix and
Helmstetter, 2011; Dammeier et al., 2011; Hibert et al., 2011]). In particular, using seismic and photogrammetric
data from rockfalls at Piton de la Fournaise volcano (Reunion Island), Hibert et al. [2011] showed that the power
law relating the seismic energy Es to the seismic duration ts

Es ∝ t𝛽s , (1)

where 𝛽 is a constant and is similar to the power law relating the loss of potential energy ΔEp of the granular
mass to the flow duration tf

ΔEp ∝ t𝛽
f
. (2)

In the following sections, we annotate 𝛽p or 𝛽s this power law coefficient in order to distinguish when it is
estimated using loss of potential energy values or seismic energy values.

As field observations have shown that ts ≈ tf , this similarity in the power laws has been used to determine
a value between 10−5 and 10−3 for the mean ratio of seismic to potential energy Rs∕p = Es∕ΔEp [Hibert et al.,
2011, 2014b], in agreement with Deparis et al. [2008] for rockfalls in the French Alps. The similarity in power
laws (1) and (2) makes it possible to calculate the volume of the granular flow from the measured seismic
energy. Numerical simulations have shown that the coefficient 𝛽 of the observed power law is highly sensitive
to the interaction of the granular flow with the underlying topography [Hibert et al., 2011].

An important question is whether or not power laws (1) and (2) and their similarity are also observed in con-
texts other than those of the Piton de la Fournaise volcano. In particular, major simplifications were made
in the seismic energy calculation, potentially affecting the quantification of the parameters involved in the
power law and thus the volume calculation. Indeed, the attenuation was not considered frequency dependent
and only vertical component seismometers were used in Hibert et al. [2011]. Furthermore, because rockfalls
at Piton de la Fournaise only involve relatively small volumes (up to thousands of cubic meters), the potential
friction weakening occurring during large volume landslides has not been identified.

By combining seismic data and granular flow modeling, our aims are (i) to show that these power laws and
energy ratios are universal and can be observed in other geological contexts and (ii) to provide evidence for
and quantify the friction weakening effect of granular flows involving large volumes from data related to the
flow dynamics.

More specifically, these aims are as follows:

1. Using 200 records of granular flows at the Soufrière Hills Volcano (Montserrat Island), we demonstrate that a
power law (1) similar to the one obtained by Hibert et al. [2011] is observed. We also find that the parameters
of this power law depend on the valley in which the granular flows occur.

2. By considering frequency-dependent parameters and all the components of the ground velocity, we
improve the calculation of the seismic energy and show that the coefficients of power law (1) are very stable,
supporting its relevance.

3. Using numerical simulations of granular flows at the Soufrière Hills Volcano, we show that power law (2)
also holds at this site. Interestingly, the power law deduced from seismic data can only be reproduced by
simulations with a friction coefficient that decreases with the volume [Lucas et al., 2014] and not with a
constant friction coefficient. The good agreement between simulations and seismic observations shows
that friction weakening is not negligible in the dynamics of large volumes.

4. Finally, we show that the scatter of the power law coefficients constrained by seismic data is essentially due
to topographic effects, which strikingly influence the mass propagation. We show that pulses in the seismic
envelope are directly related to the interaction of the flow with the topography changes.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7537



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 1. Epicenters of the 46 granular flows that could be located from records from 4 to 6 November 2009 on the
Soufrière Hills Volcano.Main figure: Epicenters of the granular flows that could be located. Estimations of the location
errors are represented by red circles. Events located within the valleys are circled in color: blue for Gages ghaut, green
for Aymer’s ghaut, purple for Gingoes ghaut, orange for White River, yellow for Dry ghaut, and red for Tuitt’s ghaut. The
event presented in Figure 3 is indicated by a star on the border of Aymer’s ghaut. The volcanic centers of the southern
part of Montserrat are represented by brown squares. The MVO seismic stations used in this study are represented as
red triangles. Insert: Location of Montserrat within the Lesser Antilles arc.

2. Characteristics and Location of Rockfalls and Pyroclastic Flows

2.1. Geological Context at the Soufrière Hills Volcano
Montserrat Island is located in the Lesser Antilles arc (insert of Figure 1) where the Atlantic plate is being
subducted beneath the Caribbean plate [Feuillet et al., 2002]. This island comprises four volcanic centers with
ages ranging from 2.5 Ma to the present [Harford et al., 2002; Le Friant et al., 2008]. The latest center in activity
is the Soufrière Hills Volcano, which dominates the southern two thirds of the island at about 1 km above sea
level. Its activity started at least 170 ka ago [Sparks and Young, 2002]. The upper structure of the Soufrière Hills
Volcano is a set of andesitic domes surrounded by a cloak of pyroclastic deposits and collapse debris [Young
et al., 1998].

After a long period of dormancy, Soufrière Hills Volcano became active in 1995 with a series of phreatic erup-
tions [Young et al., 1998]. Since 1995, the activity of this volcano has seen alternate phases of dome growth
and collapses that have produced onshore and offshore pyroclastic and debris flow deposits [e.g., Deplus et al.,
2001; Le Friant et al., 2004; Trofimovs et al., 2006; Le Friant et al., 2009; Loughlin et al., 2010]. For more details on
the setting of the Soufrière Hills Volcano, see Kokelaar [2002].

2.2. Rockfalls and Pyroclastic Flows at the Soufrière Hills Volcano
We analyze the seismic signal of a series of 200 rockfalls and pyroclastic flows occurring during a 3 day crisis
starting 4 November 2009. This period with important rockfalls and pyroclastic flows activity occurs during
the fifth phase of lava extrusion (phase 5) that lasted from 9 October 2009 to 11 February 2010. The volumes
of these granular flows are estimated to range between 103 and 106 m3 while their durations range from a few
tens of seconds to about 350 s (E. Calder, University of Edinburgh, personal communication, 2013). Both types
of events consist of a mixture of hot rock and ash. The difference lies on the size of the events, as the Montserrat
Volcano Observatory (MVO) classifies the smaller ones as rockfalls and the larger ones as pyroclastic flows.
Figure 2 presents photographs of two events occurring at the Soufrière Hills Volcano, one rockfall and one
pyroclastic flow.

Most rockfalls and pyroclastic flows originate from the dome of the Soufrière Hills Volcano (grey areas in
Figure 1) [Jolly et al., 2002]. Materials spread downward, roughly following preferential paths such as the sur-
rounding valleys. Rockfalls rarely travel far beyond the base of the talus apron, 500–800 m from the dome
summit, whereas the bigger events can propagate several kilometers and reach the sea shore.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7538



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

a) b)

Figure 2. Photos of a rockfall and a pyroclastic flow at the Soufrière Hills Volcano.(a) Rockfall dust trail on the upper part
of White River. (b) Pyroclastic flow going down Gages ghaut.

2.3. Characteristics of the Seismic Events
The seismic signals were recorded at nine digital stations belonging to the permanent seismic network of
the MVO (Figure 1). The stations were equipped with Guralp CMG-40T 3C broadband sensors, having a band-
width of 0.03–30 Hz. The sampling rate was 100 Hz. The events studied here have been identified by the MVO
with the usual procedures applied at the Soufrière Hills Volcano, as described by Miller et al. [1998]. The main
categories currently in use at MVO are VT (volcano-tectonic), hybrid, long period, and rockfalls (including both
pyroclastic flows and rockfalls). There is no automatic classification at MVO and triggered events are manu-
ally classified by a seismic analyst, mainly on the basis of waveform and with additional spectral analysis if
required.

Seismic signals associated with rockfalls are complex given that the source of the seismic waves propagates
downslope, while covering an area of increasing size. Authors like Hibert et al. [2011] and Deparis et al. [2008]
report observations that the duration of the seismic signal is about the same duration as the propagation of
materials on the slope (i.e., the duration of the seismic source). An example of a typical recorded signal is given
in Figure 3a. As in this example, these signals have emerging onsets and ragged wave train envelopes. Most of
them are characterized by dominant frequencies f between 2 and 10 Hz [Langer et al., 2006], although energy
is found over the wider range 1 < f < 20 Hz (Figure 3b).

It is often assumed that rockfall signals are mainly made of surface waves [Rousseau, 1999; Deparis et al., 2008;
Zhao et al., 2012]. However, as the source is very complex, it is difficult to identify wave packets within the
data. Figures 3c–3e present the signal filtered at increasing frequency bands, together with the correspond-
ing polarization plots (Figures 3f and 3g). The envelope of the signal at low frequency (Figure 3c) is smooth
and cigar shaped, and the wave packets are well mixed (see the polarization plots of Figure 3f ). When the
recording station is close to the seismic source, the signal has impulsive parts at high frequency (see the sig-
nal in purple on the zoom of Figure 3d), and the corresponding wave arrivals can be identified (see the zoom
on the polarization plots of Figure 3g). The elliptical motion identifies them as surface waves.

2.4. Automatic Picking Using the Kurtosis of the Signal
Seismic signals associated with rockfalls have emerging onsets. This particular feature makes the picking of
arrival times with classical seismological methods impractical and unreliable.

A specific automated method based on kurtosis has been used here [Baillard et al., 2013; Hibert et al., 2014b].
This statistical method measures the fourth moment of the amplitude distribution of the signal. The principle
is the following: kurtosis indicates the resemblance to a normal distribution. The amplitude distribution is
normal for seismic noise and events but turns leptokurtic at the transition between noise and event (i.e., the
amplitude distribution is sharper, or more peaked, than a normal distribution and the value of the kurtosis is
greater than 3). Thus, by computing the value of the kurtosis for small sliding time windows, we can find the
onset of rockfall signals, when the value of the kurtosis increases. This procedure is performed on the filtered
seismic signal for sliding windows of varying sizes. The final picking takes into account the results for the
signal filtered on multiple frequency bands. For our database, optimal parameters were sliding windows of 5 s,
10 s, and 15 s, as well as four band-pass filters with overlapping frequencies: 1–7 Hz, 5–10 Hz, 8–15 Hz, and
12–20 Hz. An example of automatic picking is illustrated in Figure 4 for stations MBFR, MBLY, and MBBY. The
estimated kurtosis function is shown at each frequency band for station MBFR in Figure 4a, and the integrated

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7539



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 3. Characteristics of the seismic signal of a granular flow. (a) Seismic signal of the granular flow represented by a
star in Figure 1 as recorded at station MBFR 1.7 km away (vertical component). This event is also used as an example in
Figure 9e and in Figures 4 and 5 and in Figure A2 in Appendix A. (b) The associated S transform. (c–e) Seismic records of
the same event for frequency bands 1–8 Hz, 8–10 Hz, and 10–20 Hz, respectively. The Hilbert envelope of the eastern
component is plotted in red. (f–h) The polarization plots for the frequency bands corresponding to Figures 3c–3e. The
identification of wave packets is difficult in this type of signal, although some transient phenomena can be spotted at
high frequencies. An example of an impulsive part of the signal is plotted in purple in the zoom of Figure 3d. The
corresponding wave arrivals can be identified on the zoom of the polarization plot of Figure 3g.

kurtosis function for all frequency bands is shown for stations MBFR, MBLY, and MBBY in Figures 4b–4d. In
this example, the signal-to-noise ratio (SNR) is reasonably high for all stations (the SNR being considered here
as the ratio between the mean amplitude of the first 10 s of the signal and the mean amplitude of the 10 s
preceding the event). When the SNR value drops below 2.5, the picking procedure becomes unreliable. The
end of each rockfall signal was considered to be the time at which the smoothed Hilbert envelope of the signal
recovers its preevent value. The Hilbert envelope is smoothed using a 5 s averaging filter.

2.5. Locating Source Areas Using the Fast Marching Method
The observed time delays were used to locate the source areas by applying the fast marching method
(FMM) [Sethian, 1996; Hibert et al., 2014b]. The FMM is a numerical technique for computing the position of

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7540



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 4. Automatic picking of an event located in Aymer’s ghaut.(a) Example of automatic picking (red vertical line)
of an event located in Aymer’s ghaut (see Figures 1 and 3) on the seismic record of station MBFR (vertical component).
The picking procedure is based on the average of the estimated kurtosis function at four different frequency bands
(see legend). (b–d) Examples of automatic picking (red vertical line) on the vertical components of stations MBFR, MBLY,
and MBBY for the same event as in Figure 4a. The averaged kurtosis function for all frequency bands is displayed as
green areas.

propagating fronts that uses an extremely fast scheme for solving the Eikonal equation. Without going into
detail, the FMM allows us to estimate grids of travel times from each MVO seismic station at a given wave veloc-
ity (see the example in Figure 5a). To solve the problem, we consider a medium with homogeneous velocity,
where rockfall signals are mainly related to surface waves. As surface waves propagate along the topography,
the computed travel times depend strongly on the topography. Computations use a digital elevation model
based on aerial lidar data acquired by the MVO in 2008. Grid spacing is 10 m.

It is now possible to compute delay maps for each station pair, corresponding to the difference in travel times
between the two stations of each pair (see the example in Figure 5b). The observed time delay for a pair of
stations can be reported on the corresponding delay map, thus appearing as a hyperbola (see the example
in Figure 5c). The width of the hyperbolas takes into account the potential errors in the estimation of time
delays. For each event, the positions of the hyperbolas corresponding to all the observed time delays are
plotted, and the result is stacked on one single map (see the example in Figure 5d). When an event can be
located, three or more hyperbolas intersect, thus delineating a source area. The possibility of delineating a
source area is considered for a realistic range of surface wave velocities. Locating is thought to be optimal for
the wave velocity at which the largest number of hyperbolas intersect with a minimum difference between
observed and computed time delays (to be more specific, it is the root-mean-square of these differences that
is minimized; see Figure 5e).

This method using hyperbolas was preferred over a random search of epicenter locations within the grid by
minimizing the RMS of the differences between observed and computed time delays. The advantages of this
method are as follows: (1) the maps are computed relatively quickly, as only simple operations are performed
on already existing grids that can be reused for different events and (2) this method excludes aberrant time
arrivals from the locating process, as the corresponding hyperbolas will not focus. The uncertainties on the
epicenter locations are directly proportional to the estimated RMS.

Locating was possible only when the signal-to-noise ratio (SNR) was sufficiently high for multiple stations
(SNR > 2.5), i.e., only for 46 events out of 200. The epicenters of the events that could be located are reported
on the map in Figure 1 with an error represented by the size of the red circles.

In most cases, this error is about the same size as the dot representing the epicenter, i.e., around 80 m.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7541



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 5. Locating of an event using the fast marching method. (a) Travel time map from station MBBY for a
homogenous surface wave velocity of 1000 m/s. Travel times were estimated by applying the FMM on a digital elevation
model (DEM) grid with 10 m spacing. (b) Delay map between stations MBBY and MBLY for a homogenous surface wave
velocity of 1000 m/s, deduced from the travel time maps of both stations. (c) The observed delay ΔTobs between arrival
times at stations MBBY and MBLY for an event in Aymer’s ghaut (presented in Figures 1 and 3) is reported on the delay
map. Values matching ΔTobs ± dt are set to 1, where dt = 0.3 s is the worst picking error possible for SNR = 2.5. This
process makes it possible to plot a hyperbola (in red) with a width that takes into account the potential errors on the
estimation of time delays. (d) The observed delays between stations MBLY, MBBY, and MBFR are reported on a single
map, thus leading to three distinct hyperbolas. The hyperbolas intersect close to the Soufrière Hills Volcano summit,
constraining the probable location of the source area (in red). (e) Number of crossing hyperbolas and RMS value as a
function of the surface wave velocity used to compute the FMM delay maps.

Most epicenters are located on the western flank of the Soufrière Hills Volcano, with a large majority of
events occurring within Gages and Aymer’s ghauts (blue and green epicenters in Figure 1). The event
distribution is consistent with the MVO report of activity for the same period, indicating a switch in the
direction of the dome growth from south-west to west in early November 2009 (2009 MVO activity report,
http://www.mvo.ms/pub/Open_File_Reports/MVO_OFR_13_04-MVO_Activity_Reports_2009.pdf).

3. Power Law Between Seismic Energy and Signal Duration
3.1. Seismic Power Law
In order to potentially identify power law (1), we estimate the seismic energy and the duration of the events.
Because the objective is to extract quantitative information from the power law in terms of landslide rheology,
it is important to reduce the uncertainties in the energy calculation as much as possible, as well as to quantify
and understand the scattering of the data.

Calculating seismic energy and signal duration is an important issue in seismology [Trifunac and Brady, 1975;
Pousse et al., 2006; Kanamori, 1977; Singh and Ordaz, 1994; Izutani and Kanamori, 2001]. We estimated the
duration of each event using the automatic method explained in the section 2.4. As the picking may differ
from one station to another, the duration ts is chosen as the average duration over all the stations. Picking

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7542

http://www.mvo.ms/pub/Open_File_Reports/MVO_OFR_13_04-MVO_Activity_Reports_2009.pdf


Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

differences are typically within the range of a few seconds (< 5 s), which is small compared to the total rockfall
durations (30 to 350 s).

To estimate the seismic energy, we assume a point-source [Kanamori and Given, 1982; Eissler and Kanamori,
1987; Dahlen, 1993; Favreau et al., 2010; Moretti et al., 2012; Zhao et al., 2012] and we consider the medium
as isotropic and homogeneous. We also assume that surface waves dominate the seismic signal and use the
same relation as [Hibert et al., 2011, 2014b]:

Es = 2𝜋r𝜌hc ∫
t2

t1

uenv(t)2e𝛼rdt (3)

with

uenv(t) =
√

u(t)2 + Ht(u(t))2, (4)

where t1 and t2 (ts = t2 − t1) are the picked onset and final times of the seismic signal, respectively, r is the
distance between the event and the recording station, h is the thickness of the layer through which surface
waves propagate, 𝜌 is the density of the ground, c is the group velocity of the seismic waves (and not the
phase velocity as in Hibert et al. [2011]), 𝛼 is a damping factor that accounts for anelastic attenuation of the
waves [Aki and Richards, 1980], and uenv(t) is the amplitude envelope of the ground velocity obtained using
the Hilbert transform (Ht).

The energy is computed using the three components of the signal (u =
√

u2
E + u2

N + u2
Z ). It is designated Es freq,

as the parameters h and 𝛼 are considered to be frequency dependent. For each event, the seismic energy
is averaged for all the stations to obtain a mean value of Es freq. The density is approximated at 2000 kg/m3.
The source-station distance r is deduced from the location of the events, and we assume an initiation at the
Soufrière Hills Volcano summit for the events that could not be located. We assume that the signal is mostly
made of Rayleigh waves, even though this is very difficult to observe on the recorded signal (see Figure 3).
Because the velocity model is insufficiently constrained, we assume that the wave velocity is constant and
slightly smaller than the shear wave velocity Vs: c ≈ 0.9Vs = 1300 m/s. In that case, the thickness of the layer
through which surface waves propagate can be approximated by h(f ) ≃ 2.5c∕f [Lay and Wallace, 1995]. By
using the seismic records of the located rockfalls at different stations, we empirically estimate 𝛼 as a function
of the frequency (see Appendix A):

𝛼 = Af 𝛾 (5)

with A = 2.4 × 10−4 ± 2.4 × 10−5 m−1 and 𝛾 = 0.4 ± 0.04. For frequency bands between 0.5 and 19.5 Hz, 𝛼
ranges from 2.2 × 10−4 to 8.7 × 10−4 m−1.

The resulting seismic energy is plotted as a function of the seismic signal duration for each of the 200 rockfalls
in Figure 6a. It ranges from 105 to 108 J, covering 4 orders of magnitude. The duration varies from 29 to 347 s.
Despite strong scattering of the data, we roughly identify the scaling law Es freq ∝ t𝛽s (grey line in Figure 6a). The
corresponding regression coefficient is 𝛽 = 2.5 with R2 = 0.54, indicating a weak relationship between Es freq

and ts. However, the p value of the model is very small compared to 0.01 (pvalue = 10−29), supporting the
hypothesis that a relationship between Es freq and ts is likely to exist in the form of the power law (1). Indeed, p
values are here designed for testing the hypothesis of no correlation between log(Es freq) and log(ts) against
the alternative that there is a nonzero correlation. As the p value is small, we can say that this correlation
is significantly different from zero [Gibbons, 1985]. The scattering of the seismic energy values around the
regression line covers 2 orders of magnitude. An empirical variance of 𝛽 = 2.5 ± 0.35 and R2 = 0.54 ± 0.05
was estimated by performing one thousand bootstrap resamplings [Bates and Watts, 1988].

3.2. Differences in Seismic Energy Calculations and Power Law Dispersions
We present here the estimation of the seismic energy Es when the parameters h, c, and 𝛼 are considered
constant (as in Hibert et al. [2011, 2014b]), as well as the estimation of the seismic energy Esz when only the
vertical component uZ(t) is used. The results are compared with the estimation of the seismic energy Es freq

presented in section 3.1.

To estimate Es, we take c = 1300 m/s, consistent with Jolly et al. [2002],𝛼 = f𝜋∕Qc, with f = 5 Hz corresponding
to the central frequency of the events, and Q = 20, as roughly estimated by Jolly et al. [2002] for the Soufrière
Hills Volcano and corresponding to a typical value for the attenuation of surface waves in andesitic volcanoes.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7543



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 6. The power law (1) for Es and Es freq. (a) Seismic energy Es freq estimated using frequency-dependent parameters (h and 𝛼) as a function of the signal
duration ts for the 200 recorded granular flows averaged over all the stations. Best fit regression lines (solid lines) and power law coefficients 𝛽 and associated R2

values are indicated in the legend. (b) Seismic energy Es as a function of the signal duration ts for the 200 recorded granular flows averaged over all the stations.
(c) Coefficient 𝛽 , R2, and p value estimated using the energy Es, Es freq, Esz , or Esz freq. The ratios of the mean (over all the events and all the stations) energy Esi to
the mean energy Es freq are plotted as percentages (Esi = [Es, Es freq, Esz , Esz freq]).

This gives 𝛼 = 6.0 × 10−4 m−1. The depth is assumed to be h = 260 m, corresponding to one wavelength of
Rayleigh waves at f = 5 Hz and c = 1300 m/s.

The resulting seismic energy Es is plotted as a function of the seismic signal duration ts in Figure 6b. The dis-
persion is very similar to that found for Es freq (Figure 6a), although slightly higher as R2 = 0.53 compared to
R2 = 0.54 for Es freq. For both computations of the seismic energy, the empirical variance of R2 (determined
using bootstrap resampling) is 0.05. Thus, the frequency dependence of the parameters does not explain the
main part of the dispersion in the data. Moreover, taking into account frequency-dependent parameters do
not change the coefficient of the power law as shown in Figure 6c. Indeed, here 𝛽 = 2.5 ± 0.33 instead of
2.5±0.35. However, the energy is slightly overestimated (110%) when the parameters of the energy calculation
are considered to be constant (Figure 6c, fourth panel).

To calculate Esz , we assume that all the components have about the same amplitude, so that u =
√

3uZ . The
objective is also to quantify the error in the energy calculation made when using only uz (as done, for example,
in Hibert et al. [2011]). Indeed, it is common to have only seismic stations with one component (vertical) and for
three component stations, the horizontal components may be very noisy [Bonnefoy-Claudet et al., 2006]. We
observe that the energy Esz calculated with the vertical component is about half of the total energy (Figure 6c).
The reason is that the horizontal components of the signal generally have a higher amplitude than the vertical
component (see Figure 3). The value of 𝛽 is, however, not very different (Figure 6c, first panel). On the other
hand, the dispersion is smaller when using only the Z component, possibly because the signal is generally less
noisy on the vertical component. As for the total energy, the energy Esz and the associated power law changes
only very slightly when frequency-dependent parameters are taken into account (Figure 6c).

3.3. Origin of the Scattering in the Power Law
To assess the relevance of the power law, the origin of the scattering of the data must be understood. This
scattering could be due to errors in picking the signal duration and to the approximations made in the energy
computation (presence of body waves, 3-D variations of the seismic velocities, and density). Another crucial
point is the topography over which the granular mass is flowing, which has been shown to change the coeffi-
cient of the power law [Hibert et al., 2011]. We focus on this aspect, as we may expect different power laws for
rockfalls flowing within different valleys. To check this, we differentiate between the rockfalls located in differ-
ent valleys, keeping only the valleys where a minimum of six events were located. Even though we have a very
small number of points (16 events in Aymer’s ghaut, 12 in Gages ghaut, and 6 in the White River valley), we
observed different coefficients for each valley, i.e., 𝛽Aymer′s = 2.9±1.6, 𝛽Gages = 3.2±1.7, and 𝛽White = 1.9±2.8.
Hibert et al. [2011] showed that the value of 𝛽 decreases with increasing rugosity of the topography. This is
compatible with our observations because the profiles of Aymer’s ghaut and Gages ghaut are very similar

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7544



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 7. Two-dimensional profiles used in the simulations of granular flows. (a) Two-dimensional profiles of the main
valleys of the Soufrière Hills Volcano (extracted from a 1 m gridded digital elevation model (DEM) and slightly smoothed
using a 20 points averaging filter), together with an exponential fit of the mean profile of these three valleys (black
dashed line), as well as this exponential fit with added rugosity (grey dashed line). The rugosity corresponds to the
mean topographical features of the profiles of Aymer’s ghaut, Gages ghaut, and White River filtered above 100 m.
Departure positions of simulated granular flows are indicated by arrows. (b) Variation of the slope angle with the
horizontal distance along the profiles.

while the fluctuation of the slope in the White River valley is greater (Figure 7). Even with these different val-
ues of 𝛽 , the scattering of the points for one valley is not smaller than for all the data together, as indicated by
the R2 values (from 0.49 to 0.53, with an empirical variance determined using bootstrap resampling ranging
from 0.2 to 0.26).

These results show that we observe the same power laws as those observed in the very different context of
rockfalls at Piton de la Fournaise with coefficients 𝛽 of the power law strongly related to the local topography
over which the masses flow [Hibert et al., 2011]. For Montserrat and Reunion Island, these coefficients range
between 1.8 < 𝛽 < 3.2.

4. Friction Weakening in Granular Flows and Topography Explain
Seismic Observations

As in Hibert et al. [2011], we perform a series of numerical simulations of granular flows over the topography
of the different valleys to test for evidence of a power law similar to (2) between the loss of potential energy
ΔEp of the mass and the duration of the flow tf . The general idea is to constrain the friction coefficient in the
granular flow model to reproduce the power law (2) with the same coefficient 𝛽 as that observed using seismic
data (power law (1)). Numerical simulations are also expected to give insight into the topography effects on
flow dynamics and into the origin of the scattering observed in the seismic power law.

4.1. Numerical Simulations
We used the SHALTOP numerical model, which describes dry granular flows over a complex 3-D topogra-
phy [Bouchut et al., 2003; Bouchut and Westdickenberg, 2004; Mangeney-Castelnau et al., 2005; Mangeney et al.,
2007]. This model is based on a depth-averaged thin-layer approximation (the flow is assumed to be thin
compared to its longitudinal extent). It describes the change of flow thickness with time h(x, y, t) in the direc-
tion normal to the topography and the depth-averaged velocity of the flow u(x, y, t) along the topography
z = b(x, y) where (x, y, z) are the coordinates in a Cartesian reference frame. This model deals with the
full tensor of terrain curvature. SHALTOP has successfully reproduced experimental granular flows
[Mangeney-Castelnau et al., 2005; Mangeney et al., 2007] and natural landslides [Lucas and Mangeney, 2007;
Kuo et al., 2009; Mangold et al., 2010; Favreau et al., 2010; Lucas et al., 2011].

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7545



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

SHALTOP uses a Coulomb-type friction law with a friction coefficient 𝜇 = tan 𝛿, where 𝛿 is the friction angle.
This friction coefficient can be considered as a measure of the effective dissipation during the flow [Mangeney
et al., 2010; Lucas et al., 2014]. We carried out simulations using either

1. a constant friction coefficient for all rockfalls or
2. a friction coefficient that varies with the landslide volume V approximately as [Lucas et al., 2014]:

𝜇 = 1
V0.0774

. (6)

This variation of the friction coefficient with the landslide volume (equation (6)) is an empirical law determined
by Lucas et al. [2014] by fitting runout distances for small to very big landslides occurring within different
contexts, including extraterrestrial landslides.

We simulated 2-D flows created by the release of a parabolic mass from rest over the along-slope profiles of
Aymer’s ghaut, Gages ghaut, and White River extracted from a 1 m gridded digital elevation model, slightly
smoothed using a 20 point averaging filter (Figure 7). These valleys correspond to the most common flow
paths during the three days of records. Simulations were also done over an exponential profile fitted on the
mean profiles of these three valleys with and without added rugosity. The added rugosity corresponds to
the mean topographical features of the three profiles filtered above 100 m. For each profile, we varied the
position of the released mass (top, middle, and bottom of the released area) and its volume (from 500 to 106

m3), according to geological observations. The potential released areas are represented in grey in Figure 1.
For each simulation, we calculated the loss of potential energy ΔEp of the mass between the initial instant
t = 0 s and the final state tf when the deposit has formed. The flow duration tf is calculated by assuming
that the mass is at rest when the maximum velocity is lower than 0.1 m/s for a flow thickness higher than
10 cm. Note that simulations using SHALTOP do not take into account erosion processes that may increase
the flow duration [Mangeney et al., 2010]. Furthermore, only single collapses are simulated here even though
successive collapses may be mixed up in one event [Moretti et al., 2015].

4.2. Conversion of Potential to Seismic Energy
Figure 8a shows ΔEp(tf ) for three series of simulations on the Aymer’s ghaut profile using three different
friction angles 𝛿 = 20∘, 𝛿 = 30∘, and 𝛿 = 35∘ (orange, blue, and red dots, respectively). These values are in the
range of the friction angles measured during laboratory experiments using samples from the Soufrière Hills
Volcano [Voight et al., 2002]. Each point cloud is fitted using a robust fit algorithm (dashed lines of slope 𝛽).
This algoritm uses iteratively reweighted least squares with a bisquare weighting function (2015 MATLAB doc-
umentation web page on robust fit, http://fr.mathworks.com/help/stats/robustfit.html). Interestingly, none
of these three series of simulations reproduces the seismic power law Es(ts) (green dots and solid line in
Figure 8a). When 𝛿 ≥ 30∘, the simulated large volumes (i.e., 106 m3) tend to stop at distances much shorter
than those observed in the field. As a result, the associated flow duration tf ≤ 80 s is too short given that large
volumes usually have signal durations ts > 200 s. For such friction angles, the simulated flow durations of small
volumes (500–3000 m3) are, however, reasonably consistent with seismic observations (20 ≤ tf ≤ 30 s). On
the other hand, when 𝛿 = 20∘, simulations of large volumes (i.e., 106 m3) lead to flow durations more con-
sistent with seismic observations while small volumes run too far and have flow durations longer than those
observed (tf ≥ 60 s). Thus, whatever the chosen friction coefficient, the simulations using a constant fric-
tion coefficient for all rockfalls cover a range of durations much narrower than the durations of the generated
seismic signal. The approximations of the SHALTOP model might be the cause of the difference between sim-
ulations and seismic observations, as nonhydrostatic forces are not completely included in the shallow water
formulation. In principle, nonhydrostatic forces will slow down the flow, as in hydrostatic, the free surface gra-
dient (which is a driving force) is always overestimated. However, we believe it unlikely, as in the literature,
all models, whether continuous or discrete, need to use very low friction coefficient to reproduce the runout
distances of large landslides [Smart et al., 2010; Lucas et al., 2014].

Instead of a constant friction coefficient for all rockfalls, we will now use the variable friction coefficient given
by relation (6). The color of each symbol in Figure 8b corresponds to the value of the friction angle 𝛿 and its
shape (circle, square, or triangle) to the departure zone (bottom, middle, and top, respectively). The use of
a variable friction coefficient gives a range of simulated durations (21 ≤ tf ≤ 219 s) similar to the seismic
durations observed in Aymer’s ghaut (29 ≤ ts ≤ 347 s). The best fit of ΔEp ∝ t𝛽

f
gives 𝛽p = 2.7 ± 1.8, very

close to 𝛽s = 2.9 ± 1.6 obtained from seismic data for the events located in Aymer’s ghaut (green dots and

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7546

http://fr.mathworks.com/help/stats/robustfit.html


Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 8. The loss of potential energy ΔEp for simulations at the Soufrière Hills Volcano. (a) ΔEp as a function of the flow
duration tf for simulations in Aymer’s ghaut compared to the seismic energy Es freq as a function of signal durations ts for
the 200 recorded granular flows. Three series of simulations were performed with friction angles 𝛿 = 20∘ , 30∘, and 35∘,
respectively. (b) Same as Figure 8a but for simulations using a friction angle 𝛿 dependent on the volume (equation (6)).
Best fit regression lines and power law coefficients 𝛽 are indicated on the plots, together with the regression statistic R2.
(c and d) Same plot as in Figure 8b but for simulations along the profiles of Gages ghaut and White River. (e) Ditto for
simulations along the exponential fit profile (triangles) and the same profile with added rugosity (squares). (f ) Ditto for a
mix of simulations along Aymer’s ghaut, Gages ghaut, and White River.

line in Figure 8b). When calculating different values of 𝛽p for simulations corresponding to different departure

zones in Aymer’s ghaut, we find 𝛽p = 3.43, 𝛽p = 2.45, and 𝛽p = 2.41 when the mass is released from the top,

middle, and bottom, respectively (Figure 8b). These results show that different departure zones may induce

scattering in the observed power law and variability in the value of 𝛽 .

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7547



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Table 1. Mean Ratio Rs∕p = Es freq∕ΔEp for the Cases of Figures 8b–8d and 8f

Rs∕p Aymer′s Rs∕p Gages Rs∕p W.River Rs∕p All

2.78 × 10−5 1.11 × 10−5 1.25 × 10−5 1.52 × 10−5

A similar series of simulations was conducted on Gages ghaut and White River profiles as shown in Figures 8c
and 8d, respectively. Despite the small number of located events (12 events in Gages ghaut and 6 events
in the White River valley), the slope 𝛽p obtained from the simulations ΔEp(tf ) when using a variable friction
coefficient is very close to the value obtained from seismic data by fitting Es(ts). For Gages ghaut, 𝛽s = 3.2±1.7
and 𝛽p = 3.4 ± 1.5 and for White River, 𝛽s = 1.9 ± 2.8 and 𝛽p = 1.3 ± 2.5.

The value of 𝛽 is thus closely related to the topography over which the granular mass flows. Indeed, for sim-
ulations over a smooth exponential profile fitted on the mean profiles of the valleys (see Figure 7), we obtain
𝛽p = 6.8, much larger than 𝛽s = 2.5 (Figure 8e). Moreover, the range of durations of the flow in the simula-
tions (18 ≤ tf ≤ 81 s) is not consistent with seismic observations (29 ≤ ts ≤ 347 s). On the other hand, when a
mean roughness is added to this exponential profile (Figure 7), the simulations are more consistent with the
seismic observations with 𝛽p = 2 and 24 < tf < 289 s.

Plotting the seismic data and numerical simulations for all the valleys together, Figure 8f shows again that
𝛽p = 2.4±0.8 and 𝛽s = 2.8±1.3 are close. Interestingly, the scattering of the seismic and simulated power laws
is also quite similar, with R2 = 0.58 ± 0.1 for seismic data and R2 = 0.59 ± 0.06 for simulations. Furthermore,
the fluctuation of the potential energy around its mean value is of 2 orders of magnitude, like the scattering of
the seismic energy. The scattering of the seismic data Es(ts) is likely due to the different valleys and departure
zones and to the scattering of the power law ΔEp(tf ) for a given valley and departure zone.

Numerical simulations of granular flow with a friction coefficient dependent on the landslide volume (see
equation (6)) allowed to reproduce better the seismic observations. In the literature, several mechanisms
have been proposed to explain the friction reduction under certain conditions for landslides, such as lubri-
cation by water or air, thermal pressurization, fragmentation, acoustic fluidization, flash heating, or erosion
[Legros, 2002; Singer et al., 2012; Davies, 1982; Tsutsumi and Shimamoto, 1997; Lucas and Mangeney, 2007; Lucas
et al., 2014].

Our models do not take into account the full complexity of natural flows. In particular, we do no model erosion
that has been shown to increase runout distances for granular flows on slopes higher than about half of their
friction angle, which is the case at Montserrat [Mangeney et al., 2010; Farin et al., 2014]. Since erosion efficiency
increases with the volume involved, part of the apparent friction weakening may be explained by erosion
effect. With the current knowledge, it is not possible to quantify this effect.

Lucas et al. [2014] showed that another physical process that could well explain this volume-dependent fric-
tion law for granular flows is an actual frictional weakening with increasing slip velocity. A micromechanical
process that can lead to dramatic velocity weakening is flash heating [Lube et al., 2004; Moretti et al., 2012]. The
real contact between two rough solid surfaces generally occurs over a small fraction of their nominal contact
area, on highly stressed microcontacts (asperities). Slip produces frictional heating at the microcontact scale.
If slip is fast enough to prevent heat dissipation by conduction, the microcontact experiences a significant
transient temperature rise that activates thermal effects such as melting, dehydration, or other phase trans-
formations. This reduces the local shear strength of the microcontact and leads to a macroscopic decrease of
the friction coefficient as a function of sliding velocity.

4.3. The Ratio Rs∕p = Es∕𝚫Ep, Comparison With Other Studies
The strong similarity of the relations Es(ts) and ΔEp(tf ) suggests a constant transfer of potential to seismic
energy for these granular flows, as observed at Piton de la Fournaise. The ratio Rs∕p = Es∕ΔEp ≃ 10−5 is almost
constant for all the valleys (Table 1), showing that only a small portion of the potential energy released is
converted into seismic energy. Our ratio Rs∕p = Es∕ΔEp at Montserrat is of the same order of magnitude as
the ratios found by Deparis et al. [2008] (essentially ranging from 10−5 to 10−4) for 10 rockfalls recorded by a
regional seismological network in the Alps with volumes 2 × 103 < V < 1.75 × 106 m3 and source-station
distances 10 < d < 100 km, by Hibert et al. [2011] (10−5 < Rs∕p < 10−3) for granular flows at the Piton de la
Fournaise volcano with 1 < V < 5 × 103 m3 and 0.3 < d < 2.1 km and by Berrocal et al. [1978] (Rs∕p ≃ 10−6)
for the 1974 Mantaro slide in Peru with d > 80 km.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7548



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

When comparing in detail with Piton de la Fournaise, for a given seismic (or flow) duration, the seismic energy
is 1 order of magnitude higher in Montserrat and the loss of potential energy is more than 2 orders of magni-
tude higher, leading to a smaller ratio Rs∕p. The higher seismic energy obtained in Montserrat could be partly
explained (i.e., a factor of 2) by the fact that all the components of the ground velocity are taken into account.

4.4. Flow Dynamics and Seismic Energy
Numerical simulations provide the details of the interaction between the flow and the topography. They also
make it possible to calculate the force applied by the flow to the ground that generate seismic waves [Favreau
et al., 2010; Moretti et al., 2012, 2015]. Figures 9a–9d present the details of a typical simulation along the
Aymer’s ghaut profile when using a volume-dependent friction coefficient (6). An initial mass of V =50, 701 m3

(blue area in Figure 9a) with a corresponding friction angle 𝛿 = 23.4∘ is released about 4300 m from the
sea shore. The mass accelerates down the slope until it reaches a position where the slope angle 𝜃 is smaller
than the friction angle 𝛿, at time t = 10 s (Figures 9a–9c). A second acceleration phase starting at t = 15 s is
observed when the mass again reaches a region where 𝜃 > 𝛿 until the maximum velocity is reached at t = 20 s.
The mass then slows down as it reaches another portion of the slope with 𝜃 < 𝛿, until it stops. At 36 s, most
of the mass has already stopped, while a small fraction of the material is still flowing down until it forms a
second bulge 1.2 km from the initial position (grey area in Figure 9a). The force F that the mass applies to the
ground reflects the acceleration/deceleration phases as illustrated in Figure 9c where the norm of this force F
calculated with the SHALTOP model is represented (for details on the force calculation, see Favreau et al. [2010]
and Moretti et al. [2012, 2015]). Essentially, the force has a direction opposite to the flow when it accelerates
and the same direction than the flow when it decelerates [see, e.g., Zhao et al., 2012]. The direction of the force
could be observed on its horizontal component Fh represented on Figure 10b. Figures 10a and 10b show Fh

and the angle of the slope 𝛼center at the position of the center of mass for the simulation shown in Figure 9a.
Fh decreases with increasing 𝛼center. Each time 𝛼center becomes lower than the friction angle (red dashed
lines in Figure 10a), Fh becomes positive, with a time delay of a few seconds. Thus, the mass decelerates.
Conversely, each time 𝛼center becomes higher than the friction angle, Fh becomes negative, with a time delay
of a few seconds. Thus, the mass accelerates.

Figure 9d shows the instantaneous frictional power W(t) = F ⋅u that was found to be strongly correlated with
the seismic power envelope filtered between 0.1 and 0.05 Hz [Schneider et al., 2010]. One of the objectives is
to observe potential correlations with the higher frequency seismic signal.

In the frequency range investigated here (f > 1 Hz), it is not possible to invert the seismic signal to determine
the force applied by the flow to the ground or to calculate the seismic signal generated by the force F as was
done at lower frequencies by Brodsky et al. [2003], Favreau et al. [2010], Moretti et al. [2012], Ekström and Stark
[2013], Yamada et al. [2013], Allstadt [2013], and Moretti et al. [2015]. At these frequencies, seismic waves are
sensitive to small inhomogeneities within the Earth, which seismic velocity models cannot capture. Another
way of comparing the numerical flow model with seismic data in this frequency range is to compare the
simulations with the observed seismic power (i.e., seismic energy per second). To do this, we choose an event
with a seismic duration close to the flow duration calculated with the model that was located in Aymer’s Ghaut.
An example is represented in Figures 9e and 9f, where the seismic duration of the event ts = 87 s is close to
the flow duration tf = 97 s. The calculated force F correlates well with the seismic power estimated using 1 s
sliding windows. Both exhibit two peaks around t = 3–5 s and t =36–40 s. The frictional power is less peaked
and less clearly related to the seismic power estimated using sliding windows. Most of the seismic energy
ranges within 0–6 Hz. Note that almost no energy is recorded between 6 and 20 Hz during the first peak,
while 25% of the energy in the second peak is between 6 and 20 Hz.

Qualitatively, the same observations can be made for the 5 times smaller event of Figures 9g–9l. The maximum
force (first peak) is 5 times smaller than that calculated in Figure 9c, suggesting that this maximum force is
proportional to the volume involved as verified by the third event (Figures 9m–9r) and in agreement with
Ekström and Stark [2013] and Zhao et al. [2012]. The second peak of the force is much smaller than for the
larger event, as is also observed for the seismic power and the frictional power. The force peaks are, however,
more correlated over time with those of the seismic power than the frictional power peaks.

The same observations can be made once again for the even smaller simulated rockfall (V =3906 m3

and 𝛿 = 27.8∘) of Figures 9m and 9n. For this simulated rockfall, most of the mass stops upslope in the region

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7549



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 9. Simulations of granular flow and seismic energy of events located at Aymer’s ghaut. (a–d) Outputs of a 2-D
simulation of granular flow along the Aymer’s ghaut profile, with an initial mass of 50,701 m3 and friction angle
𝛿 = 23.4∘ . (a) Thickness of the simulated mass above the ground at the initial state (blue mass), at 3 s, 15 s, and 36 s after
the start of the simulation and at the end of the simulation (grey mass). (b) Change of the slope angle with the horizontal
distance along the profile. (c) Maximum velocity of the simulated mass (blue line) and norm of the basal force applied
by the flow to the ground (red line) as a function of time. Simulation end time tf is indicated. (d) The frictional power at
the base of the flowing mass (basal force multiplied by the velocity). (e) Seismic energy Es estimated using 1 s sliding
windows of the event of Figure 3 located at Aymer’s ghaut. End time ts is indicated. (f ) Seismic energy Es estimated
using 1 s sliding windows of the event in Figure 3 for high frequencies only (6–20 Hz). (g–j) outputs of a 2-D simulation
of granular flow along the Aymer’s ghaut profile, with an initial mass of 11,267 m3 and friction angle 𝛿 = 25.9∘ . (k, l)
Same as Figures 9e and 9f for an event located at Aymer’s ghaut. (m–p) Outputs of a 2-D simulation of granular flow
along the Aymer’s ghaut profile, with an initial mass of 4 × 103 m3 and friction angle 𝛿 = 25.9∘ . (q, r) Same as Figures 9e
and 9f for another event located at Aymer’s ghaut. Note the exaggeration of the vertical axis in Figures 9a, 9g, and 9m.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7550



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 10. Simulations of granular flows at Aymer’s ghaut. (a) Angle of the slope 𝛼center at the position of the center of
mass for the simulation of Figure 9a. The horizontal red line indicates the value of the friction angle 𝛿. The vertical red
line indicates the time tf , i.e., the sliding duration. (b) The horizontal component Fh and the norm of the basal force F
applied by the flow to the ground as a function of time for the simulation of Figure 9a. (c and d) Same as Figures 10a
and 10b for the simulation of Figure 9q.

where 𝜃 < 𝛿 for the first time (see Figures 9m and 9n). A small portion of the mass is able to propagate farther
down toward a new zone where 𝜃 > 𝛿 and then stops when 𝜃 < 𝛿 again. The second velocity peak observed in
Figure 9o is likely related to the acceleration of this small propagating mass. A small force and frictional power
are associated with this motion for tf > 40 s. Similarly, low seismic power is also observed for ts > 40 s support-
ing a very strong correlation between the seismic power and the accelerations/decelerations of the flowing
mass. Figures 10c and 10d show the angle of the slope 𝛼center at the position of the center of mass and the
horizontal component of the basal force Fh, respectively, for the simulation shown in Figure 9k. Fh decreases
with increasing 𝛼center. Each time 𝛼center becomes lower (higher) than the friction angle (red dashed lines in
Figure 10c), Fh becomes positive (negative), with a time delay of a few seconds. Thus, the mass decelerates.

Note that the two force peaks for the small event of Figures 9m and 9n are not as well pronounced as for the
larger events (Figures 9a and 9g) and that the frictional power exhibits only one peak. Interestingly, the seismic
power between 0 and 6 Hz exhibits one peak while two peaks are observed at higher frequencies (Figure 9r).
The peaks in the seismic power are, however, not well in phase with those of the force and here the frictional
power correlates slightly better with the seismic power.

For most of the events located in Aymer’s ghaut, the simulated basal force and the frictional power corre-
late very well with the seismic power at f > 1 Hz for flows with tf ≃ ts. Only four seismic events were poorly
correlated with the simulations, two of them being mixed with other seismic sources.

Very similar seismic signals and flow dynamics are observed for the events located in Gages ghaut (Figure 11),
also exhibiting a strong correlation between force F, frictional power, and seismic power. Two simulations
along the profile of Gages ghaut (flow thickness, velocity, force applied to the ground, and frictional power)
and two seismic events (seismic power), located at Gages ghaut, of similar durations (tf ≃ ts) are compared.
The simulations are performed with a friction coefficient that decreases with the volume (equation (6)). As for
most simulations in Aymer’s ghaut, the force Fh reflects the acceleration and deceleration of the flow when
it reaches slope angles larger and smaller than the friction angle, respectively (see Figures 11b and 11h).
The forces exhibit two peaks that correlate with the peaks observed for the seismic power estimated slid-
ing windows for seismic events of the same duration as the flow duration (see Figures 11e and 11f, and 11k
and 11l).

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7551



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure 11. Simulations of granular flow and seismic energy of events located at Gages ghaut.(a–d) Outputs of a 2-D
simulation of granular flow along the Gages ghaut profile, with an initial mass of 1,064,392 m3 and friction angle
𝛿 = 18.9∘ . (a) Thickness of the simulated mass above the ground, at t = 0 s (blue mass), 3 s, 27 s, and 63 s and at the end
of the simulation (grey deposit). (b) Change of the slope angle with the horizontal distance along the profile. (c)
Maximum velocity of the simulated mass (blue line) and norm of the basal force applied by the flowing mass to the
ground (red line) as a function of time. Simulation end time tf is indicated. (d) The frictional flow rate (the basal force
multiplied by the velocity). (e) Seismic energy Es freq estimated 1 s sliding windows for an event located in Gages ghaut.
End time ts is indicated. (f ) Seismic energy Es freq for the same event for high frequencies only (6–20 Hz). (g–j) Outputs
of a 2-D simulation of granular flow along the Gages ghaut profile, with an initial mass of 267,092 m3 and friction angle
𝛿 = 20.8∘ . (k, l) Same as Figures 11e and 11f for another event located at Gages ghaut. Note the exaggeration of the
vertical axis of Figures 11a and 11g.

5. Conclusion

Using 200 records of granular and pyroclastic flows at the Soufrière Hills Volcano (Montserrat Island), we
have shown that the seismic energy Es follows a power law Es ∝ t𝛽s

s , where ts is the seismic event duration.
Using numerical simulations of these flows, we have shown that a similar power law is observed for the loss
of potential energy ΔEp ∝ t

𝛽p

f
, where tf is the numerical flow duration. The coefficients are very similar, i.e.,

𝛽s ≃ 𝛽p ≃ 2.4 ± 1 and are shown to mostly depend on the topography where the granular mass flows. These
two power laws have very similar scattering. As ts ≃ tf , these strong similarities between the power laws sug-
gest that the ratio of the seismic energy to the loss of potential energy Rs∕p is almost constant, as observed
for rockfalls at Piton de la Fournaise. Here Rs∕p ≃ 10−5, corresponding to the smallest ratios observed by
Hibert et al. [2011] and Deparis et al. [2008], possibly due to the larger source-station distances in our case.
These power laws, observed both at Montserrat and Reunion Island, make it possible to calculate the

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7552



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

conversion of potential to seismic energy. Given that they seem to hold in these very different contexts, they
may possibly apply more generally.

Despite the low efficiency of the potential to seismic energy conversion, numerical simulation of the flows
shows that the seismic energy correlates well with the force applied by the granular flow to the ground. In
particular, we show that the peaks in seismic power can be explained by the interaction of the flow with the
topography that generate phases of acceleration/deceleration of the granular mass related to the relative
values of the slope angle and the friction angle of the material involved. While the link between the seismic
signal and the flow dynamics has been widely demonstrated at low frequencies (f < 0.1 Hz), our results
obtained with high-frequency data (f > 1 Hz) suggest that information on the change of the force applied by
the flow to the ground with time can be obtained by analysis of the high-frequency seismic power. As a result,
numerical models of natural granular flows can also be constrained by comparing the calculated force with
the high-frequency seismic power over time.

Note that the similarity in the power laws and the correlation between the flow dynamics and the seismic
energy could only be found when using a friction coefficient in the granular flow model that decreased with
the volume of the released mass. Indeed, using a constant friction coefficient for all rockfalls and pyroclastic
flows makes it impossible to determine the duration of small to large events or their runout distances. The
difference between the simulated results using a constant friction coefficient and the seismic data may result
from physical processes that are not taken into account in the model such as erosion processes and to the
approximations made in the thin-layer depth-averaged model used here. However, the very good agreement
between the simulation and seismic data using a volume-dependent friction coefficient derived from com-
pletely different data (data on the deposits of a wide range of small to large landslides in different geological
context) supports the hypothesis that these seismic data may carry the signature of this friction weakening
effect. This is the first time that analysis of seismic data suggests that a friction coefficient depending on the
volume can be used to simulate granular flows with volumes spanning between 500 and 106 m3, such as
those observed in Montserrat. More precisely, the empirical volume-dependent friction coefficient proposed
by Lucas et al. [2014] is found to reproduce well the seismic power laws and the signature of the flow dynamics
on the seismic energy power. As a result, the friction coefficient varies here within the range 0.34 < 𝜇 < 0.62
(i.e., 18∘< 𝛿 < 32∘) for volumes of 500 < V < 106 m3.

The small values of the friction coefficient needed to reproduce large pyroclastic flows in Montserrat and more
generally large landslides worldwide is in agreement with previous studies [Mangeney et al., 2000; Heinrich
et al., 2001; Bayarri et al., 2009; Nikolkina et al., 2011; Pirulli et al., 2007; Lucas et al., 2014]. As shown by Lucas
et al. [2014], the friction weakening with volume can be the result of a friction decrease with velocity or other
dynamic parameters as also suggested for earthquakes. Defining friction weakening in landslides more pre-
cisely is one of the most important challenges in the domain, and we show here that high-frequency seismic
data may provide an ideal tool to go farther in this direction.

Appendix A: An Attenuation Model
We used the seismic records of the rockfalls that could be located to estimate the attenuation factor 𝛼 as a
function of the frequency. Indeed, the amplitude A of seismic records decreases with increasing distance r to
the source according to the following equation: A ∝ r−ne−𝛼r , with n = 0.5 for surface waves and n = 1 for
body waves (consistently with the geometrical spreading of these waves). The linear relationship derived by
Kanai et al. [1984] follows from this property:

loge

[
a1

a2

(
Δ1

Δ2

n)]
= loge

(
𝛾1

𝛾2

)
− 𝛼

(
Δ1 − Δ2

)
(A1)

with a1 and a2 the maximum amplitudes at two stations, Δ1 and Δ2 the corresponding epicentral distances,
and 𝛾1 and 𝛾2 constants relating to the ground and installation conditions of each sensor.

The slope 𝛼 of the linear relationship is estimated at various frequencies using the filtered seismic signals of
all the located rockfalls, for all station pairs. The signals are filtered several times between 0.5 and 19.5 Hz
using 1 Hz band-pass filters. As the seismic source is complex, the maximum amplitudes a1 and a2 are selected
within the first 20 s of the filtered signals. Figure A1 presents the corresponding point clouds for n = 0.5 and
n = 1 at 4.5–5.5 Hz. The linear fit of the point clouds gives an estimation of 𝛼 at 5 Hz when considering that
the signal is dominated by body waves (n = 1) or by surface waves (n = 0.5). For all frequency bands, the

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7553



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure A1. Estimation of the attenuation of high frequencies at the Soufrière Hills Volcano.(a) Estimation of the
attenuation at 5 Hz using the linear relationship of Kanai et al. [1984]. (b) Estimation of the attenuation assuming that
the seismic signal of rockfalls is mostly composed of surface waves (n = 0.5).

linear regressions have higher corresponding R2 statistics when considering that the signal is dominated by
surface waves (n = 0.5). This result is in agreement with our assumption that the seismic signal of rockfalls is
mostly made of surface waves. An empirical curve of the evolution of 𝛼 with frequency is deduced for n = 0.5
(Figure A1).

If we assume that the quality factor Q is frequency dependent and that Q(f ) = Qof n, with n ∈ [0 1] [Erickson
et al., 2004; Campbell, 2009], then 𝛼 would vary with frequency given the following expression: 𝛼 = 𝜋

Qoc
f 1−n.

Our results indicate a possible value of Qo that is very small (Qo < 50 and Qo ≈ 10 if we consider c = 1300 m/s),
consistent with the estimation made by Jolly et al. [2002] for the Soufrière Hills Volcano.

The attenuation 𝛼 is used to compute the seismic energy Es at each frequency. Figure A2 shows the estimation
of Es freq (solid lines) and of the integrated seismic energy ∫ fx

f1
Es(f )df (dashed lines) at each frequency for an

event located in Aymer’s ghaut (see Figures 1 and 3 of the article) and for signals recorded at stations MBFR,
MBLY, MBWH, MBGB, and MBGH. Even if a correction is applied to take the attenuation 𝛼 into account, the
energy at high frequency is somewhat less for station MBFR compared to station MBLY. As a result, the seismic
energy calculated from the seismic signal at MBFR is about 20% smaller than that calculated from the signal

Figure A2. Frequency-dependent Es . (a) The seismic energy Es estimated at different frequencies for stations MBFR,
MBLY, and MBWH (at distances of 2.0, 2.2, and 3.9 km from the source area, respectively) for the event presented in
Figures 1 and 3 (solid lines) and the integrated seismic energy ∫ fx

f1
Es(f )df (dashed lines). (b) Ditto for stations MBGB and

MBGH at distances of 6.0 and 3.7 km from the source area.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7554



Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Figure A3. The seismic energy Es freq as a function of the signal duration ts for different stations. Es freq is estimated for
stations MBFR, MBGH, MBLY, MBGB, and MBWH at distances of 2.0, 3.7, 2.2, 6.0, and 3.9 km from the source area,
respectively. The power law Es ∝ t𝛽s is fitted for each station (solid lines) and the associated 𝛽 values are indicated on
the plot.

at MBLY, even though the source-station distances are almost the same. This could be due to site effects or to
the different path of the waves. In the paper, we averaged Es for all stations in order to average discrepancies
due to the changes of frequency content related to different seismic paths or site effects. Figure A3 represents
the dispersion of the estimation of the seismic energy as a function of the signal duration ts for stations MBFR,
MBGH, MBLY, MBGB, and MBWH. This shows that the average seismic energy at more distant stations is lower,
even when taking attenuation into account. At depths inferior to 500 m, the frequency response of the soil
is highly nonlinear. Our homogenous model for attenuation is too simple to reproduce such phenomena, as
well as path effects due to lateral variations of the ground. Both effects will partly explain the dispersion of
the estimation of the seismic energy for the various seismic stations. The impact of these effects increases
with distance.

References
Aki, K., and P. G. Richards (1980), Quantitative Seismology, vol. 1424, Freeman, San Francisco, Calif.
Allstadt, K. (2013), Extracting source characteristics and dynamics of the August 2010 mount meager landslide from broadband

seismograms, J. Geophys. Res. Earth Surf., 118(3), 1472–1490, doi:10.1002/jgrf.20110.
Baillard, C., W. C. Crawford, V. Ballu, C. Hibert, and A. Mangeney (2013), An automatic kurtosis-based p-and s-phase picker designed for local

seismic networks, Bull. Seismol. Soc. Am., 104, 394–409.
Bates, D. M., and D. G. Watts (1988), Nonlinear Regression Analysis and Its Applications, John Wiley, New York.
Bayarri, M., J. O. Berger, E. S. Calder, K. Dalbey, S. Lunagomez, A. K. Patra, E. B. Pitman, E. T. Spiller, and R. L. Wolpert (2009), Using statistical

and computer models to quantify volcanic hazards, Technometrics, 51(4), 402–413.
Berrocal, J., A. Espinosa, and J. Galdos (1978), Seismological and geological aspects of the Mantaro landslide in Peru, Nature, 275(5680),

533–536.
Bonnefoy-Claudet, S., F. Cotton, and P.-Y. Bard (2006), The nature of noise wavefield and its applications for site effects studies: A literature

review, Earth Sci. Rev., 79(3), 205–227.
Bouchut, F., and M. Westdickenberg (2004), Gravity driven shallow water models for arbitrary topography, Commun. Math. Sci., 2(3),

359–389.
Bouchut, F., A. Mangeney-Castelnau, B. Perthame, and J.-P. Vilotte (2003), A new model of Saint Venant and savage-hutter type for gravity

driven shallow water flows, C. R. Mathematique, 336(6), 531–536.
Brodsky, E. E., E. Gordeev, and H. Kanamori (2003), Landslide basal friction as measured by seismic waves, Geophys. Res. Lett., 30(24), 2236,

doi:10.1029/2003GL018485.
Campbell, K. W. (2009), Estimates of shear-wave q and 𝜅0 for unconsolidated and semiconsolidated sediments in eastern North America,

Bull. Seismol. Soc. Am., 99(4), 2365–2392.
Dade, W. B., and H. E. Huppert (1998), Long-runout rockfalls, Geology, 26(9), 803–806.
Dahlen, F. (1993), Single-force representation of shallow landslide sources, Bull. Seismol. Soc. Am., 83(1), 130–143.
Dammeier, F., J. R. Moore, F. Haslinger, and S. Loew (2011), Characterization of Alpine rockslides using statistical analysis of seismic signals,

J. Geophys. Res., 116, F04024, doi:10.1029/2011JF002037.
Davies, T. R. (1982), Spreading of rock avalanche debris by mechanical fluidization, Rock Mech., 15(1), 9–24.
Deparis, J., D. Jongmans, F. Cotton, L. Baillet, F. Thouvenot, and D. Hantz (2008), Analysis of rock-fall and rock-fall avalanche seismograms in

the French Alps, Bull. Seismol. Soc. Am., 98(4), 1781–1796.
Deplus, C., A. Le Friant, G. Boudon, J.-C. Komorowski, B. Villemant, C. Harford, J. Ségoufin, and J.-L. Cheminée (2001), Submarine evidence

for large-scale debris avalanches in the Lesser Antilles arc, Earth Planet. Sci. Lett., 192(2), 145–157.
Eissler, H. K., and H. Kanamori (1987), A single-force model for the 1975 Kalapana, Hawaii, earthquake, J. Geophys. Res., 92(B6), 4827–4836.
Ekström, G., and C. P. Stark (2013), Simple scaling of catastrophic landslide dynamics, Science, 339(6126), 1416–1419.
Erickson, D., D. E. McNamara, and H. M. Benz (2004), Frequency-dependent lg q within the continental United States, Bull. Seismol. Soc. Am.,

94(5), 1630–1643.

Acknowledgments
The research was funded by Institut
de Physique du Globe de Paris,
the French ANR LANDQUAKES and
SYNAPS@ programmes, and by the
ERC Consolidator SLIDEQUAKES grant.
Classified seismic data were provided
by the Montserrat Volcano Observa-
tory. The numerical model SHALOP
is being developed by the Institut de
Physique du Globe.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7555

http://dx.doi.org/10.1002/jgrf.20110
http://dx.doi.org/10.1029/2003GL018485
http://dx.doi.org/10.1029/2011JF002037


Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Farin, M., A. Mangeney, and O. Roche (2014), Fundamental changes of granular flow dynamics, deposition, and erosion processes at high
slope angles: Insights from laboratory experiments, J. Geophys. Res. Earth Surf., 119(3), 504–532, doi:10.1002/2013JF002750.

Favreau, P., A. Mangeney, A. Lucas, G. Crosta, and F. Bouchut (2010), Numerical modeling of landquakes, Geophys. Res. Lett., 37, L15305,
doi:10.1029/2010GL043512.

Feuillet, N., I. Manighetti, P. Tapponnier, and E. Jacques (2002), Arc parallel extension and localization of volcanic complexes in Quadeloupe,
Lesser Antilles, J. Geophys. Res., 107(B12), 2331, doi:10.1029/2001JB000308.

Gibbons, J. (1985), Nonparametric Statistical Inference, 2nd ed., M. Dekker, New York.
Harford, C., M. Pringle, R. Sparks, and S. Young (2002), The volcanic evolution of Montserrat using 40Ar/39Ar geochronology, in The Eruption

of Soufrière Hills Volcano, Montserrat, from 1995 to 1999, vol. 21, edited by T. H. Druitt and B. P. Kokelaar, pp. 93–113, London.
Heim, A. (1932), Bergsturz und Menschenleben, Fretz and Wasmuth, Indiana Univ.
Heinrich, P., G. Boudon, J. Komorowski, R. Sparks, R. Herd, and B. Voight (2001), Numerical simulation of the December 1997 debris

avalanche in Montserrat, Lesser Antilles, Geophys. Res. Lett., 28(13), 2529–2532.
Helmstetter, A., and S. Garambois (2010), Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their

correlation with rainfalls, J. Geophys. Res., 115, F03016, doi:10.1029/2009JF001532.
Hibert, C., A. Mangeney, G. Grandjean, and N. Shapiro (2011), Slope instabilities in Dolomieu crater, Réunion Island: From seismic signals to

rockfall characteristics, J. Geophys. Res., 116, F04032, doi:10.1029/2011JF002038.
Hibert, C., G. Ekström, and C. P. Stark (2014a), Dynamics of the Bingham canyon mine landslides from seismic signal analysis, Geophys. Res.

Lett., 41, 4535–4541, doi:10.1002/2014GL060592.
Hibert, C., et al. (2014b), Automated identification, location, and volume estimation of rockfalls at Piton de la Fournaise volcano, J. Geophys.

Res. Earth Surf., 119, 1082–1105, doi:10.1002/2013JF002970.
Izutani, Y., and H. Kanamori (2001), Scale-dependence of seismic energy-to-moment ratio for strike-slip earthquakes in Japan, Geophys. Res.

Lett., 28(20), 4007–4010.
Jolly, A., G. Thompson, and G. Norton (2002), Locating pyroclastic flows on Soufriere Hills Volcano, Montserrat, West Indies, using amplitude

signals from high dynamic range instruments, J. Volcanol. Geotherm. Res., 118(3), 299–317.
Kanai, K., K. Yamabe, and A. Habasaki (1984), Study of the attenuation of seismic waves, in Proceedings of the Eighth World Conference on

Earthquake Engineering, vol. 2, pp. 273–280, Prentice-Hall, San Francisco, Calif.
Kanamori, H. (1977), The energy release in great earthquakes, J. Geophys. Res., 82(20), 2981–2987.
Kanamori, H., and J. W. Given (1982), Analysis of long-period seismic waves excited by the May 18, 1980, eruption of Mount St. Helens—A

terrestrial monopole?, J. Geophys. Res., 87(B7), 5422–5432.
Kokelaar, B. (2002), Setting, chronology and consequences of the eruption of Soufrière Hills Volcano, Montserrat (1995–1999),

in The Eruption of Soufrière Hills Volcano, Montserrat From 1995 to 1999, vol. 21, pp. 1–43, Geol. Soc. Mem., London.
Kuo, C., Y. Tai, F. Bouchut, A. Mangeney, M. Pelanti, R. Chen, and K. Chang (2009), Simulation of Tsaoling landslide, Taiwan, based on Saint

Venant equations over general topography, Eng. Geol., 104(3), 181–189.
Lacroix, P., and A. Helmstetter (2011), Location of seismic signals associated with microearthquakes and rockfalls on the Séchilienne

Landslide, French Alps, Bull. Seismol. Soc. Am., 101(1), 341–353.
Langer, H., S. Falsaperla, T. Powell, and G. Thompson (2006), Automatic classification and a-posteriori analysis of seismic event identification

at Soufrière Hills Volcano, Montserrat, J. Volcanol. Geotherm. Res., 153(1), 1–10.
Lay, T., and T. C. Wallace (1995), Modern Global Seismology, vol. 58, Academic Press, San Diego, Calif.
Le Friant, A., C. Harford, C. Deplus, G. Boudon, R. Sparks, R. Herd, and J. Komorowski (2004), Geomorphological evolution of Montserrat

(West Indies): Importance of flank collapse and erosional processes, J. Geol. Soc., 161(1), 147–160.
Le Friant, A., E. Lock, M. Hart, G. Boudon, R. Sparks, M. Leng, C. W. Smart, J.-C. Komorowski, C. Deplus, and J. Fisher (2008), Late pleistocene

tephrochronology of marine sediments adjacent to Montserrat, Lesser Antilles volcanic arc, J. Geol. Soc., 165(1), 279–289.
Le Friant, A., C. Deplus, G. Boudon, R. Sparks, J. Trofimovs, and P. Talling (2009), Submarine deposition of volcaniclastic material from the

1995–2005 eruptions of Soufrière Hills Volcano, Montserrat, J. Geol. Soc., 166(1), 171–182.
Legros, F. (2002), The mobility of long-runout landslides, Eng. Geol., 63(3), 301–331.
Loughlin, S., R. Luckett, G. Ryan, T. Christopher, V. Hards, S. De Angelis, L. Jones, and M. Strutt (2010), An overview of lava dome

evolution, dome collapse and cyclicity at Soufrière Hills Volcano, Montserrat, 2005–2007, Geophys. Res. Lett., 37, L00E16,
doi:10.1029/2010GL042547.

Lube, G., H. E. Huppert, R. S. J. Sparks, and M. A. Hallworth (2004), Axisymmetric collapses of granular columns, J. Fluid Mech., 508, 175–199.
Lucas, A., and A. Mangeney (2007), Mobility and topographic effects for large Valles Marineris landslides on Mars, Geophys. Res. Lett., 34,

L10201, doi:10.1029/2007GL029835.
Lucas, A., A. Mangeney, D. Mège, and F. Bouchut (2011), Influence of the scar geometry on landslide dynamics and deposits: Application to

Martian landslides, J. Geophys. Res., 116, E10001, doi:10.1029/2011JE003803.
Lucas, A., A. Mangeney, and J. P. Ampuero (2014), Frictional velocity-weakening in landslides on Earth and on other planetary bodies,

Nat. Commun., 5, 3417, doi:10.1038/ncomms4417.
Mangeney, A., P. Heinrich, R. Roche, G. Boudon, and J. Cheminée (2000), Modeling of debris avalanche and generated water waves:

Application to real and potential events in Montserrat, Phys. Chem. Earth Part A, 25(9), 741–745.
Mangeney, A., L. Tsimring, D. Volfson, I. Aranson, and F. Bouchut (2007), Avalanche mobility induced by the presence of an erodible bed and

associated entrainment, Geophys. Res. Lett., 34, L22401, doi:10.1029/2007GL031348.
Mangeney, A., O. Roche, O. Hungr, N. Mangold, G. Faccanoni, and A. Lucas (2010), Erosion and mobility in granular collapse over sloping

beds, J. Geophys. Res., 115, F03040, doi:10.1029/2009JF001462.
Mangeney-Castelnau, A., F. Bouchut, J. Vilotte, E. Lajeunesse, A. Aubertin, and M. Pirulli (2005), On the use of Saint Venant equations to

simulate the spreading of a granular mass, J. Geophys. Res., 110, B09103, doi:10.1029/2004JB003161.
Mangold, N., A. Mangeney, V. Migeon, V. Ansan, A. Lucas, D. Baratoux, and F. Bouchut (2010), Sinuous gullies on Mars: Frequency,

distribution, and implications for flow properties, J. Geophys. Res., 115, E11001, doi:10.1029/2009JE003540.
Miller, A., R. Stewart, R. White, R. Luckett, B. Baptie, W. Aspinall, J. Latchman, L. Lynch, and B. Voight (1998), Seismicity associated with dome

growth and collapse at the Soufriere Hills Volcano, Montserrat, Geophys. Res. Lett., 25(18), 3401–3404.
Moretti, L., A. Mangeney, Y. Capdeville, E. Stutzmann, C. Huggel, D. Schneider, and F. Bouchut (2012), Numerical modeling of the Mount

Steller landslide flow history and of the generated long period seismic waves, Geophys. Res. Lett., 39, L16402, doi:10.1029/2012GL052511.
Moretti, L., K. Allstadt, A. Mangeney, Y. Capdeville, E. Stutzmann, and F. Bouchut (2015), Numerical modeling of the mount meager landslide

constrained by its force history derived from seismic data, J. Geophys. Res. Solid Earth, 120, 2579–2599, doi:10.1002/2014JB011426.
Nikolkina, I., N. Zahibo, T. Talipova, and E. Pelinovsky (2011), Pyroclastic flow from Soufrière Hills Volcano, Montserrat: Solid block model,

Int. J. Geosci., 2(3), 326–335.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7556

http://dx.doi.org/10.1002/2013JF002750
http://dx.doi.org/10.1029/2010GL043512
http://dx.doi.org/10.1029/2001JB000308
http://dx.doi.org/10.1029/2009JF001532
http://dx.doi.org/10.1029/2011JF002038
http://dx.doi.org/10.1002/2014GL060592
http://dx.doi.org/10.1002/2013JF002970
http://dx.doi.org/10.1029/2010GL042547
http://dx.doi.org/10.1029/2007GL029835
http://dx.doi.org/10.1029/2011JE003803
http://dx.doi.org/10.1038/ncomms4417
http://dx.doi.org/10.1029/2007GL031348
http://dx.doi.org/10.1029/2009JF001462
http://dx.doi.org/10.1029/2004JB003161
http://dx.doi.org/10.1029/2009JE003540
http://dx.doi.org/10.1029/2012GL052511
http://dx.doi.org/10.1002/2014JB011426


Journal of Geophysical Research: Solid Earth 10.1002/2015JB012151

Norris, R. D. (1994), Seismicity of rockfalls and avalanches at three cascade range volcanoes: Implications for seismic detection of hazardous
mass movements, Bull. Seismol. Soc. Am., 84(6), 1925–1939.

Pirulli, M., M.-O. Bristeau, A. Mangeney, and C. Scavia (2007), The effect of the Earth pressure coefficients on the runout of granular material,
Environ. Model. Software, 22(10), 1437–1454.

Pouliquen, O. (1999), Scaling laws in granular flows down rough inclined planes, Phys. Fluids, 11(3), 542–548.
Pousse, G., L. F. Bonilla, F. Cotton, and L. Margerin (2006), Nonstationary stochastic simulation of strong ground motion time histories

including natural variability: Application to the k-net Japanese database, Bull. Seismol. Soc. Am., 96(6), 2103–2117.
Rousseau, N. (1999), Study of seismic signal associated with rockfalls at 2 sites on the Réunion island (Indian Ocean): Mahavel cascade and

Soufrière cavity, PhD thesis, Université de Paris, Paris.
Schneider, D., P. Bartelt, J. Caplan-Auerbach, M. Christen, C. Huggel, and B. W. McArdell (2010), Insights into rock-ice avalanche dynamics by

combined analysis of seismic recordings and a numerical avalanche model, J. Geophys. Res., 115, F04026, doi:10.1029/2010JF001734.
Sethian, J. A. (1996), A fast marching level set method for monotonically advancing fronts, Proc. Natl. Acad. Sci., 93(4), 1591–1595.
Singer, K. N., W. B. McKinnon, P. M. Schenk, and J. M. Moore (2012), Massive ice avalanches on Iapetus mobilized by friction reduction during

flash heating, Nat. Geosci., 5(8), 574–578.
Singh, S., and M. Ordaz (1994), Seismic energy release in Mexican subduction zone earthquakes, Bull. Seismol. Soc. Am., 84(5), 1533–1550.
Smart, K., D. Hooper, and D. Sims (2010), Discrete Element Modeling of Landslides in Valles Marineris, Mars, AGU, San Francisco, Calif.,

13–17 Dec.
Sparks, R., and S. Young (2002), The eruption of Soufrière Hills Volcano, Montserrat (1995–1999): Overview of scientific results, in The

Eruption of Soufrière Hills Volcano, Montserrat, From 1995–1999, vol. 21, edited by T. H. Druitt and B. P. Kokelaar, pp. 45–69, Geol. Soc.
Mem., London.

Suriñach, E., et al. (2005), Seismic detection and characterization of landslides and other mass movements, Nat. Hazards Earth Syst. Sci., 5(6),
791–798.

Trifunac, M. D., and A. G. Brady (1975), A study on the duration of strong earthquake ground motion, Bull. Seismol. Soc. Am., 65(3), 581–626.
Trofimovs, J., et al. (2006), Submarine pyroclastic deposits formed at the Soufrière Hills Volcano, Montserrat (1995–2003): What happens

when pyroclastic flows enter the ocean?, Geology, 34(7), 549–552.
Tsutsumi, A., and T. Shimamoto (1997), High-velocity frictional properties of gabbro, Geophys. Res. Lett., 24(6), 699–702.
Vilajosana, I., E. Surinach, A. Abellan, G. Khazaradze, D. Garcia, and J. Llosa (2008), Rockfall induced seismic signals: Case study in Montserrat,

Catalonia, Nat. Hazards Earth Syst. Sci., 8(4), 805–812.
Voight, B., et al. (2002), The 26 December (boxing day) 1997 sector collapse and debris avalanche at Soufrière Hills Volcano, Montserrat, 21,

363–408
Yamada, M., H. Kumagai, Y. Matsushi, and T. Matsuzawa (2013), Dynamic landslide processes revealed by broadband seismic records,

Geophys. Res. Lett., 40, 2998–3002, doi:10.1002/grl.50437.
Young, S. R., R. S. J. Sparks, W. P. Aspinall, L. L. Lynch, A. D. Miller, R. E. Robertson, and J. B. Shepherd (1998), Overview of the eruption of

Soufriere Hills Volcano, Montserrat, 18 July 1995 to December 1997, Geophys. Res. Lett., 25(18), 3389–3392.
Zhao, J., et al. (2012), Model space exploration for determining landslide source history from long-period seismic data, Pure Appl. Geophys.,

172, 389–413.

LEVY ET AL. FRICTION WEAKENING IN GRANULAR FLOWS 7557

http://dx.doi.org/10.1029/2010JF001734
http://dx.doi.org/10.1002/grl.50437

	Abstract
	References


<<
  /ASCII85EncodePages false
  /AllowTransparency false
  /AutoPositionEPSFiles true
  /AutoRotatePages /All
  /Binding /Left
  /CalGrayProfile (None)
  /CalRGBProfile (ECI-RGB.icc)
  /CalCMYKProfile (Photoshop 5 Default CMYK)
  /sRGBProfile (sRGB IEC61966-2.1)
  /CannotEmbedFontPolicy /Warning
  /CompatibilityLevel 1.6
  /CompressObjects /Off
  /CompressPages true
  /ConvertImagesToIndexed true
  /PassThroughJPEGImages false
  /CreateJobTicket false
  /DefaultRenderingIntent /Default
  /DetectBlends false
  /DetectCurves 0.1000
  /ColorConversionStrategy /sRGB
  /DoThumbnails false
  /EmbedAllFonts true
  /EmbedOpenType false
  /ParseICCProfilesInComments true
  /EmbedJobOptions true
  /DSCReportingLevel 0
  /EmitDSCWarnings false
  /EndPage -1
  /ImageMemory 524288
  /LockDistillerParams false
  /MaxSubsetPct 100
  /Optimize true
  /OPM 1
  /ParseDSCComments true
  /ParseDSCCommentsForDocInfo true
  /PreserveCopyPage true
  /PreserveDICMYKValues true
  /PreserveEPSInfo false
  /PreserveFlatness false
  /PreserveHalftoneInfo false
  /PreserveOPIComments false
  /PreserveOverprintSettings true
  /StartPage 1
  /SubsetFonts true
  /TransferFunctionInfo /Preserve
  /UCRandBGInfo /Remove
  /UsePrologue false
  /ColorSettingsFile ()
  /AlwaysEmbed [ true
  ]
  /NeverEmbed [ true
    /Courier
    /Courier-Bold
    /Courier-BoldOblique
    /Courier-Oblique
    /Helvetica
    /Helvetica-Bold
    /Helvetica-BoldOblique
    /Helvetica-Oblique
    /Symbol
    /Times-Bold
    /Times-BoldItalic
    /Times-Italic
    /Times-Roman
    /ZapfDingbats
  ]
  /AntiAliasColorImages false
  /CropColorImages false
  /ColorImageMinResolution 300
  /ColorImageMinResolutionPolicy /OK
  /DownsampleColorImages true
  /ColorImageDownsampleType /Bicubic
  /ColorImageResolution 300
  /ColorImageDepth -1
  /ColorImageMinDownsampleDepth 1
  /ColorImageDownsampleThreshold 1.00000
  /EncodeColorImages true
  /ColorImageFilter /DCTEncode
  /AutoFilterColorImages true
  /ColorImageAutoFilterStrategy /JPEG
  /ColorACSImageDict <<
    /QFactor 0.76
    /HSamples [2 1 1 2] /VSamples [2 1 1 2]
  >>
  /ColorImageDict <<
    /QFactor 0.76
    /HSamples [2 1 1 2] /VSamples [2 1 1 2]
  >>
  /JPEG2000ColorACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 15
  >>
  /JPEG2000ColorImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 15
  >>
  /AntiAliasGrayImages false
  /CropGrayImages false
  /GrayImageMinResolution 300
  /GrayImageMinResolutionPolicy /OK
  /DownsampleGrayImages true
  /GrayImageDownsampleType /Bicubic
  /GrayImageResolution 300
  /GrayImageDepth -1
  /GrayImageMinDownsampleDepth 2
  /GrayImageDownsampleThreshold 1.00000
  /EncodeGrayImages true
  /GrayImageFilter /DCTEncode
  /AutoFilterGrayImages true
  /GrayImageAutoFilterStrategy /JPEG
  /GrayACSImageDict <<
    /QFactor 0.76
    /HSamples [2 1 1 2] /VSamples [2 1 1 2]
  >>
  /GrayImageDict <<
    /QFactor 0.76
    /HSamples [2 1 1 2] /VSamples [2 1 1 2]
  >>
  /JPEG2000GrayACSImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 15
  >>
  /JPEG2000GrayImageDict <<
    /TileWidth 256
    /TileHeight 256
    /Quality 15
  >>
  /AntiAliasMonoImages false
  /CropMonoImages false
  /MonoImageMinResolution 1200
  /MonoImageMinResolutionPolicy /OK
  /DownsampleMonoImages true
  /MonoImageDownsampleType /Bicubic
  /MonoImageResolution 400
  /MonoImageDepth -1
  /MonoImageDownsampleThreshold 1.00000
  /EncodeMonoImages true
  /MonoImageFilter /CCITTFaxEncode
  /MonoImageDict <<
    /K -1
  >>
  /AllowPSXObjects true
  /CheckCompliance [
    /None
  ]
  /PDFX1aCheck false
  /PDFX3Check false
  /PDFXCompliantPDFOnly false
  /PDFXNoTrimBoxError true
  /PDFXTrimBoxToMediaBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  ]
  /PDFXSetBleedBoxToMediaBox true
  /PDFXBleedBoxToTrimBoxOffset [
    0.00000
    0.00000
    0.00000
    0.00000
  ]
  /PDFXOutputIntentProfile (None)
  /PDFXOutputConditionIdentifier ()
  /PDFXOutputCondition ()
  /PDFXRegistryName ()
  /PDFXTrapped /False

  /CreateJDFFile false
  /Description <<
    /ENU ()
  >>
  /Namespace [
    (Adobe)
    (Common)
    (1.0)
  ]
  /OtherNamespaces [
    <<
      /AsReaderSpreads false
      /CropImagesToFrames true
      /ErrorControl /WarnAndContinue
      /FlattenerIgnoreSpreadOverrides false
      /IncludeGuidesGrids false
      /IncludeNonPrinting false
      /IncludeSlug false
      /Namespace [
        (Adobe)
        (InDesign)
        (4.0)
      ]
      /OmitPlacedBitmaps false
      /OmitPlacedEPS false
      /OmitPlacedPDF false
      /SimulateOverprint /Legacy
    >>
    <<
      /AllowImageBreaks true
      /AllowTableBreaks true
      /ExpandPage false
      /HonorBaseURL true
      /HonorRolloverEffect false
      /IgnoreHTMLPageBreaks false
      /IncludeHeaderFooter false
      /MarginOffset [
        0
        0
        0
        0
      ]
      /MetadataAuthor ()
      /MetadataKeywords ()
      /MetadataSubject ()
      /MetadataTitle ()
      /MetricPageSize [
        0
        0
      ]
      /MetricUnit /inch
      /MobileCompatible 0
      /Namespace [
        (Adobe)
        (GoLive)
        (8.0)
      ]
      /OpenZoomToHTMLFontSize false
      /PageOrientation /Portrait
      /RemoveBackground false
      /ShrinkContent true
      /TreatColorsAs /MainMonitorColors
      /UseEmbeddedProfiles false
      /UseHTMLTitleAsMetadata true
    >>
    <<
      /AddBleedMarks false
      /AddColorBars false
      /AddCropMarks false
      /AddPageInfo false
      /AddRegMarks false
      /BleedOffset [
        0
        0
        0
        0
      ]
      /ConvertColors /ConvertToRGB
      /DestinationProfileName (sRGB IEC61966-2.1)
      /DestinationProfileSelector /UseName
      /Downsample16BitImages true
      /FlattenerPreset <<
        /PresetSelector /MediumResolution
      >>
      /FormElements true
      /GenerateStructure false
      /IncludeBookmarks false
      /IncludeHyperlinks false
      /IncludeInteractive false
      /IncludeLayers false
      /IncludeProfiles true
      /MarksOffset 6
      /MarksWeight 0.250000
      /MultimediaHandling /UseObjectSettings
      /Namespace [
        (Adobe)
        (CreativeSuite)
        (2.0)
      ]
      /PDFXOutputIntentProfileSelector /DocumentCMYK
      /PageMarksFile /RomanDefault
      /PreserveEditing true
      /UntaggedCMYKHandling /UseDocumentProfile
      /UntaggedRGBHandling /UseDocumentProfile
      /UseDocumentBleed false
    >>
  ]
>> setdistillerparams
<<
  /HWResolution [600 600]
  /PageSize [612.000 792.000]
>> setpagedevice