Measuring Basal Force Fluctuations of Debris Flows Using Seismic Recordings and Empirical Green's Functions

We present a novel method for measuring the fluctuating basal normal and shear stresses of debris flows by using along‐channel seismic recordings. Our method couples a simple parameterization of a debris flow as a seismic source with direct measurements of seismic path effects using empirical Green's functions generated with a force hammer. We test this method using two large‐scale (8 and 10 m3) experimental flows at the U.S. Geological Survey debris‐flow flume that were recorded by dozens of three‐component seismic sensors. The seismically derived basal stress fluctuations compare well in amplitude and timing to independent force plate measurements within the valid frequency range (15–50 Hz). We show that although the high‐frequency seismic signals provide band‐limited forcing information, there are systematic relations between the fluctuating stresses and independently measured flow properties, especially mean basal shear stress and flow thickness. However, none of the relationships are simple, and since the flow properties also correlate with one another, we cannot isolate a single factor that relates in a simple way to the fluctuating forces. Nevertheless, our observations, most notably the gradually declining ratio of fluctuating to mean basal stresses during flow passage and the distinctive behavior of the coarse, unsaturated flow front, imply that flow style may be a primary control on the conversion of translational to vibrational kinetic energy. This conversion ultimately controls the radiation of high‐frequency seismic waves. Thus, flow style may provide the key to revealing the nature of the relationship between fluctuating forces and other flow properties.


Introduction
Debris-flow and lahar monitoring sites commonly make use of along-channel seismic instrumentation for event detection, monitoring, flow characterization, and scientific research. Applications span many scales, from debris flows in small alpine catchments (e.g., Arattano & Moia, 1999;Badoux et al., 2009;Hürlimann et al., 2019;Kean et al., 2015;McArdell et al., 2007) to great volcanic lahars (e.g., Hadley & LaHusen, 1991;Lavigne et al., 2000;Leonard et al., 2008;Marcial et al., 1996;Suwa et al., 2011). The limitation of such seismic monitoring systems is that they are qualitative or site specific. Many investigators have observed correlations between seismic amplitude and flow properties such as flow height, sediment concentration, kinetic energy, momentum, and wetted area (e.g., Arattano & Moia, 1999;Coviello et al., 2019;Doyle et al., 2010;Hadley & LaHusen, 1991;Hikida et al., 2007;Suwa et al., 2000), and several have successfully made order-of-magnitude discharge estimates based on seismic data alone (e.g., Almeida, 2016;Andrade et al., 2006;Lavigne et al., 2000;Suwa et al., 2000). However, these estimates have all been empirical, calibrated against seismic recordings of past flows with known discharges and nearly always specific to one monitoring location. This is not an option for new systems where flows have not yet been recorded or are infrequent but potentially very large and deadly. Universally applicable relations between the seismic signal and flow characteristics that could be applied to any new monitoring site have been elusive.
Few monitoring sites rely solely on seismic methods, in part due to the lack of clear quantitative links between seismic signals and flow characteristics. By necessity, most sites combine seismic sensors with nonseismic monitoring methods such as tripwires, pendulums, cameras, flow depth sensors, and manned watch stations (e.g., Hürlimann et al., 2019;Lockhart & Murray, 2004;Marcial et al., 1996;Tungol & Regalado, 1996). Yet each alternative monitoring method presents its own logistical and practical challenges. For example, depth sensors only provide point measurements, video methods only work in daylight, tripwires require manual resets, and many systems have high power requirements (Arattano et al., 2008;Itakura et al., 2005;Suwa et al., 2011). Seismic methods are not subject to many of these limitations and have the potential to give information on distributed flow characteristics and basal processes that are difficult to measure by other means. Thus, a better way to link seismic recordings quantitatively to flow characteristics is clearly needed.
Limitations of currently available seismic methods result, in part, from a lack of theoretical understanding of how debris flows generate seismic waves. In turn, this understanding is limited by the wide variability in attenuation and scattering of seismic waves as they travel through the subsurface. Path effects are frequency dependent and are especially strong, spatially variable, and hard to quantify for the high-frequency (1 to hundreds of Hz) seismic waves that are preferentially observed for debris flows and commonly used in along-channel debris flow monitoring (Allstadt et al., 2018). Furthermore, debris flows are a complex seismic source. They are elongated in space, so many parts of the flow contribute to the seismic signal observed at a given site simultaneously, but each contribution to the signal is subject to different path effects. Flow characteristics also vary in space and time as do the characteristics of the bed material, which may also be subject to deformation, erosion, or deposition. All of these factors affect the characteristics of the forcing that generates the seismic signals (e.g., Allstadt et al., 2018;Cole et al., 2009;Doyle et al., 2010;Kean et al., 2015). Thus, parameterizing debris flows as seismic sources requires accounting for these complexities within a simplified framework.
In this paper, we present a framework for removing one of the biggest barriers to progress: the effect of the path between the source and seismic station. To do so, we combine direct measurements of path effects with a parameterization of spatially and temporally variable debris flows as seismic sources. To make direct measurements of path effects, we extend the empirical Green's function (EGF) technique originally used for earthquake applications (e.g., Aki, 1967;Hartzell, 1978) to along-channel debris flow monitoring. Green's functions relate a unit impulse in one location to the response at another location, providing information about the path effects between those two locations (Aki & Richards, 2002).
Whereas most EGF studies use small natural earthquakes as Green's functions to study rupture propagation and other details of larger earthquakes (e.g., Hartzell, 1978;Hough, 1997;Kim et al., 2016;Mueller, 1985) or use smaller landslide signals to study the dynamics of larger landslides (Zhang & He, 2019), we instead use a controlled active source in the form of a force hammer. Hammer sources have been used in active source seismology to obtain EGFs and image the subsurface for many decades (e.g., Reynolds, 2011). However, we are unaware of other studies that use hammer-based EGFs specifically for source studies and that use the time series of the forcing exerted by the hammer as we do in this investigation. The most conceptually similar method is that of Jolly et al. (2014) and Walsh et al. (2016), who dropped bags of rocks from helicopters to calibrate seismic amplitudes and energies for eruption, debris avalanche, and lahar monitoring.
In the following sections, we demonstrate and validate our technique by using large-scale laboratory experiments conducted at the U.S. Geological Survey debris-flow flume, in which 8-10 m 3 of wet debris flowed over a concrete substrate. Much work remains to extend these methods to larger scales and natural settings, but here we establish a foundation that serves as a starting point.

Seismic Wave Generation by Debris Flows
Debris flows generate seismic waves by exerting time-varying forces on the surface of the Earth. Debris flows are spatially elongated, and flow characteristics are variable in space and time over many different scales, implying that the basal forces are similarly variable. The observed high-frequency seismic signal represents the sum of signals from all contributing source areas.
Debris-flow basal normal and shear forces can be measured directly using force plates that provide information averaged over the area of the force plate. Fluctuating forces that can exceed the mean forces by more than an order of magnitude have been observed in force plate measurements of many natural and experimental debris flows (e.g., Hsu et al., 2014;Iverson, 1997;McArdell et al., 2007;McCoy et al., 2013). While any given debris flow likely exerts basal forces that vary over a wide range of timescales, only those forces that change over relatively short time scales (the fluctuating forces) will generate high-frequency seismic energy. Fluctuating boundary forces, F f (t), are defined as where F(t) is the instantaneous force and F t ð Þ is the mean time-averaged force. Dividing both sides by area of forcing (e.g., force plate area) in Equation (1 replaces force (F) with basal stress (σ). The larger the size of the force plate, the smaller the measured basal stress fluctuations since a larger force plate averages out the effects of more simultaneous individual basal interactions (Iverson, 1997).
The instantaneous force fluctuates because particles at the bed of a debris flow are agitated and oscillate as translational momentum is converted to fluctuation momentum due to interactions of particles and clusters of particles with adjacent particles and bed roughness (Iverson, 1997). As this shear work generates kinetic energy of random particle motions, the agitation of the flow increases, and higher agitation results in higher rates of energy dissipation through grain collisions (Iverson, 1997;Jenkins & Askari, 1999). Variations in granular flow dynamics result in differing rates of conversion of kinetic energy to seismic energy and are reflected in the resulting seismic signal (Farin et al., 2018).
The three simple theoretical debris-flow seismicity models that have been published (Farin et al., 2019;Kean et al., 2015;Lai et al., 2018) all parameterize the source entirely as a series of individual random particle impacts whose statistical distributions relate to flow characteristics such as flow width, length, and mean velocity. Though these models have not yet been thoroughly verified against independent data, this treatment of the fluctuating forcing as a stochastic process is in line with the findings of experimental and modeling studies that have found that the statistics of the distribution of fluctuating forces correlate, in some cases, with other bulk flow properties such as effective grain-size diameter, mean flow velocity, mean time-averaged basal forces, flow depth, and collisional energy (Ferguson et al., 2004;Gardel et al., 2009;Hsu et al., 2014;McCoy et al., 2013).
Individual particles in a granular flow respond to and transmit stresses from the surrounding flow (Iverson, 1997;Iverson & Vallance, 2001). Under some conditions, including within denser parts of flows with more enduring contacts between grains (conditions that may be met at the base of even rapid granular flows), the translation of forces from the body of the flow to the bed may occur through force chains (Campbell, 2006). Force chains are discrete, transient, filamentary structures that can locally focus basal forces many times greater than the weight of the material on a minority of basal particles. These force chains are expected to rapidly rearrange during a flow due to shearing, particle collisions, vibrations, and evolution of enduring grain contacts (Campbell, 2006;Furbish et al., 2008), thus causing rapid basal force fluctuations (Estep & Dufek, 2012). These fluctuations could also be a source of high-frequency seismic waves. The boundary between force-chain-dominated basal fluctuations and those caused by random individual particle fluctuations (and whether they are actually distinct) is not well understood but likely lies on a continuum (e.g., Campbell, 2006).
Basal force fluctuations that could generate seismic waves can also result from other aspects of unsteady flow such as interactions of clusters of particles with the bed, rapid larger-scale accelerations or decelerations in the flow such as a surge-front moving over larger-scale bumps and bends in the path (e.g., Iverson et al., 2016), and dynamic flow features such as standing waves and jets (Drake, 1990;Hsu et al., 2014;Iverson, 1997;Iverson et al., 2010;McCoy et al., 2013). The mean forces, whether time or spatially averaged, vary more slowly due to changes in the overlying thickness and density of the material as well as bulk flow accelerations and decelerations (Hsu et al., 2014;Iverson, 1997;Iverson et al., 2010;McCoy et al., 2013). Whether or not these slower changes generate observable seismic waves depends on (1) their timescale, which controls the frequency of resulting seismic waves, (2) the amplitude of the force changes and the resulting seismic waves relative to the ambient seismic noise of the area, and (3) the bandwidth of the seismic instrument being used. The most likely scenario is that the fluctuating forces, and thus the seismic methods, are fundamentally missing some information because they are bandlimited. The high-frequency seismic wavefield will have no direct sensitivity to the mean of the forcing except as reflected in the basal particle interaction variations. While the changes to mean forces may at times be too slow to themselves generate observable seismic waves (i.e., we are only directly sensitive to F f (t) in Equation (1), this does not mean that measurements of the fluctuating stresses have no relation to other useful flow properties. For example, changes to the mean forces or other bulk properties may result in changes in the statistical properties of the fluctuating forces that in turn may alter the observed high-frequency seismic wavefield over time.
In summary, a diversity of potential sources of basal force fluctuations could generate seismic waves at different frequencies at different points in time and space for a given debris flow. Additionally, parts of the forcing signal may not be directly reflected in the seismic signal. Therefore, in this study, we remain agnostic about which processes cause or dominate the fluctuations and focus on inferring the force fluctuations regardless of the exact source. We also investigate the statistical relationship between fluctuating stresses and other flow characteristics that may not directly be reflected in the seismic signal.

Seismic Path Effects
As signals travel from the source to the recording station, they are altered by the structure of the Earth by phenomena such as scattering, attenuation, and dispersion. We lump all these effects into the term "path effects." The Green's function represents the path effects. If it is known, then we can obtain information about the source function from the seismic data by removing (deconvolving) the Green's function and the seismometer instrument response from the seismic data. For simplicity, considering just one source location and one point force and assuming a linear system, the seismic signal may be expressed as where _ u p t; x ð Þis the ground velocity recorded at time t on the p component of a station located at x, which is equal to the sum over each direction j of the force F j (t 0 ), applied at location x 0 at time t 0 convolved with the Green's function _ G pj t; x; t 0 ; x 0 ð Þ . _ G pj relates a unit impulse force in direction j applied at x 0 at time t 0 to the observed seismogram at location x in the p direction at time t, and the superscript dot indicates we are considering ground velocity instead of the traditionally considered ground displacement (both sides of the original equation have been differentiated with regard to t). If we know _ G pj and _ u p , we can invert the equation to obtain the source, F j . Within this framework, the difficulty lies in obtaining sufficiently accurate estimates of _ G pj . Moreover, as alluded to earlier, in the case of a debris flow as a seismic source, F j is also not so simple, because debris flows are distributed sources that vary in space and time, not a single force at one location, as is assumed in Equation 2.
There are methods of modeling Green's functions if the velocity and attenuation structure of the subsurface are known or can be estimated (e.g., Herrmann, 2013;Shearer, 2009). Such an approach is feasible when considering low-frequency waves where small-scale details of the earth structure are unimportant, allowing simple 1-D models to be used. This approach is straightforward for the very low-frequency signals sometimes generated by large, rapid, landslides (Allstadt, 2013;Ekström & Stark, 2013;Moretti et al., 2012) but is not as feasible for debris-flow seismic signals because, with the exception of extremely large and rapid lahars, we observe primarily higher-frequency energy (from~5 Hz to hundreds of Hz) from debris flows (e.g., Allstadt et al., 2018, and references within). Path effects are especially potent for higher-frequency energy because (1) higher frequencies are more sensitive to small-scale heterogeneities and (2) amplitudes drop off with frequency at exponential rates (Shearer, 2009). Furthermore, seismic signals from surface sources such as debris flows are commonly dominated by surface waves (Sánchez-Sesma et al., 2011). This dominance further exacerbates the path effects problem because higher-frequency surface waves with their shorter wavelengths propagate in and are most sensitive to the shallowest and typically most attenuating and heterogeneous layers. Attenuation is especially difficult to measure and can be highly spatially variable and frequency dependent. Yet knowing the attenuation between source and site is critical if we want to use signal amplitudes to obtain information about source amplitudes with any degree of certainty. Therefore, estimates of _ G pj that are sufficiently accurate in the frequency range of interest are especially necessary for debris flow seismicity studies.
In this study, we directly measure _ G pj using controlled active sources, such that it constitutes an EGF. In the following sections, we describe our methods for generating the EGFs. We next construct a framework for using these EGFs for debris flows by parameterizing debris flows as a seismic source that allows us to obtain estimates of basal fluctuating forces. We then demonstrate and validate the method using data from large-scale flume experiments and discuss our results and their implications.

Experimental Setup
We conducted our study at the U.S. Geological Survey debris-flow flume. The flume is a 95-m-long, 2-mwide, 1.2-m-deep concrete channel that has been used for decades to study debris-flow dynamics (Iverson et al., 2010). The flume has a constant slope of 31°until 74 m downslope from the gate where it begins to curve and finally flatten to about 2°where it adjoins a large, nearly planar concrete runout pad ( Figure 1). Downslope distances along the flume are referenced to the position of the gate, where the gate is located at x = 0 m. Below the gate, the surface of the flume bed is bumpy concrete, as detailed by Iverson et al. (2010). Above the gate and on the runout pad, the surface is smooth brushed concrete.
In this study, we focus on two gate release experiments conducted in May 2017. Sediment was loaded into a hopper behind the 2-m-high headgate. Water was added using subsurface conduits and surface sprinklers until the sediment was nearly saturated, and then the gates were opened. Both experiments used debris that consisted of a 50/50 mix of wet sand and gravel with about 0.2 m 3 of cobbles added randomly throughout the mass ( Table 1). The 23 May experiment used 8 m 3 of debris compared to 10 m 3 for the 25 May experiment. The 10-m 3 sediment mass had a steeper surface slope (17°, compared to 9°), and a pond containing about 0.3 m 3 of water formed at the bluntly tapered upslope end of the debris. Iverson and George (2019) provide further details regarding these initial debris geometries, grain size distributions, and other debris physical properties.

Instrumentation
Instrumented flume cross sections are located in the hopper and at about x = 32, 66, 80, and 90 m (Figure 1), though only data from x = 32 and 66 m were used in this study. The instrumentation at these two cross sections consists of 1. An infrared laser distance meter mounted 1.2 m above the bed that measures flow thicknesses normal to the slope over an area of about 1 cm 2 in the center of the flume with~0.5-cm accuracy. Lasers are located above the normal force plates. 2. Force plates that measure basal normal and shear forces that are mounted flush to the flume bed and have a surface roughness matching that of their surroundings (bumps or smooth). Each plate has a circular load-contact surface of 500 cm 2 . At each cross section, the shear plate is 0.6 m downslope from the normal plate. A slight eastward tilt of the lower part of the flume means that the full flow thickness does not always pass under the laser or over the force plates at the x = 66-m cross section. 3. Three to four basal pore pressure sensors collocated with each normal force plate, mounted around its perimeter.
More instrumentation details can be found in Iverson et al. (2010) and Iverson and George (2019). All data from the instrumented cross sections were digitized by the standard flume data acquisition system (DAS), a DATAQ DI-722-EN, at a sample rate of 1000 Hz. To sync data recorded by the flume DAS precisely with the seismic data, we used the arrival time of the gate seismic signal at the seismic station closest to the gate.
To ensure we do not misinterpret sensor resonance as real signal when we compare seismically derived force fluctuations with direct force plate measurements, we investigated the frequency response of the force plates using impulsive ping tests and controlled-oscillation tests (supporting information section S1). We found that the force plate response is roughly flat and thus directly usable at lower frequencies (less thañ 50 Hz), but the response has strong resonances at higher frequencies that vary in detail among plates and under different loading conditions. This resonance results in "ringy" signals with unphysical fluctuations to negative values. We observed these effects in both individual ping tests and unfiltered data from our flume experiments ( Figures S2 and S3). Such characteristics have also been observed in other studies that use force plates (e.g., Hsu et al., 2014;McCoy et al., 2013). Therefore, we use only <50-Hz signals for comparison of force plate data with seismically derived fluctuating force estimates.
Photos and videos of each experiment were collected from multiple viewpoints. In this study, we use the data from video cameras mounted on a tower at the distal end of the runout pad, at x =~107 m, which had a view of the entire experiment. Videos were synced to the flume DAS using timing LED lights that illuminate instantaneously when the gate opens. Video recordings of these experiments (and many others) are available for viewing at https://pubs.usgs.gov/of/2007/1315/ (Logan et al., 2018).
Seismic data in our experiments were recorded by three-component Fairfield Zland nodes, which are small, portable, self-contained 5-Hz seismometers. We oriented the nodal seismometers to North and buried them directly in the soil to reduce noise (e.g., Farrell et al., 2018), covering them with just 2-4 in. of soil in order to maintain GPS connection for timing. The nodes were installed every~3 m along the length of the flume,~2-4 m from the flume wall ( Figure 1). Seismic data were digitized at 1000 Hz with a gain of 12 dB. We used a total station survey to obtain high-precision relative locations of seismic stations as well as other instrumentation locations and flume features because distances between features were as small as 2-3 m. Relative locations were corrected to UTM coordinates using block rotation around the UTM coordinates of the lower corner of the flume wall.

EGF Generation
In order to relate the observed seismic signals to the amplitude of forces applied at specific source locations, we generated EGFs by directly recording the response of the seismic network to point forces of known magnitude applied at 27 hammer locations, spaced approximately every 3 m along the flume (Figure 1b). To produce point forces, we used a 5.4-kg sledge force hammer (PCB Piezotronics Model 086D50) that outputs a digital time series of the forces measured by the hammer head that we digitized at 1000 Hz on a seismic digitizer (RefTek RT130) to ensure precise timing. More details on how we set up the hammer system and practical considerations for its use are found in section S2.
To protect the brittle flume tiles from damage, we struck a steel-reinforced concrete block that was crafted with a pattern on the bottom that mated with the bump pattern of the flume tiles to promote mechanical coupling. To generate EGFs for forces normal to the bed, we struck the center of the block with the hammer at each location ≥5 times and struck a side panel on the block ≥10 times for shear forces (due to the lower amplitudes of the shear forces). We conducted tests to confirm that the concrete block did not significantly affect energy transfer into the ground (section S2). Hammer impact locations were measured using the total station.
To cancel out noise and augment the signal during the processing stage, we used median stacking for each EGF direction and location after normalizing to a 1 N·s impulse by either deconvolving the forcing time series (simulating a Dirac delta impulse, which we refer to as true EGFs) or dividing by the impulse (pseudo-EGFs). The main difference is that the pseudo-EGFs are smoother ( Figure S7). True EGFs are more valid at higher frequencies because the hammer blow duration (~5 ms) is too long to be considered purely impulsive. We use the true EGFs in this study, but there is very little difference between the two for the frequencies we consider. See section S3 for full data processing details and EGF preparation methods. Examples of true EGFs for vertical components of Stations 631 and 670 are shown in Figure 2. The amplitudes are highest for the closest hammer locations to each station but remain elevated for hammer locations about 10 m in either direction. Also notable is that the waveforms vary quite significantly. For Station 631, the amplitudes are much smaller for locations downslope of x = 40 m, suggesting strong lateral variation in the subsurface.
We examine the frequency content of the EGFs between each hammer location and its closest seismic station (section S4) and find that while the spectral range varies from location to location, the central 80% of the energy is typically between 12 and 69 Hz ( Figure S8). Even though the stations are located just meters away from the hammer blow locations, there is very little energy above 70 Hz. This absence can be explained by attenuation and the characteristics of the hammer source. There is also very little energy below 12 Hz, possibly because the hammer source does not efficiently generate those frequencies. We avoid using EGFs at frequencies where they lack energy because when we use them in the inverse sense to estimate forces from a seismic signal, we divide by the EGF spectrum. In frequency ranges where amplitudes are low, we would divide by very small numbers and risk amplifying noise. Since we cannot use the force plate data above 50 Hz (due to sensor resonance) to validate whether the resulting amplified high-frequency fluctuating forces are real or simply amplified noise, we focus on the 15-to 50-Hz frequency band in this study but also briefly discuss results for higher frequencies.

Station Sensitivity
We used the EGFs to quantify the sensitivity of each seismic station to different parts of the flume bed, allowing us to understand what areas of the flume each station "sees" the best. To compute the sensitivity maps, we first removed the station response (using the nominal responses from the manufacturer) to correct all the EGFs to ground velocity, so the results would be directly comparable from station to station. Then, we computed the energy spectral density, which gives a relation between forcing and signal energy, of each EGF using the multitaper method (Thomson, 1982). We interpolated these spectra to an evenly spaced grid of frequency and distance along the flume using 2-D linear interpolation. The result is a map of the sensitivity of a given station to a 1 N·s Dirac delta impulse in space and frequency, shown for the station closest to the x = 32-m cross section in Figure 3. Since all EGFs are normalized to 1 N·s and corrected to physical units, we can directly compare the response of a given station to forcing at different locations and different frequencies, and we can also compare the relative sensitivity of different stations to forcing applied at different parts of the flume. Figure 3 shows that, as expected, sensitivity to lower frequencies is higher at greater distances along-flume than to higher frequencies. It is important to note that areas of lower sensitivity do not imply no sensitivity to forcing but merely proportionally less sensitivity.
We compute sensitivity maps for each seismic station and summarize the sensitivity in the frequency band of interest by integrating the sensitivity maps over 15-50 Hz and computing the ratio of the contribution from each individual location to the total energy summed across all locations. These sensitivity functions, which show the sensitivity of each station location x to each hammer location x 0 , are shown in Figure 4. To then quantify the span of primary sensitivity, A H , of each seismic station, we integrate the sensitivity functions outward from the peak contributing until we reach 75% of the total energy (gray lines). The majority of the energy within this frequency band comes from an~15-m section of flume closest to each station and is dominated by energy coming from the closest part of the flume, but there is variability.

Single Station Inversion
In this section, we describe conceptually how we parameterize the flow as a seismic source as well as the inversion procedure we used to obtain the fluctuating force time series from the seismic records and EGFs. The full mathematical details can be found in section S5. Note that we do not fully adopt any of the published debris flow seismicity models (Farin et al., 2019;Kean et al., 2015;Lai et al., 2018) because they focus on forcing from individual particle impacts and steady flow conditions. We also do not adapt models from studies of basal force fluctuation statistics (Ferguson et al., 2004;Gardel et al., 2009;Hsu et al., 2014;McCoy et al., 2013) because the results are largely empirical and specific to only the model, experiment, or natural flows being investigated. These past studies do give important context to our findings, however, so we revisit them in the discussion.
In this study, we use a single three-component seismic station to solve for the normal and shear fluctuating forces applied to a single subarea of the bed. More advanced approaches, such as using all seismic stations and all hammer locations to simultaneously solve for forcing over the entire flume, are possible with our experimental setup. However, such approaches require more constraints and complex inversion techniques and can only be used in places where dense station coverage is possible. In contrast, the single-station, single subarea approach is the most common configuration and assumption for along-channel monitoring. Therefore, we leave the multistation approach for future studies and here focus on the simplest approach.
To set up the problem, we start with one seismic station. That station is sensitive primarily to forcing applied over a subarea of the bed, A H (derived in section 3.2.1). We assume the EGF between that station and any hammer location, x 0 , that falls in A H is sufficiently representative of forcing anywhere in A H , though the hammer location closest to the station provides the best estimate. With this setup, the simplest application of the EGF approach is to solve for the equivalent fluctuating forces exerted at x 0 that would generate the observed seismic signal _ u p t; x ð Þ. The equivalent forcing represents the aggregate contribution of all the forcing occurring over many spatial scales at a given moment in time. A simplified way of thinking of this is as a system of particles in a volume above A H that are all moving around in different ways, resulting in the center of mass of that volume also moving around as the particles move relative to each other. The equivalent forcing felt by that subarea is the total mass of that volume times the acceleration of the center of mass of the system of particles.
We assume forces applied outside of A H have a small enough contribution that they can be neglected. Note that this is a significant simplification. It is most valid when the flow is passing over all parts of the flume with relatively equal energy, and it is not valid at all when the flow front has not yet reached the station's area of sensitivity because the station will nevertheless detect energy from the flow farther upslope, resulting in the estimated forcing at A H starting too early. Our simplification is also not valid if the source is not the flow at all but some other external seismic source such as an earthquake or human noise. We will revisit these and other limitations in detail in section 5.3.3.
We compute the fluctuating forces separately for each station using each hammer location that is within A H and take the median of all these solutions as our best estimate. This averaging provides a more stable The corresponding closest distance to the flume (in m) for each station is shown in parentheses. Thin gray lines span the extent that contributes 75% of the total energy for each station. The circles indicate the distance along the flume to the closest hammer location to each station. estimate than using just the closest hammer location, and the variability of the individual solutions also provides an estimate of uncertainty due to the assumption that spatially distributed forcing can be represented by forcing applied at one hammer location. This division of the flow into volumes of particles over bed area A H can be repeated independently for all seismic stations to observe space and time variations in fluctuating forcing (with overlaps in sensitivity; see Figure 4).
To implement this inversion mathematically, we replace F j (t 0 ) in Equation (2 with F fS (t 0 ) and F fN (t 0 ), where the subscript f indicates fluctuating forces and where j = 1 is parallel to the slope (shear), so we replace it with subscript S for simplicity, and j = 2 is normal, so we replace it with N, and we neglect forcing in the cross-channel direction (j = 3). These changes convert Equation (2 (for the vertical component) to We can rewrite Equation (3 in the frequency domain to remove convolution and include all station components as Þ . This is essentially a set of three linear equations for two unknowns (for each frequency and location) that we can solve using least squares to obtain the best fitting values of In practice, we solve for all frequencies at once. See section S4 for full details of the inversion process.
As discussed earlier, the fluctuating forces represent a stochastic process for which the small-scale statistical properties relate to macroscopic flow properties that are not directly generating high-frequency seismic waves. Therefore, the fluctuating forces are best considered in time windows with a length long enough to provide a statistically representative sample of variance but short enough so that the macroscopic flow properties do not change significantly within the time considered. Autocorrelation of the force plate data ( Figure S10) shows that the force amplitudes decorrelate slowly with time. This observation implies that there is no clear best choice for window length, but that any choice requires balancing a trade-off between temporal resolution and statistical stability. We choose sliding time windows of 0.5 s overlapping by 80% and explore the effect of different window lengths in section 5.3.1.
The result of the procedure described above is the frequency domain representation of fluctuating forces in both shear F f S f ; t win ð Þ and normal F f N f ; t win ð Þ directions at each time window, t win . This can be plotted directly as a 2-D image (i.e., a spectrogram of fluctuating forces) to visualize frequency changes of the fluctuating forces with time. However, because the fluctuating force spectrum of each time window should be treated as a statistical representation of the forcing over a given time period, we apply further processing to compute the root-mean-squared (RMS) fluctuating force F rms,j for each t win . We compute this in the time domain from the spectrum directly for each t win and each direction j by where W is the total number of samples in t win and M is the number of samples in the FFT. Since RMS is equal to the standard deviation for a signal with zero mean (like fluctuating forces), this equation represents the frequency domain equation for the standard deviation of the signal with normalization factors for the discrete FFT (see section S4 for full derivation). Note that because we are limiting ourselves to the 15-to 50-Hz band in this study, when applying Equation 5, we only sum the components of F fj that lie within that frequency range (including the corresponding negative frequencies). Thus, what we obtain is the RMS of only the 15-to 50-Hz fluctuating forces (which is approximately equal to computing the RMS of the band-pass-filtered fluctuating force signal).

Journal of Geophysical Research: Earth Surface
To validate this seismically derived F rms,j against the force plate data directly, we need to apply further processing. Recall that F rms,j is the equivalent forcing, which is the point source that, when applied at the corresponding hammer location, would return the observed seismic signal representative of the integral of forcing over A H . In reality, however, the basal forcing consists of a sum of spatially distributed shear and normal forces occurring over a range of areas and spatial scales that are not well correlated and that contribute in different amounts depending on how sensitive the station is to a given location. The forces between two close locations may be similar statistically. That is, they may have similar spectral amplitudes, but they are likely poorly correlated in time (random in phase). Since we do not have information about the true spatial correlation of the fluctuating forces and we want to validate against force plate data, then we actually want to convert our result to forcing at the scale of the force plate, which is what the force plates measure. To do this, we make the assumption that forcing is uncorrelated at scales larger than the force plates. The sum of k identical individual forcings that occur randomly in time scales in the frequency domain by the square root of the total number of forcing locations (e.g., Tsai et al., 2012). Thus, we can use this approach to scale F fj f ð Þ down to the fluctuating forces experienced by a small area of bed with the same areas as the force plate F fpj f ð Þ by and fluctuating boundary stress in the j direction, as measured at the scale of A plate by where A plate is the area of the force plate. We use Equation (7 to compute σ fj from F fj and modify Equation (5 to compute σ rms,j for each sliding time window t win separately for each seismic station. We show σ rms,j in most of the figures presented in subsequent sections of this paper.

Processing Nonseismic Data for Comparison
We process data from force plates, lasers, and pore pressure sensors at the instrumented cross sections to enable direct comparison to the seismic results by computing the mean time-averaged basal stresses, flow thicknesses, and pore pressures over the same sliding windows (0.5 s with an overlap of 80%) as described in section 3.3. As with the seismic data, we associate these sensor data to the midpoint time of each window. We then use these directly measured parameters to compute effective stress σ e t ð Þat time t by σ e t ð Þ ¼ σ N t ð Þ − p t ð Þ , where σ N t ð Þis the mean time-averaged basal normal stress and p(t) is the basal pore pressure (median of the three to four sensors located at each force plate). We also estimate time-dependent debris bulk density ρ(t) following Iverson et al. (2010) by where h(t) is the flow thickness normal to the bed, g is Earth's gravitational acceleration, and θ is the slope of the bed (31°). Equation 8 assumes that σ N t ð Þ balances the slope-normal component of the weight of the debris. For comparison with the seismically derived results, we compute the RMS fluctuating stresses based on observed force plate time series data in the 15-to 50-Hz frequency band using Equation 5. We also compute the ratio of the fluctuating stresses and the mean stresses over time.
To obtain the location of the flow front and major surge fronts with time and to estimate flow velocities, we use the methods of Rengers et al. (2019) to create a 2-D space-time snapshot of the flow from the experiment video recordings ( Figure 5). We extract the color of the centerline of the flume at each moment in time in the video, correct it for perspective using marked intervals on the flume bed, and repeat over each frame of the video. We stack these time slices in a figure such that the x axis represents distance and the y axis represents time so that we can see a single snapshot of changes in color over space and time during each experiment representing details of the debris-flow structure. The space-time snapshots clearly depict the contact between the flow front and the bare flume due to the strong color contrasts, and the saltating front is clearly differentiated by its more granular character. We extract along-flume flow front velocities with time by taking the inverse of slopes of the flow front and saltating front locations over time. Additionally, we estimate the instantaneous surface velocity at the instrumented cross sections by taking the slope of linear features crossing the location from the space-time flow snapshot ( Figure 5) and interpolating and smoothing between measurements (see section S6 for details). These linear features represent surge fronts and roll waves. Following Iverson et al. (2010), we use the term roll wave to describe the sequences of trailing smaller secondary surges resulting from flow instabilities.
We also investigate two additional quantities that we compute from combinations of the aforementioned measurements: bulk flow kinetic energy per unit bed area (KE A ) and downslope linear momentum (mom A ) per unit bed area, computed at the instrumented cross sections by where ρ is flow bulk density, h is flow thickness normal to the bed, and v is the instantaneous, depth-averaged downslope velocity. We have only surface velocity estimates, so we assume plug flow with a uniform flow velocity profile for these computations.

Experiment Observations
We first summarize overall visual and seismic observations of the entirety of each experiment and then focus in on the instrumented cross-section data, primarily at x = 32 m, to gain a more detailed picture of each flow. From the space-time flow snapshots (Figures 5a and 5b), extracted velocities (Figures 5c and 5d), videos (Logan et al., 2018), and seismic data (Figure 6), we found that the two experimental debris flows behaved similarly at first but quickly diverged in their behavior, a consequence of their different starting conditions (Table 1). In each case, as the gates unlatched at t = 0 s and the flow pushed them open against the flume sidewalls (~1 s later), the two largest seismic signals were generated. The first was due to the force drop associated with the flume gate latch opening, and the latter was due to the gates hitting the padded flume sidewalls. The amplitudes of these signals far exceeded the amplitudes of the signals from the debris flows themselves, especially at stations close to x = 0 m. As a result, we exclude the first 1.8 s of data from our analysis.

Journal of Geophysical Research: Earth Surface
Upon leaving the hopper, the debris flows in both experiments reached peak frontal velocities at t =~2.5 s. The flow fronts began to decelerate and stall around x = 60 m, but the stalled fronts were eventually overtaken and merged with subsequent surges (white arrows, Figure 5). However, the first surge to overtake the front in the 23 May experiment did so more than a second after the front stalled, and then the overtaking surge also stalled and was overtaken by yet another major surge. In comparison, the first surge in the 25 May experiment immediately overtook the stalled flow front and continued down to the bottom of the flume without slowing. The behavior resulted from the additional debris and the extra pulse of trailing water from the aforementioned source area pond on 25 May, as occurred in past similar flume experiments (Iverson et al., 2010;Logan et al., 2018). Though both debris flows reached similar peak velocities, the 25 May flow maintained higher overall velocities and reached the bottom of the flume much sooner than the 23 May flow did, and it was therefore briefer in duration (Figures 5c and 5d). This difference is reflected in the seismic signal, especially beyond x =~60 m: The 23 May signals are more diffuse and generally have a longer duration, whereas the 25 May signals have higher and more distinct peaks and shorter overall durations (Figures 6a and 6b).
All the debris-flow seismic signals have emergent waveforms that build up for about 0.5-1 s prior to the flow front passing each station. This pattern indicates that a given seismic station is dramatically more sensitive to the flow as it passes nearby than to more distant parts of the flow (Figures 6a and 6b), confirming what we inferred from the sensitivity maps ( Figure 3). The signals tend to peak~1-2 s after the flow front first passes closest to the station, suggesting that the peak amplitudes occur when the entire area of sensitivity is occupied by the passing flow rather than when the leading edge of the flow front is closest to the station. The power spectral densities (Figures 6c and 6d) show that even though stations are mere meters away from the flow, the seismic energy is still concentrated at frequencies <100 Hz, with peaks for most stations occuring at~20-30 Hz. In each experiment, a cloud of saltating gravel and rolling cobbles outran the main debris-flow front. The cloud emerged at~x = 45 m and continued down the flume without stalling. Once the saltating cloud broke away from the debris-flow front, the seismic signals for nearby stations became spikier and contained more higher-frequency energy as the saltating grains proceeded downslope ( Figure 6). The seismic manifestation of the saltating cloud is more apparent for the 23 May experiment (between dashed and dotted line >x = 60 m; Figure 6a) because due to their differing flow front velocities, the saltating cloud reached the bottom of the flume~7 s prior to the main debris-flow front on 23 May but only~1.5 s prior to the flow front in the 25 May experiment. The biggest seismic energy spikes, such as those at Stations 670 and 901 during the 23 May experiment, are caused by scattered cobbles flying out of the flume and either directly hitting or landing very close to the seismometers (as confirmed by video recordings) and contaminating the signal. These impacts must be excluded from the analysis, so they do not get misinterpreted as high amplitude debris-flow signals.

Instrumented Cross-Section Observations
The data from the force plates, lasers, and pore pressure sensors at the instrumented cross sections give us more detailed insight into flow dynamics but only at fixed locations. Since these types of data have been thoroughly explored in other studies at the debris-flow flume (e.g., Iverson et al., 2010), we avoid a detailed analysis and focus on the aspects most relevant to the fluctuating basal forces and primarily on the data from the x = 32-m cross section (Figure 7). Results for the x = 66-m cross section are shown in Figure S13, but these data are of limited use for this investigation for reasons described below.
The 15-to 50-Hz basal stress fluctuations measured by the force plates at x = 32 m show a peak when the main flow front first arrives, which is also when flow bulk density is lowest (~500-1,000 kg/m 3 ) and prior to when the pore pressures start to rise dramatically. These low bulk densities reflect the agitated, dilated, and unsaturated state of the flow front. Basal pore pressures do not rise appreciably until after the flow front passes. Arrival and passage of the front also coincides with the time of the peak downslope flow velocities and with the highest seismic amplitudes. The formation of the flow front is due to grain size segregation (e.g., Iverson et al., 2010;Johnson et al., 2012), and therefore, the characteristic grain size within the front is coarser than within the subsequent flow. When the main flow front first arrives, the ratio of the fluctuating stresses to the mean stresses is~1, but this ratio dramatically decreases nonlinearly as the rest of the flow passes. The ratio eventually levels out around 0.01 for normal stresses and 0.05 for shear stresses (Figures 7e and 7f) when the tail of the debris flow is dominated by a series of trailing, water-saturated roll waves (t > 10 s). The 25 May experiment has more numerous, smaller amplitude roll waves at x = 32 m than the 23 May experiment, which is instead dominated by fewer larger surges. This distinction is also clear in Figures 5a and 5b. The small roll waves during the 25 May experiment are barely perceptible in the force plate fluctuating stresses.
In between the main front and the smaller, trailing roll waves, larger secondary surges pass the x = 32-m cross section (dotted lines, Figure 7). Relative to the main flow front, these are accompanied by smaller peaks in flow thickness and fluctuating stresses and also higher pore pressures, with no obvious decrease in flow bulk density. During passage of these surges, the ratio of fluctuating stress to mean stress does not noticeably rise from its descent except slightly for one surge on 23 May. This behavior suggests that these secondary surges are wetter and less agitated and/or consist of smaller particles than the main flow front. Thus, they are also likely less efficient at generating seismic waves. However, even though the passage of these secondary surges is not obvious in the raw seismic data (Figure 6), the surges may be more seismically important than the smaller, trailing roll waves because they do correspond with sharp rises in force plate fluctuating stresses.
The x = 66-m data ( Figure S13) show that a prominent saltating flow front has clearly formed and outrun the main debris-flow front at these locations. As a result, the structure of the fluctuating stress time series is less clear and does not have such a distinctive peak as when the main flow front arrives as at x = 32 m. As is apparent in the video recordings (Logan et al., 2018), the saltating gravel and rolling cobbles that precede the debris-flow front only sporadically hit the force plate or pass beneath the laser at x = 66 m, yielding an incomplete record. Furthermore, the slight eastward tilt of the lower half of the flume directs the flow away from the force plates, such that the leading edge of the main flow does not pass over the force plate. Thus, we cannot rely on the data from x = 66 m until the main flow front fills the channel (solid green line in Figure S13). After this point, patterns similar to those we described for x = 32 m are evident.

Inversion for Basal Fluctuating Forces
Here, we apply the methods outlined in section 3.3 to solve for the best fitting fluctuating forces and stresses at each seismic station using the closest hammer location (see Figure 1 for map) and compare against the force plate-derived measurements. Note some stations are excluded from the fluctuating force results. Stations close to the gate area are excluded because they are dominated by gate noise, and we replaced Station 901 with nearby 402 for the 23 May experiment because 901 was hit by a cobble. Results are presented in Figures 8-10 and discussed below, starting with the instrumented cross sections for validation and then to all locations to examine the spatiotemporal evolution of fluctuating stresses.
In the spectrograms of the fluctuating normal force power ( F fN f ð Þ 2 =t win Þ at x = 32 m and x = 66 m (Figure 8), after the gate signal at t = 0 s, we see that the debris flow signal is broad in frequency, with some bands of more concentrated energy at specific frequency ranges. The relative proportions of energy at each frequency do not change substantially during the main progression of the flow, but amplitudes peak when the flow front and surges pass close to the station (white vertical lines). There is substantial energy at >50 Hz, suggesting that we likely can obtain useful information at even higher frequencies than we use in this study, but we cannot verify against force plate data. There is also energy below 15 Hz, even prior to the gate opening. Recall that we hypothesized that the amplification of noise could occur at frequencies where EGF

10.1029/2020JF005590
Journal of Geophysical Research: Earth Surface energy is low, and that is likely the case here. The EGFs shown in Figure 8 (thick white lines along y axis) have very low amplitudes below 15 Hz. There are also some bands of apparent amplification, especially at higher frequencies than shown in Figure 8 (see Figures S14 and S15 for the full spectrograms up to 500 Hz). The bands of amplification may be related to persistent noise signals that have been observed in data recorded by the seismometer model used in this study (Farrell et al., 2018). Farrell et al. (2018) found that the frequency of this noise varies by instrument and with temperature. The EGFs may not properly account for this instrument noise during the experiments if the frequency of the noise shifted due to differing temperatures.
We next compare the seismically derived RMS fluctuating stresses (σ rms,j ) against the force plate data (Figure 9). We find that the seismic solution using the closest hammer location to the station and the median of all solutions for hammer locations within the area of sensitivity are very closely matched. They both compare well in amplitude and timing to the fluctuating stresses from the force plates in the 15-to 50-Hz band at x = 32 m. The main difference is that the seismically derived signal starts earlier than the force plate signal. This time difference occurs because the station is sensitive to areas upstream of the force plate (Figure 3), so the seismometer "sees" the flow earlier and the assumptions established in 3.3 are less valid when the flow is not entirely covering A H . The seismically derived and force plate-derived fluctuating forces also compare well at x = 66 m but only after the main flow front arrives ( Figure S16). They do not compare well during the passage of the saltating cloud, as expected, because the force plates do not measure the behavior of this part of the flow well. The span of the individual solutions (gray lines, Figure 9) for each hammer location that is within A H gives an estimate of the uncertainty due to the hammer location selection. For Station 631 (the closest station to the x = 32-m force plate), the range is small, with individual amplitudes being, on average, a factor of 0.8 to 1.3 times the median solution. The mean for all stations is a factor of 0.7 to 1.6 times the median ( Figure S11). The median solution is more stable because it is not as prone to errors in any individual EGF.
The seismic method does appear to detect the arrival of the main secondary surge of the 25 May experiment at the x = 32-m cross section at around~7 s, even though the surge is just starting to form at that point (Figures 9d and 9e, dotted line). By the time that secondary surge reaches the x = 66-m cross section, it has essentially merged with the initial flow front, so it is not seen distinctly there in results for either experiment ( Figure S16). The seismic method does not reproduce the isolated spikes in the fluctuating stresses later in the flow that occur when cobbles roll over the force plates or roll waves pass (especially in Figures 9a and  9b). This insensitivity is expected because force plates effectively provide point measurements and the seismic method provides a distributed measurement. Roll waves are numerous and closely spaced, so it is unlikely they would be distinguishable from each other in the seismic data. Also, as noted in section 4.2, the roll waves do not cause significant increases in fluctuating stresses.
The seismically derived fluctuating stress estimates demonstrate how the fluctuating stresses change in time and space due to the density of our station spacing ( Figure 10). The figure shows the median of the fluctuating stresses computed using all hammer locations for each seismic station. Note that the close (~3 m) spacing of these station-hammer pairs results in spatial overlap in the area of sensitivity of the seismic stations used for each hammer location. Therefore, the same fluctuating energy is mapped to more than one location at a time. However, the result is informative because it shows that the timing of the peak RMS fluctuating stresses roughly corresponds to (but lags slightly) the flow front arrival until about x = 50 m for both experiments. After this point, the fluctuating forces instead track the saltating front, but the shape of the time series broadens as the saltating and main fronts separate in time. Interestingly, the arrival of the main front in the 23 May experiment after x = 60 m is not seismically perceptible at this frequency range. It is masked by the signal generated by the saltating material arriving in front of it. The lack of signal also likely reflects the low speed of the main front by this time. The first secondary surge is not clearly distinguishable, but it appears as a broadening and, in some cases, a weak second peak in the stresses on stations between x = 31 m and x = 66 m for both experiments. The arrival of the next secondary surge for the 23 May experiment is not noticeable in the seismic results except possibly close to where it overtakes the main front.

Correlations Between Fluctuating Stresses and Flow Characteristics
In this section, we investigate the relation between the fluctuating stresses and other flow characteristics. We use the fluctuating stresses measured directly by the force plates at the x = 32-m cross section rather than the seismically derived result to ensure precise timing of the data as the flow front arrives at the location of the other cross-section instruments. We plot these data as scatterplots with fluctuating stresses on the x axis and  (c) and (f). Prior to taking the ratios of the force plate-derived data, the relative timing is approximately corrected for the spatial offset between the normal and shear plates (0.6 m) using the instantaneous velocities extracted from Figure 5 at the x = 32-m cross section. Dashed horizontal lines indicate the median ratios for the seismic (black) and force plate (blue) fluctuating stresses.

10.1029/2020JF005590
Journal of Geophysical Research: Earth Surface each different flow characteristic on the y axis, color the dots by time, and look for systematic patterns such as linear relations or changes in slope that correspond to changes in flow characteristics ( Figure 11).
We find systematic relationships between the fluctuating stresses and all flow characteristics investigated, but the relationships are neither linear nor well described by any other simple functional form, especially if the entire flow is considered. The mean stresses, flow thickness, and surface velocity all have positive correlations with fluctuating forces. In contrast, bulk density has a high negative correlation, indicating a strong inverse relationship (Table 2), except during the passage of the main flow front, when it has a positive or no correlation. Mean normal stress has the lowest correlation coefficient, just 0.4, compared to 0.8 for mean shear stress, which has the highest correlation. Flow thickness, velocity, and bulk density all also have high correlations, exceeding 0.8, though the correlation is negative for bulk density.
The coarse, unsaturated flow front has a distinct relationship with the investigated flow characteristics compared to after pore pressures rise. The scatterplot data form outlier loops and arcs on most panels of Figure 11 during the time period when the coarse flow front arrives (e.g., see panel h for an example "loop" and panel w for an example "arc"). During this time, the correlation between fluctuating stress and the flow parameter is sometimes in the opposite sense to the rest of the time period. The pattern these initial outliers follow is similar between the two experiments (see Figure 5) and is similar for mean and shear fluctuating stresses except for a slightly earlier peak for the shear fluctuating stresses (despite the shear force plate being downhill from the normal force plate). Correlations between the stress fluctuations and mean normal stresses and effective stresses are much higher if we exclude data from before the pore pressures start to rise significantly (this transition point, shown by the cyan circles in Figure 11, also coincides with the timing of peak  (Table 2). There is another inflection point in the relation with mean normal stress when the peak pore pressures are reached (green circle on Figure 11). All of the relations are most straightforward after the rate of change of the bulk density has slowed upon reaching~1,800 kg/m 3 (gray circle). The fact that fluctuating stress is Figure 11. Scatterplots showing the relationship between normal and shear fluctuating stresses (from force plates) and measured and computed flow properties for each of the two experiments (labeled at top) measured at the x = 32-m instrumented cross section. Subplots compare the RMS of fluctuating normal and shear stress (σ rms,N and σ rms,S respectively) on the y axis against mean normal stress (a-d), mean shear stress (e-h), mean effective stress (i-l), flow thickness (m-p), flow velocity (q-t), bulk density (u-x), kinetic energy per unit area (y, z, 0, 1), and momentum per unit area (2-5). Colors indicate time from gate opening. Cyan outlines represent the time when pore pressures start to rise significantly, green outlines indicate the approximate time of the peak pore pressures, and gray outlines indicate when the bulk density reaches 1,800 kg/m 3 . Solid lines indicate power law fits and dotted lines indicate linear fits from Table 2. The fit lines shown are the same across rows (same for all components and all experiments for a given flow property); apparent differences are related to axis extent differences.

10.1029/2020JF005590
Journal of Geophysical Research: Earth Surface less correlated with mean normal stress than flow thickness, and bulk density has a negative correlation, suggests that a more dilated, and thus agitated, flow generates higher fluctuating forces for a given mass.
The correlation coefficient quantifies the linear correlation. However, most of the curves in Figure 11 are nonlinear, especially after the initial arrival of the coarse flow front. We test whether a nonlinear relation better explains the trends by comparing a linear fit to a power law fit to the relation between each flow parameter and the combined shear and normal fluctuating stresses for both experiments (Table 2). We use nonlinear least squares and exclude the coarse flow front. We find that a power law fit is significantly better than a linear fit for mean normal stresses (r 2 = 0.71 vs. r 2 = 0.55), slightly worse than linear for velocity and bulk density, and about the same for the other factors. For mean shear stresses and flow thickness, the two parameters best fit by a power law (r 2 > 0.8), the best fit powers are 0.6 and 0.5, respectively. This suggests that the fluctuating stresses might be roughly proportional to the square of the mean shear stresses and flow thickness, at least for the part of the flow with elevated pore pressures. However, linear fits are nearly as good for both factors.

Ratio of Shear to Normal Basal Stress Fluctuations
Interestingly, the seismically derived shear stress fluctuations are higher in amplitude than the normal stress fluctuations by a factor of~2 (Figures 9c and 9f). In fact, this ratio is nearly constant for all hammer/station pairs, which have a mean ratio of 2.1 (std 0.5), despite their independence. A constant ratio is not unreasonable, given that the slope angle and bump pattern is constant. The force plates are offset by 0.6 m, so to obtain a rough estimate of the measured ratio for comparison, we correct the timing of the shear force plate data using the instantaneous velocity estimates presented in Figure 5 for the x = 32-m plates. We find that the ratio varies quite a bit more than the seismic result, possibly due in part to an imperfect correction for the spatial offset, but the mean (and median) ratio of shear to normal stress fluctuations is also~2 (Figure 9). However, in contrast to the near-constant ratio of the seismically derived result, the force plate-derived ratios generally transition from slightly lower ratios during the passage of the main flow front to slightly higher as later parts of the flow pass.
We can infer that in general the ratio of shear to normal stress fluctuations relates to the mean angle of fluctuating forcing; a higher ratio implies shallower mean forcing angles and a lower ratio implies steeper forcing angles (closer to normal). The evolution from lower ratios (mean 1.2-1.4) to higher ratios (~2) during the passage of the main flow front over the force plate suggests the main flow front has a higher angle of forcing than the tail. Though neither data source is without fault, if this evolution is real, it is not detectable by our seismic method. One possible explanation is that because the seismic method is sensitive to a larger area of forcing than the force plates, spatial averaging may reduce its sensitivity to localized changes in the stress ratio.

Relation of Fluctuating Forces to Flow Characteristics
We have demonstrated that we can estimate spatiotemporal variations in debris flows' basal fluctuating forces by using the observed along-channel seismic signals. That in itself is a step forward, but it leaves  open the question of what we can then learn from the fluctuating forces. Basal fluctuating forces themselves are sometimes of interest, for example, for understanding debris flow erosion (e.g., Hsu et al., 2014;McArdell et al., 2007;McCoy et al., 2013). Unfortunately, few physics-based studies relate fluctuating forces for dense, realistic, debris flows to other flow characteristics, such as those more closely related to hazards (e.g., flow speed and depth). In this section, we compare our findings from section 4.4 to the handful of empirical studies that have given us insight into the problem. We do so while keeping in mind that comparing to one flow parameter at a time may be a flawed approach given that many properties are also correlated with each other and multiple properties may be changing simultaneously. This is true for our results and also for most of the studies we mention below, so we focus on general trends.
The relationships shown in Figure 11 and quantified in Table 2 are generally consistent with observational and experimental findings from past studies. Numerous observational case studies have shown that within narrow size limits and under similar flow conditions, the amplitude of the seismic signal relates closely (but not necessarily linearly) with flow depth and/or discharge (Doyle et al., 2010;Lavigne et al., 2000;Marchi et al., 2002;Marcial et al., 1996;Tungol & Regalado, 1996). By proxy, the amplitude of the seismic signal also correlates with the mean normal and shear stresses because they largely reflect the weight of the flow if significant local accelerations are absent. Given that these past studies all focused on high-frequency alongchannel seismic recordings which are generated by fluctuating basal forces, we would expect the fluctuating basal forces to also relate to flow depth or discharge. McCoy et al. (2013) and Hsu et al. (2014) directly measured the basal fluctuating stresses of small debris flows using force plates and found a positive correlation between mean normal stress and fluctuating stress in force plate measurements within the limited flow size ranges of their investigations. Since mean normal stress is directly proportional to flow thickness for uniform, steady flows, their results are consistent with the findings of the aforementioned observational seismic studies.
The negative correlation with bulk density we observe is consistent with the past observation that denser, plug-like flows were seismically quieter than less dense flows (Cole et al., 2009) and that fluctuating forces are higher at coarser, erosive flow fronts than wetter tails (Berger et al., 2011, and this study; Figure 7). Unsaturated flow fronts are generally more dilated and have lower bulk density than trailing parts of debris flows (Iverson et al., 2010). The observation that the unsaturated, coarse, initial flow front behaves differently than the rest of the flow in relation to most of the flow properties investigated is not only scientifically interesting but also noteworthy because the initial front is thought to be one of the strongest generators of seismic energy (e.g., Arattano & Moia, 1999;Farin et al., 2019;Lai et al., 2018) and of the most interest for monitoring.
The reasonably good linear correlation of fluctuating stresses with flow thickness (0.8) and the inverse relationship with bulk density (−0.8) can be understood by considering granular temperature, T. Granular temperature is a measure of the degree of agitation (i.e., grain-scale, randomly directed kinetic energy) of a granular flow. It is equal to twice the fluctuating kinetic energy per unit mass of solid grains (Iverson, 1997;Ogawa, 1978), which is the sum of the squared fluctuating grain velocities normalized by mass. Fluctuating velocity is the velocity of the grain minus the mean translational flow velocity at the grain's location. A higher T means higher average fluctuating velocities, which also means higher fluctuating forces because impact force depends on velocity at impact, as well as the mass of the particle (e.g., Farin et al., 2015;Hertz, 1882). Therefore, faster particle fluctuations and larger average particle diameters would result in higher fluctuating forces. Larger particles are more common in larger, faster flows, and higher particle fluctuations are also more likely in larger flows: The higher the translational flow velocity, the more translational momentum is available to be converted into fluctuation momentum (Iverson, 1997). Granular temperature is difficult to measure directly (we cannot compute it without keeping track of motions of every grain). However, for steady, uniform flows of identical, nearly elastic spheres with a high solid volume fraction, according to Jenkins and Askari (1999), granular temperature can be approximated as . The conditions required for Equation (11 are not valid for natural debris flows, but the general concept may explain the strong (but imperfect) relation between flow thickness (or by proxy, mean stresses) and fluctuating basal stresses observed in this study (and observed indirectly via the seismic wavefield in many other studies).
Despite the expected dominance of velocity on the fluctuating stresses (Farin et al., 2019;Hsu et al., 2014;Lai et al., 2018), the fit of a linear or power law model to the relationship between the two is not very high (r 2 = 0.63 and 0.58) compared to other factors. This difference may arise because flow is not in steady state as the aforementioned studies require. Moreover, our measurements of velocity are indirect and uncertain, plus our study does not control for grain diameter, which is estimated to exert just as strong an influence as flow velocity in the theoretical models (in which both particle diameter and velocity scale to the third power in their influence on seismic amplitudes; Farin et al., 2019;Lai et al., 2018). Though we know the grain-size distributions of the starting conditions with great confidence, we do not have control or the ability to measure in detail how the grains segregated throughout the flow. It is also possible that the effect of flow velocity may not be as dominant as the theoretical models and past studies suggest.
With our experimental setup, we are able to estimate bulk properties that depend on velocity as well as flow mass. While the kinetic energy and momentum per unit area are both linearly correlated (~0.8) with fluctuating stresses if the flow front is excluded, they are about as equally well fit by a power law (Table 2). This finding is consistent with the findings of Coviello et al. (2019) who noted a roughly linear trend between seismic amplitude (which is linearly proportional to force) and kinetic energy per unit area for natural flows. However, because mean flow velocities do not vary dramatically in our experiments, the fit to momentum per unit area is just as plausibly linear as the fit to kinetic energy for the time period after peak pore pressures occur. A linear fit with momentum is in line with findings of Hibert et al. (2017), who found a roughly linear relation between high-frequency amplitudes and momentum of large landslides. Therefore, with the currently available data, we cannot say that one variable is more closely related to fluctuating forces than the other; the relationship between both factors appears to be complicated and indirect. One possibility is that the flow profile changes substantially over the course of the flow, but we have only surface velocities and assume plug flow for all time periods. Another possible explanation is that the amount of bulk translational kinetic energy that is converted into vibrational kinetic energy actually depends on the characteristics of the flow and also of the bed.

Practical Considerations 5.3.1. Time and Frequency Choices
In section 3.3, we chose to use a frequency band of 15-50 Hz. However, the energy we observe in the solution above 50 Hz (Figure 8) may still provide useful information about the flow, we are just unable to validate it against the force plate data because of force plate resonances at higher frequencies (section S1). In spite of this, to explore what we might be excluding through our frequency choice, we estimate the RMS fluctuating forces via seismic inversion at x = 32 m for the 25 May experiment over a range of frequency bands and compare to the force plate RMS for the same frequency bands (Figure 12). We found that the frequency band does matter. At the highest frequency ranges (>100 Hz), the result is dominated by very sharp peaks in both the force plate result and the seismic result, but the frequency of the sharp peaks is not the same. Due to the short duration and high amplitudes, these peaks probably reflect cobble impacts very close to the seismic station and on the force plate at different times. Such impacts are amplified in the seismic result (especially >200 Hz) due to the low EGF amplitudes at this frequency. At the other extreme, the force plate has significant fluctuating force energy in the 5to 15-Hz band, whereas there is very little energy in the seismic result in that range. This finding reflects the bandlimitedness of our seismological technique. To further investigate, we split the frequency range we used in this study of 15-50 Hz into two bins, 15-30 and 30-50 Hz, and find that the frequencies at which the fluctuating forces are concentrated are slightly offset between the force plate and the seismic result when the flow front first arrives. The seismic result has higher amplitudes than the force plate at 30-50 Hz earlier in the flow and vice versa for 15-30 Hz. Some of this mismatch may be due to energy at other frequencies leaking into the range of interest (especially from lower frequencies because the force plate data have substantial energy at lower frequencies), but much of it is likely because the seismic method is sensitive to a larger area of forcing than the force plate and the assumptions of the method are not well met when the flow first approaches the station.
In section 3.3, we also chose a somewhat arbitrary sliding time window of 0.5 s for all our processing. We investigate the influence of this choice by redoing the inversion for the 25 May experiment at x = 32 m with different window lengths. We found that the selection of the time window length is much less important than frequency range selection (Figure 12). At most window lengths, the force plate and seismic results are similar in both shape and amplitude; the results just progressively become smoother. The exception is the 0.1-s time window, which gives a seismic result that fluctuates much more and to higher amplitudes than the force plate result, suggesting that a 0.1-s window length is too short to statistically represent the fluctuating forces. Note that the force plate data are likely affected by resonance at frequencies greater than 50 Hz and do not serve as a validation at those frequency ranges.

Station Spacing
One question that may arise in future applications of this technique is "what station spacing is ideal to capture the behavior of an entire debris flow?" In this study, we used a very dense spacing of 3 m. We found that there was significant overlap in the sensitivity of adjacent stations (Figure 4) once we did our hammer tests, and we probably could have obtained similar sensitivity to the flume with~15-m spacing if we had done this testing ahead of time. In future studies in other locations, we suggest that the station spacing be optimized to minimize overlaps by doing hammer tests ahead of time using a temporary deployment of geophones (or by moving geophones and repeating hammer tests along the channel) and computing the sensitivity maps before a permanent installation to ensure the final station design has sensitivity to the entire area of interest without significant overlaps.

Limitations of the Single-Station Approach
The single station, single subarea approach is useful in settings where we want to measure debris-flow behavior at that specific subarea location, although we also need to be tolerant of some spillover of seismic energy from adjacent areas of the bed. The downside is that all energy recorded by a given seismic station gets mapped to the subarea even if it is coming from somewhere else. It is almost always the case that the seismic signal will simultaneously record contributions from many different areas when we have a spatially distributed source and the source to station distance is less than the length of the source area. This effect was obviously present in our experiments, as demonstrated by the sensitivity maps ( Figure 3) and by the prominence of the gate signal on all stations ( Figure 6). One can imagine in a natural, uncontrolled settingespecially one in a seismically noisy environment such as a volcano-that we run the risk of interpreting seismic signals from earthquakes, explosions, or human activity as local debris flows if the method is blindly applied to all recorded signals. One way to overcome this problem would be to use a dense array of at least three stations and to run the analysis only on the signal of the "beam" coming from the direction of the channel segment of interest at the expected surface wave velocity range. Another possibility to distinguish local sources from other more distant sources is to look for the presence of higher frequencies in the signal even if they are not directly used because high frequencies (hundreds of Hz) are present for near-station sources but not usually for distant sources due to attenuation (e.g., Lavigne et al., 2000;Marcial et al., 1996).
Another major limitation is the influence of the channel substrate on the observed seismic signal. Channel substrate characteristics vary from site to site and even change over time at a given site. Observations have revealed that softer substrates consisting of loose sediment significantly dampen the seismic signal of a debris flow as opposed to harder substrates such as compacted sediment and bedrock (e.g., Cole et al., 2009;Kean et al., 2015), likely because the softer sediment dissipates the energy of impacts and lengthens the contact time. The method presented in this study could be used to better understand this effect in future research, but it would also have to be modified and validated for deformable substrate conditions.

Conclusions
We have demonstrated a technique that allows us to obtain seismic estimates of basal fluctuating forces of a debris flow that compare well with independent basal force plate measurements. Our technique uses a single three-component seismic station and normal and shear EGFs that can be generated with a force hammer. A more stable solution can be obtained by taking the median of the solution computed independently for several hammer locations within the station's area of sensitivity. There are limitations and drawbacks to the single station approach we use: (1) Other seismic sources can be incorrectly mapped as flow energy at our channel segment; (2) we need three-component seismic data and both shear and normal EGFs; (3) we are limited in the accuracy of what we can attain to the frequency band in which the seismic data and EGFs both have sufficient energy; and (4) we have validated the method only for a rectangular concrete channel. In spite of these limitations, the method we demonstrate represents a step toward advancing seismic debris flow monitoring into a more quantitative domain.
The next step in this approach entails validating the technique in a natural laboratory where flows and channel substrates are relatively simple and well studied, and independent measurements from force plates, flow depth sensors, and other instrumentation are synchronously collected with seismic data (e.g., Kean et al., 2015;McCoy et al., 2013). Once validated in such a setting, we can start to explore how to scale the method up or modify it for larger channels such as the wide, braided, sediment-choked channels common in volcanic settings. Another possibility would be to expand the method to solve for multiple hammer locations simultaneously to obtain a spatiotemporal distribution of fluctuating forces that does not include overlaps in sensitivity, as our present method does.
Beyond further validation and scaling of the method, some basic scientific questions arose in our study that need further investigation if such quantitative methods are to be widely useful. Though our parameterization of the debris flow seems sufficiently valid in this case because it compares well to independent force plate measurements, our assumption that the equivalent forcing could be decomposed into smaller subareas where forcing has the same statistical characteristics but is completely random in time and space within a given time window likely would break down for flows that experience coherent accelerations of packets of material due to variations in topography or flow structure. We need to better understand the spatiotemporal coherency of basal forcing in natural debris flows. We currently lack a complete understanding of the actual physical processes that generate the fluctuating forces, whether it be individual impacts or force chains or something in between.
We also showed here that measuring the fluctuating basal forces does not directly constrain the flow depth or normal stress or any other bulk flow property, and therefore, the seismic wavefield also does not directly measure such properties. Within the framework of the current experimental setup and available data, we cannot say anything conclusive about the physical relationship between fluctuating stresses and any of the examined bulk characteristics other than that there are definitely systematic trends among many flow factors of interest. However, the relationships are complex and likely involve interrelated factors. Moreover, all flow properties also correlate and relate physically with each other, so a simple scaling of fluctuating stresses with any one measured flow property is bound to fail because many interrelated factors are clearly at play. A combination of properties is likely required. The interdependence of multiple flow properties explains why it has been difficult to find universal relations between along-channel seismic amplitudes and specific flow properties despite the obvious correlations that have frequently been observed. This is the case even in our relatively simple experimental setup. In nature, even more factors can vary simultaneously. A multiparametric inversion of many flow properties at once may shed some light on the problem, as could numerical modeling. And although not definitive, the relationships we observed are still encouraging because they mean with further study, there may be a way to refine and interpret the relationships, enabling us to use seismically derived fluctuating stresses to quantify bulk flow properties.
The distinct behavior of the unsaturated flow front, combined with the observed negative correlation with bulk density, and the observation that the ratio of the fluctuating forces to the mean forces drifts systematically downward as the flow becomes muddier and finer grained all suggest that the key to the relationship of fluctuating stresses (and thus the seismic wavefield) with other bulk flow properties lies in the style of the debris flow. In this context, flow style includes the degree of agitation, how the grains and fluids interact with each other, the distribution of grain sizes, and how flows interact with channel roughness. Flow style is likely a primary control on the conversion of bulk translational kinetic energy, which contains information on many flow properties of primary interest (thickness, density, and velocity), to vibrational kinetic energy (granular temperature and larger-scale flow structures) that generates seismic waves.