Passive acoustic measurement of bedload grain size distribution using self-generated noise

Monitoring sediment transport processes in rivers is of particular interest to engineers and scientists to assess the stability of rivers and hydraulic structures. Various methods for sediment transport process description were proposed using conventional or surrogate measurement techniques. This paper addresses the topic of the passive acoustic monitoring of bedload transport in rivers and especially the estimation of the bedload grain size distribution from self-generated noise. It discusses the feasibility of linking the acoustic signal spectrum shape to bedload grain sizes involved in elastic impacts with the river bed treated as a massive slab. Bedload grain size distribution is estimated by a regularized algebraic inversion scheme fed with the power spectrum density of river noise estimated from one hydrophone. The inversion methodology relies upon a physical model that predicts the acoustic field generated by the collision between rigid bodies. Here we proposed an analytic model of the acoustic energy spectrum generated by the impacts between a sphere and a slab. The proposed model computes the power spectral density of bedload noise using a linear system of analytic energy spectra weighted by the grain size distribution. The algebraic system of equations is then solved by least square optimization and solution regularization methods. The result of inversion leads directly to the estimation of the bedload grain size distribution. The inversion method was applied to real acoustic data from passive acoustics experiments realized on the Isère River, in France. The inversion of in situ measured spectra reveals good estimations of grain size distribution, fairly close to what was estimated by physical sampling instruments. These results illustrate the potential of the hydrophone technique to be used as a standalone method that could ensure high spatial and temporal resolution measurements for sediment transport in rivers.

1 Introduction

Sediment transport analyses in river catchments are one of the key activities stipulated by the European water framework directive (European Commission, 2001) and also applied in French environmental policies. Climate changes and anthropological actions impact the sediment transport in rivers such that it produces changes in the river morphology and may put at risk ecosystems and hydraulic structures, eventually. One of the major concerns of sediment transport in rivers is determining the total discharge of bedload transport (Gray et al., 2010). Bedload transport models are highly sensitive to incipient motion, which is directly related to river bed grain size distribution (GSD). Bedload GSD is linked to both surface and substrate GSD. In his paper, Parker (1990) constructed a two-size fraction transport model, assuming that the bedload GSD is identical to substrate GSD, for stable armored river beds, and becomes identical to surface GSD whenever the armor is destroyed. The development of surface-based and mixed-size transport models has received considerable attention (Heimann et al., 2015; Kuhnle, 1993; Parker, 1990; Recking, 2016; Wilcock and Kenworthy, 2002; Wilcock and McArdell, 1993). Knowing the bedload GSD solves the problem of initiation of motion and, therefore, enhances the accuracy of transport rate prediction. Therefore, measuring bedload leads to not only transport rates, but also to bedload GSD to calibrate models (Parker, 2002; Wilcock et al., 2009). However, obtaining bedload samples during exceptional hydraulic events may be difficult by using traditional bedload sampling techniques (e.g., pressure-difference samplers) (Bunte et al., 2010). To measure a wide range of discharge flows, the scientific community has been interested in developing indirect, or surrogate, methods that achieve continuous measurements no matter the hydraulic conditions (Gray et al., 2010; Hubbell, 1964). This paper is dedicated to the monitoring of bedload GSD using the acoustic noise naturally generated by bedload transport in rivers, the so-called bedload self-generated noise (SGN).

Acoustics surrogate methods are divided into two categories: active and passive methods (Gray et al., 2010; Hubbell, 1964). Examples of active methods are the acoustic Doppler current profiler, aDcp (Rennie and Millar, 2004), or the acoustic mapping velocity technique (Muste et al., 2016). Active methods use emissions of well-known signals but, actually, to the best of our knowledge, no active instrument was conceived to estimate bedload GSD. Besides, the major problem of the active instruments is that they do not properly behave during high flow discharges. This is why the passive instruments are preferred instead of the former. These instruments use seismic or acoustic signals generated by bedload particle impacts. Recorded signals contain information on both sediment impact rate and bedload particle sizes. One of the most widely used techniques consists in recording the signal of particle impacts on steel objects like plates (Rickenmann et al., 2014; Wyss et al., 2016a), pipes (Mao et al., 2016; Mizuyama et al., 2010) or column pipes (Papanicolaou et al., 2009). Other passive instruments consist in directly recording bedload SGN by using passive acoustic monitoring (PAM) (Barton, 2006; Bedeus and Ivicsis, 1963; Geay, 2013; Geay et al., 2017a; Thorne, 1986a, b) or seismic monitoring (Gimbert et al., 2014; Roth et al., 2016; Tsai et al., 2012). Measuring bedload GSD with passive methods has been achieved using plates (Barrière et al., 2015; Krein et al., 2014; Rickenmann et al., 2014; Wyss et al., 2016b) or pipes (Dell'Agnese et al., 2014; Mizuyama et al., 2010; Papanicolaou et al., 2009), and SGN (Geay et al., 2017a; Johnson and Muir, 1969; Jonys, 1976; Thorne, 1986a, b), by using experimental laws of calibration. Concerning seismic methods, bedload GSD measurements were not yet proposed as a direct application.

The existence of a link between the GSD and the features of vibrational signals has been demonstrated in several experiments (Belleudy et al., 2010; Bogen and Møen, 2001; Krein et al., 2008; Turowski et al., 2011). Coupling geophones with steel plates (Barrière et al., 2015; Wyss et al., 2016a) produced composite power laws by linking both peak amplitude and peak frequency to the grain size. Using the Japanese pipe, Mao et al. (2016) proposed an empirical model based on multi-channel recorded amplitude ratios to estimate different percentiles of grain diameters (D₁₆, D₅₀ and D₈₄). The only metric exploited in this kind of measurement is the amplitude of shocks on steel structures. Thus, these passive techniques involving shocks on steel structures offer a high-quality signal, or signal-to-noise ratios (SNRs). The analyzed physics is the same as in the case of SGN measurements by PAM, which is the rigid body radiation caused by Hertzian impacts between sediments. In the case of SGN measurements, unlike the steel structure impact measurements, the SGN signal amplitudes are not usable for grain-size inversion because of the issues concerning the sound propagation throughout the reach (the amplitudes depend on the distance between the shocks and the hydrophone). This makes the amplitude a futile metric to infer grain-size information from SGN signals.

Several studies in the field highlighted that the frequency content (i.e., spectrum shape) of SGN signals is heavily dominated by grain sizes. For example, Jonys (1976) showed by laboratory experiments with ceramic spheres that spectral peak frequency is linked to sphere diameter. The author found a peak frequency at about 4 kHz for 19 mm diameter particles, at 2.2 kHz for the 38 mm diameter and at 1 kHz for the 75 mm diameter. This means that a doubling of grain size is almost equivalent to halving of peak frequency. Extensive research on GSD estimation by SGN recordings was done by Thorne (1986b), where he presented two strategies for inversion of acoustic spectra to estimate GSD. Results were encouraging as GSD was roughly estimated. These techniques are based on experimental measurements that have been made in a rotating drum with specific conditions that are different from the conditions found in rivers (e.g., impact velocities, acoustic propagation). Besides, his inversion techniques raise issues because of the broadband nature (shape) of spectra, even for uniform sediments. The author himself assumed that this was the major cause of inaccurate estimations of GSD from composite spectra.

This paper proposes an inversion method that solves the issue of spectrum shape and which accurately estimates the entire bedload GSD curve. This proposed method is conceived as being transferable to a large set of operational contexts. The procedure of inversion is based on a physical direct model which is presented in the first part of this paper. In the second part, the inversion algorithm is presented in the form of a technique for solving least square (LS) problems with a regularization condition about the positivity of the GSD curve. Simulated acoustic spectra and their inversion are used to test the robustness of LS methods to measurement uncertainties. In the third part, the LS inversion algorithm is applied to field measurements made in the large gravel-bed Isère River, France. GSDs estimated with our method are compared to GSDs measured with a pressure-difference sampler. Additionally, the cross-sectional variability of bedload GSD is analyzed using both acoustic and direct measurements. Finally, results are discussed to give a technical overview of the proposed inversion method.

2 SGN model

2.1 Analytic model of Hertzian impact between a sphere and a slab

This section deals with spectral modeling of the impact between a sphere and a slab (Akay and Hodgson, 1978; Hunter, 1957), because the main assumption of this study is that the acoustics of gravel is described by impacts between bedload sediments and the river. To prove the validity of our model, the study includes some comparative facts with the sphere–sphere spectral model of Thorne and Foden (1988).

As a brief introduction, the collision between bed particles radiates energy. Such a rigid body radiation phenomenon is due to both vibrations and accelerations. These processes are very well separated with respect to their dominant frequencies, such that the spherical mode vibrations generate much higher frequencies than the acceleration-based sound (Barton, 2006; Thorne and Foden, 1988). The acoustic effect of accelerating rigid bodies is physically modeled by Kirchhoff (1883). A framework was constructed by Goldsmith (2003), Hertz (1882) and Hunter (1957) to model acceleration profiles from elastic impacts between two solid rigid bodies like two spheres or a sphere and a slab. In a mathematical sense, the acoustic pressure field generated from the acceleration of a rigid body is evaluated by the integral convolution from Eq. (1) (Akay and Hodgson, 1978; Koss and Alfredson, 1973; Thorne and Foden, 1988). The integral consists of the convolution between Kirchhoff's impulse response p_I and an acceleration profile A. In the case of elastic (Hertzian) impacts the acceleration occurs during the impact, and so the integral is evaluated by intervals with respect to a contact duration T_c. The contact duration T_c is modeled by Hertz's law and it is put in a simplified form in Eq. (1), for both sphere–sphere and sphere–slab impact models.

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }p\left(t\right)=\underset{\mathrm{0}}{\overset{\mathit{\chi }}{\int }}{p}_{\mathrm{I}}\left(t^{\prime }\right)\cdot A(t-t^{\prime })\mathrm{d}{t}^{\prime },\text{(2)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{c}}={\mathit{\vartheta }}^{\left(\mathrm{1}\right)}{\left(\mathit{\zeta }{\mathit{\rho }}_s\right)}^{\mathrm{0.4}}a{U}_{\text{imp}}^{-\mathrm{0.2}},\end{array}$$

where χ is the time interval of convolution, with χ=t, if 0 $\le \mathit{\tau }\le T_{\mathrm{c}}$, and χ=T_c, if τ>T_c, with τ a delayed time due to sphere geometry, $\mathit{\tau }=t-(r-a)/c$, r is the distance between the observation point and the impact (see also Fig. 1a and b), a is the radius of the sphere, c is the sound celerity, ρ_s is material density and U_imp is the impact velocity. The parameter ϑ⁽¹⁾ is a constant, ϑ⁽¹⁾=9.229, for the impact between two spheres of the same radii and ϑ⁽¹⁾=10.601, for the impact of a slab and a sphere. The parameter $\mathit{\zeta }=(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})/\left(\mathit{\pi }{E}_{\text{long}}\right)$ is a material parameter, and it contains Young's modulus (E_long) and the Poisson ratio (ν).

Figure 1 (a) Setup for the impact between two spheres, here of the same radius a; the acoustic dipole source is illustratively depicted by the gray patch; (b) setup for the impact between a sphere of radius a and a semi-infinite rigid plane; to be noted is the boundary condition of the hard bottom (reflector) assumed in the framework of the “method of images”; thus, the impacting sphere is mirrored in the slab, so the acoustic fields are subtracted; the acoustic dipole source is illustratively depicted by the gray patch; (c) the elementary acoustic process of bedload noise in the river: the particle of equivalent diameter D=2a impacts the armored river bed (a massive slab) which generates a transient recorded by a hydrophone.

Figure 2 (a) Analytical waveform of sound from impact between a granite sphere of diameter D=20 mm and a granite slab, where the impact velocity U_imp=1 m s⁻¹, the directivity angle θ=0^∘ and the sensor is r=1 m from the impact; the arrow indicates the contact duration T_d; (b) the analytical spectrum modeled with Eq. (7) using the same parameters as in (a); the spectrum is an energy spectral density and is measured in Pa² s Hz⁻¹; (c) analytical spectra of sphere–slab impacts modeled by Eq. (7) as a function of diameter, $D=\mathit{\{}\mathrm{1},\mathrm{10},\mathrm{20},\mathrm{30}\mathrm{\dots }\mathrm{150}\mathit{\}}$ mm; impact velocity U_imp=1 m s⁻¹, the directivity angle θ=0^∘ and the sensor is r=1 m from the impact. (d) Peak frequency f_peak and power peak variations, from spectra modeled by Eq. (7), with the diameter and the sphere's diameters; the diameters are coded by colors. The power law f_peak=aD^b is given, where the sphere–slab impact tests consider three impact velocities ($U_{\text{imp}}=\mathit{\{}\mathrm{0.01},\mathrm{0.1},\mathrm{1}\mathit{\}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$) and the law of sphere–sphere impact is underlined by the dotted line; the material is granite. From bottom to top, the regression laws of sphere–slab impact vary from U_imp=0.01 (bottom) to U_imp=1 (top) m s⁻¹. The sphere–sphere impact tests are done using U_imp=1 m s⁻¹ and the same other parameters as with sphere–slab impacts. (e) Detail where the two vertical dotted lines locate the f_peak of the impact spectrum from 150 mm diameter particles for both impact models.

Table 1 Parameters used to model analytical spectra of sediment size mixtures, Eq. (6b), and the typical values adapted for the underwater environment. The typical singular values are used in inversion further in this paper. The ranges of values are used in the global sensitivity analysis. K is the number of size classes used to inverse acoustic spectra.

Table 2 First-order sensitivity indices S_i computed by the FAST method, assuming the peak frequency as the output of the model, f_peak: 0 % means no influence, and 100 % means total influence on the model output.

Figure 3 GSA flowchart to compute the first-order sensitivity indices by the FAST method; the spectrum is simulated with Eqs. (7) and (6b), with a log-normal distribution generated using a diameter in the range 1 to 150 mm and SD σ in the range 0.01 to 10. The rest of the input parameters are defined in Table 1. From the simulated spectra, the f_peak are computed, and finally the first-order sensitivity indices S_i are calculated using the FAST method. The results are shown in Table 2.

Figure 4 (a) Simulated PSD from the uniform PMF γ_m of sediments, 10 kg per 1 mm size class, from 10 to 50 mm, where the impact velocity is U_imp=1 m s⁻¹; the other input parameters are defined in Table 1; (b) the PMF solution obtained by the classical LS inversion, Eq. (12). The parameters used to simulate the PSD (grain size and impact velocity) shown in (a) are exactly the same as those used in modeling the dictionary Δ; (c) the PMF solutions obtained from the inversion of the spectrum shown in (a) using the two algebraic methods: the classical LS and the NNLS algorithm. The impact velocity used in modeling is U_imp=0.1 m s⁻¹, whereas in the simulation it is U_imp=1 m s⁻¹ (the other input parameters remain the same as in the simulation); the solution γ_m is post-processed by smoothing with a Gaussian moving window of 5 mm; (d) the simulated PSD shown in (a) with added variance (see text for the noise simulation procedure); (e) the cumulative GSD obtained from the inversion of the noised spectrum by the NNNLS algorithm. The estimated solution γ is used to reconstruct the spectrum, shown in green solid line in (d).

Figure 5 Experimental setup: (a) Isère River basin geographical location (http://histgeo.ac-aix-marseille.fr) and Google Earth© picture showing the river morphology near the bridge where measurements were taken; (b) instruments used during the trials; from left to right: Toutle TR sampler and the floating river-board with hydrophone; (c) the bridge from where acoustic drifts and sediment physical samplings were realized.

The general form of the acceleration profile is provided by Goldsmith (2003) and it is rewritten in a unified form for both sphere–sphere and sphere–slab impact models; see Eq. (1).

$$\begin{array}{}\text{(3)}& {\displaystyle }{\displaystyle }A\left(t\right)=\left\{\begin{array}{ll}{\mathit{\vartheta }}^{\left(\mathrm{2}\right)}\frac{U_{\text{imp}}}{T_{\mathrm{c}}}\sin \left(\mathit{\pi }\frac{t}{T_{\mathrm{c}}}\right),& \mathrm{0}\le t\le T_{\mathrm{c}}\\ \mathrm{0},& t>T_{\mathrm{c}}\end{array}\right.\end{array}$$

where the constant ϑ⁽²⁾=1.5708, for sphere–sphere impact, and ϑ⁽²⁾=3.353 for sphere–slab impact.

The first important observation from Eq. (1) is the half-period sinusoidal form of the Hertzian acceleration. The two modeled acceleration laws show close frequencies as the constant ϑ⁽¹⁾ from T_c's formula is not dramatically different from one case to another. If the frequency of acceleration of sphere–sphere impact is 1000 Hz, then the frequency of acceleration of sphere–slab impact is 909 Hz, which is almost only 10 % of the deviation. The maximum amplitude of acceleration for the impact between two spheres of radius a is almost 2 times less than the impact between a sphere of radius a and a slab, considering the same U_imp and T_c.

The integral convolution in Eq. (1) is transformed into multiplication in the complex Fourier space. Thus, the analytical magnitude spectrum of the noise from the rigid body acceleration, F_acc, is given in Eq. (4a).

$$\begin{array}{}\text{(4a)}& {\displaystyle }{\displaystyle }{F}_{\text{acc}}\left(\mathit{\omega }\right)=F\left(p_{\mathrm{I}}\right(t\left)\right)\cdot F\left(A\right(t\left)\right),\end{array}$$

where F(p_I) is the Fourier transform (FT) of Kirchhoff's impulse response p_I (Koss and Alfredson, 1973), for a sphere of radius a, defined in Eq. (4b), F(A) the FT of Hertzian acceleration due to elastic impact between two identical radii and the same material spheres, defined in Eq. (4c), and ω is the angular frequency which is a measure of rotation rate, in radians per seconds, and it is equal to 2πf; f is the linear frequency, a measure of number of occurrences per second.

$$\begin{array}{}\text{(4b)}& {\displaystyle }& {\displaystyle }F\left(p_{\mathrm{I}}\right)={\displaystyle \frac{\mathit{\rho }{a}^{\mathrm{3}}c}{r^{\mathrm{2}}}}\left(c+j\mathit{\omega }r\right){\displaystyle \frac{\mathrm{2}{c}^{\mathrm{2}}-(\mathit{\omega }a{)}^{\mathrm{2}}-j\mathrm{2}\mathit{\omega }ac}{{\left[\mathrm{2}{c}^{\mathrm{2}}-(\mathit{\omega }a{)}^{\mathrm{2}}\right]}^{\mathrm{2}}+{\left(\mathrm{2}\mathit{\omega }ac\right)}^{\mathrm{2}}}}\cos \mathit{\theta }\text{(4c)}& {\displaystyle }& {\displaystyle }F\left(A\right)={\mathit{\vartheta }}^{\left(\mathrm{3}\right)}{\displaystyle \frac{U}{{\mathit{\pi }}^{\mathrm{2}}-{\left(\mathit{\omega }{T}_{\mathrm{c}}\right)}^{\mathrm{2}}}}\left(e^{-j\mathit{\omega }{T}_{\mathrm{c}}}+\mathrm{1}\right)\end{array}$$

j is an imaginary unit and ${\mathit{\vartheta }}^{\left(\mathrm{3}\right)}=\pm {\mathit{\pi }}^{\mathrm{2}}/\mathrm{2}$ for sphere–sphere impact, and ϑ⁽³⁾=1.067π² for sphere–slab impact.

As we know, the nature of Hertzian sound is the oscillation of a rigid solid along a particular direction, and so each of the two objects in collision will give rise to a dipole acoustical source, shown in Fig. 1a. The case of sphere–slab impact is treated below. Hence, the amplitude of an oscillating sphere is dependent on the cos θ term and the phase of the acoustic pressure field changes by 180^∘ at θ=90^∘; i.e., the rarefaction wave changes into a compression wave or vice versa. In the case of the sphere–slab impact shown in Fig. 1b, the total pressure field is modeled as the addition between the compression wave and the slab-reflected rarefaction wave of the acoustic dipole. Thus, the addition becomes a subtraction as the reflected rarefaction wave keeps its sign (does not shift in phase), so there are two waves (compression and rarefaction) arriving at the sensor almost at the same time (Akay and Hodgson, 1978). This acoustic process is modeled by the so-called method of images by which one considers a mirrored sphere replacing the slab and being responsible for the rarefaction wave generation.

The same subtraction is applied in the case of complex spectra to obtain the total spectrum F_im, Eq. (5). In this formula, the first term of the right member is attributed to the impacting sphere, whereas the second term pertains to the mirror. The time delay T_d of sound arrival due to distance of measurement and the sphere's geometry make the two terms not perfectly cancel out or not arrive at the same time at the sensor.

$$\begin{array}{}\text{(5)}& {\displaystyle }{F}_{\text{im}}\left(\mathit{\omega }\right)=F_{\text{acc}}\left(\mathit{\omega }\right)-F_{\text{acc}}\left(\mathit{\omega }\right)\cdot e^{-j\mathit{\omega }{T}_{\mathrm{d}}}\end{array}$$

Introducing Eqs. (1), (1) and (4a–c) into Eq. (5), one obtains the complex magnitude spectrum of the impact between a sphere and a slab. The spectrum contains complex numbers, so one applies the multiplication of the spectrum and its conjugate to compute the magnitudes of the energy spectrum, Eq. (6).

$$\begin{array}{}\text{(6)}& {\displaystyle }|F_{\text{im}}{|}^{\mathrm{2}}=F_{\text{im}}\cdot F_{\text{im}}^{\ast },\end{array}$$

where $F_{\text{im}}^{\ast }$ is the complex conjugate of F_im.

The quantity $|F_{\text{im}}{|}^{\mathrm{2}}$ from Eq. (6) is noted as in Eq. (7) by E and its unit of measurement is Pa² s Hz⁻¹. Thus, the analytical model of impact used in this paper is an energy spectral density and it will be used to inverse acoustic spectra measured in the field.

$$\begin{array}{}\text{(7)}& {\displaystyle }\mathit{E}\left(\mathit{\omega }\right)=|F_{\text{im}}{|}^{\mathrm{2}}\end{array}$$

An example of an analytical model computed in time by Akay and Hodgson (1978) and reformulated in Appendix B is presented in Fig. 2a. The impacting sphere is 20 mm in diameter, the material is granite and the impacting velocity is 1 m s⁻¹. The shape of the waveform is an approximately one and a half period sinusoid. The subtraction of the two pressure fields, the rarefaction and the reflected compression wave fields are observed. It is also important to notice that the first arrival at the sensor is the compression wave. Thereafter, the other part of the acoustic dipole (the rarefaction wave) arrives with delay T_d at the sensor. The power spectrum density modeled by Eq. (7) is shown in Fig. 2b. Here, the spectrum has a principal lobe and numerous side lobes. The principal lobe has a peak at the frequency of approximatively $\mathrm{1}/(\mathrm{1.1}\cdot T_{\mathrm{c}})$ and the side lobes are approximatively associated with the term cos (ωT_c) also observed by Thorne and Foden (1988).

In Fig. 2c it is shown that the frequency peaks of spectra from both types of impact model decrease with the sphere's diameters (from 1 to 150 mm), as experimentally observed by Thorne (1986b). Frequency peak as a function of diameter D, in the case of sphere–slab impact, $f_{\text{peak}}\left(D\right)=a\cdot D^b$, is given in the case of three impact velocities, $U_{\text{imp}}=\mathit{\{}\mathrm{0.01};\mathrm{0.1};\mathrm{1}\mathit{\}}\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$. The exponents of the regression laws prove the exact inverse proportionality between f_peak and D. Besides, the power peaks and peak frequencies increase, for a certain diameter, when the impact velocity increases. There is only a doubling of f_peak when U_imp changes by an order of magnitude. This is also proved by the formula of Eq. (1) of T_c (almost the reciprocal of f_peak) where the parameter U_imp is raised to a weak exponent of −0.2.

The f_peak in the case of sphere–sphere impact, modeled for impact velocity U_imp=1 m s⁻¹, is higher than in the case of the sphere–slab impact. Here, the analytical model of the sphere–sphere spectrum was computed using Eqs. (4a) and (5), with the two pressure fields auditioned instead subtracted. This gives the same results as the spectral model reported by Thorne and Foden (1988). To give an idea, a 150 mm diameter particle in sphere–sphere impact has a spectrum of f_peak=1700 Hz (see detail of Fig. 2c), whereas sphere–slab impact has f_peak=1500 Hz, so the 200 Hz represents cca. 15 % of variation between the cases.

It is worth mentioning that U_imp greatly influences the power peak, if the former is changed by 1 order of magnitude. On the other hand, the power peak of the sphere–sphere impact is slightly weaker than the sphere–slab impact. In this paper, we choose to use a slab model to model bedload SGN as it simplifies the inverse problem. Indeed, the task of determining the dimensions of impacted particles is skipped. Therefore, we consider that the riverbed could be modeled as a slab. This hypothesis could be supported when the riverbed is armored or paved, but may be false when the riverbed is totally mobile and when the impacts between particles of different diameters are very common.

2.2 PSD model of the SGN generated by a mixture of sediments

In the previous section, the analytic energy spectral density (ESD) was defined for the impact between a sphere and a slab. In this section, we model the power spectral density (PSD) of a sediment mixture using these analytic ESDs and the impact rate of each class of diameter, or the number of impacts per second. Assuming that particle collisions are random and independent noise sources, the model of the PSD of a mixture, denoted by P, can be expressed as a linear summation of the elementary ESD, denoted by E_i (Johnson and Muir, 1969; Jonys, 1976; Thorne, 2014) weighted by the impact rate I_i. The acoustic bedload model under discussion is defined in the scalar form in Eq. (6) and the matrix form in Eq. (6a).

$$\begin{array}{ll}\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{P}=\sum _{i=\mathrm{1}}^K{\mathit{E}}_i\cdot I_i{\displaystyle }& {\displaystyle }\underset{\mathit{P}}{\underbrace{\left(\begin{array}{c}{P}_{\mathrm{1}}\\ P_{\mathrm{2}}\\ \mathrm{⋮}\\ P_{N_{\text{FFT}}}\end{array}\right)}}=\underset{\mathbf{\Delta }}{\underbrace{\left(\begin{array}{cccc}\stackrel{{\mathit{E}}_{\mathrm{1}}}{\overbrace{E_{\mathrm{1},\mathrm{1}}}}& \stackrel{E_{\mathrm{2}}}{\overbrace{E_{\mathrm{1},\mathrm{2}}}}& \mathrm{\dots }& \stackrel{E_K}{\overbrace{E_{\mathrm{1},K}}}\\ E_{\mathrm{2},\mathrm{1}}& E_{\mathrm{2},\mathrm{2}}& \mathrm{\cdots }& E_{\mathrm{2},K}\\ \mathrm{⋮}& \mathrm{⋮}& \mathrm{\ddots }& \mathrm{⋮}\\ E_{N_{\text{FFT}},\mathrm{1}}& E_{N_{\text{FFT}},\mathrm{2}}& \mathrm{\dots }& E_{N_{\text{FFT}},K}\end{array}\right)}}\cdot \underset{\mathit{I}}{\underbrace{\left(\begin{array}{c}{I}_{\mathrm{1}}\\ I_{\mathrm{2}}\\ \mathrm{⋮}\\ I_K\end{array}\right)}}\end{array}$$

$$\begin{array}{}\text{(9a)}& \Rightarrow \mathit{P}=\mathbf{\Delta }\cdot \mathit{I}\end{array}$$

where Δ is the dictionary of elementary ESD of impacts between spheres and slab and I is the vector of impact rates per diameter class or, basically, a histogram. Class i takes integer values, from the lowest limit, 1 mm, to the highest one, K mm, where K is the largest diameter considered in modeling. Here, we consider K equal to 150 mm. The parameter NFFT is the number of values contained in the spectrum or the number of Fourier transform points on which the spectrum is modeled.

The histogram I can be transformed in the probability mass function γ by normalizing it by its sum of elements. The cumulative form of γ will be denoted by Γ. Thus, one of the main assumptions is that Eq. (6a) can be written in terms of probabilities γ, as in Eq. (6b):

$$\begin{array}{}\text{(9b)}& {\displaystyle }{\displaystyle }\mathit{P}=\mathbf{\Delta }\cdot \mathit{\gamma },\end{array}$$

where γ_i represents the probability of having a number of impacts of particles per second for size class i.

Therefore, the random variable here is I and γ is the probability of impacts, and so the quantity γ (I=I_i) is a discrete probability, given that we operate on size classes of 1 mm in diameter. This probability is computed from a histogram of the number of impacts per second, so one needs to transform it into a histogram in mass of sediments M, to be compatible with the measured GSD by physical sampling. In consequence, γ (I=I_i) will be scaled by $D_i^{\mathrm{3}}$, as in Eq. (7), in order to obtain γ_m(M=m_i). Finally, the grain size distribution (GSD), or the cumulative distribution form of γ_m, will be Γ_m(M≤m_i), expressing the probability of sediments finer than D_i, as defined in Eq. (7).

$$\begin{array}{}\text{(10)}& {\displaystyle }& {\displaystyle }{\mathit{\gamma }}_m(M=m_i)=\mathit{\kappa }\cdot \mathit{\gamma }(\mathit{I}=I_i)\cdot D_i^{\mathrm{3}},i=\mathrm{1},\mathrm{\dots },K\text{(11)}& {\displaystyle }& {\displaystyle }{\mathbf{\Gamma }}_{m_j}(M\le m_i)=\sum _{i=\mathrm{1}}^j{\mathit{\gamma }}_{m_i}(M=m_i),j=\mathrm{1},\mathrm{\dots },K\end{array}$$

where D_i is the diameter in meters and κ is a constant which is the rest of the mass-volume formula coefficient including the material density.

This section gave a formal definition of the PSD of a bedload size mixture defined by its GSD. The proportions are considered to be a probability mass function (PMF) of the rate of impacts. The size classes concerned in this study are integer numbers, from 1 to K, with a resolution of 1 mm per size class.

2.3 Global sensitivity analysis of the spectrum generated by a mixture of sediments

This analysis was done to determine the importance of input parameters for the shape of the PSD modeled with Eq. (6b). The parameters are defined in Table 1. Global sensitivity analysis (GSA) is done to assess the impact of input parameters on the model output, which in our case is the peak frequency f_peak of spectra modeled by Eq.(6b). We use the Fourier amplitude sensitivity test (FAST) (Cukier et al., 1973) to compute the first order indices of sensitivity S_i for each input parameter. The coded version of the FAST algorithm is presented in Cannavó (2012).

The flowchart of the GSA is presented in Fig. 3. The FAST analysis uses the typical range of parameters found in rivers, defined in the second column of Table 1. The input log-normal distributions (GSD) have median diameter values in the range from 1 to 150 mm and SDs σ from 0.01 to 10. All other input parameters needed for the model in Eq. (6b) are given in Table 1. As the output model analyzed is the peak frequency f_peak of the simulated PSD curves, the analysis does not claim to completely describe the model, but pertinent ideas could be drawn on the model's behavior.

The results in terms of first-order indices are presented in Table 2. The SD σ of log-normal GSD has the greatest influence on the PSD shape. This is because σ affects the values of all percentiles of the GSD curve. The median D₅₀ is almost 2 times less important than σ. The third-greatest parameter as a degree of influence on output is surprisingly Young's modulus, but this is due to a very wide range of values (here $E_{\text{long}}={\mathrm{10}}^{\mathrm{10}}\mathrm{\dots }\mathrm{7}\times {\mathrm{10}}^{\mathrm{10}}$ Pa, i.e., from quartz to granite materials). Such variation is not possible at the reach scale, where the sediments are of the same material. The impact velocity U_imp comes shortly after Young's modulus and confirms the conclusion of the previous local analysis, that U_imp has ca. 10 % of importance in the output. Other relatively important parameters are Poisson's ratio and the density of sediments, which means that the type of material also plays a role in the dynamics of the f_peak. The distance of measurement, r, also plays also a role in the f_peak variation. The angle of the point of observation with respect to the impact, θ, and the propagation medium properties, ρ and c, are considered of little influence on the values of f_peak.

In conclusion, the first-order global sensitivity analysis of the peak frequency shows a comprehensive view of its dynamics with input parameter variation. It is found that the peak frequency is mainly affected by two parameters, the distribution's SD and the median diameter, together making up ca. 65 % of output variation, whereas the material properties (i.e., density, Poisson's ratio, Young's modulus) have almost ca. 20 %, and the impact velocity U_imp has ca. 10 %. In conclusion, the acoustic model is quite complex and care must be taken regarding the recording of the power spectra on the field, as their shape heavily affects the estimation of GSD. Also, the impact velocity is regarded as a minor factor of uncertainty and, because it is almost impossible to be measured for each grain size class, the 10 % uncertainty in peak frequency is almost unavoidable. The material properties should not be a problem with the condition that the sediments are the same. For a complete GSA, the computation of high-order sensitivity indices can be made using Sobol's methodology (Sobol, 2001), but this type of analysis is beyond the scope of this article.

2.4 Assumptions on the proposed SGN spectrum model

Modeling of single impacts requires definition of parameters typical for river environments, in Table 1. Using the global sensitivity model, it has been shown that PSD shapes are essentially influenced by four parameters: the shape of the GSD curve, the median diameters of the colliding particles, the impact velocities and the material. Grain sizes are estimated later using the inversion algorithm presented in Sect. 3. Concerning the other model parameters, as they do not affect the PSD shape, they will be fixed for the inversion process, using realistic values. These parameters are listed in the third column of Table 1. The main assumptions of the SGN spectrum model are the following.

• i.

  The geometry of the channel and of the material: the river bed is considered a massive slab and moving particles are considered spherical.

• ii.

  Sediment transport assumptions: impact velocities are assumed to be invariant with grain size. This assumption is supported by the relative size effects on bedload transport (Einstein, 1950; Recking, 2016; Wilcock and McArdell, 1993) referring to the mobility of finer and coarser particles.

• iii.

  Acoustic propagation:

  • as the bedload GSD is assumed to be homogeneous everywhere in the space, the propagation effects like the attenuation with distance (geometrical spreading models) will not impact the spectrum shape;

  • the attenuation due to diffraction from bed and water surface roughness or from the suspended sediments is not considered. The issue of the nonlinear propagation will be detailed in the discussion part of this paper.

3 Inverse model to estimate the GSD of bedload particles

The inversion uses least square (LS) optimization methods to compute the inverse of dictionary Δ. Normally K<N_FFT, so Δ is a non-square matrix. Moreover, the matrix Δ is possibly rank deficient because the spectra generated by impacts of coarser particle sizes show very similar shapes, that is, the coarser the particle, the more similar the produced sound. This is also shown in Fig. 2b, where one could observe that the points on the graph agglomerate as the grain diameters increase. In this case, the pseudo-inverse algorithm is used to solve the algebraic system of Eq. (6b). The optimization problem is defined as in Eq. (12). The least square solution to this problem is the PMF of the rate of impacts γ. The estimated PMF is further transformed into the final GSD of mass of sediments according to Eqs. (7)–(7).

$$\begin{array}{}\text{(12)}& {\displaystyle }\widehat{\mathit{\gamma }}(I=I_i)=\text{ minimize}({\mathbf{\Delta }}^{+}\cdot \mathit{P}-\mathit{\gamma })\end{array}$$

where ${\mathbf{\Delta }}^{+}={\left({\mathbf{\Delta }}^t\cdot \mathbf{\Delta }\right)}^{-\mathrm{1}}\cdot {\mathbf{\Delta }}^t$ is the pseudo-inverse, and Δ^t means the transpose of the matrix Δ.

Eq. (12) conveys the idea of minimizing the error between the model and the measurement. This minimization operation is realized in the sense of the least square optimization.

3.1 Numerical test of the LS method

A simulation case is proposed here to test the robustness of the LS inverse method. The simulated PMF, or grain size distribution, γ_m is uniformly distributed between 10 and 50 mm. The uniform distribution means that 1 kg of D=10 mm has the same probability of producing impact noise as 1 kg of D=11 mm, and so on. To obtain γ, the simulated PMF γ_m is converted back to impact rates by dividing by D³, D in meters. Using an impact velocity of 1 m s⁻¹ and the rest of the input parameters defined in Table 1, the simulated PSD P is shown in Fig. 4a. Here, the dictionary Δ contains spectra from 1 to 150 mm and the grain size distribution has 1 mm resolution. Applying Eq. (12) on the simulated spectrum and considering exactly the same parameters in modeling and in simulation, it is found that the estimated γ_m is exactly the same as the simulated γ_m, as is expected; see Fig. 4b.

However, if the impact velocity used in modeling the dictionary is set to a value (U_imp=0.1 m s⁻¹) which is different than the one used in simulation (U_imp=1 m s⁻¹), then high instabilities are observed on the estimated γ_m; see Fig. 4c. This is explained by the fact that there is a high similarity between the elementary spectra E, especially for the larger size classes. Thus, the matrix Δ is ill-conditioned and the problem is ill-posed. Ill-conditioning is linked to the high condition number of the normal matrix (Δ^t.Δ). It is defined as the ratio between the largest and smallest eigenvalues of a matrix. A well-conditioned algebraic system requires that the normal matrix should have a condition number as close as possible to 1 (Strang, 2006). In these tests, Δ's condition number reaches huge values on the order of 10¹²–10²⁰. In consequence, the similar spectra from the matrix Δ produce high instability in solution.

To avoid the instability in the LS solution, the non-negative least squares (NNLS) algorithm (Lawson and Hanson, 1974) is proposed to solve the LS problem. This optimization algorithm, Eq. (13), casts non-negative constraints on solution γ. The non-negative factorization is widely used, for example, in various domains like image processing or chemometrics. The side-effect of using this algorithm is the strong regularization of solutions. The regularization aims to keep the sum of components in γ constant. The solution of the NNLS algorithm (see Fig. 4c) shows that the instabilities are completely removed off. Besides, it is important to note that the estimated diameters are inside the simulated interval of diameters.

$$\begin{array}{}\text{(13)}& {\displaystyle }\widehat{\mathit{\gamma }}(I=I_i)=\text{ minimize}({\mathbf{\Delta }}^{+}\cdot \mathit{P}-\mathit{\gamma }),\widehat{\mathit{\gamma }}\ge \mathrm{0}\end{array}$$

3.2 Robustness of the NNLS algorithm to PSD noise

The signal processing tools in this paper refer to using the PSD as the method of spectral representation of bedload signal. The use of PSD is worthwhile because the type of bedload signal is a stationary random one. Random stationary signals are signals varying in time but whose average and SD of amplitude values over some fixed periods are constant.

A particular concern for the signal processing of random processes is the minimization of the variance on the PSD. This work makes use of the periodogram algorithm for the PSD estimation, which means the Fourier transform is applied on local portions (windows) of the random signal, with an overlap of 50 %, and then the local results are averaged in narrow bandwidths (Oppenheim and Verghese, 2010). The averaging is useful because it mitigates the variance on the PSD. In this work, the quality of spectra is vital for accuracy of estimations. The uncertainty principle tells us that the smaller the temporal window, the greater the uncertainty in locating two very close frequencies on the spectrum, so a trade must be made between the PSD variance and its spectral resolution. If the bedload signal is too short, the quality of spectra toward the low-frequency bands is worsened because in one single bandwidth of the Fourier transform there is spectral information of impacts from multiples grain sizes. Finally, the longer the signal the better the spectral resolution and the less the variance on the PSD curve.

The NNLS algorithm will be tested on three simulated spectra which have different degrees of variance. The simulated γ used is identical to that in Sect. 3.1. The simulated signal is obtained by convolving a realization of a white noise with a transfer function being the modeled spectrum shown in Fig. 4a. The simulated noised PSD is shown in Fig. 4d. The results of the inversion using the NNLS algorithm show that, even for the worst scenario of variance on a spectrum, the inversion method correctly reconstructs the simulated GSD; see Fig. 4e.

Finally, we conclude that the NNLS algorithm is robust with respect to PSD noise and fits to this kind of inversion problem. The inversion procedure will now be tested on in situ measurements.

4 Application to real data

4.1 Isère River and experimental setup

The Isère River is a piedmont gravel river bed located in southeastern France, and it is one of the main tributaries of the Rhône River, which reaches the Mediterranean Sea. The monitoring section is located in the city of Grenoble (45^∘11^′52.8^′′ N, 5^∘46^′14.88^′′ E); see Fig. 5a. In this reach, the mean slope is about 0.06 %, the area of the watershed is 5500 km² and the annual average flow rate is 180 m³ s⁻¹. At the time of experiments, 29–30 June 2016, the monitored discharge was on average 300 m³ s⁻¹. The measurement section has a rifle-pool morphology with riprap-protected embankments. Two different types of instrument were used: SGN measurements using hydrophone and direct sampling using a pressure-difference sampler, shown in Fig. 5b. All these measurements were carried out from a suspension bridge (Fig. 5c).

Figure 6 (a) Positions of the floating board during drift experiments, with essential positions marked on the bridge, X= {14, 35, 58} m across the river; (b) the PSD estimated from the 12 drifts, in units of Pa² Hz⁻¹; to be noted is the change in peak frequencies: the leftmost position (Drift no. 12) has the highest frequency, meaning that the finer size fractions are transported, and the particles are getting coarser up to the right bank; (c) the measured SPL map from the 12 drifts, in units of dB re 1 µPa; the maximum values are found in the middle of the Isère River's cross section.

Figure 7 Signal representations of the SGN recorded during hydrophone experiments on the Isère River (France): (a) temporal signal in units of Pa; (b) time–frequency representation (spectrogram), with the color code normalized with respect to power values, in Pa² Hz⁻¹; the specific frequency bandwidth of the bedload acoustic effects and of the hydrodynamic noise agitation (extraneous sources) are indicated; (c) the PSD curve, also in Pa² Hz⁻¹, estimated using either the average or the median power values, in time, from the spectrogram in (b).

4.1.1 SGN measurements

SGN measurements were made using a HTI99 hydrophone (High Tech, Inc., http://www.hightechincusa.com/) with a sensibility of −160 dB re 1 V µPa⁻¹ ± 3 dB from 10 Hz to 125 kHz. The hydrophone was connected to an autonomous-waterproof autonomous recorder – SDA14 (RTSYS©, http://www.rtsys.eu). The gain of the recorder was set to 15 dB. Signals were sampled at a 312 kHz frequency with a resolution of 24 bits and saved as wav files. The scope of these field experiments was to trace maps of the SGN on the local reach. The hydrophone and the recorder were attached to a free-floating river-board. The hydrophone position was about 1 m below the water surface and 1.5 m on average above the river bed. The SGN map consists in launching 12 drift measurements from the bridge which are located due to a GPS device connected to the acoustic recorder. Each drift consists of recordings of about 30 to 40 s, or in terms of distance, between 50 and 100 m. The river-board positions during the drifts are shown in Fig. 6a. The recorded signals were processed to compute acoustic spectra. The 12 acoustic spectra recorded across the river (Fig. 6b) are inversed to estimate the bedload. The river cross section is about 60 m. Also, the 12 drift measurements are synchronized with GPS data to compute the SGN map in terms of sound pressure level (SPL), as is shown in Fig. 6c. The variability of SGN noise from the left to right banks can be observed from both the spectra and SPL map.

4.1.2 Definition of the SGN spectrum

SGN signals are measurements of bedload transport noise propagating in the river environment. Several representations of the acoustic signal are presented hereby, computed on the signal recorded in the middle of the Isère River (X=34 m): (a) the temporal waveform, in Fig. 7a; (b) the spectrogram, in Fig. 7b, as the scaled squared magnitude of short-time Fourier transform, in Pa² Hz⁻¹; and (c) the PSD, also expressed in Pa² Hz⁻¹, computed by either averaging or medianizing the PSD spectrogram, in Fig. 7c. Two main sources of noise can be distinguished in the recordings: below and above 400 Hz (Fig. 7). Bedload impacts can clearly be heard in the higher frequency band, sounding like the crackling of flames. Sounds occurring below 400 Hz are non-propagating sounds as they are localized below the cutoff frequency of the river waveguide (Geay et al., 2017b; Rigby et al., 2016). They are related to turbulence-induced noise around the sensor and to mechanical movements of the structure sharing the hydrophone. In the Isère River experiment, the SGN signal measured by drifts is almost free of hydrodynamic noise, which is proved by the typical median spectrum presented in Fig. 7c. In this study the inversion will be applied on such high signal-to-noise ratio PSD curves.

The median procedure is used to provide better smoothing as it filters more efficiently the unwanted low-frequency noises (Geay et al., 2017a). As in Fig. 7c, the suppression of the lower-frequency spikes can be noticed, attributed to the hydrodynamic noise, when median PSD is used instead of the average one.

4.1.3 Pressure-difference sampling

A Toutle River (TR) sampler, depicted in Fig. 5b, has been used to sample bedload particles (entrance width of 305 mm by 152 mm). There were two mesh sizes used for sampling: 0.2 and 1.3 mm. Sample durations were between 4 and 8 min. Finally, each bedload sample was dried, weighted and sieved in the laboratory. The sampled sediments were classified into six size classes: K= {<0.5; 0.5–2; 2–8; 8–16; 16–32; 32–64} mm. The TR sampler has been deployed in three cross-sectional positions (at X=27 m, X=35 m and X=44 m, marked on the bridge from the left to right river banks). The number of repetitions for each cross-sectional position is indicated in Table 3. Bedload fluxes ($\mathrm{g}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$) have been averaged for each position of the sampler. GSDs have been computed for each position and for each mesh size used.

Figure 8 Measured bedload flux in three positions across the Isère River, $X=\mathit{\{}\mathrm{27},\mathrm{35},\mathrm{44}\mathit{\}}\mathrm{m}$ and (b) measured GSD curves in these positions, using the TR sampler with two mesh sizes, 0.2 and 1.3 mm.

Figure 9 (a) Estimated GSD by the NNLS algorithm in the center of the Isère River (X=34 m), using different values of impact velocities U_imp= {0.01, 0.1, 1, 5} m s⁻¹. Measured GSD by TR sampler (X=35 m) is represented by the yellow envelope for the two mesh sizes (see Fig. 9b for fraction sizes finer than 1 mm); (b) The D₁₆, D₅₀, and D₈₄ estimated by NNLS across the Isère River compared to the regression laws of Thorne (1985, 1986b) for estimating the equivalent diameter D_eq; we also indicate by arrows the range of D₅₀ measured by the Toutle River sampler (in positions $X=\mathit{\{}\mathrm{27},\mathrm{35},\mathrm{44}\mathit{\}}\mathrm{m}$). The impact velocity used in the inversion is U_imp=1 m s⁻¹.

Table 3 Number of repetitions for each measurement.

4.2 Results

4.2.1 Direct measurements of bedload

Results of TR sampler measurements are shown in Fig. 8a–b. A maximum of bedload flux was found in the middle of the cross section (X=35 m) (Fig. 8a). A value of 100 $\mathrm{g}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$ has been measured. In side positions, the flux was found to be 5 times smaller, around 20 $\mathrm{g}{\mathrm{s}}^{-\mathrm{1}}{\mathrm{m}}^{-\mathrm{1}}$. Concerning grain size distributions, most of the measurements indicate a D₅₀ between 7 and 20 mm. Notice that measurements made with the 0.2 mm mesh size towards the left bank (X=27 m) indicate a GSD toward much finer sediments (D₅₀ of about 0.3 mm) (Fig. 8b). Bedload samples closest to the left bank were indeed constituted of huge amounts of fine sediment mixed with vegetable debris (about 60 % of the total mass sampled). In the central and right positions, neither vegetable debris nor silts were sampled. TR sampler measurements showed grain size sorting along the river cross section, varying from silts, near the left bank, to gravel, near the right bank.

In the following, the GSD measured in the central position (X=35 m) will be considered. Its flux was indeed the largest measured, and it is considered to be the principal source of bedload noise throughout the river.

4.2.2 SGN spectra inversion

All the median PSDs of SGN signals recorded across the Isère River have been presented in Fig. 6b. The seventh drift will be studied, the one positioned in the center of the cross section at X=34 m, which is the closest to the middle position of TR sampling measurements. In this position, it can be observed that a maximum bedload acoustic energy has been recorded. Additionally, a maximum flux of sediments was sampled in this position. The results of spectrum inversion, using a modeled dictionary Δ with size classes from 1 mm to K=100 mm (100 size classes), are shown in Fig. 9a. The results are compared to the GSD measured by the TR sampler in position X=35 m. Four different values of the impact velocity U_imp are tested (from 0.01 to 5 m s⁻¹) and it is noticed that the impact velocity U_imp between 0.01 and 0.1 m s⁻¹ leads to a very good match between estimation and TR sampling measurements, except for the very small size classes from 1 to 5 mm. The value of impact velocity U_imp=0.1 m s⁻¹ will be used in the inversion of all other spectra measured across the Isère River.

Secondly, the GSD variations, represented by the percentiles D₁₆, D₅₀ and D₈₄, are estimated by the inversion of 12 drift measurements taken across the Isère River. The model uses the impact velocity of 0.1 m s⁻¹ and the rest of the parameters defined in Table 1. The estimated percentiles are compared to equivalent diameters D_eq computed by regression laws found by Thorne (1985, 1986b) and redefined below in Eq. (9) and Eq. (9). The equivalent diameter D_eq is a measure of particle size, and it is the diameter of the circle with the center as the centroid mass. The D_eq is computed using f_peak and, respectively, the centroid frequency f_centr. They are also compared to the TR sampler measurements, in the three positions across the Isère River.

$$\begin{array}{}\text{(14)}& {\displaystyle }& {\displaystyle }{f}_{\text{peak}}={\displaystyle \frac{\mathrm{224}}{D_{\text{eq}}^{\mathrm{0.9}}}},\text{(15)}& {\displaystyle }& {\displaystyle }{f}_{\text{centr}}={\displaystyle \frac{\mathrm{209}}{D_{\text{eq}}^{\mathrm{0.88}}}},\text{(16)}& {\displaystyle }& {\displaystyle }\underset{f_{\mathrm{1}}}{\overset{f_{\text{centr}}}{\int }}\mathit{P}\mathrm{d}f=\underset{f_{\text{centr}}}{\overset{f_{\mathrm{2}}}{\int }}\mathit{P}\mathrm{d}f,\end{array}$$

where P is the PSD and (f₁,f₂) is the frequency band defined by a value of 10 dB below the power peak. It is observed that the estimated D₅₀ by the NNLS algorithm is 10–14 mm, which is in the upper limit of the D₅₀ measured by the TR sampler (ca. 7 mm), in the middle of the river; X=35 m. On the one hand, the percentile D₁₆ almost matches the equivalent diameter D_eq estimated by Thorne's regression law f_centr(D_eq), Eq. (9), which is on average 50 % below the D₅₀ measured by the TR sampler. On the other hand, the percentile D₈₄ is closer to the equivalent diameter D_eq estimated by using the peak frequency regression law f_peak(D_eq), Eq. (9), overestimating the measurements of the TR sampler.

Figure 10 (a) Modeled spectra using a log-normal GSD, $D=\mathrm{1},\mathrm{2},\mathrm{\dots },\mathrm{150}\mathrm{mm}$, where D₈₄=2 D₅₀, D₅₀=10 mm (see medallion); typical input parameters are given in Table 1 and U_imp=1 m s⁻¹. Concerning the sphere–sphere impact, the impactor has the same size as the impactee; (b) inversion using the NNLS algorithm by acoustic spectra shown in (a).

5 Discussion on real data results

This work deals with the development of a novel estimation strategy of bedload GSD from acoustic PSD. The spectrum inversion used the model based on sphere–slab impact, where the impacting sphere diameters range from K=1 mm to K=100 mm. The inversion of field experiments on the Isère River have shown in Fig. 10a interesting results in conformity with the assumptions enounced in Sect. 2.4.

The inversion considered four values of impact velocity $U_{\text{imp}}=\mathit{\{}\mathrm{0.01};\mathrm{0.1};\mathrm{1};\mathrm{5}\mathit{\}}$m s⁻¹. The best fit to the measured GSD by the TR sampler is when the impact velocity U_imp is between 0.01 and 0.1 m s⁻¹, which could be possible for a large gravel bedded river like the Isère. To verify this, the apparent velocity of the bed material (see Rennie and Miller, 2004, for a definition) was measured by an aDcp at the moment of hydrophone experiments. This estimated value was at a maximum around 0.01–0.02 m s⁻¹, which can be in accordance with the impact velocity modeling the best NNLS estimates.

The cross-sectional variation of the estimated D₁₆, D₅₀ and D₈₄ by the NNLS algorithm follows the same trend of increasing values from the left to right banks as the bedload D₅₀ measured by the TR sampler (Fig. 9b). However, the cross-sectional variability of sampled diameters is higher than the estimated one. This is explained by the fact that the hydrophone has the spatial integrative characteristic (Geay et al., 2017b). The phenomenon of signal integration is typical for rivers like the Isère, where high fluxes of bedload transport are concentrated only in a small portion across the section, i.e., in its center. In this case spatial homogeneity as stated in Sect. 2.4 is no longer valid. However, the powerful acoustic source makes noise all over the cross section, causing the sound sources to appear ubiquitous. This may be the reason that the inversion of acoustic PSD measured in the center (X=34 m), for U_imp=0.1 m s⁻¹, still shows a good match to the sampling measurements in that position, only because of the high powerful acoustic source localized in this position.

Despite the consistent variation of the GSD across the river bed, measured by the sampler, the acoustic spectrum shapes shown in Fig. 6b are relatively stable, in the interval $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$ to $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{3}}$ Pa Hz⁻¹. This suggests that measurements by hydrophone installed from one of the banks are not dramatically different from measurements by free floating hydrophones along the watercourse.

The propagation of sound throughout the local reach also raises some concerns about the quality of measured acoustic spectra. The proposed model Eq. (6b) has been elaborated by assuming a simple geometrical spreading model of the acoustic waves in the river. Bedload SGN spectra monitored by a hydrophone are not only dependent on bedload sizes, but are also affected by propagation effects. For example, an alpine river has been modeled as a Pekeris waveguide (Geay et al., 2017b). Consequently, it has been shown that the monitored spectra were slightly dependent on the hydrophone position in the lower frequency band. Another propagation effect concerns the frequency cutoff phenomena, due to acoustic propagation in waveguides (Geay, 2013; Geay et al., 2017b; Jensen et al., 2011; Rigby et al., 2016). In our case, the Isère River has enough large depth that the bandwidth of bedload is not being impacted. The pebble-sized particles that are up to 64 mm give SGN of dominating frequencies well above 1000 Hz, whereas the channel's depth of 2.5 m fixes the cutoff frequency to about 148 Hz, assuming a perfectly rigid bottom. Therefore, the bandwidth of bedload is far superior to the frequency cutoff in the Isère River, so there are no risks of inversion. However, SGN monitoring and the inversion technique for GSD determination are particularly adapted to large rivers. Generally, propagation effects are frequency-dependent and higher frequency ranges are more affected by attenuation or scattering effects. A solution to the nonlinear effects of acoustic propagation would be to determine the river's transfer function by active acoustic experiments (Rigby et al., 2016) and to construct laws of attenuation that will compensate for the loss (Wren et al., 2015).

At first sight, our comparison with Thorne (1985, 1986b)'s regression laws would be very naïve due to the nature of theories: we considered the sphere–slab impact, whereas the regression laws are from sphere–sphere impact phenomena. Therefore, the inversion is put into discussion when the river bed is no longer armored, and so, the model of impact between sphere and slab is debatable. Here, we target the large gravel rivers. The dictionaries Δ for both impact models use an impact velocity U_imp=1 m s⁻¹, the material is granite and a GSD is simulated according to Recking's procedure (Recking, 2013), where D₈₄=2 D₅₀. When comparing the shapes of both simulated PSDs, shown in Fig. 10a, their respective frequency peaks f_peak are nearly identical. Likewise, the slopes of the spectra are found to be quite similar. Fig. 10b shows that the two solutions show no difference, except for a little disparity in the region of small grains. This proves that sphere–slab framework modeling of the collision between sediments and the river bed could work not only for stable conditions, but also for hydraulic events.

Another strong assumption used in modeling the PSD model of mixed impacts is that the particles are of spherical shapes. It is intuitively reasoned that the particle sphericity, shape factor and roundness also affect the acoustics of impacts. There are multiple possible ways of reckoning the equivalent diameter of a non-spherical particle. The particle's radius may be computed with respect to the curvature of the region of contact (see Chadwick et al., 2012; Goldsmith, 2003), with respect to the particle's mass centroid (Thorne, 1986b), which is in fact the a axis of the particle, or with respect to the b axis of the particle (Wyss et al., 2016b). Laboratory tests were conducted at the GIPSA laboratory, during which two pebbles of a size in the range 32–46 mm were impacted in a water pool along the three ellipsoid axes a, b, and c. The methodology of measuring the ellipsoid axes is found in Bunte and Abt (2001). It was found that the measured centroid frequencies take values from 3000 to 8000 Hz. If the regression law Eq. (9) is used, then the estimated diameters span the range from 23 to 73 mm, which is the repartition of all possible radii of curvature of the respective zones of contact. If the mode of sediment transport by sliding is the most frequent, then the particle c axis could be used to infer an equivalent diameter. If the rolling mode is more frequent, then the b axis would be more appropriate to work with. Finally, if the saltation is concerned, which makes the point of this work, then axes a and b are equally probable to be taken into account in modeling impacts.

6 Conclusion

A new strategy has been presented for data processing on hydrophone measurements for monitoring the bedload GSD in a gravel river bed. This strategy defines a forward model and a spectrum inversion approach. Firstly, the forward model combines generated spectra from collisions between a sphere and a slab. Secondly, the inversion procedure treats the forward model as a linear system of equations and uses algebraic methods of solving least square problems to obtain the GSD.

The forward model is based on a weighted sum of analytical energy spectral densities modeling the physics impact between a sphere and a slab. The weighting coefficients of the model represent a probability mass function which gives in the end the grain size distribution of bedload particles. The global sensitivity analysis of the PSD model of mixed impacts determined that the shape of the GSD has the biggest influence on the shape of the acoustic spectrum computed by Eq. (6b). Other important parameters are the median diameter and the impact velocity. However, the influences are from mixed interactions of parameters, and it is very hard, if not impossible, to obtain a complete analysis of the sensitivity of the analytical model of Eq. (6b).

The PSD model of mixed impacts works under the following strong assumptions: (1) the GSD is distributed everywhere in space and, in the same way, (2) the acoustic propagation is not frequency-dependent and, so, the spectrum shape is not affected by propagation in the river, (3) the impact velocity is invariant with the grain size, and (4) the impacting particles are of a spherical shape. The in situ experimentations showed that the integrative sound from all over the reach could render the first assumption verified (or true). In the case of the Isère River, the concentration of high transport rates in the middle of the cross section permits reliable measurements of bedload GSD by hydrophone from river banks.

The inversion method is a non-negative least square algorithm and it eliminates the negative solutions caused by ill-conditioned matrices. Concerning the least square approach for inversion, it is robust to noise.

The inversion of spectra from field trials on the Isère River proved that the method is highly reliable with no consideration of a priori information on bedform morphology of hydrological conditions. Surrogate methods for sediment transport in rivers were conceived with the idea of having access to information all over the reach and real time. In contrast to geophones and the Japanese pipe, the hydrophone technique does not require particular efforts to be installed in the watercourse.

Appendix A: Notations.

$a,b,c$	length of sediment (ellipses) axis	mm
A	Hertzian acceleration	m s−2
c	sound celerity in water	m s−1
Δ	modeled dictionary of individual energy spectra
Ei	Energy spectral density of the impact of the size class i, Eq. (7)
Ei,x	energy of collision in a narrow frequency bandwidth x	Pa2 s Hz−1
Elong	elastic modulus (Young's modulus) of rigid body	Pa
ESD	Energy spectral density	Pa2 s Hz−1
D	generic notation for the grain diameter	mm
Deq	equivalent diameter (with respect to the grain's mass center)	mm
Di	grain size for i from 1 to K	mm
DTR	grain size class measured by Toutle River TR sampler	mm
D16, D50, D84	the 16th, 50th and 84th percentiles of the grain size distribution	mm
f	linear frequency	Hz
fcentr	centroid frequency	Hz
fpeak	peak frequency	Hz
FAST	Fourier amplitude sensitivity test
FT	Fourier transform
Fs	Sampling frequency	Hz
F	Fourier transform operator
Fimp	linear complex magnitude spectrum of the elastic impact, Eq. (5)	Pa
$|F_{\text{imp}}{|}^{\mathrm{2}}$	energy spectral density of the elastic impact, Eq. (6)	Pa2 s Hz−1
GSA	Global sensitivity analysis
GSD	Grain size distribution
γ	solution of the inversion written as a probability mass function
γm	solution of the inversion written as a probability mass function, computed from the mass histogram of sediments (Eq. 10)
Γ	solution of inversion (GSD) in the cumulative form
Γm	solution of inversion (GSD) in the cumulative form, computed from ${\mathit{\gamma }}_m^{}$ (Eq. 11)

I	Histogram of rate of impacts
Ii	Rate of impact of the size class i	no. imp s−1
j	imaginary unit
K	number of grain sizes classes
LS	Least square problem
ν	Poisson's ratio of rigid body
NFFT	number of points for FT computation
NNLS	Non-negative least squares
ω	Angular frequency	rad s−1
P	Power spectrum density of the noise from an elastic impact, Eq. (9a–b)	Pa2 Hz−1
PMF	Probability mass function
PSD	Power spectral density	Pa2 Hz−1
r	reference measurement distance between the sensor and the center of the impact (see Fig. 1a and b)	m
ρs	density of sediment	kg m−3
ρ	density of water	kg m−3
SGN	Self-generated noise (noise generated by the transported sediments in collision)
Si	sensitivity indices from the first-order global sensitivity analysis
σ	the SD of a normal distribution (used in sensitivity analysis)
Td	phase shift between the signals from the two objects in collision	s
t	time	s
τ	delayed time	s
t′	time variable used in the convolution Eq. (1)	s
θ	angle of directivity acoustic sources – sensor	∘
Tc	duration of Hertzian contact	s
Td	delayed time (delayed propagation due to the geometry of particles)	s
Uimp	impact velocity	m s−1
X	position on the cross section of the Isère River (marked on the bridge from the left to right banks)	m

Appendix B

Table B1 Coefficients C_1…6 used in the analytical model of impact of Akay and Hodgson (1978).

The acoustic pressure field generated by the Hertzian impact between sphere and slab is used to model elementary spectra contained in the dictionary Δ. The analytical temporal solutions, obtained from the integral convolution of Eq. (1) and using the geometric setup of Fig. 1b, are rewritten below from Akay and Hodgson (1978)'s paper. Thus, equations 6a-b and 7 from the paper of Akay and Hodgson (1978) are reformulated here in Eqs. (B1)–(B2) and, respectively, (B3). This analytical solutions model is a two-branch function, depending on the duration contact T_d. Thereafter, the total acoustic pressure field, during and after the impact, is obtained by subtracting the individual pressure fields. Such a resulting waveform was shown in Fig. 2a and was modeled using Eq. (B3). It is important to note that another way to compute the energy spectral density of the impact is to numerically compute the Fourier transform on this equation.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\to \mathrm{0}\le \mathit{\tau }\le T_{\mathrm{d}}\\ {\displaystyle }& {\displaystyle }{p}_i^{\left(\mathrm{1}\right)}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathit{\rho }ca}{\mathrm{2}{r}_i}}\cdot {\displaystyle \frac{C_{\mathrm{1}}}{C_{\mathrm{2}}}}\cdot \cos {\mathit{\theta }}_i\\ \text{(B1)}& {\displaystyle }& {\displaystyle }\times \left(\begin{array}{c}{C}_{\mathrm{3}}\sin \left(\frac{\mathit{\pi }\mathit{\tau }}{T_{\mathrm{c}}}\right)-C_{\mathrm{5}}\cos \left(\frac{\mathit{\pi }\mathit{\tau }}{T_{\mathrm{c}}}\right)\\ -\left[C_{\mathrm{6}}\sin \left(\frac{c\mathit{\tau }}{a}\right)-C_{\mathrm{5}}\cos \left(\frac{c\mathit{\tau }}{a}\right)\right]e^{-\frac{c\mathit{\tau }}{a}}\\ +\left(\mathrm{1}-\frac{a}{r_i}\right)\left[C_{\mathrm{6}}\sin \left(\frac{\mathit{\pi }\mathit{\tau }}{T_{\mathrm{c}}}\right)-C_{\mathrm{4}}\sin \left(\frac{\mathit{\pi }\mathit{\tau }}{T_{\mathrm{c}}}\right)\right]\\ -\left(\mathrm{1}-\frac{a}{r_i}\right)\left[C_{\mathrm{6}}\cos \left(\frac{c\mathit{\tau }}{a}\right)+C_{\mathrm{5}}\sin \left(\frac{c\mathit{\tau }}{a}\right)\right]e^{-\frac{c\mathit{\tau }}{a}}\end{array}\right),{\displaystyle }& {\displaystyle }\to \mathit{\tau }>T_{\mathrm{d}}\\ {\displaystyle }& {\displaystyle }{p}_i^{\left(\mathrm{2}\right)}\left(\mathit{\tau }\right)={\displaystyle \frac{\mathit{\rho }ca}{\mathrm{2}{r}_i}}\cdot {\displaystyle \frac{C_{\mathrm{1}}}{C_{\mathrm{2}}}}\cdot \cos {\mathit{\theta }}_i\times \\ \text{(B2)}& {\displaystyle }& {\displaystyle }\left(\begin{array}{c}-\left\{C_{\mathrm{6}}\sin \left[\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}\right]+C_{\mathrm{5}}\cos \left[\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}\right]\right\}{e}^{-\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}}\\ -\left[C_{\mathrm{6}}\sin \left(\frac{c\mathit{\tau }}{a}\right)-C_{\mathrm{5}}\cos \left(\frac{c\mathit{\tau }}{a}\right)\right]e^{-\frac{c\mathit{\tau }}{a}}\\ -\left(\mathrm{1}-\frac{a}{r_i}\right)\left\{C_{\mathrm{6}}\cos \left[\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}\right]-C_{\mathrm{5}}\sin \left[\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}\right]\right\}{e}^{-\frac{c(\mathit{\tau }-T_{\mathrm{c}})}{a}}\\ -\left(\mathrm{1}-\frac{a}{r_i}\right)\left[C_{\mathrm{6}}\cos \left(\frac{c\mathit{\tau }}{a}\right)+C_{\mathrm{5}}\sin \left(\frac{c\mathit{\tau }}{a}\right)\right]e^{-\frac{c\mathit{\tau }}{a}}\end{array}\right),\text{(B3)}& {\displaystyle }& {\displaystyle }p\left(\mathit{\tau }\right)=p_{\mathrm{1}}^{\left(\mathrm{1}\right)\cup \left(\mathrm{2}\right)}\left(\mathit{\tau }\right)-p_{\mathrm{2}}^{\left(\mathrm{1}\right)\cup \left(\mathrm{2}\right)}(\mathit{\tau }-T_{\mathrm{d}}),\end{array}$$

where i=1 (the impacting sphere) or 2 (the mirrored sphere), both with the same radius a, $\mathit{\tau }=t-\frac{\mathrm{2}{r}_i-D}{\mathrm{2}c}$ is the delayed time, D is the sphere's diameter, c is the the velocity of speed, r is the reference distance, θ is the angle between source and sensor (see Fig. 1b) and the constants C_1…6 are

$${\displaystyle }{\mathit{\alpha }}_m={\left[{\displaystyle \frac{\mathrm{15}}{\mathrm{16}}}{\displaystyle \frac{\mathit{\pi }{U}_{\text{imp}}^{\mathrm{2}}({\mathit{\zeta }}_{\mathrm{2}}+{\mathit{\zeta }}_{\mathrm{1}})m_{\mathrm{2}}}{\sqrt{a}}}\right]}^{\mathrm{0.4}},{\mathit{\zeta }}_{\mathrm{1},\mathrm{2}}={\displaystyle \frac{\mathrm{1}-{\mathit{\nu }}_{\mathrm{1},\mathrm{2}}^{\mathrm{2}}}{\mathit{\pi }{E}_{lon{g}_{\mathrm{1},\mathrm{2}}}}}.$$