Joint Sensing of Bedload Flux and Water Depth by Seismic Data Inversion

Rivers are the fluvial conveyor belts routing sediment across the landscape. While there are proper techniques for continuous estimates of the flux of suspended solids, constraining bedload flux is much more challenging, typically involving extensive measurement infrastructure or labor‐intensive manual measurements. Seismometers are potentially valuable alternatives to in‐stream devices, delivering continuous data with high temporal resolution on the average behavior of a reach. Two models exist to predict the seismic spectra generated by river turbulence and bedload flux. However, these models require estimating a large number of parameters and the spectra usually overlap significantly, which hinders straightforward inversion. We provide three functions contained in the R package “eseis” that allow generic modeling of hydraulic and bedload transport dynamics from seismic data using these models. The underlying Monte Carlo approach creates lookup tables of potential spectra, which are compared against the empirical spectra to identify the best fitting solutions. The method is validated against synthetic data sets and independently measured metrics from the Nahal Eshtemoa, Israel, a flash flood‐dominated ephemeral gravel bed river. Our approach reproduces the synthetic time series with average absolute deviations of 0.01–0.04 m (water depth, ranging between 0 and 1 m) and 0.00–0.04 kg/sm (bedload flux, ranging between 0 and 4 kg/sm). The example flash flood water depths and bedload fluxes are reproduced with respective average deviations of 0.10 m and 0.02 kg/sm. Our approach thus provides generic, testable, and reproducible routines for a quantitative description of key metrics, hard to collect by other techniques in a continuous and representative manner.


Introduction
Understanding the boundary conditions and nonlinear dynamics of bedload transport by streams is essential for understanding process geomorphology and long-term landscape evolution and also from an engineering and hazard perspective, especially if the transport happens under episodic flood conditions. Accordingly, there has been significant effort in collecting instrumental data on important parameters determining flow characteristics and boundary conditions. Classic approaches involve either labor-intensive manual sampling (e.g., Bunte & Abt, 2005;King et al., 2004) or the permanent construction of monitoring infrastructure in the stream bed (e.g., Habersack et al., 2016). Any sensors within the stream need to be sufficiently resilient to maintain operation under the harsh conditions during flood events (Geay et al., 2017). Typical in-stream observatories include pressure gauges, temperature sensors, and turbidity sensors. Bedload dynamics are monitored with time-resolving slot samplers (Cohen et al., 2010) and acoustic impact sensors, such as pipe microphones, geophones and accelerometers, or plate geophones (Mizuyama et al., 2010;Rickenmann, 2017). All acoustic bedload sensors, with the exception of hydrophones deployed in the water column (Geay et al., 2019), deliver direct and indirect data on the target parameters, provide point measurements, or can at best be installed along a line crossing the channel (e.g., Hilldale et al., 2014), whereas interest is often directed to the average dynamics of a given reach.
In recent years, a complementary approach has gained increasing attention: streamside instrumentation with seismic sensors (Barrière et al., 2015;Burtin et al., 2008;Roth et al., 2016;Schmandt et al., 2017). Such sensors, typically off-the-shelf seismometers or geophones, are installed at a safe distance from the inundated channel and record the ground motion due to in-stream processes. A sensor can be deployed within less than an hour, records data continuously and autonomously for several months, and is, in principle, able to transmit the data in near real time to processing and evaluation facilities. Hence, seismic 10.1029/2019WR026072 monitoring shows potential for recording bedload flux, which has recently been demonstrated under laboratory and field conditions Schmandt et al., 2017). However, unlike signals derived from bedload impact sensors and similar to the soundscape of rivers recorded by in-stream hydrophones (Geay et al., 2017), seismic signals derive from a multitude of sources (e.g., Roth et al., 2017), and therefore, the identification, extraction, and processing of signals to determine bedload flux are challenging.
Physical models have been suggested to predict the seismic frequency spectra due to bedload transport (Tsai et al., 2012) and due to hydraulic processes within a channel (Gimbert et al., 2014). Dietze (2018) has shown how to use such models to infer water depth for creeks. This involved computing a lookup table of potential spectra that differ only due to changes in river depth and identification of the best reference data fits to the time series of empirical spectra. Here, we expand this approach to bedload flux, based on the notion that the spectra generated by turbulence and bedload transport should be sufficiently distinct (cf. Dietze et al., 2019;Gimbert et al., 2014). In our approach, fits of the empirical data with precalculated reference spectra are optimized based on random combinations of the target parameters. Applying the approach to a case study at the Nahal Eshtemoa, Israel, we show how seismic stations can be used to continuously estimate key hydraulic and bedload transport parameters. We explore the validity of the approach based on synthetic data and by comparing the model output against independent measurements of target parameters. We show the value of seismic stations to gather insight on the anatomy of bedload transporting floods and discuss the potential and the limitations of the technique.

Study Site and Instrumentation
The Nahal (river) Eshtemoa is an ephemeral, flash flood-dominated gravel bed river in the semiarid northern Negev Desert, Israel, draining the southern Hebron mountains in a catchment of about 112 km 2 . Close to the town of As-Samu, the stream crosses a gently undulating landscape in an alluvial valley. A straight, 5-m-wide reach with 1-m-high banks and carbon-cemented gravels is instrumented by a comprehensive in-stream observatory (Laronne et al., 1992), including Reid-type slot samplers, plate geophones, a pipe microphone, water quality sensors and sampler, as well as pressure transducers for the determination of water depth and water surface slope (Powell et al., 2001;Reid & Laronne, 1995). Since 2016, a Nanometrics TC120s broadband seismometer has been installed in the right bank (Figures 1a and 1b). It is sampled by a Nanometrics Centaur data logger at a recording frequency of 200 Hz and a gain of 2.

Computational Environment
The R package "eseis" (Dietze, 2018) is a free and open source toolbox for handling the work flow of generic environmental seismology. With the latest developer version (0.5.0), it contains models to predict the seismic spectrum due to turbulent channel flow (Gimbert et al., 2014) and impacting bedload particles (Tsai et al., 2012). Both models can be explored in an interactive graphical user interface (GUI; Figure 2). Three additional functions, denoted by the prefix fmi, are devoted to the approach of fluvial model inversion presented in this study. Data preparation, processing, analysis, and visualization steps were performed with R v. 3.5.3 (RCoreTeam, 2015). The R functions, data sets, and utilized scripts as commented markdown files are available in the supporting information (https://osf.io/5uzpw/) to reproduce the presented results.

Data Processing
Flood water depth and bedload flux time series were recorded at 1-min resolution. The bedload flux time series starts when at least 4 kg of sediment has been collected in the slot samplers during an event, which represents the sensitivity threshold of the sensors. We used the median of the values measured by the three bedload samplers to generate a representative bedload flux per unit stream width. The recorded seismic files were converted to hourly SAC files (IRIS, 2017) and organized in the consistent structure as used by the functions of the eseis package. For the relevant part of the flood (05:40 to 11:00 UTC; cf. Figure 1b) we calculated a spectrogram from the vertical component of the seismic time series using the method of Welch (1967) with 10-s-long, nonoverlapping windows, created by averaging 5-s-long and 80% overlapping subwindows ( Figure 1c).

Model Approach
Our approach assumes that the recorded seismic spectrum is dominated by channel activity, that is, a combination of turbulent flow and sedimentary particles impacting the channel floor during bedload transport, whereas other sources such as the effects of wind and rain or anthropogenic activity are of subordinate importance. Under these conditions, we can exploit the combination of the seismic models of Tsai et al. (2012) and Gimbert et al. (2014).
The seismic spectrum due to particles impacting the bed is based on the Hertzian impact theory. Impact energies are calculated for a mobile sediment layer, composed of particles whose sizes are drawn from a raised cosine distribution function (Tsai et al., 2012). The layer is further characterized by its velocity, height, and settling velocity. The seismic power is calculated for each frequency increment of the output spectrum, accounting for seismic ground properties that determine frequency-dependent amplitude changes. The seismic spectrum due to the interaction of turbulent fluid flow with bed particles (Gimbert et al., 2014) is described by drag, lift, and cross-stream force fluctuations. The force fluctuations are converted to ground velocity and, subsequently, to a seismic power spectrum representative of the vertical spatial component. We used a Monte Carlo approach to randomly vary the two parameters of interest, water depth, and bedload flux, to generate 5,000 different potential seismic conditions that serve as a lookup table. In addition, to account for flow without bedload transport, we calculated another 1,000 realizations where bedload flux was set to zero and only water depth was varied. In the Nahal Eshtemoa case, we allowed water depth (h w ) to range from 0.01 m (minimum value required to allow model evaluation) to 1.20 m (120% of bankfull depth). Bedload flux q s was varied between 0 and 15 kg/sm (200% of the range reported for other floods, (cf. Cohen et al., 2010). The selected boundaries are arbitrary and can be extended, if needed-for example, when the model output yields values that clearly undershoot the expected empirical data. For each parameter combination, we generated a seismic reference spectrum and calculated root-mean-square errors with the corresponding observed spectrum. For each time step, we then selected the values for water depth and bedload flux corresponding to the reference spectrum with the smallest root-mean-square error. To account for short-term variability of the seismic record, the model results were smoothed with a running average (R package caTools v. 1.17.1.2, Tuszynski, 2014) using a window size of 18 samples, that is, 180 s.

Estimation of Unknown Model Parameters
Both the turbulence (Gimbert et al., 2014) and the bedload (Tsai et al., 2012) models require constraints on a set of 17 parameters ( Table 1). Some of these parameters can be determined from field measurements, namely, the median gain size D 50 (d s ), logarithmic grain-size standard deviation (s s ), channel width (w w ), channel bed gradient (a w ), and the distance between the center line of the river and the seismic station (r 0 ). Other parameters can be estimated at reasonable accuracy based on prior measurements, such as the specific density of the fluid (r w ) and of the bedload material (r s ). And yet others are simply set according to computational needs and convention, such as the reference frequency ( f 0 ), frequency range (f), and resolution (res) for which the model yields results. Several parameters describe the seismic ground characteristics due to the site properties. This set of parameters (material quality factor q 0 and its increase with frequency e 0 , Rayleigh wave phase velocity at the reference frequency v 0 and its variation coefficient p 0 , and the Greens function Greens function displacement amplitude coefficients -n 0 0.5, 0.8 Note. Values in parentheses give initial parameter ranges used to find the most plausible values for seismic model parameters (see section 2.5). a Cohen et al. (2010). displacement amplitude coefficients n 0 ) can be constrained by performing an active seismic survey (e.g., Foti et al., 2018). However, when that is not possible (as is the case in this study), they must be estimated.
In a first step we make use of the interactive GUI provided with the R package eseis ( Figure 2). This application allows changing all relevant model parameters and instantaneously visualizes the updated model outputs, together with an optionally provided empirical spectrum. We used this tool to explore the meaningful parameter space, which is able to create model spectra that match the overall shape and amplitude of a series of empirical spectra. We focused on empirical spectra at the beginning of the flood event, where sharp rises of broadband seismic signals (Schmandt et al., 2017;Tsai et al., 2012) indicate pulses of bedload movement close the seismic sensor and later stages of the flood when the bedload signal is no longer visible in the seismic spectrogram and most of the seismic signal is presumably generated by turbulence. We adjusted the parameters q 0 , v 0 , p 0 , e 0 , and n 0 to roughly match the shape of the resulting fluvial, bedload, and joint spectra to the empirical example spectra mentioned above. Thereafter, we changed the parameters water depth and bedload flux to adjust the seismic power of the model spectra until they visually matched the empirical spectra. The quality of the match was subsequently quantified and optimized by minimizing the root-mean-square error. From this set of combinations optimized to first order we started changing the seismic parameters toward lower and higher values, respectively, until the match of empirical and model spectra obviously diverged. We defined these parameter ranges as the limits for the subsequent step of parameter range optimization. In a second step we performed the inversion of the example flood data set in an extended Monte Carlo experiment. Since the lower and upper Greens function parameters n 0 did not have significant impact on the model spectra shape when changing them between 0.4 and 0.8 and 0.5 and 0.9, respectively, we set them arbitrarily to 0.5 and 0.8. We created 10 5 random parameter combinations of the most sensitive parameters (q 0 , v 0 , p 0 , and e 0 ) and the target parameters (h w and q s ), exploring the range of the former set of parameters to identify the most likely values throughout the event (i.e., the medians of the distributions).

Model Validation
In order to infer the ability of the model approach to estimate water depth and bedload flux, we created several synthetic data sets, inverted them, and compared the resulting model time series of the target parameters to the input data ( Figure 3). All except the varied parameters are the same as denoted in Table 1. Synthetic Data Set 1 imposes a constant water depth of 0.5 m. The bedload is injected after 2 hr of the modeled time period, resulting in an instantaneous rise to 5 kg/sm, which is held constant for another 2 hr Four synthetic data sets were tested and are organized by columns (a-d).
Each panel shows the resulting synthetic spectrogram (PSD), the fit deviation matrix depicting the root-mean-square error between empiric spectra and best fit reference spectra, the input (blue line shows water depth and orange line bedload flux) and modeled time series (black lines), and the distribution of model errors (residuals) in target parameter units.
until it is reduced linearly to zero for the rest of the time. This data set is mainly used to test the sensitivity of the model to fluctuations in a parameter when the other is changed. Synthetic Data Set 2 assumes synchronously rising and falling water depth and bedload flux, both of which are modeled as lognormal distribution curves. This scenario reflects a river where water depth and bedload flux do not show a hysteresis effect and where the seismic signal overlap is constant through time. Synthetic Data Set 3 features a lognormal bedload time series that rises steeper and narrower than the lognormal water depth time series, thus generating a bedload wave traveling in front of the flood wave. This scenario inherits a clockwise hysteresis pattern. Synthetic Data Set 4 uses the empirically measured water depth and bedload flux values to generate a seismic spectrogram. It is used to explore how precisely the target variables can be estimated by the model approach under ideal conditions: all panel shows the resulting synthetic spectrogram (PSD) water and bedload flux.
Model quality is assessed by the absolute difference between synthetic spectra and the best fit modeled reference spectra. This error can be studied both in time and frequency space. Another measure of model quality is the error (residual) between water depth or bedload flux and the respective model estimates.

Characteristics of the Flood
The flood hydrograph (data available under https://osf.io/5uzpw/) shows a rapid rise of water depth although the actual onset of the event is not shown here because we define the event by the onset of the bedload sampler records, that is, at 05.40 UTC. After the flood's double peak occurred (0.84 and 0.83 m), water depth dropped logarithmically for at least 13 hr ( Figure 1b). The three bedload samplers monitored a maximum average value of 4.29 kg/sm. The highest bedload fluxes were recorded within the first 2 min. Thereafter, values declined progressively to almost zero around 05:55 UTC, when two further smaller bedload waves (peak flux 1.08 kg/sm) emerged for 30 min. Bedload transport ceased at 07:10 UTC. With the onset of the flood, the seismic spectrogram shows a broadband (10-90 Hz) increase in seismic power up to −100 dB, which progressively grades into the background for about 1 hr. At about 07.50 UTC, a period of broadband spike appearance occurs that lasts for at least 2.5 hr.

Model Validation With Synthetic Flood Time Series
The ability of the model to reconstruct the synthetic time series of target parameters (which were used to generate noise-free spectrograms that were inverted) provides the accuracy baseline for the actual inversion of the empirical data set. Synthetic Data Set 1 (Figure 3a) yielded absolute differences between best fit model and input spectrogram of less than 0.5 dB and target parameter errors of 0.02 ± 0.04 m (water depth) and −0.03 ± 0.06 kg/sm (bedload flux). The modeled time series resemble the onset of changes and are only slightly affected by changes in the corresponding parameter. Synthetic Data Set 2 (Figure 3b) has only minor spectral differences (less than 0.26 dB) and model errors (0.02 ± 0.04 m and −0.06 ± 0.13 kg/sm, respectively). The concurrent changes in water depth and bedload flux are captured well. However, during the second half of the synthetic event the model produced increasingly larger deviations. Synthetic Data Set 3 (Figure 3c) has the largest spectral differences (up to 1.75 dB) but yielded the smallest target parameter errors (−0.01 ± 0.03 m and −0.001 ± 0.03 kg/sm, respectively). These errors mainly appear toward the end of the synthetic data set, when the continuously declining water depth curve is represented by stepwise model results. The synthetic data set produced by the real world time series of water depth and bedload flux (Figure 1b) produced spectral differences of up to 0.47 dB and target parameter errors for water depth and bedload flux of −0.04 ± 0.03 m and −0.001 ± 0.02 kg/sm, respectively. The water depth is thus overestimated, especially when bedload transport ceases.

Model Parameter Estimation
Explorative model parameter adjustments (Figure 2) revealed that the shape of the fluvial and bedload model spectra can vary significantly. In turn, the parameter range that lets the models and their summed effect converge in shape to those of the empirical spectra during the peak water depth and the falling limb of the flood is small. Thus, we defined the limits within which q 0 was allowed to vary to 15-20, for v 0 to 800-900 m/s, for p 0 to 0.4-0.7 and for e 0 to 0.01-0.25 (cf. Table 1). As expected, changes in the input parameters water depth and bedload flux result in amplitude shifts with no visible effects on the shape of the spectrum (Figure 4a). In contrast, higher ground quality factors (Figure 4b) lead to systematic increase of spectral power especially for higher frequencies, an observation which is not visible in the empirical data ( Figure 1c). A contrary effect occurs for the Rayleigh wave phase velocity v 0 (Figure 4c), where lower frequencies decrease in power with increasing parameter values. The wave velocity variation coefficient p 0 (Figure 4d) mainly affects the amplitude of the bedload spectrum and the convexity of the turbulence spectrum. The parameter describing quality factor increase with frequency e 0 (Figure 4e) results in declining seismic power for higher parameter values. This parameter is not included in the turbulence model and has therefore no effect on the latter.
Running the Monte Carlo approach with the range of seismic parameters as defined in Table 1 yielded convergent results with median values and quartiles of the distributions well within the defined parameter range (Figure 4f). The effect of the parameters is independent of each other. Thus, the best fitting combination of parameters for each of the 10-s-long empirical spectra can in principle be anywhere within that imposed range. Since this is not the case the parameter distribution is assumed to be unimodal and adequately represented by the median as a most likely value. Therefore, we chose the medians (q 0 = 16.77, v 0 = 859, p 0 = 0.62, e 0 = 0.07) for the subsequent Monte Carlo run to estimate the actual target parameters.

Model Results for the Empirical Data Set
The seismic data of the example flood event (Figure 1c) show contribution of the expected frequency bands between 5 and 70 Hz (Gimbert et al., 2014;Schmandt et al., 2017;Tsai et al., 2012). However, above 70 Hz there is increased seismic energy. That pattern appears to be a horizontally flipped version of the <70-Hz signals and cannot be physically explained. Therefore, and to avoid introducing a systematic bias, we truncated the spectrogram to the frequency range 10-70 Hz, an interval to which the seismic models are most sensitive. Furthermore, to reduce scatter in the frequency domain and to improve computational speed (the frequency vector of the raw spectrogram had 1,000 values), we spline interpolated the frequency vectors of the spectrogram to 100 values between 5 and 70 Hz, corresponding to the modeled spectra (cf. supporting information S1).
The best fit spectra deviations (Figure 5b) range between 0 and 15 dB. The highest deviations appear at the continuous narrowband signals (23, 47 Hz) as well as during the period with numerous short-term, broadband signals (7:50-10:10 UTC). Smaller deviations, up to 10 dB occur during the early stage of the flood (5:40-7:50 UTC). They affect the upper and lower frequencies of the modeled spectra as well as the central bands (30-45 Hz).
The modeled water depth (Figure 5c) is in general agreement with the independent water depth measurements, although the falling limb of the flood is underestimated by 0.10 m on average (i.e., median of the absolute deviations). During 7:50 and 10:10 UTC (gray polygon in Figure 1), when the spectrogram (Figure 5a) exhibits several broadband spikes, the model shows significant overestimation effects. Overall, the seismic results are more variable than the 1-min resolution control data (180-s running standard deviations of 0.041 versus 0.029 m). Results of seismic bedload flux are also in the same range as the slot sampler data (0.02 kg/sm average deviation), and most of the short excursions of increasing and decreasing bedload flux values in the slot sampler time series are coincident with the seismic model results, both in terms of timing and amplitude.

Model Quality
The synthetic data sets (Figure 3) allow insight on three different dimensions. First, they show the general applicability and validity of the Monte Carlo-based inversion approach. Second, they provide the baseline of accuracy, that is, the minimum deviations to expect when modeling an empirical data set. Third, the scenarios allow insight as to how different combinations of flood and bedload flux evolutions appear in seismic spectrograms.
In all cases, the input time series were depicted by the model, with average deviations of less than 0.04 m for water depth and 0.04 kg/sm for bedload flux. Thus, for inversions of empirical data sets one should anticipate at least these ranges of model deviations. Under the ideal conditions of synthetic data sets generated without any noise or contribution of additional seismic sources, the deviations of the best fit reference spectra from the synthetic spectra time series (deviation matrices in Figure 3) are negligible. An exception is test Data Set 3 (Figure 3c), which shows misfits of up to 2-dB coincident with the step-like evolution of the modeled water depth during times of virtually zero bedload flux. This step-like behavior disappears when more than the 1,000 Monte Carlo cycles are used to generate the reference spectra (not shown). Thus, it is important to provide a sufficiently large number of potential parameter combinations for the reference spectra, especially when better model fits for the falling water depth limb are of interest. A similar effect is visible in the fourth synthetic data set ( Figure 3d) and, to a lesser degree also in Data Set 2, where the falling water depth curve is systematically overestimated as soon as bedload flux fades.
The imposed time series of synthetic Data Set 1, constant water depth, and a step-like onset of constant bedload movement throughout an event are far from what one would expect in natural systems. However, this scenario shows that model deviation is systematically higher for times when only one of the two expected seismic sources is active (i.e., water depth is overestimated when no bedload is transported). In the case of synchronous evolution of flood stage and bedload flux (Data Set 2) the model results maintain this synchronicity. This is encouraging because when a seismically derived data set exhibits such a pattern, it is difficult to judge merely from the properties of the spectra, whether there indeed are two seismic sources present. In the case of a bedload wave traveling in front of a flood (Figure 3c), that is, a clockwise hysteresis pattern in the h w -q s relationship, the combined effect of turbulence and bedload movement result in a spectrogram with a trend of rising dominant frequency with time. Such patterns were observed in natural settings, such as a flash flood observatory in New Mexico (Dietze et al., 2019) but could not be attributed to a likely cause. Here, we can provide this cause, which is simply the combination of two seismic sources with different time evolution paths. The trend toward higher frequencies remains visible without any hysteresis effect, albeit weaker (e.g., Figure 3b). Since even water depths up to 0.5 m only contribute as much as −135 dB to the total seismic signal (Figure 3a), it appears that most of the seismic energy is contributed by the bedload part at this distance of channel to sensor.

Evolution of the Flood Event
The example event shows the typical features of flash floods in the Nahal Eshtemoa (Halfi et al., 2018): a suddenly rising water depth that remains unstable due to the high turbulence and a bedload bore at the front of the flood, occasionally followed by further bedload bores. The passing of these bores is recorded in the example flood by both the slot samplers and the seismic sensor (Figures 5d and 5a), the latter showing this as broadband spikes of seismic energy after the onset of the flood. With the end of the bedload transport period the spectrogram only shows noticeable seismic energy between 20 and 50 Hz that gradually decreases in amplitude with time.
This trend is interrupted between 07:50 and 10:10 UTC (gray shaded area in ( Figure 5) by recurring broadband seismic pulses. We interpret these pulses as the effects of persons working at the observatory for data collection and station maintenance reasons. These activities included documented walking and operating at close proximity to the seismic sensor and a car idling at the bank. A further set of seismic signals, the temporally constant, narrowband horizontal lines in the spectrogram (Figure 1c) around 22 and 44, 30 and 60 Hz, are interpreted as the signature of measurement devices operating at the observatory. When excluding this period of signal contamination, the temporal variation of the seismic signal-derived bedload flux shows three important components of the average-channel bedload flux: (i) a first very large wave in bedload flux up to 5 kg/sm, which drastically recedes within minutes of the arrival of the flood bore, (ii) a second multiple-rise peaking at about 1 kg/sm, and (iii) and a third, smaller rise (0.5 kg/sm) with a long, 1-hr recession. Given the 120-s averaging required with respect to the sensitivity of the slot sampler bedload monitoring equipment, it is remarkable that a single sensor deployed on the bank of a river can determine the main and relative features of bedload flux.

Benefits, Limitations, and Outlook
In comparison to classic approaches to constraining hydraulic and sediment transport parameters in fluvial systems, the seismic method introduced here shows several advantages. The sensors can be deployed easily and quickly (placement of sensors in small hand dug pits, connection to rugged loggers, and power supply by off the shelf batteries, (cf. Dietze et al., 2017), at a safe distance from the hazardous conditions of flood-prone streams. Modern seismic stations can record ground motion data at high frequency, under harsh conditions and even transmit the data in near real time to analysis facilities where they can be automatically analyzed.
In the case shown here, and once the data set of reference spectra is precalculated, inverting an empirical spectrum requires less than 1-s computation time on a single CPU. Thus, efficient and near-real-time information about floods and the potentially hazardous bedload they transport can be provided also for remote locations in a continuous manner.
While classic approaches, such as slot samplers, are only able to measure the bedload flux at discrete cross-sectional intervals (the slot aperture in the Nahal Eshtemoa observatory is 11 cm and the devices are spaced about 1 m) and a representative estimate of bedload flux must be based on averaging data of several samplers, the seismic approach implicitly provides an integrated estimate of the full cross section. Thereby, due to inelastic attenuation and geometric spreading effects, the amplitudes of the seismic wave field decay with distance to the source. Thus, for small distances between river and seismometer (like in our case), the seismic approach will emphasize bedload flux values closer to the sensor. The magnitude of that effect may be approximated by changing the parameter distance to river r 0 in the interactive GUI ( Figure 2). The supporting information contains another synthetic test data set that shows how bedload fluxes of 1 and 4 kg/sm (and vice versa) for the left and right channel half, respectively, result in spectral power differences of 0.5-3.5 dB, which would yield different model inversion results. However, a more robust field experiment with actively moved pebbles in the channel would be more appropriate (cf. Schmandt et al., 2017) to shed light onto the actual effects of nonuniform bedload flux across a channel.
None of the available bedload formulae can replicate such natural fluvial sediment wave phenomena as presented here (e.g., Cudden & Hoey, 2003;Gomez et al., 1989), even though theoreticians, notably Einstein (1950) and experimentalists (e.g., Aberle et al., 2012;Dhont & Ancey, 2018;Ghilardi et al., 2014;Iseya & Ikeda, 15, 1987;Lisle et al., 2001) have long been aware of their presence. Indeed, based on a century of geomorphological research, it is known that fluvial systems are complex (Schumm, 1991(Schumm, , 2005; they do not transport bedload at certain time scales as simply as an "efficient" machine (Bagnold, 1966) nor merely determined by average reach shear stress (Parker, 1990). Instead, the fluvial system responds in complex manners, as in this case one sensor and the respective technique demonstrate. With the seismic approach we are able to provide robust data with high temporal resolution, which are crucial to determine river activity, river stability, river change, and the transport of bedload to various ecologically sensitive reaches, to reservoirs, and to the oceans.
However, in comparison to classic methods, the seismic approach also has drawbacks. First, the recorded signals represent measurements of ground velocity due to a multitude of sources, which are inverted for the parameters of interest using a combination of physical models. These models are formulated under a series of assumptions (cf. Gimbert et al., 2014;Tsai et al., 2012) and require information about a large number of parameters. Although the model output is in appropriate physical units (m and kg/sm) that does not require development of a further transfer function, they are not direct measurements of the parameters of interest. This point also needs to be considered in the light that the seismic approach does not necessarily reflect the same process as, for example, the Reid-type sampler, which records all particles that fall through the 11-cm-wide slot while omitting all particle that pass between two such slots. The seismic record is an amalgam of the impacts of all bedload particles in a given reach and therefore provides a spatially integrated result, which may differ from spatially discrete direct measurements due to cross-sectional nonuniform bedload fluxes. Likewise, the bedload model assumes a sediment-covered bedrock channel. Apparently, the cemented bed structure resembles this bedrock behavior well enough.
The selection of seismic model parameters is crucial for the inversion results. Thus, at best one performs an active seismic survey to independently constrain these parameters (e.g., Foti et al., 2018). Since this was not possible in this study, we introduced a stepwise approach as an alternative: (i) visual exploration of parameter effects on model output with respect to empirical seismic observations under partly known flood conditions (Figure 2), (ii) long Monte Carlo chains to identify the parameter combination that best explains the empirical data set (Figure 4f), before (iii) actually inverting the data with the most plausible set of parameters ( Figure 5) along with relevant metrics for model errors.
Seismic sensors are not only subject to the seismic sources of interest but also record a range of further processes, as the period of maintenance activities shows. Atmospheric processes such as wind and rain (Dietze et al., 2017;Roth et al., 2017) generate seismic signals in a similar frequency range. Burtin et al. (2008) and Cook et al. (2018) showed that the seismic footprint of rivers and the bedload they transport can be detected over tens of kilometers. Thus, nearby trunk streams may also add their seismic signature to the signals recorded at the stream of interest. Therefore, the deployment site for a seismic station intended to record water depth and bedload flux must be chosen with care. They should be out of the range of unwanted seismic sources such as roads and railroads, industrial buildings with running machines, should be shielded from the signals of wind and rain (at best by burying the sensor several decimeters to meters in the ground) and be installed several kilometers from other neighboring streams. If the latter is not possible, the Monte Carlo-based inversion must include the other stream as an additional source of water turbulence and bedload transport.
The approach is vulnerable to transgressive or sudden changes in one or more of the seismic model parameters, for example, if soil moisture changes drastically or frozen ground thaws during the summer period, both of which cause changes in the seismic wave velocity and quality factor (James et al., 2019). Likewise, reorganization of the channel bed by mobilization, redeposition, and injection of material, form example, from bank failures, can change some of the parameters assumed to be stable. Finally, floods beyond bank full depth will result in a sudden and significant change in parameters such as width and depth. Mathematically, the models might be calculated for the different cross sections of the suprabank new river, but this would require setting up more extensive synthetic data sets and exploring the quality of the results of combined model spectra from multiple independent river cross sections.
Future applications of the seismic approach introduced here could be near-real-time warning systems or continuous observation devices for streams otherwise hard to instrument, for example, due to conservation requirements or steep topography. In principle, it is also possible to survey large, navigable rivers with high bedload fluxes during floods, as long as the position of the sensor(s) is chosen carefully to minimize the overlap of spectral components and recording of other seismic sources. A continuous record of bedload transport in combination with time series of suspended sediment load opens the perspective for the holistic view on catchment-wide sediment dynamics. Finally, installation of a series of sensors along a stream over a greater distance allows for tracking and detailed insight into flood waves, as recently highlighted for a lake outburst flood in Nepal (Cook et al., 2018). The generic layout of the inversion approach, as illustrated during the seismic parameter range estimation, can in principle be used to invert for parameters other than water depth and bedload transport, as well. Given that all model parameters are well constrained, one can explore reorganization of the bed by comparing model fits with respect to grain-size distribution parameters (s d and s s ) from data before and after a flood event.

Conclusions
The seismic method is a valid approach to quantifying key hydraulic and bedload transport parameters, not merely as proxy data in its own data dimension and unit space (i.e., dB) but as estimates of the target parameters in the respective units: water depth in meters and bedload flux in cubic metres per second and metre or kilogram per second and metre. However, this is only possible if (i) one or more stations are placed at appropriate distances from the river as seismic source, so that the signals of both sources are powerful and distinct enough to be mapped out in the spectra, (ii) the empirical data are free of (or cleaned from ; e.g., Bottelin et al., 2013) unwanted signal components, and (iii) the relevant model parameters are sufficiently well constrained, either by independent measurements or at least by optimizing free parameters with respect to the target parameters during a control period. The approach yields a continuous output with average deviations of 0.10 m (water depth) and 0.02 kg/sm (bedload flux), respectively.
The comparably uncomplicated and quick installation, potential of almost real-time data transmission and quick processing, renders the seismic approach a complementary source of data otherwise difficult to obtain. This opens up perspectives such as exploring the boundary conditions that control the onset of motion in episodically active river systems, investigating the coupling of processes that shape different landscape elements such as rock walls, debris flows, bank failures, and migrating rivers and deliver high-resolution field data to long-standing concepts of fluvial geomorphology.
The model code has been implemented using a user-driven, free, and open software environment. Sensors and data loggers are becoming more and more affordable. The density of existing seismic networks along with the availability of their measurement data increases progressively. These three tendencies provide the base for other scientists to engage with the method, develop their own measurement systems, or make use of the large amount of existing data to pursue their research hypotheses. Seismic and stream observatory data, as well as analysis scripts and R functions (intended for transparency and reproducibility reasons of the introduced routine), are provided in the supporting information (DOI 10.31223/osf.io/n5gcm; https://osf.io/ 5uzpw/). The latest version of the R package eseis, containing the functions used in this study is available on GitHub (https://github.com/ coffeemuggler/eseis/). The study was funded by the Israel Science Foundation grant 832/14 and the BSF grant 2018619 both to J. B. L. M. D. is funded through the H2020 Marie Curie action ITN SUBITOP (Grant 674899). All field work was supported by GFZ internal funds. We thank the Editor and three reviewers for their encouraging and constructive input that helped improve the quality and impact of the paper.