A physical model of the high‐frequency seismic signal generated by debris flows

We propose a physical model for the high‐frequency (>1 Hz) spectral distribution of seismic power generated by debris flows. The modeled debris flow is assumed to have four regions where the impact rate and impulses are controlled by different mechanisms: the flow body, a coarser‐grained snout, a snout lip where particles fall from the snout on the bed, and a dilute front composed of saltating particles. We calculate the seismic power produced by this impact model in two end‐member scenarios, a thin‐flow and thick‐flow limit, which assume that the ratio of grain sizes to flow thicknesses are either near unity or much less than unity. The thin‐flow limit is more appropriate for boulder‐rich flows that are most likely to generate large seismic signals. As a flow passes a seismic station, the rise phase of the seismic amplitude is generated primarily by the snout while the decay phase is generated first by the snout and then the main flow body. The lip and saltating front generate a negligible seismic signal. When ground properties are known, seismic power depends most strongly on both particle diameter and average flow speed cubed, and also depends on length and width of the flow. The effective particle diameter for producing seismic power is substantially higher than the median grain size and close to the 73rd percentile for a realistic grain size distribution. We discuss how the model can be used to estimate effective particle diameter and average flow speed from an integrated measure of seismic power. © 2019 The Authors. Earth Surface Processes and Landforms Published by John Wiley & Sons Ltd. © 2019 The Authors. Earth Surface Processes and Landforms Published by John Wiley & Sons Ltd.


Introduction
With the development of global and local seismic networks, seismology has become a promising tool to obtain quantitative information about surface processes that are difficult to observe, such as landslides and debris flows (e.g., Kanamori and Given, 1982;Hibert et al., 2011;Ekström and Stark, 2013;Kean et al., 2015), bedload flux and turbulent flow in rivers (e.g., Burtin et al., 2011;Tsai et al., 2012;Gimbert et al., 2014) and crack formation in granular media (Michlmayr et al., 2012). For instance, effort has been made to understand why the amplitude of the envelope of the high-frequency (> 1 Hz) seismic signal generated by gravitational flows typically has a cigar shape with an emergent arrival, a maximum, and a long decay (e.g., Suriñach et al., 2005;Schneider et al., 2010). Fluctuating basal forces, thought to be primarily due to multiple impacts of particles with the bed (e.g., Iverson, 1997;Yohannes et al., 2012;McCoy et al., 2013;Hsu et al., 2014), are likely the source of these high-frequency seismic signals produced by granular flows. Suriñach et al. (2005) first suggested that the maximum seismic amplitude is reached when the landslide is closest to the seismic station because high-frequency waves attenuate rapidly with distance (Aki and Richards, 2002). Schneider et al. (2010) and Hibert et al. (2017), however, showed that the envelope amplitude corre-lates with the frictional work rate, i.e. the energy dissipated due to friction (Schneider et al., 2010) and with the momentum of the center of mass (Hibert et al., 2017). From coupled measurements of debris flow stage and seismic amplitude, Kean et al. (2015) observed that the snout of debris flows is thicker and generates stronger high-frequency seismic amplitudes than the flow body. In addition, high-frequency seismic waves are affected by complex bed topography (e.g., turns, change of bed slope) (Favreau et al., 2010;Moretti et al., 2012;Allstadt, 2013), by strong energy dissipation in heterogenous, fragmented media (Aki and Richards, 2002) and by the presence of an erodible bed that can lower the radiated high-frequency seismic energy (Kean et al., 2015). All of these studies show that seismic signals generated by mass failures are difficult to interpret because they depend on multiple parameters that cannot easily be separated from each other.
To date, only two physical models have been proposed to explain high-frequency seismic signals caused by granular flows, like debris flows, in terms of flow characteristics such as flow thickness, speed, particle size and slope angle (Kean et al., 2015;Lai et al., 2018). Both models modify a previous model for the high-frequency seismic signals caused by bedload transport (Tsai et al., 2012), and both models therefore assume stochastic impacts with the channel bed produce high-frequency seismic surface waves that travel to a seismic station nearby. The Kean et al. (2015) model does not predict absolute amplitudes for the seismic signal (or impact rates) and instead assumes that a reference flow is measured for which flow properties are known. In contrast, the Lai et al. (2018) model assumes that impacts are from a boulder-rich snout impacting on a rough channel bed and that this signal dominates the observed record.
Despite the models of Kean et al. (2015) and Lai et al. (2018), several important questions remain unanswered: (1) How much of the high-frequency seismic signal does the snout of a debris flow generate compared to the main body of the flow? (2) Under what conditions of downslope length, average particle size and distance between the debris flow and the seismic station does the snout or front dominate? (3) More generally, what characteristics of a debris flow have the strongest influence on the seismic amplitude spectra and which of these characteristics can be inferred from measurements of these spectra? For example, does the complex geometry of impacts not considered by previous models affect the predicted seismic signal significantly?
Our goal is to build on Lai et al. (2018) to more fully describe a physical model for the high-frequency (> 1 Hz) seismic signal generated by particle impacts in debris flows. The seismic amplitude generated by a debris flow is not simply the sum of the seismic amplitudes caused by multiple particle impacts because these signals can interfere constructively or destructively with each other. However, Tsai et al. (2012) showed that the power spectral density (PSD), i.e. the distribution of seismic power as a function of frequency, generated by multiple incoherent particle impacts on the bed of rivers could be evaluated from the sum of the seismic energies emitted by the impacts. This approach requires one to know, for a given particle diameter, the basal force per particle impact and the rate at which the particles impact the bed. We follow the same approach here and express the basal force per impact and rate of particle impacts in the different regions of debris flows. Our approach makes many simplifying assumptions not meant to represent the full physics of debris flows. Instead, our approach is to simplify the physics to only include the pieces that are most important for quantifying the magnitude of high-frequency ground motions. In future work, it may be interesting to try to relax some of the simplifying assumptions made, but we believe it is useful to start with the simplified model described here.
The model is presented in the next section, including a conceptual model, a model for the seismic wave generation, and models for how impact rates and basal impulses are calculated. In the Results Section, we determine which region of a debris flow dominates the generated seismic signal and the conditions for which it does so, and we discuss how to infer debris flow parameters from a measure of the emitted seismic power.

Conceptual model of debris flow structure
We focus on the high-frequency seismic signal of debris flows, a dense and agitated channelized flowing mixture of particles saturated with water, that are often triggered by rainfall or snowmelt, and which propagates downslope rapidly (e.g., Iverson, 1997). We propose the conceptual model shown in Figure 1a, in which a debris flow is assumed to be composed of four successive regions: (1) the body, (2) snout, (3) snout lip, and (4) saltating front, each of which generate seismic signals through incoherent impacts on the channel bed but may do so in different manners. Each region is assumed to have a statistically constant thickness, h, and statistically constant speed, u x .z/, in the downslope, x, direction that varies with depth (z direction). While debris flows often accelerate and gain mass in steep areas, and decelerate and deposit in lowland areas, assuming a velocity that is statistically constant but which varies stochastically is a simple but reasonable starting point. This assumption allows for excitation of high-frequency ground motion while not requiring detailed complex topography to calculate and is therefore perhaps the simplest situation in which high-frequency ground motions can be calculated. The flow is also assumed to have a constant bulk density, , and to propagate through a constant width channel at constant slope angle, Â , with respect to the horizontal. While assuming constant flow characteristics is not completely realistic, making more complex assumptions would not be expected to significantly affect the ground motions that are the focus of this work.
We assume that the stochastic component of velocity is described by an additional fluctuating component of velocity, u 0 , in a random direction (see Figure 1a), with root mean square (RMS) fluctuating speed given by ıu D p < u 02 >, a measure of the intensity of particle agitation in a granular flow, which is also sometimes described as a granular temperature (e.g., Ogawa et al., 1980). Without a fluctuating component of speed, the downslope flow would theoretically be accommodated purely in steady-state shear, and therefore not generate high-frequency ground motions. The rate of impact of each particle on the channel bed is determined by the characteristic fluctuating speed near the bed, ıu, divided by a characteristic length scale between bed impact events, s. We initially distinguish ıu, which is caused by impacts, from the average downslope speed, u x .z/, but later show how they could be closely related.
In our model, we distinguish the snout from the main flow because granular flow experiments and field observations have shown that the snout of debris flows contains a higher concentration of coarse particles than the body (e.g., Iverson, 1997;Stock and Dietrich, 2006;Johnson et al., 2012;Hsu et al., 2014;de Hass et al., 2015). During flow propagation, coarse particles are pushed upwards, because only fine particles can fall into open voids, and are transported towards the flow front because flow speeds are higher close to the free surface. Therefore, particle impacts in the flow snout impart larger forces on the bed, causing greater bedrock erosion (e.g., Stock and Dietrich, 2006;Hsu et al., 2008;). In our model, the snout is assumed to have a larger effective particle diameter, D snout , and a potentially different thickness, h snout , compared to those of the main flow body (D body and h body , respectively), but it is otherwise indistinguishable from the body.
At the downslope boundary of the snout, hereafter called the "snout lip," particles can fall off the snout, through air, and impact the bed. These particles impact the bed with large momentum because their downslope speed is close to the front speed and they fall a distance comparable to the total thickness of the flow. The downslope length of the snout lip is assumed to be a characteristic particle diameter, D snout , in the snout. After impacting the bed, particles at the lip are either (1) deposited at the flow base, in side levees, or re-entrained (e.g., Iverson, 1997;Thornton and Gray, 2008;Johnson et al., 2012), or (2) may bounce and roll downslope, in front of the snout, as reported by Iverson (1997) in debris-flow experiments with a relatively smooth bed. These saltating particles, which we call the "saltating front," can in some cases rebound higher than the flow thickness. Consequently, the saltating front could be responsible for some of the high-frequency radiated seismic energy. However, neither the snout lip nor the saltating front The snout has different characteristic thickness h snout and particle diameter D snout than the body (h body and D body ). Particles in the flow have an average downslope speed u x .z/ and an average fluctuating speed ıu.z/, and the volumetric particle fraction is . In the lip, particles fall in front of the snout and impact the bed. The saltating front is constituted of particles that rebound on the bed independently of each other in front of the coherent flow. (b) Schematic of the relevant angles during the impact of a particle on the bed roughness. The T-z plane is the plane of impact. The impact angles are polar angle˛with respect to the z axis, and azimuthal angle˛T with respect to the downslope x direction (see Appendix A). During the impact, the normal impact speed u N changes sign and is reduced by a factor of e b , the coefficient of restitution at the bed. (c) Schematic showing the debris flow path used for PSD computation. The position of particles in the channel is noted .x, y/ and x D 0 corresponds to the closest distance r 0 from the seismic station. [Colour figure can be viewed at wileyonlinelibrary.com] are expected to cause significant seismic power for thin flows, where D is similar to h, since the energy of impacts occurring in these regions is only due to particles free-falling from significant heights.
We focus on the high-frequency signals generated by particle impacts only. Other collision-related interactions at the scale of coherent grain clusters (Lee et al., 1974;Iverson and LaHusen, 1989) or at the flow scale, control the bulk static stresses the flow applies on the ground (Iverson, 1997;Estep and Dufek, 2012) and fluctuations in these stresses can generate high-frequency (> 1 Hz) seismic signals. Fluctuations in bulk static stresses can be caused by the development of secondary roll waves, i.e. coarse-grained wavefronts with various amplitudes, speeds and periods, on the flow surface (Zanuttigh and Lamberti, 2007;Iverson et al., 2010), or by variations in normal flow acceleration due to changes in bed topography (see supplementary materials). Because we consider a debris flow propagating on a stochastically rough bed with no turns and known large bumps, we do not model fluctuations of deterministic stresses caused by flow acceleration. We expect our conceptual model to describe the first-order high-frequency seismic spectrum generated by a propagating debris-flow surge even with roll waves, as long as the variability in flow thickness statistically averages out over time. We also do not discuss the low-frequency seismic signals expected to be caused by large scale mass acceleration. However, the statistical approach as described here could be coupled with inversions of low-frequency (< 0.1 Hz) signals generated by large-scale landslide movement over a complex topography (e.g., Allstadt, 2013).
Debris flows are often dense (Iverson, 1997) and thus particles closer to the bed have a higher chance of impacting the bed, and we assume that only particles directly above the bed can impact the bed. Particle collisions within the flow are assumed to only transmit their momentum to other particles, which can in turn impact the bed. We consider a debris flow that propagates directly on rigid bedrock because the transmission of impact forces through a layer of sediments is complex and still not well understood. For example, the presence of an underlying sediment bed can dramatically reduce the high-frequency seismic amplitude generated by a debris flow, by several orders of magnitude, compared to the same flow on bedrock (McCoy et al., 2013;Kean et al., 2015;Bachelet et al., 2018). We thus expect our model to overestimate the amplitude for debris flows over a highly attenuating sediment bed.
Overall, the aim of the present model is to give an order-of-magnitude accurate prediction of the high-frequency (> 1 Hz) seismic ground motion generated by a debris flow, provided that the flow occurs on bedrock, the elastic parameters are known, and debris flow characteristics (e.g., speed and thickness) do not vary significantly during the time period that the seismic spectrum is evaluated. While our assumptions are an idealization of real debris flows, we believe that additional complexity is currently unnecessary to model the main high-frequency seismic signals generated and could be incorporated later within a similar framework. With this simple model, we can understand how the amplitude of the seismic signal depends on the first-order flow characteristics, and it provides a base upon which future models could be built.
Seismic power spectral density generated by a debris flow Our goal is to predict the seismic ground motion, or equivalently to predict the frequency-domain power spectral density (PSD) of the seismic ground motion. With the conceptual model described above, the total PSD generated by a debris flow, PSD tot , contains contributions from the four successive regions of the flow PSD tot D PSD body C PSD snout C PSD lip C PSD saltating . (1) Within each region, we follow Tsai et al. (2012) and assume that the amplitude of the impact force, F j , in direction j is the same for every impact of particles of the same diameter, D, and that the time between each particle impact is random. We later modify this assumption to account for the complexity of particle motions with impacts at different orientations and speeds in Section 3. The Fourier transform of the force Q F j .f / is unchanged by the sum of the impacts and its amplitude increases as the square root of the number of impacts N as (Tsai et al., 2012). The PSD of the vertical component of ground velocity generated at distance, r, by each region of the flow per unit particle diameter, denoted here PSD.f , r, D/ (or PSD n .f , r, D/ for each specific region n), is then computed from the sum of the squared ground speeds radiated by the impacts over the surface area of the flow region, where R impact is the rate of particle impact per unit surface area of the debris flow per unit diameter D (in units of m 3 s 1 ). Q V z1 .f , r/ is the Fourier transform of the vertical component of the ground speed (in units of m s 1 /Hz) generated at distance r D p x 2 C y 2 by each impact, where x and y are the downslope and transverse directions, respectively. L and W are the debris flow downslope extent and average width, respectively. We note that our definition of R impact differs from the R i defined by Lai et al. (2018), which is the total rate of impact (number per unit time). The total PSD generated by the debris flow region, denoted PSD n .f , r/ (in units of m 2 s 2 /Hz), is obtained by integrating the PSD n .f , r, D/ in Equation (2) over the particle size distribution, The vertical ground speed generated by a single particle impact is given as a function of frequency, f , by where Q F j1 is the Fourier transform of the single force in direction j D x, y, z, and Q G jz are the Fourier transforms of the vertical displacement Green's functions due to a force in direction j (Aki and Richards, 2002). The Green's function can either be measured (e.g., with hammer tests) or approximated. Since gravitational flows on the Earth's surface such as debris flows are expected to mostly generate surface waves (Dammeier et al., 2011), a simple approximation is that of surface-wave propagation in a one-dimensional (1D) structure (e.g., Tsai et al., 2012;Gimbert et al., 2014) with k, the wave number, 0 , the rock density at the surface, v c and v u , the phase and group speeds, respectively, and Q, the quality factor. N jz are dimensionless numbers determined by the velocity structure that control the relative amplitude of the components Q G jz and have typical values of N zz D 0.6, N xz D 0.8 cos ' and N yz D 0.8 sin ', where ' is the source-station azimuth (Gimbert et al., 2014). Note that the horizontal components of ground motion PSD can just as easily be modeled within the same framework by replacing the Green's functions in Equation (4) with those for the component of ground motion of interest. In cases where a velocity model is available, Equation (5) could be replaced by a numerical calculation of the Green's function.
Using an idealized straight channel with a debris flow located at distance r 0 from a seismic station ( Figure 1c) leads to our final model. The position in the channel is denoted .x, y/ and x D 0 is at the closest distance from the seismic station. The downslope boundary, x max .t/, of a specific region of the debris flow at time t coincides with the upslope boundary, x min .t/, of the next region in the downslope direction. The lip is assumed to have a negligible downslope extent and is located at the downslope boundary of the flow snout x snout max .t/. The total seismic power recorded at the seismic station at time t is then where index n stands for body, snout, lip or saltating front.
Given this framework, in order to evaluate the total PSD, the main challenge is to determine the rate of particle impacts R impact and the Fourier transform of the forces Q F j1 within each region of the flow. After determining these in the following subsection, and comparing them to each other (see Results), we then evaluate the effect of various debris flow characteristics on the PSD generated (see Results).

The body and the snout
We assume that the impact rates and impulses are governed by the same mechanisms in the body and the snout of the debris flow. The only difference between the two regions is that we assume different particle diameters D body and D snout , respectively. We recognize that the average grain size difference may imply differences in flow regimes, and this possibility is considered in the next section.
Rate of particle impact in the body and snout As described earlier, we assume that the body and snout of the flow has an RMS fluctuating speed given by ıu, and we assume impacts are caused by these fluctuating departures from the steady state flow. If is defined as the volumetric particle fraction (the remainder being assumed to be fluid), and p.D/ is the particle size probability distribution (in units of m 1 ), then there are about p.D/=D 2 particles per meter square of the bed, and so the total rate of particle impacts per unit surface area of the debris flow, per unit particle diameter, is given by where s is the characteristic length scale between bed impact events. If it were assumed that all particles in the flow have the same diameter, D 0 , the distribution p.D/ would be a delta function, p.D/ D ı.D D 0 /, and integration of R impact over the grain size distribution would yield R R impact dD D ıu=.sD 2 0 /. To fully determine R impact through Equation (7), we still require an estimate of ıu and s, which depend on the dynamics of the flow. Next, we consider two end-member cases for the flow dynamics, with (1) relatively thin plug-like flows (e.g., Turcotte and Schubert, 2002) for which grain sizes, D, are comparable to the flow thickness, h, i.e. h 1 10D, which is common for debris flow snouts, and (2) relatively thick flows with small grains for which D h. For the thin flow case (h 1 10D), we envision large clasts being pushed and dragged along a rough channel bed in a 'washboard' fashion, equivalent to Lai et al. (2018). In this case, the main cause of fluctuations in velocity is that the clasts impact roughness elements of the channel bed, forcing the velocity to depart from the average flow speed ( N u x ). If bed roughness has a length scale D b , then most impact directions will result in significant perturbations to N u x , implying that ıu N u x . Thus, the impact rate of each grain can be estimated as (i.e. ıu D N u x and s D D b in Equation (7)). Lai et al. (2018) modeled the high-frequency seismic signal from the major debris flows that occurred in Montecito, California, on 9 January 2018, by using Equation (8) applied to the boulder snout in which boulder sizes were several meters in diameter. Boulder-rich debris flows (e.g., Iverson, 1997;Stock and Dietrich, 2006) fall into this thin-flow category, whereas mudflows and the often water-rich tail-end of debris flows likely do not.
For the thick flow limit (D h), we assume that u x .z/ has a convex velocity profile of the form with C a function of the flow parameters and p 1. This form for u x .z/ includes the possibilities of a Bagnold profile (p D 3=2), a Newtonian viscous profile (p D 2), and a linear profile (p D 1) as possible special cases, which have all been proposed for granular flows (e.g., Silbert et al., 2001;Cassar et al., 2005;Johnson et al., 2012). With this velocity profile, a particle of size D near the bed would be expected to have an average velocity given approximately by the profile velocity at its center of mass, i.e.
The fluctuating speed of particles at the bed is caused by interactions with the bed roughness, and likely scale with the magnitude of the particle velocity, or ıu u x , and like the thin-flow limit, s D D b . The impact rate for each grain then is If the average flow speed N u x can be measured, an uncertainty on the shape of the flow profile from linear (p D 1) to viscous (p D 2) leads to an uncertainty of a factor of 1.5 on the basal speed u x .D=2/ and thus on the fluctuating speed ıu. Since C does not enter Equation (11) explicitly, for our purposes it is not important how C depends on the flow dynamics.
There are potentially other causes for ıu other than direct interactions with the bed roughness. For example, the kinetic theory of dense granular gases (e.g., Hutter, 1993;Rao and Nott, 2008;Andreotti et al., 2013) predicts that if local production of kinetic energy is balanced by local dissipation of energy through shearing, then ıu near the bed can be estimated as where g 0 D . s f /g= s is the reduced gravitational acceleration, with s and f the particle and interstitial fluid densities, respectively, and e is the coefficient of restitution of the particles within the flow (e < 1 is assumed in this expression, following conservation of energy). (See Supplementary Material for details, including discussion of the appropriate s for this model (Lun et al., 1984;Jenkins and Askari, 1999;Lee and Huang, 2012).) For typical debris flows, we expect the bed roughness to be sufficiently large that it dominates over particle agitation given by Equation (12), but debris flows with very smooth channels may be in such a different regime.
To summarize, Equation (8) or Equation (11) can now be inserted for ıu=s in Equation (7) to determine both R snout impact and R body impact . Basal impulses If a particle impacts the bed with a normal speed u N forming an angle of impact˛with respect to the normal to the bed and an angle˛T with respect to the downslope direction (Figure 1b), and if we assume that the impact is instantaneous (compared to the frequencies of interest, see Tsai et al. (2012)), we can follow Tsai et al. (2012) and write the Fourier transforms of the tangential, Q F x1 and Q F y1 , and normal, Q F z1 , impact forces on the bed as impulses where m is the particle mass and e b is the basal coefficient of restitution, being 0 for a fully inelastic impact and 1 for a fully elastic impact. We assume that the impacts are instantaneous because seismometers are typically sensitive to waves of periods much longer than the duration of a particle impact, which is from 1 s to 1 ms, for particle diameters ranging from 1 mm to 1 m (Hertz, 1882;Johnson, 1985). We can neglect true viscous damping of the impact forces by the fluid because it is expected to only affect the finest particles (Lamb et al., 2008), which do not radiate significant seismic energy. The microscopic impacts that we model cause high apparent viscosities in particle-water mixture flows, for which damping would be higher. This damping from high apparent viscosities therefore should not be separately considered. Our approach to compute the seismic power PSD.f , r, D/ in Equation (2) supposes that the Fourier transform of the forces, i.e. the impulses I j , are the same for each impact of particles of the same diameter. This is not true in reality because particles have a complex motion and geometry and could impact the bed with various possible speeds u N and orientations, i.e. angles˛and˛T . Therefore, we evaluate characteristic impulses N I j per particle impact, defined as the root mean square (RMS) of the impulses I j in Eqs. (13)-(15), for all possibles values of the particle speed, u N , and impact angles,ą nd˛T . We expect that the predicted basal impulses N I j should be within a factor of two of the real impulses even if friction (which we ignore) is significant. If the friction coefficient is large, e.g. > 0.3, Maw et al. (1976) showed that the friction force changes sign during the impact due to the inertia of the impactor. The integral of the friction force over the impact duration (the tangential impulse) is therefore negligible because the positive and negative parts of the force cancel. In contrast, if friction is small, e.g. < 0.3, then the additional friction impulse I T that would need to be added to the normal impulse I N D .1 C e b /mu N in the expressions of the impact impulses I x , I y and I z is, at maximum, equal to I N and thus at most 30% of I N . In the following, we tolerate an uncertainty of 30% on the impulses and we only consider the normal impulse I N , as in Equations (13)-(15).
With the characteristic basal impulses N I j defined as above, it remains only to relate the average impact speed, u N , with other flow characteristics. Since u N is assumed to be due to an average speed plus a fluctuating speed, and all of these speeds scale with average flow speed N u x (Equations (8) and (11)), then u N should also scale with N u x . A more precise estimate can be obtained by directly calculating the impact speeds due to a distribution of impacts with a given average particle speed u x plus a fluctuating speed ıu (initially assumed independent of u x ) that is uniformly distributed over all possible angular directions (see Figure 2). After performing the RMS averaging (see Appendix A), the result can be summarized as with f j a dimensionless function of the speed ratio ıu=u x , obtained by computing numerically N I j =.1 C e b /mu x (Figure 3a; see also Appendix A). Good fits to these dimensionless functions are f x 0.090 p .ıu=u x / 2 =0.38 C 1, f y 0.052 p .ıu=u x / 2 =0.13 C 1, and f z 0.23 p .ıu=u x / 2 =0.18 C 1 (black dashed lines in Figure 3a). When ıu u x (Equations (8) or (11)) then u x f x =ıu 0.17, u x f y =ıu 0.15 and u x f z =ıu 0.59 (similar to the analytically determined limiting values u x f x =ıu D u x f y =ıu D 0.146 and u x f z =ıu D 0.539 as discussed in Appendix A for when u x D 0). Thus, the RMS vertical impulse is approximately 3-4 times the horizontal impulses, and about 40% smaller than when approximating u N D ıu D u x and ignoring the geometrical complexities.

The snout lip
In the snout lip region, particles may fall from the front of the snout and impact the bed (Figure 1a). For a particle that falls freely from height z, the normal speed before impact is u z D p 2gz cos Â (neglecting air resistance). If we further assume that all values of z from 0 to h are equally likely, then 2gz cos Â dz D q 2g 4h 9 cos Â . We can therefore use an average height given by z D 4h=9 0.5h so that an average normal impact speed is u z p gh cos Â . The final downslope speed is equal to the sum of the average flow speed, N u x , and the speed gained during free fall, i.e. u x D N u x C p gh sin Â . The total speed of impact in the lip is u lip impact D p u 2 x C u 2 z . The free fall height is negligible when D h so that one would expect PSD lip 0 in the thin flow limit, and this contribution can then be ignored (discussed in more detail in the Results Section).  this is plotted as a function of u z =u x , where u x and u z are the particle downslope and normal impact speeds, computed from ballistic trajectories of the particle rebounding on the bed in the lip and front, respectively. The particle does not have a fluctuating component of speed ıu in these last two regions. Black dashed lines represent the fits of these functions given in the main text. The diameter D b of the bed bumps is assumed to be the same as the particle diameter D in the debris flow. These dimensionless functions are geometrical computations based on the possible angles of impact (see Appendix A) and depend only on the speed ratios, and not on flow parameters. [Colour figure can be viewed at wileyonlinelibrary.com] Rate of impact in the lip A particle impacts the bed only one time in the lip and is then entrained by the flow or moves into the saltating front ( Figure 1a). Because the flow snout is dense, we assume the falling particles follow each other closely so that the time between two impacts of successive particles is D=u lip impact , where u lip impact is the particle fall speed in the lip. The rate of impact per unit surface area of the lip, per unit particle diameter, is then, according to equation (7 Basal impulses in the lip The characteristic basal impulses per particle impact in the snout lip is found from the RMS average of impulses, as in Equation (16)

The saltating front
For thick flows, a saltating front may contribute seismic energy. A saltating front is not dense (Iverson et al., 2010) so that ballistic trajectories of particles can be considered independent of each other. Once a particle has impacted the bed in the snout lip with speed u impact , it rebounds with a speed u reb (Figure 1a). Because of bed roughness, the particle can rebound in multiple possible directions. Applying standard Newtonian mechanics to a particle in free fall with initial speed w reb and acceleration g cos Â , the time between the two impacts of the particle is t.˛/ D 2w reb =.g cos Â/ and the height of rebound above the bed is h reb .˛/ D w 2 reb =.2g cos Â/, where w reb D u reb O z is the vertical component of u reb . Characteristic time of flight t and rebound height N h reb of a rebounding particle can be evaluated as the RMS of t.˛/ and h reb .˛/, respectively, averaged over all values of˛.

Rate of impact in the saltating front
We deduce the rate of impact for a given particle rebounding on the bed as the inverse of the RMS of the successive times of flight t until the rebound height N h reb is smaller than the typical size D b =2 of the bed roughness. Therefore, the rate of particle impacts per unit surface of the saltating front, per unit particle diameter, is where s is the particle fraction in the saltating front. product of R impact and the squared impulse N I 2 z for different particle diameter D in the body, the snout, the lip and the saltating front of a debris flow composed of particles of the same diameter D. The equations to compute the impact rate (Equation (7)) and impulses (Equation (16)) in the flow body and the snout are the same. Results for the thin-flow model (solid line) and thick-flow model (dashed line) are indicated. (d) Normalized PSDs integrated over frequency f for the flow body, snout, lip and saltating front as a function of the particle diameter D. The PSDs are computed using a log-'raised cosine' particle size probability distribution p.D/ (Tsai et al., 2012) (black line) with median diameter D 50 D 0.3 m and standard deviation g D 0.5. Inset: Normalized cumulated particle size distribution that allows us to determine the Xth percentile particle diameter D X corresponding to the maximum of the PSDs in each region of the debris flow. The flow parameters are u x D 10 m s 1 , h D 1 m, s D 2500 kg m 3 , D 0.6, e D 0.5, Â D 10 ı and W D 10 m. The particle fraction in the saltating front is assumed to be s D 0.1 and the diameter of the particle forming bed roughness is assumed to be D b D D. [Colour figure can be viewed at wileyonlinelibrary.com] Basal impulses in the saltating front Calculating numerically the RMS of the successive impulses until N h reb is smaller than D b =2, we obtain the characteristic impulse per particle impact in the saltating front with j D x, y, z, and f saltating j is a function of the speed ratio u z =u x , where u x and u z are the downslope and normal particle speeds given in the previous subsection. f

Results of the Model
With expressions for the rate of particle impact R n impact and characteristic basal impact impulses N I n j determined for each of the different regions n of the debris flow, we now compare them in the next section. Specifically, we combine Eqs. (7), (17) and (19) for R n impact with Eqs. (16), (18) and (20) for Q F n j1 N I n j in Equation (6) to compute the PSDs generated by each region of a debris flow and investigate how the flow parameters control the generated PSDs.
We show model results for a simplified example debris flow, i.e. with parameters N u x D 10 m s 1 , h D 1 m, W D 10 m, s D 2500 kg m 3 , D 0.6, e D 0.5, and Â D 10 ı (Figures 4a and 4b). The volume particle fraction in the saltating front is assumed to be s D 0.1, and the diameter of the particles forming the bed roughness is assumed to be D b D D, as might be expected for debris flows that propagate over previous debris flow deposits of similar grain size. We believe the thin-flow model to generally be more appropriate for boulder-rich debris flows that generate large seismic signals, and henceforth we concentrate on predictions for this thin-flow case except when 'thick-flow' is noted.

Comparison of impact rates and basal impulses in a debris flow
The rate of impacts (per unit area per unit grain size) in the body/snout (where ıu D N u x , s D D b D D) scales with N u x =.D b D 2 /, whereas in the thick-flow model it scales with N u x =.hD b D/, and the rate of impact in the lip scales with N u x =D 3 . Thus, the thin-flow model and lip have an order of magnitude higher rate of impact than the thick-flow model over the range of D plotted in Figure 4a. The rate of impacts is much smaller in the saltating front than in the other regions due to long times of flight of the saltating particles between two impacts (see Figure 4a).
Since the impact speeds are all on the order of N u x (except for that of the thick-flow model, which scales as N u x D=h), impact impulses per impact, N I z , all scale with D 3 (except for that of the thick-flow model, which scales as D 4 and is an order of magnitude smaller) (see Figure 4b). The characteristic particle diameter is typically larger in the snout, lip and saltating front regions than in the main flow body and if, for example, D body D 0.1 m and D snout D D lip D D saltating D 0.5 m, N I z would be 125 times larger in the snout than in the flow body.

Effective particle diameter
The seismic power PSD.f , r, D/ (Equation (2)) varies as the product of R impact N I 2 z . In the flow body and snout, the impact rate R impact scales with particle diameter as D 3 if D b D. The squared impulse N I 2 z scales as m 2 ıu 2 D 6 to first order because m / D 3 and ıu is independent of D. Thus, for a constant mass flux rate, the seismic power PSD.f , r, D/ radiated per unit area of debris flow should be roughly proportional to D 3 (Figure 4c). Because of this strong dependence of the PSD on the third power of the particle diameter D, similar to the model of Tsai et al. (2012), the seismic signal generated by a debris flow should be dominated by the largest particles in the flow. In the thick-flow case, the same arguments suggest an even higher dependence of PSD on D 6 .
The particle diameter contributing the most to the seismic signal is determined by a tradeoff: the PSD is proportional to D 3 , whereas the fraction of these largest particles which cause a stronger signal is generally small. The sensitivity to the high-end tail of the grain size distribution implies that errors in this tail result in magnified errors in predicting the PSD. To avoid the disproportional contribution of unrealistically large particles on the PSD, we follow Tsai et al. (2012) and use a log-'raised cosine' probability distribution for p.D/, which has similar characteristics as a log-normal distribution but with truncated tails (see Figure 4d). The precise grain size with maximum PSD is strongly affected by the standard deviation of the truncated normal distribution, g , though this maximum always occurs significantly higher than the median grain size D 50 . For g D 0.5, similar to Tsai et al. (2012), the peak in PSD occurs at D 89 , the 89th percentile grain size. This differs slightly from the D 94 of Tsai et al. (2012) due to the model of Tsai et al. (2012) being for fluvial bedload transport, not debris flows. In the thick-flow case, the peak in PSD occurs at D 98.6 for the same g . In addition, we can also calculate the grain size for which the predicted power would be equivalent if all grains were of that size. With g D 0.5 fixed, R D 3 p.D/dD D 3 73 , so that a grain size distribution p.D/ D ı.D D 73 / would yield the same predicted power as using the actual p.D/ in the thin-flow limit, or similarly for p.D/ D ı.D D 86 / in the thick-flow limit. These are the grain sizes that should be used in theory that utilizes a single "average" D as in Lai et al. (2018). (Using D 89 instead of D 73 would result in a factor of 2.4 overprediction of the PSD.) We now evaluate what region of the debris flow dominates the seismic signal. We compute the total PSD generated by the debris flow as a function of frequency for various times during the debris flow, corresponding to different downslope positions of the flow body (see spectrograms in Figures 5 and 6). Separate contributions from each debris flow region are also shown as a function of frequency when the boundary between the body and the snout is closest to the seismic station (for x D 0, Figures 5a and 5d) and in plots of the PSD integrated over all frequencies (Figures 5c, 5f, 6a and 6d). We assume a debris flow composed of particles of the same diameter D body in the flow body and D snout > D body in the snout, lip and saltating front, and initially assume that L body =L snout D 10 and L saltating D L snout .
The PSDs decrease as the distance r 0 between the seismic station and the debris flow increases (compare Figures 5a and  5d), just as in Tsai et al. (2012) and Lai et al. (2018). Peak power is at lower frequencies as r 0 increases because high frequencies attenuate faster with distance (Tsai et al., 2012). For example, the maximum of the PSD generated by the flow body is -71.6 dB (i.e. max .PSD.f // 6.9 10 8 m 2 s 2 Hz 1 ) at f p 40 Hz for r 0 D 50 m ( Figure 5agh) and is -107.2 dB (i.e. max .PSD.f // 1.9 10 11 m 2 s 2 Hz 1 ) at f p 8 Hz for r 0 D 500 m (Figure 5bgh) given the parameters used. The relation between the peak frequency f p and distance r 0 for the flow snout is similar to that observed by Lai et al. (2018) for the same choice of flow parameters (Figure 5h). Note that the PSDs generated by the flow body and snout cover a wider frequency range than the PSDs generated by the lip and the saltating front when the station is close to the debris flow path (e.g., for r 0 D 50 m, Figure 5a), but all regions are within the same frequency range at farther distances (e.g., for r 0 D 500 m, Figure 5d). Figure 5 demonstrates that the PSD generated by the lip is always many orders of magnitude lower than that of the snout, so the lip can be ignored as an unimportant contributor. The main reason for this is that the area of the lip is small, but impact rates and speeds are similar to those in the snout. Similarly, the PSD generated by the saltating front is also orders of magnitude lower than the PSD generated by the snout, even when assuming the saltating front length is similar to the snout length (which is likely an overestimate for most natural debris flows), primarily because the particle fraction and impact rates are smaller. Even when the saltating front is closer to the seismic station, the snout PSD is still larger (Figure 5c), so we can also comfortably ignore the saltating front in the generation of seismic signals, except when the saltating front is much longer than the snout or for the thick-flow case (for which the snout signal is lower), or in the rare case that the average height of bounces in the saltating front would be significantly larger than the snout thickness. For the remainder of this work, we therefore concentrate solely on the signals produced by the snout and main body of the flow. While our model suggests that the snout always dominates over the main body of the flow when the station is close to the snout (e.g., for x D 0, Figure 5a and Figure 5d), Figure 5 also shows that the body dominates once the snout is sufficiently far downstream. Thus, the rise phase of the seismic amplitude, i.e. when the flow is propagating towards the seismic station, is always dominated by the snout whereas the decay phase, i.e. when the flow is propagating away from the seismic station, can sometimes be dominated by the flow body (unless the D snout =D body 1 or L body =L snout 1, in which case the snout always dominates). These findings are affected by the relative length ratio L body =L snout (Figure 6abc), the relative grain size ratio D snout =D body (Figure 6def) and the source-station distance (Figure 6g). The snout can be better discerned in the total debris flow spectrogram as the ratios L snout =L body and D snout =D body increase (Figures 6bc and 6ef). These findings are consistent with both observations of natural flows where coarser surge fronts generate higher amplitudes than later parts of the flows (e.g., Arattano and Moia, 1999;Kean et al., 2015) and supports an assumption central to the model of Lai et al. (2018) that modeling the signal from the snout models most of the signal.
How to retrieve debris flow parameters from seismic data Because the snout or the body of a debris flow likely dominate the generated seismic signal, we can generally ignore the signal from the lip and saltating front, and simply use the impact rate R impact D N u x p.D/=.D b D 2 / and impulse N I j D .1Ce b /m N u x f j derived for the body and snout (i.e. from Equations (7) and (16), respectively, assuming a thin flow) in the first two terms of the sum in Equation (6) This derived expression (21) of the PSD is similar to that obtained by Tsai et al. (2012) for the PSD generated by particle impacts in the bedload of rivers. The factor N 2 jz f 2 j (for j D x, y, z, see Equations (5) and (16)) is a function of the bed roughness and the impact directions. It equals 0.6 2 0.59 2 C 0.8 2 0.15 2 D 0.14 in the thin-flow limit (for which ıu D u x , Figure 3a). This factor was not present in the original model of Tsai et al. (2012) for simplicity, whereas Lai et al. (2018) use 0.6 2 to account for the vertical impact force. Since the vertical force does indeed dominate, accounting for geometrical complexities thus introduces approximately a factor of 0.59 2 D 0.35 on the PSD amplitude. Lai et al. (2018) also assumes the debris flow is dominated by the snout and that this snout is short compared to its distance from the instrument, preventing one from making a prediction for how the PSD changes when the main flow begins to dominate the signal, as analyzed in the present paper. We further note that if one were to use a single average source-station distance, r avg , as done in Lai et al. (2018), the appropriate way of defining that average would be as in Equation (22) where S is the area over which the average distance is to be computed, making r avg frequency dependent. The predictions of Lai et al. (2018) can thus be made equivalent to those here by defining r avg by Equation (23) and setting their w i D 0.59.1 C e b / N u x . The integrated PSD over frequency f and diameter D can be expressed as a function of the flow parameters by integrating equation (21) as iPSD.r 0 , t/ D N u 3 x D 3 73 W.1 C e b / 2 g.r 0 /, where g.r 0 / expresses the distance dependence caused by wave propagation (see Figure 7a). g.r 0 / depends on ground properties 0 , v c , v u , Q as The seismic ground motion thus depends significantly on the elastic parameters of the ground, as 2 0 v 3 c v 2 u , and has an additional exponential dependence on v u and Q coming from .r 0 , f / (see Equation (22), Tsai et al. (2012)). Therefore, accurate determination of these parameters is important even if one only wants to make PSD predictions that are order-of-magnitude accurate. The D 73 dependence would be different if the standard deviation of the grain size distribution were different (though it would always be larger than D 50 ). We also note that g.r 0 / depends approximately linearly on the length L n (i.e. Equation (22) is linear with dy if y constant or dy=y 1) so that the seismic signal scales linearly with the total areal extent W L n . Finally, a similar integrated expression can be obtained in the thick-flow limit iPSD.r 0 , t/ thick D 1 with the characteristic diameter D 73 replaced by D 86 in iPSD.r 0 , t/ thin . The pre-factor 1 8 .p C 1/ 3 D 86 h 3 D 86 h 3 1 in the thick-flow limit, for which h D 86 . Therefore, the high-frequency seismic signal generated by a debris flow in the thick-flow limit is expected to be much smaller than that generated by a debris flow in the thin-flow limit, for equivalent N u x , W and ground properties.
Equation (24) shows that the debris flow parameters that most affect the generated PSD are the average speed, N u x , and the effective particle diameter, D 73 (see Figure 7b), which both enter the PSD to the third power. The width, W, and length, L n , (either for the main flow or the snout) affect Equation (24) linearly. It should be noted, though, that the variables discussed are not independent of each other; for example, narrowing of the flow in a canyon potentially causes faster flow and larger debris to be entrained. Since only a small range of 0.1 < e b < 0.5 and 0.4 < < 0.6 are realistic (Iverson, 1997), e b and only have second-order effects on the integrated PSD (see Figure 7cd). If the effective particle diameter D 73 is known, the uncertainty on the average speed N u x from an observation of iPSD.r 0 , t/ is about 20% when the e b and are in the range quoted.
The particle size distribution, and thus D 73 , may be determined from sampling of previous debris flow deposits (e.g., Stock and Dietrich, 2006;Johnson et al., 2012), which then allow one to directly infer the average speed N u x from a measure of the average PSD if the estimates of D are accurate enough (Doyle et al., 2010). Alternatively, Lai et al. (2018) showed that the peak frequency of the observed high-frequency seismic signal is expected to depend on the average distance of the snout from the station, so that the shift in peak frequency with time could be used to indirectly infer the average speed N u x , though uncertainty exists regarding whether other physics can cause similar systematic shifts in frequency within this frequency band (Doyle et al., 2010). Allstadt (2013) also showed that the trajectory of the flow center of mass and thus its average speed N u x could potentially be determined from low-frequency (< 0.1 Hz) signals, though successful results so far are only for much more massive events. Yet another technique using several seismic stations could be used to determine N u x . For example, if two stations are located a few meters away along the debris flow path but at same distance r 0 from the path, the cross-correlation of the seismic signals recorded by these stations (e.g., Burtin et al., 2010) can constrain the time the debris flow front spends traveling from one station to the other and thus the average flow speed N u x (Arattano and Marchi, 2005). Debris flow speed can also be measured with non-seismic methods, for example with cameras, multiple laser altimeters, or centrifugal force estimates (e.g., McArdell et al., 2007;Takahashi, 2007;Arattano and Marchi, 2008). In such cases where the average flow speed is determined, a measure of the integrated seismic power iPSD.r 0 , t/ then allow us to determine the effective particle diameter D 73 in the snout or the body of the debris flow, depending on which region dominates.

Conclusions
We have proposed a physical model for the high-frequency (> 1 Hz) seismic signal generated by particle impacts in a debris flow in the inertial regime that extends the model of Lai et al. (2018). As in this earlier work, the spectral distribution of seismic power at a given distance from a debris flow depends on the rate of particle impacts per unit surface area of the flow and on the squared basal impact impulses. Unlike Lai et al. (2018), who considered just the snout of the debris flow, we model the entire debris flow structure and investigate their respective contributions to the expected seismic signal. We separate the debris flow into four successive regions, the body, the snout, the snout lip and the saltating front, where the impact rate and impulses per impact are controlled by different physical mechanisms. We also account for geometrical complexity of the bed in the computation of seismic power and predict a factor of 2 lower PSD for the same parameters. We confirm that the thin-flow limit (D=h 1) may generally be a more realistic assumption than the thick-flow limit (D=h 1) for natural debris flows (i.e., with meter-scale debris), unless grain sizes are extremely small and yet not suspended (e.g., for some mudflows). In this thin-flow limit, we find that the seismic signal is always dominated by the snout for incoming flows and dominated by either the snout or the body as the snout moves away from the station. We therefore support the previous assumption of Lai et al. (2018) that the early seismic signal due to debris flows is dominated by the snout.
Model results further support that the most important parameters for the seismic power generated by a debris flow are the particle diameter D and the average flow speed N u x , both to the third power (in the thin-flow limit), as in the model of Lai et al. (2018). Seismic power depends linearly on width W and length L of the snout and body. Other parameters, such as the particle fraction, coefficient of restitution and slope angle, have smaller second-order influences on the seismic amplitude ( 20% uncertainty), except indirectly, if they affect the size, speed or particle sizes of a debris flow. Only in the thick-flow limit does flow thickness h have any direct effect on seismic power, and in this case the seismic power has an even stronger dependence on particle diameter to the sixth power.
Though this model simplifies some of the well-known complexities of natural debris flows, it provides a starting point upon which more complex models can be built and provides the theoretical framework that is necessary if we ever hope to use seismic recordings of debris flows to obtain quantitative estimates of flow parameters. We also showed that the observed PSD has a strong dependence on the elastic parameters of the subsurface through which the waves are traveling (path effects). These need to be measured or controlled for in order to obtain information about the source. Even considering the path effects are known, our model demonstrates that several interrelated flow characteristics (flow width, flow velocity, particle size) all contribute to the signal generated simultaneously. This explains why past studies that have tried to fit simple models relating seismic amplitude to a single flow parameter such as discharge, were not able to produce good model fits (e.g., Lavigne et al., 2000). Progressing towards more quantitative applications of seismic debris flow monitoring requires that multiple contributing factors be taken into account simultaneously.