Wavelet and index methods for the identification of pool – riffle sequences

The accuracy of hydraulic models depends on the quality of the bathymetric data they are based on, whatever the scale at which they are applied (e.g., 2D or 3D reach-scale modeling for local flood studies or 1D modeling for network-scale flood routing). The along-stream (longitudinal) and cross-sectional geometry of natural rivers is known to vary at the scale of the hydrographic network (e.g., generally decreasing slope, increasing width, etc.), allowing parameterizations of main cross-sectional parameters with proxy such as drainage area or a reference discharge quantile (an approach coined downstream 5 hydraulic geometry, DHG). However, higher-frequency morphological variability is known to occur for many stream types, associated with varying flow conditions along a given reach: alternate bars, pool-riffle sequences, meanders, etc. To better take this high-frequency variability in bedforms into account in hydraulic models, a first step is to design robust methods to characterize the scales at which it occurs. In this paper, we propose and benchmark several methods to identify bedform sequences in pool-riffle morphology, for six small French rivers: the first one called the index method, based on three morphological and 10 hydraulic descriptors; the second one called wavelet ridge extraction, performed on the continuous wavelet transform (CWT) of bed elevation. Finally, these new methods are compared with the bedform differencing technique (BDT, O’Neill and Abrahams (1984)), compared by computing a score that gives a percentage of agreement along the total surveyed length and by calculating the number of bedforms and the pool spacings for each method. The three methods were found to give similar results on average for wavelength estimation, with agreement from 64% to 84% and a similar number of bedforms identified. 15 The filter-like behavior of the wavelet ridge analysis tends to give more robust results for the estimation of mean bedform amplitude, which varies from 0.30 to 0.81 with an SNR (signal-to-noise ratio) from 2.68 to 7.91. Otherwise, BDT gives higher mean bedform amplitude but lower SNR values from 0.85 to 1.73.


Introduction
Flood forecasting models are based on the description of river morphology (cross-sectional geometry), and this is their essential input despite its scarcity and price (Saleh et al., 2013;Grimaldi et al., 2018).In fact, the most important aspect to know is the river bathymetric data at the local scale, detailed and specific to the site and local conditions (Alfieri et al., 2016), for the accurate modeling of river hydraulics, which is essential for predicting floodplain flooding (Neal et al., 2015;Trigg et al., 2009).
Many researchers are working on determining the best simplified representation of channel geometry (Saleh et al., 2013;Grimaldi et al., 2018;Neal et al., 2015;Orlandini and Rosso, 1998), based on the variability of cross sections but without the knowledge of the bed elevation variability on a small scale.This longitudinal variability in the river geometry has greater impact on the simulation of the water level than the cross-sectional shapes (Saleh et al., 2013) and it must be taken into account in the hydraulic models designed to improve flood forecasting.This longitudinal variability is related to the channel morphology types.
Several classifications of river channel morphology were suggested (e.g., Montgomery and Buffington (1997); Rosgen (1994)).These classifications are broken down into five major categories, based on variables such as channel patterns and channel slope (Robert, 2014).In this paper, we retain the following five main alluvial reach morphologies: cascade, step-pool, plane bed, pool-riffle, and dune-ripple (Montgomery and Buffington, 1997).We summarize all these types in Fig. 1 where each morphology is defined with its characteristics.
In this study, we focus mainly on alluvial pool-riffle sequences, even though the methods presented here could be used to iden-Figure 1.Five essential river morphologies with their typical bed materials, the reach profile corresponding to each type, bed slope, typical pool spacing λ * = λ w bf with, λ the reach wavelength in meters and w bf the average of the bankfull width in meters.These characteristics are taken from the Montgomery and Buffington (1997) study.
tify sequences of various bedforms.The objective is to identify the pool (or riffle) spacing or dimensionless reach's wavelength in the river channel, which is the distance between pools (or riffles) divided by average channel width (Richards, 1976a;Keller and Melhorn, 1978;Carling and Orr, 2000) or bankfull width (Leopold et al., 1964).In this paper, we use a normalization by the bankfull width (w bf ).
with λ the reach wavelength in meters and w bf the reach average bankfull width in meters.The purpose of this identification is to produce a high-frequency variability of cross sections in order to apply it to hydraulic models (see section 2).
From a morphological point of view, pool-riffle sequences are defined as rhythmic undulations in bed topography (O'Neill and Abrahams, 1984).According to Richards (1976a), that definition is better than those based on hydrodynamic variables (e.g., Froude number, velocity) because it changes less with discharge.In fact, pools are defined as topographic lows along a longitudinal stream profile (Fig. 2 (B) and (D)) and research has shown that they generally have an asymmetrical cross section shape.Conversely, riffles are topographic highs (Fig. 2 (C) and (D)) that have symmetrical cross section shapes (O'Neill and Abrahams, 1984;Knighton, 1981).
All these identification methods, which we present in detail in the following section, have shown limits in calculating the wavelengths of pool-riffle sequences and in most cases their results are often difficult to interpret in terms of bedform amplitude.
This amplitude, which varies according to each bedform, goes with the notion of the pseudo-period.We therefore choose to work with wavelet analysis, a new method in the hydraulic and morphological field, that estimates the local variability strength of a signal and allows one to extract the signal amplitude and wavelength.The wavelet transform has been used for numerous studies in several domains, especially geophysics and electronics (Torrence and Compo, 1998), and to analyze time series that contain nonstationary power at many different frequencies (Daubechies, 1990).In this study we apply continuous wavelet transform (CWT) to spatial series instead of time series, to calculate the dimensionless pool (or riffle) spacing λ * between ri−1 and ri intervals and between pi−1 and pi intervals.
We start with the state of the art of pool-riffle identification and pool (or riffle) spacings in section 2. We present our reaches and data set in section 3. Also, we present two methods to identify pool-riffle sequences in several reaches in France.Each method is studied separately in sections 4 and 5.The first one, called the index method, combines hydraulic and morphological descriptors and splits pools and riffles using a threshold on a criterion index.The second one, called wavelet method, uses ridge curve analysis to determine the wavelength and amplitude of the signal that contains the pool-riffle sequences.Section 6 compares these methods and the bedform differencing technique (BDT) developed by O'Neill and Abrahams (1984) to determine if they yield the same results.To accomplish this, a score is defined to compare the pool-pool and riffle-riffle conformity in the index and wavelet methods.In addition, we calculate the wavelength λ, the pool spacing λ * , the number of bedforms for each method, and the results of the hydraulic variables (F r, R h , and u y ) for each method, and compare their In fact, the choice of the BDT is explained in the paper of Krueger and Frothingham (2007), which finds that the BDT is a more appropriate technique to use to identify these bedforms than the Fr method, especially during low flow.
2 Pool-riffle and pool spacing identification 5 Among the morphological methods to identify pool-riffle sequences, the regression analysis by Richards (1976a)  residuals.The limitation of this method is that small-scale undulations in the bed topography are incorrectly identified as separate bedforms and can be classified as pools or riffles.This is why O' Neill and Abrahams (1984) developed the BDT as a refinement of Richards' methodology.Other researchers have investigated the accuracy and agreement of pool-riffle identification methods (e.g., Frothingham and Brown (2002) compared two geomorphological methods of pool-riffle identification: BDT by O' Neill and Abrahams (1984) and the areal difference asymmetry index by Knighton (1981)).
The BDT uses bed elevations measured at a fixed interval along the channel.Using this data, we calculate the bed elevation difference series between local extrema (maximum and minimum) of the bed profile.The BDT introduces a tolerance value (T ), which is the minimum absolute value of the cumulative elevation change required for the identification of a pool or riffle.
The value of T is based on the standard deviation (S D ) of the bed elevation difference series and eliminates the erroneous classification of small undulations in the bed profile.Several values of T should be tested, and for odd-numbered series that follow each bedform, the value of the cumulative is updated and the absolute value of the cumulative is compared with T to ascertain whether the series culminates in a pool or riffle (O'Neill and Abrahams, 1984).
However, what a hydrologist would classify as a pool or riffle would not necessarily agree with the classification of a geomorphologist (Krueger and Frothingham, 2007).In fact, hydrologists define pools and riffles by flow parameters such as water depths, velocities (Allen, 1951), water surface slope (Yang, 1971) and Froude number (Jowett, 1993).Yang (1971) affirmed that the fundamental difference between riffles and pools is the difference in energy gradients, and in a complete cycle of a pool-riffle sequence, the riffle is defined as the portion that has an energy gradient (water surface slope) steeper than the average energy gradient of the complete cycle, whereas the pool is the portion that has an energy gradient milder than the cycle average.Jowett (1993) proposed and compared multiple hydrological identification methods (velocity/depth ratio, Froude number, and water surface slope).In fact, the Froude number is the dimensionless velocity/depth ratio F r = u m √ gy , where u m is the mean water column velocity, y the water depth, and g the acceleration due to gravity (9.81m.s−2 ).It relates inertia forces to gravity forces and is important wherever gravity dominates.In addition, the Froude number has been recognized as a criterion to distinguish between pools and riffles (Wolman, 1955;Bhowmik and Demissie, 1982).
The study of Jowett (1993) suggested that pool and riffle habitats are relatively easy to distinguish, but that run habitats (transitional facies) are difficult to distinguish from riffle habitats.Pools occur only where the local stream gradient is low and riffles or rapids in areas where the gradient is high.Moreover, Allen's classification gave pools a Froude number less than 0.15 and a velocity/depth ratio less than about 0.8, and riffles a Froude number greater than 0.25 or a velocity/depth ratio greater than about 1.8.
Another study by Krueger and Frothingham (2007) applied a common geomorphological and hydrological method to identify pool-riffle sequences.It compared the identification agreement between the two methods on several reaches on Ransom Creek, New York.Other studies concerning pool-riffle sequences occurred in these decades, for instance de Almeida and Rodríguez (2011)  Furthermore, pool spacings have several definitions.Some authors have defined the wavelength as the distance between riffle crests (e.g., Harvey (1975); Hogan et al. (1986)), others have chosen the distance from the bottom of successive pools (e.g., Keller andMelhorn (1973, 1978)).Some have also chosen channel width (e.g., Richards (1976a, b); Dury (1983)), otherwise others have used bankfull channel width (e.g., Leopold et al. (1964)).These differences raise questions about the certainty of these ratios (dimensionless pool spacings λ * ) and their dependence on geometric or hydraulic parameters.Moreover, the majority of researchers use the average channel width instead of the bankfull width because both give a similar pool-riffle spacing interval.Here, we are working with w bf and with two wavelength calculation methods (Wavelet & index) that identify the part where the pool-riffle sequences are located and their repetitions on the reach.
A great deal of research has investigated the variability of pool spacing in relation to geometric or hydraulic parameters.Rosgen ( 2001) developed an empirical relationship between the ratio of pool-to-pool spacing/bankfull width and the channel slope expressed as a percentage based on a negative power function λ * = 8.2513S −0.9799 .In addition, Montgomery et al. (1995) showed that there is an influence of large woody debris (LWD) on channel morphology that leads to a relation between LWD and pool spacing in a pool-riffle sequence, and found that 82% of pools were formed by LWD or other obstructions, and increased numbers of obstructions led to a decrease in pool-riffle spacing.Moreover, research has linked variation in spacing to channel characteristics including gradient (Gregory et al., 1994).Also, Harvey (1975) showed that pool-riffle spacing correlated strongly with discharges between the mean-annual flood and a 5 year recurrence interval (Thompson, 2001).
Therefore, the definition of the characteristics and the measurement methods allowed us to expect some variation from one study to another in the estimated relationship between pool spacing and bankfull width (Richards, 1976a;O'Neill and Abrahams, 1984;Gregory et al., 1994;Knighton, 1998).Aside from the interval [5w bf , 7w bf ] defined by Leopold et al. (1964) and the interval [2w bf , 4w bf ] defined by Montgomery et al. (1995) in forested streams, other values of the spacing pool exist, such as the Carling and Orr (2000) interval, which is [3w bf , 7.5w bf ] and decreases to [3w bf , 6w bf ] as sinuosity increases (Clifford, 1993;Carling and Orr, 2000).

Data set and study reaches
Six reaches of small French rivers are used in this study (Navratil, 2005;Navratil et al., 2006): the Graulade in St Sylvain Montaigut (1), the Semme in Droux (2), the Olivet in Beaumont Village (3), the Ozanne in Tirzay lès Bonneval (4), the Avennelles in Boissy-le-Châtel Les Avennelles (5), and the Orgeval in Boissy-le-Châtel Le Theil (6) (Fig. 3).These reaches were chosen in such a way that we have pool-riffle sequences, i.e., all reaches have slopes less than 0.015, mobile gravel beds, stable banks, and well-defined floodplains along at least one side of the channel (Navratil et al., 2006).These reaches are located in the Loire River Basin (four reaches) and the Seine River Basin (two reaches), their length ranges from 155 to 1034 m, and their drainage area is from 19 to 268 km 2 (Table 1).All reaches are located at or near the stream gauging stations of the French national hydrometric network.Long-term (about 20 years) hydrological records are available for most reaches (Navratil Many cross sections are surveyed along the river reaches at the level of hydraulic controls and morphological breaks in order to describe the major variations in terms of width, height, and slope in the main channel and the floodplain and at the level of pool-riffle sequences.The hydraulic data set is also surveyed at least for two low discharges and for one almost bankfull discharge (z ws , u, ...).
These data are modeled between Q min and Q max (Table 1) with Fluvia (Baume and Poirson, 1984) 4 Index method Bedforms are defined by both morphological and hydrological characteristics.In this paper, we develop a method that combines hydraulic and geomorphological variables into a single index.It takes into account several descriptors that are selected using the PCA technique (principal component analysis) during the minimum discharge Q min (Table 1).First, we choose nine variables namely: bed elevation z, water surface elevation z ws , wetted surface A, wetted perimeter P , hydraulic radius R h , maximum depth y max , mean depth y m , velocity u, and Froude number F r. PCA is a method that represents these dependent and intercorrelated variables in order to extract the important information and reproduce it in new orthogonal variables (principal components), and its goal is to find similarities and eliminate redundant variables (Abdi and Williams, 2010).
Before processing, it is necessary to normalize the data and remove the trend from z and z ws ; the detrended variables are then normalized by: with p the normalized variable, p i the detrended variable, p m its mean, and σ(p i ) its standard deviation.
The PCA assumes that the directions that contain the high variance are the most important (called principal).In Fig. 4 (A), we find two dimensions that contain the most explained variances, the first and the second ones, so we will work on these two axes.Figure 4 (B) illustrates the results and shows the representation quality (cos 2 ) of the variables on the PCA graphs, which is the square of the coordinate of the variable and it is also represented by the arrow's length.A high cos 2 indicates a good representation of the variable on the main axes, and in this case, the variable is positioned near the correlation circumference circle.A low cos 2 indicates that the variable is not perfectly represented by the main axes, and the variable is close to the circle's center.Consequently, to eliminate the redundancy and accurately characterize the bedform variability, according to Fig. 4 and the results of all the reaches, three descriptors are chosen to define the index for all reaches: bed elevation z, hydraulic radius R h , and Froude number F r.The identifying technique of pool-riffle sequences is based on an index that is defined as a function: with N the number of the p j descriptors chosen and 1 the Heaviside step function, defined as: The function 1 returns to 1 if the normalized descriptor is positive, otherwise it returns to 0 (e.g., bed elevation z in Fig. 5 (B)).This method is applied to the three descriptors (the opposite of hydraulic radius −R h , bed elevation z, and Froude number F r) for the six reaches.The index I, in Fig. 6, can take four different values (0, 1 3 , 2 3 , and 1).In fact, we observe that these results are depicted in stairs, and to make the identification of the bedform easier, we introduced colored bars in Fig. 6  By ignoring the intermediate morphologies, we present the discrete pool-riffle identification of the Orgeval reach in Table 2 with 16 bedforms, eight pools, eight riffles, an average wavelength of 38.5 m for riffles and 39.7 m for pools, and λ * ≈ 3.9m.Interval abscissa (m) [105.5,143.5] [143.5,158.5] [158.5,182.5] [182.5,210.5] [210.5,227.5]Index method Riffle Pool Riffle Pool Riffle Interval abscissa (m) [227.5,248.5] [248.5,261.5] [261.5,269.5] [269.5,290.5] [290.5,317.5

Wavelet method
Classical mathematical methods, such as Fourier analysis, extract the wavelengths in the frequency domain only for stationary signals.Wavelet transforms can do the the same for nonstationary signals and find the localized wavelength.Analyzing a signal basically consists of looking for similarity between the signal and well-known mathematical functions.In this paper, we use the continuous wavelet transform with the Morlet wavelet (Gabor, 1946) applied to spatial series instead of time series, so periods and frequencies in time series are replaced by wavelengths and wavenumbers in spatial series.In this section, we introduced the wavelet analysis for spatial series in section 5.1, the wavelet ridge analysis and steps to extract amplitude and wavelength of the pool-riffle channel in section 5.2, and finally the pool-riffle identification for the Orgeval reach in section 5.3.

Definitions
Wavelets are used because of their aptitude to space-wavenumber and space-scale analysis.The advantage of analyzing a signal with these functions, is that it makes it possible to study the features of the signal locally with the level of detail matched to their scale (Kumar and Foufoula-Georgiou, 1997).They are mathematical functions that cut up data into different frequency components and then study each component with a resolution matched to its scale; however, the structure of the calculations remains the same (Graps, 1995).
Wavelet analysis is very popular in many fields such as fluid mechanics (e.g., Schneider and Vasilyev (2010); Higuchi et al. (1994); Katul et al. (1994); Katul and Parlange (1995b, a)), meteorology (e.g., Kumar and Foufoula-Georgiou (1993); Kumar where ψ is the mother wavelet function that depends on the dimensionless "position" parameter η and β is the dimensionless frequency, here taken to be 6 as recommended by Torrence and Compo (1998).Starting with this wavelet mother, a family ψ s,x called wavelet daughters is obtained by translating and scaling ψ.
where x is the translation factor, which represents a position at which the signal is analyzed, and s the dilation or scale factor.
If s > 1, the daughter wavelet has a frequency lower than the mother wavelet, whereas if s < 1, a wavelet with a frequency higher than the mother wavelet is generated.
Given a spatial series f (η), its continuous wavelet transform W f (x, s) with respect to the wavelet ψ is a function of two  variables where: (*) indicates the complex conjugate.This complex function can also be written as: where R is the absolute value (modulus) and φ the phase (argument) at position x with the scale s.
To respect the nomenclature in the spatial definition and facilitate the extraction, we choose the angular wavenumber (in rad.m −1 ) k = 2π λ instead of the scale factor.We associate a wavelength λ = 2παs with the scale parameter s, where α is the Fourier factor associated with the wavelet, with k ψ = 1 α the peak wavenumber of the mother wavelet, and Thus, the wavelet transform of the function f (x) is defined in the space-wavenumber as:

Wavelet ridge analysis
Using the wavelet ridge analysis, we extract the oscillatory components, characterized by pool-riffle sequences along the river, and their properties from the space-frequency representations of wavelet transform.This analysis is based on the existence of spatial space/wavenumber curves, called wavelet ridge curves or simply ridges (Lilly and Olhede, 2010), where the signal concentrates most of its energy (Carmona et al., 1999) and where R 2 is high.A ridge curve is a group of points called ridge points.
The objective is to extract these ridge points that contain not only the wavelength but also the amplitude of the signal.This signal represents the pool-riffle shape.Therefore, the wavelet ridge is used to identify frequency-modulated sinusoidal signal components using a curve at the space-k plane along which the energy is locally maximum (Ozkurt and Savaci, 2005).Indeed, this method is a filtering tool that evacuates all details below the wavelength resolution (or k-resolution) and which keeps only its pattern.In this paragraph, we detail the ridge points extraction for Morlet wavelets.
For the wavelet ψ(η), the ridge point of W f (x, k) is a space/wavenumber pair (x, k) satisfying the amplitude and phase ridge point conditions (Lilly and Olhede, 2010): These conditions state that for each fixed position point x, a phase ridge point φ (or an amplitude ridge point) corresponds to the wavenumber at which a local maximum in the transform magnitude occurs (Lilly and Olhede, 2010).
According to the mathematical development of the partial derivatives in Appendix B and since we are limited only to extracting the phase ridge points φ(x, k), we have Im ln W f (x) (x, k) = φ(x, k), and the equation to solve becomes: Therefore, we search locally in position x for the wavelet daughter function whose wavenumber K(x) takes the maximum variance of the signal, and represents the pool-riffle sequences.
To determine the function K(x), we will work in the space of two variables (x, k) using the phase function φ(x, k) of the wavelet transform, or more precisely its partial derivatives in x (see Appendix B).In fact, the function K(x), that we are looking for could be seen as a parametric curve (x, K(x)) defined on R × R + of the wavelet transform.It is therefore implicitly defined by: The curve K(x) belongs to the level curves that verify Eq. ( 17).This property is illustrated in Fig. 8 (A): the set of points of K(x) (in thick white) verifying Eq. ( 18).
There may be several curves that verify this property; in practice we choose the curve that crosses continuously (without interruption) the domain of the wavelet transform (from one cone of influence to another).It also represents the local wavenumber, which is defined on a support < L, with L the total reach length.
The phase function Φ is then obtained by evaluating the function φ(x, k) along the curve (x, K(x)), in thick black in Fig. B1 (A) in Appendix B.
In the end, we can extract the wavelength function of pool-riffle sequences, which corresponds to a pseudo-period function of the signal f , and which is: Also, the shape's amplitude, with which pools and riffles vary, is corrected by a coefficient k ψ K(x) :

Pool-riffle identification
In this part, we apply the ridge wavelet analysis, explained in the previous subsection, to the spatial function z (bed elevation) and we limit the study only with univariate analysis of the pool-riffle sequences of the Orgeval reach (6).We use the toolbox of Torrence and Compo (1998) [atoc.colorado.edu/research/wavelets/],and transform it into Scilab and add modifications to extract the wavelength λ and amplitude A m .
To obtain the pool-riffle sequences, we define a function I W , which is extracted from the daughter wavelet: In fact, we present a wavelet power plot (Fig. 8 (A)) applied to bed elevation z in the wavenumber/distance space to show the higher energy region and the extracting ridge curve that contains all the information about the reach (wavelength, phase, and amplitude).The function I W translates the ridge curve variability in the spatial domain, it is a frequency-and amplitudemodulated sinusoid, and as a result, it gives the bedform variability represented in Fig. 8 (B).This plot shows great conformity between the bed elevation and the signal identified, and the advantage of this method is that, at different scales, it not only gives the wavelength of the pool-riffle repetitions, but it gives the amplitude of each bedform shape associated with the wavelength.
Moreover, wavelet ridge analysis can be extended to the multivariate case (Lilly and Olhede, 2009).
If we analyze the results of this method on the Orgeval reach (6), we find that the wavelet method eliminates part of the bed elevation details and takes the part that contains the high power (Fig. 8 (A) and (B)) on which it identifies the bedforms graphically.Figure 8 (C) shows us the pool-riffle sequences in spectral representation as in the index method section.Therefore, we consider that the parts that are in black are pools, white are riffles, and gray represents intermediate morphologies (runs, flat, or rapids).3 with seven bedforms, four pools, three riffles, a wavelength of 53.4m, and λ * ≈ 5.3m.

Results and discussion
In this section, we present the results, the comparison between methods and the discussion.First, we introduce the comparison method based on a score technique that compares the two methods index and wavelet.Second, we present the application of this comparison on the index and wavelet methods, and we introduce the BDT results.Finally, these results are discussed according to data presented in section 3.

Comparison method
We compare the wavelet method and the index method using the score technique.The idea is to define a symmetric agreement score that quantifies to what extent the methods conform, it measures at which point we could have the same results, and finally it gives a percentage score that represents the common parts of pools and riffles of the total length.Let us assume for example that the wavelet method identifies n pools and m riffles, and the index method identifies n pools and m riffles.Consequently, for the wavelet method, the pool intervals of this type of reach are PW = , where pW i , rW i , pI j , and rI j are successively pool and riffle interval abscissa number i identified by the wavelet method, and pool and riffle interval abscissa number j identified by the index method.Thus, the score is defined as:

Application and results
Each method scans the bed elevations of the reach and identifies pool-riffle sequences.The number of these bedforms varies from one reach to another and goes from 4 to 49 bedforms.Before presenting the results, we assume that pools and riffles are intervals that are different from what O'Neill and Abrahams (1984) considered in their identification.In fact, BDT assumes it works on points and it is only used to check the wavelength and the number of bedforms of the reach (Table 5).
Another issue is that the work is done on a length hundreds of meters long (which is the K(x) support in Fig. 8) in which we represent the results of the three methods at the same time; this length is calculated between the centers of two intervals.In addition, the other intermediate morphologies are ignored in this comparison.
By studying the case of the Orgeval reach ( 6), which has an identification interval between 81 and 241 m, and using the results of the two methods illustrated in Fig. 9 below and Tables 2 and 3 that we have 82% accuracy to find the same bedform using the two methods.If we analyze this reach, we have seven intervals; on each of them we identify the corresponding bedform for each method (Tables 2 and 3) and calculate the score.As a result, we find high percentages for most of the reaches, from 64% to 84%, showing that these two methods are accurate in this paper (Table 4 presents the scores of all the reaches).To validate and discuss these two methods, the technique of O'Neill and Abrahams (1984) (BDT) is chosen.This technique uses a tolerance value (T), which defines the minimum absolute value needed to identify a pool or a riffle (Krueger and Frothingham, 2007).It is calculated using the standard deviation (S D ) of the series of bed elevation differences from upstream to downstream for each reach and corrected by a coefficient chosen according to the reach.For this, we test several tolerance values, and for the Graulade (1), Ozanne (4), and Orgeval (6) reaches we find the same results.We choose to check two tolerance values for each reach with T = S D and T = 2S D .
For the wavelet method (Fig. 8), the wavelength extraction is among its objectives, while for the index and BDT, the computation of the wavelength is done by first averaging the pool-to-pool spacings, then averaging the riffle-to-riffle spacings, and finally calculating their mean.This kind of calculation is used to prevent overestimation or underestimation of the wavelength.
For example, for the Orgeval reach (6), we find a wavelength of 53.4 m for the wavelet method, 50.8 m for the index and 49.8 for BDT.However, we choose to compare the number of bedforms instead of the number of pools and riffles for each method, and the dimensionless pool spacings instead of the reach wavelengths.

Discussion
To provide an accurate assessment of each method, the variability of these results (Tables 5 and 6) must be discussed with some of the parameters in Table 1 and with the literature.In fact, the computation of the wavelength of each reach and the calculation of the numbers of bedforms remain the best morphological solution.
Moreover, we compare not only the dimensionless pool spacings λ * and number of bedforms but also the Froude number F r, the hydraulic radius R h , the velocity/depth ratio, and the signal-to-noise ratio of each method.In this section, we choose the BDT results, which are closer to the other methods and to reality.
First, in the case of shorter lengths (the Avennelles reach), the index method underestimates the number of bedforms (four bedforms) and for longer lengths (the Olivet reach) the wavelet method also underestimates the number of bedforms ( 30bedforms instead of 49 and 46 bedforms).Second, by examining the relation between the number of bedforms and λ * , we find that when there is the same number of bedforms for the three methods (the Graulade and Orgeval reaches), with a deviation of Table 5. Summary of results: n is the number of pools and m is the number of riffles, λIndex and λBDT are the average between the wavelengths calculated between two pools and between two riffles, BDT1 means BDT for T = SD and BDT2 means BDT for T = 2SD, λ * is calculated by λ w bf , and w bf is taken from Table 1.The amplitude is calculated in the case of BDT for T = SD and it is the difference between two successive extremums, µ(Am)BDT is their mean, σ(Am)BDT their standard deviation, and SNRBDT is the signal-to-noise ratio, which is the ratio of the signal's mean to standard deviation.However, the mean amplitude of the wavelet method is calculated using Eq. ( 21). one bedform at most, we find that the dimensionless pool spacings are very close as well, with a deviation of 0.3 at most.In the case of the Semme (2) and Ozanne (4) reaches, the number of bedforms are different for the three methods with a deviation of two bedforms at most, we find that dimensionless pool spacings vary from one method to another with a deviation of 0.6 at most.This explains why the three methods are almost equivalent and deviations represent only noise, which corresponds to parts ignored by a method and taken into account by another.
On the other hand, the Olivet (3) and Avennelles (5) reaches produce a large deviation in the wavelet and the index methods 5 successively, and this is due to the methods themselves: the index method sweeps all the bedforms present in a studied length and mainly the higher ones while maintaining an accurate wavelength.However, the Avennelles reach shows that this method fails for small reaches where there are few bedforms and subsequently it does not give a good wavelength.Otherwise, the On the other hand, we have:

Figure 2 .
Figure2.Different views of pool-riffle sequences.(A) Plan view pattern that includes bankfull width w bf , floodplain extent, talweg line, velocity u, pool interval abscissa (pi−1, pi, and pi+1), and riffle interval abscissa (ri−1 and ri); (B) cross-sectional view of a pool in the interval pi with a section width w and a steeper water depth y calculated from the talweg elevation, which is the deepest part of the bottom, and y = ymax = zws − z with z is the bed elevation and zws the water surface elevation, and ym the mean water depth; (C) cross-sectional view of a riffle in the interval ri with a shallower water depth y, higher bed elevation z and high bankfull width w bf ; (D) longitudinal profile that makes it possible to see the water surface, the bed slope, the pools and riffles, and the wavelength λ calculated between two successive riffles or pools analyzed pool-riffle morphodynamics using unsteady 1-D flow morphology and the bed-sorting model operated on a continuous basis.Parker et al. (2008) studied the dynamics and pool-riffle evolution.Moreover, Rodríguez et al. (2013) established detailed laboratory flow measurements that allow for the reconstruction of 3-D flow patterns in a pool-riffle design Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-381Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.that aims to provide a level of flow variability similar to what would be expected in a natural stream.

Figure 3 .
Figure 3. Location of the study reaches in France.

Figure 4 .
Figure 4. Results of the PCA for the Orgeval reach.(A) Bar plot that represents the percentage of explained variances within the dimensions axis; it limits the number of axes to a number that represents high variance.We keep the first two main components (Dimension 1 and 2) which contain an acceptable percentage of 82% of information (variances) in data; (B) this graph is known as the circle of correlations represented on the two chosen axes 1 & 2. It shows the relationships between all the variables (bed elevation z, water surface elevation zws, wetted surface A, wetted perimeter P , hydraulic radius R h , maximum depth ymax, mean depth ym, velocity u, and Froude number F r), and therefore those that are positively correlated are grouped together.Negatively correlated variables are positioned on opposite sides of the circle's origin.The distance between the variables and the origin measures how representative the variables are (cos 2 ).Variables that are far from the origin are well represented by the PCA.

Figure 5 .
Figure 5. (A) The sampled and normalized descriptors for the Orgeval reach.z bed elevation, F r Froude number, and −R h opposite of the hydraulic radius, represented on the x axis, which is the distance in meters; (B) the Heaviside step function 1 applied to bed elevation z of the Orgeval reach and represented at the same figure, the objective is to calculate 1 of each descriptor and to take their average.

Figure 7 .
Figure 7. Wavelet Morlet mother function, the plot gives the real part and the imaginary part of the wavelets in the space domain (distance).

Figure 8 .
Figure 8. Steps to identify pool-riffle sequences for the Orgeval reach using bed elevation z with the wavelet method.(A) Power of the wavelet, which is represented in a wavenumber/distance space, with the region where there is maximum variability depicted by the white curve (ridge curve); (B) variation of IW function compared to the bed elevation z. (C) Identification of pools (black), riffles (white), and intermediate morphologies (gray).
Figure 9. Plot showing the identification of pools (black) and riffles (white) ignoring the intermediate morphologies for the Orgeval reach, for the index method (A) and the wavelet method (B).

::::
Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-381Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.extract the pool-riffle amplitudes.It could be extended in a multivariate case with a quantitative approach by crossing several variables.On the other hand, the index method suggests identifying all the major and minor bed profile undulations, but it cannot extrapolate a wavelength to all the reaches based on only a small part of this reach.Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-381Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.that represents alternating pools and riffles k : Wavenumber (rad.m −1 ) K : Local wavenumber that corresponds to the maximum variance of the signal (rad.m −1 ) K s : Strickler coefficient k ψ : Peak wavenumber, which is the spatial frequency domain that has a maximum amplitude (rad.m −1 ) : K(x) support (m) L : Total reach length (m) ([a, b]) : Length of [a, b] (m) N : Number of selected descriptors n and m : Number of pools and riffles identified, respectively, by the wavelet method n and m : Number of pools and riffles identified, respectively, by the index method p Interval of the abscissa of the ith pool identified by the wavelet method pI j Interval of the abscissa of the ith pool identified by the index method p m : Descriptor's mean P : Wetted perimeter (m) PI : Pool interval identified by the index method PW : Pool interval identified by the wavelet method Q : Reach discharge (m 3 .s−1 ) Q min : Minimum discharge modeled (m 3 .s−1 ) Q max : Maximum discharge modeled (m 3 .s−1 ) rW i Interval of the abscissa of the ith riffle identified by the wavelet method rI j Interval of the abscissa of the ith riffle identified by the index method R(x, s) or R(x, k) or R : Absolute value or modulus of the wavelet transform at a position x and with a scale s or wavenumber k R h : Hydraulic radius (m) Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-381Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.RI : Pool interval identified by the index method RW : Riffle interval identified by the wavelet method s : Dilation or scale factor S : Reach slope (m/m) S D : Standard deviation of the bed elevation diffrence series (m) T : Tolerance value, which is the minimum absolute value of the cumulative elevation change required for the identification of pools and riffles (m) u : Velocity (m.s −1 ) u m : Mean velocity (m.s −1 ) w : Reach width (m) w bf : Reach bankfull width (m) x : Translation factor in the wavelet transform or the abscissa position (m) y : Water depth measured from the the talweg elevation y = z ws − z (m) or talweg elevation (m) z ws : Water surface elevation measured from the 0NGF (m) α : Fourier factor associated with the wavelet (m.rad −1 ) β : Dimensionless frequency taken to be 6 recommended by Torrence and Compo (1998) λ : Reach wavelength (m) λ * : Typical pool (riffle) spacing or dimensionless reach wavelength µ() : Variable mean σ(w) : Standard deviation of the width along the reach (m) σ(p i ) : Descriptor's standard deviation Φ : Corresponding phase at the position x and the wavenumber K with Φ(x) = φ(x, K(x)) (rad) φ(x, s) or φ(x, k) or φ : Phase or argument at a position x and with a scale s or wavenumber k (rad) W f (x, s) : Continuous wavelet transform of f (x) with the wavelet ψ +∞ −∞ (iβ − αkx)f (η) + αkηf (η) • ψ * αk(η − x) dη =(αk)(iβ − αkx) W f (x) (x, k) + (αk) 2 W xf (x) (x, k) Hydrol.Earth Syst.Sci.Discuss., https://doi.org/10.5194/hess-2018-381Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 11 September 2018 c Author(s) 2018.CC BY 4.0 License.Finally:∂φ ∂x = Im (αk)(iβ − αkx) + (αk) 2 W xf (x) (x, k) W f (x) (x, k) (B4) ∂φ ∂x = (αk) β + (αk) 2 Im W xf (x) (x, k) W f (x) (x, k)The previous expression numerically avoids the derivative of the function φ(x, k), which varies quickly for large wave numbers.

Figure B1 .
Figure B1.Steps of determining the local wavenumber K(x) using the wavelet transform of the bed elevation z, represented in the four panels.(A) The phase function φ(x, k); (B) the power's cone of influence of the wavelet to characterize the region where there is a maximum variability; (C) the function ∂φ ∂x ; (D) the function k
, a one-dimensional openchannel steady and step backwater model.It is based on interpolations rather than extrapolations, it uses simplified shallow water equations, and it assumes that the friction forces are well represented by the Manning-Strickler formulation.The Strickler coefficient is the model's calibration parameter, adjusted visually or with a minimization function to model water depths that correspond to the depths observed..This code provides us a series of multi-section flows as output for the following variables: bed elevation z, water surface elevation z ws , discharge Q, width w, wetted surface A, wetted perimeter P , hydraulic radius R h , maximum depth y max , mean depth y m , velocity u, Strickler coefficient K s and Froude number F r, where: Figure 6.Index function I, applied to the Orgeval reach, varies between (0, It is depicted on the x axis, which represents the longitudinal distance in meters.The pool-riffle sequences are identified by a spectrum that represents pools in black, riffles in white, and intermediate morphologies (runs, flat, or rapids) in gray.

Table 4 .
The score technique results