the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Automatic identification of alternating morphological units in river channels using wavelet analysis and ridge extraction
Nicolas Le Moine
Roger Moussa
Oldrich Navratil
Pierre Ribstein
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. The alongstream (longitudinal) and crosssectional geometry of natural rivers is known to vary at the scale of the hydrographic network (e.g., generally decreasing slope, increasing width in the downstream direction), allowing parameterizations of main crosssectional parameters with largescale proxies such as drainage area or bankfull discharge (an approach coined downstream hydraulic geometry, DHG). However, higherfrequency morphological variability (i.e., at river reach scale) is known to occur for many stream types, associated with varying flow conditions along a given reach, such as the alternate bars or the pool–riffle sequences and meanders. To consider this highfrequency variability of the geometry in the hydraulic models, a first step is to design robust methods to characterize the scales at which it occurs. In this paper, we introduce new wavelet analysis tools in the field of geomorphic analysis (namely, wavelet ridge extraction) to identify the pseudoperiodicity of alternating morphological units from a general point of view (focusing on pool–riffle sequences) for six small French rivers. This analysis can be performed on a single variable (univariate case) but also on multiple variables (multivariate case). In this study, we choose a set of four variables describing the flow degrees of freedom: velocity, hydraulic radius, bed shear stress, and a planform descriptor that quantifies the local deviation of the channel from its mean direction. Finally, this method is compared with the bedform differencing technique (BDT), by computing the mean, median, and standard deviation of their longitudinal spacings. The two methods show agreement in the estimation of the wavelength in all reaches except one. The method aims to extract a pseudoperiodicity of the alternating bedforms that allow objective identification of morphological units in a continuous approach with the maintenance of correlations between variables (i.e., at many station hydraulic geometry, AMHG) without the need to define a prior threshold for each variable to characterize the transition from one unit to another.
 Article
(4995 KB) 
Supplement
(306 KB)  BibTeX
 EndNote
Hydraulic modeling is based on the description of river morphology (crosssectional geometry), and this is the essential input of models despite its scarcity and cost of acquisition. 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). This component is essential for accurate modeling of river hydraulics such as flood modeling (e.g., Neal et al., 2015; Trigg et al., 2009), river restoration (e.g., Wheaton et al., 2004), ecohydraulics (e.g., Pasternack and Brown, 2013), environmental modeling, and fluvial process (e.g., Rodríguez et al., 2013). Longitudinal variability in river geometry may have a greater impact on the simulation of the water level than the crosssectional shapes (Saleh et al., 2013), and it must be taken into account in the hydraulic models. This variability of river geometry at a small scale in the longitudinal and crosssectional direction yields a variation in flow parameters and is known to occur for many morphological channel types, each type being characterized by typical morphological units (MUs), e.g., pools, riffles, steps, point bars, and meanders.
1.1 Stateoftheart methods for a quantitative assessment of morphological variability within a reach
Morphological units are topographic forms that shape the river corridor (Wadeson, 1994; Wyrick and Pasternack, 2014). They form alternating and rhythmic undulations continuously varying along the river (Thompson, 2001). This continuity is challenging to represent; for this reason, most of the methods that model these patterns divide the topography into discrete units to analyze them (Kondolf, 1995; Wyrick et al., 2014).
Among the most frequently observed alternating MUs, pools and riffles have been recognized as fundamental geomorphological elements of meandering streams (Krueger and Frothingham, 2007). Pools are located in the outer edge of each meander loop and defined as topographic lows along a longitudinal stream profile with high depth and low velocity (Fig. 1a, b and d), and research has shown that they generally have an asymmetrical crosssection shape. Conversely, riffles are topographic highs with shallow depths and moderate to high velocities located in the straight parts of the reach between adjacent loops (Fig. 1a, c and d), and they have symmetrical crosssection shapes (O'Neill and Abrahams, 1984; Knighton, 1981).
For many years, many researchers have been trying to develop techniques to identify MUs and especially pools and riffles using hydraulic variables or topographic ones, or both (Table 1). In 1D identification, some studies used bed topography only to determine the characteristics of MUs. Richards (1976a) proposed the zerocrossing method, which fits a regression line to the longitudinal profile of the bed elevation and defines pools as points that have negative residuals and riffles as points with positive residuals. O'Neill and Abrahams (1984) developed the bedform differencing technique (BDT) as a refinement of Richards' methodology. This one uses bed elevations measured at a fixed interval along the channel to calculate the bed elevation difference series between local extrema (maximum and minimum) of the bed profile. The BDT introduces a tolerance value (T); it 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 of the bed elevation difference series and eliminates the erroneous classification of small undulations in the bed profile. Knighton (1981) proposed the Areal Difference Asymmetry Index to identify the location of pools and riffles by their symmetrical or asymmetrical areas. This index is defined as the ratio of the difference between the area of the right and the left of the channel centerline on the total crosssectional area.
Yang (1971)Richards (1976a)Milne (1982)Knighton (1981)Nordin (1971)Box and Jenkins (1976)O'Neill and Abrahams (1984)Jowett (1993)Milne and Sear (1997)Schweizer et al. (2007)Hauer et al. (2009, 2011)Wyrick and Pasternack (2014)Wyrick et al. (2014)Brown and Pasternack (2017)On the other hand, some studies focused only on hydraulic parameters to identify MUs. For example, Yang (1971) proposed an identification of pools and riffles using the energy gradient and affirmed that the fundamental difference between riffles and pools is the difference in energy gradients. Also, Jowett (1993) proposed a classification criterion with Froude number and velocity–depth ratio to distinguish between pools, runs, and riffles.
All these methods handle topographic or hydraulic parameters separately. Recently, however, several researchers have improved MU identification through the use of the covariance of several parameters in a multidimensional approach. Schweizer et al. (2007) used a joint depth and velocity distribution to predict pools, runs, and riffles without the knowledge of the river bathymetry. Hauer et al. (2009) used a functional linkage between depthaveraged velocity, water depth, and bottom shear stress to describe and quantify six different hydromorphological units (riffle, fast run, run, pool, backwater, and shallow water) using a conceptual Mesohabitat Evaluation Model (MEM) under various flow conditions. These methods use digital elevation models (DEMs) to extract more information about MUs. For this purpose, Milne and Sear (1997) began with depth to define pool–riffle sequences using ArcGis tools and DEMs to model the geometry of river channels based on fieldsurveyed cross sections on a 3D basis. But, by choosing depth alone, the difference between two bedforms with the same depth becomes difficult to know. By contrast, it is easy with different bed slopes and bed roughness that yield different velocities and shear stresses (Wyrick et al., 2014). So to overcome this and take into account the lateral variation of rivers, Wyrick et al. (2014) proposed a new method for the objective identification and mapping of landforms at the morphological unit scale. They used spatial grids of depth and velocity at low flow estimated using a 2D hydrodynamic model and an expert classification scheme that determine the number and the nomenclature of MUs and the range of base flow depth and velocity of each type.
Brown and Pasternack (2017) chose two variables: the minimum bed elevation and the channel top width across several flow discharges. They calculated the geomorphic covariance structure (GCS); it is a bivariate spatial relationship amongst or between standardized and possibly detrended variables along a river corridor. They found that there is a positive correlation between these two variables. Also, they used an autocorrelation function and power spectral density to prove a quasiperiodic pattern of wide and shallow or narrow and deep cross sections along the river. This pioneering work and other studies (e.g., Richards, 1976b; Carling and Orr, 2000) proved that a single longitudinal cycle might contain a pool with a narrow and deep cross section, a riffle with a wide and shallow cross section, in addition to transitional forms. The work that we introduce in this paper aims to present a spectral method that extracts this pseudoperiodicity from a river to characterize the alternating MUs, and especially pool–riffle sequences, and to identify the key parameter (the wavelength) that characterizes the scale of variability of the river topography. This information can be further used to build a synthetic river such as the RiverBuilder (Pasternack and Zhang, 2020) or the channel builder for simulating the river morphology of Legleiter (2014).
Some of the methods presented in the literature have shown limits in calculating the wavelengths of pool–riffle sequences. Others have given results that are often difficult to interpret in terms of bedform amplitude. This amplitude, which varies according to each bedform, involves the use of the pseudoperiod. A few methods are developed to extract this pseudoperiod from alternating MU rivers. We, therefore, choose to work with wavelet analysis that estimates the local variability strength of a signal and extracts the signal amplitude and wavelength. In this study, we apply continuous wavelet transform (CWT) to calculate the wavelength λ and the dimensionless wavelength spacing λ^{*} (longitudinal spacing) which is
where w_{bf} is bankfull width.
In reality, the longitudinal spacing λ^{*} has several definitions. Some authors have defined the wavelength λ as the distance between riffle crests (e.g., Harvey, 1975; Hogan, 1986) or the distance from the bottom of successive pools (e.g., Keller and Melhorn, 1973, 1978). Other authors have chosen channel width w (e.g., Richards, 1976a, b; Dury, 1983) instead of bankfull channel width w_{bf} (e.g., Leopold et al., 1964). These differences raise questions about the selection of these ratios 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 a new automatic wavelength calculation method that uses the whole covariance structure of a set of hydraulically independent variables without the need for ad hoc thresholding of these variables.
Some researchers have investigated the variability of longitudinal spacing depending on geometric or hydraulic parameters. Rosgen (2001) developed an empirical relationship between the ratio of pooltopool spacing/bankfull width and the channel slope expressed as a percentage based on a negative power function of slope S:
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 longitudinal spacing in a pool–riffle sequence, and they found that 82% of pools were formed by LWD or other obstructions, and increased numbers of obstructions led to a decrease in the 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 meanannual flood and a 5year recurrence interval (Thompson, 2001). Recently, Wyrick and Pasternack (2014) measured the spacing of six different morphological units using a tool in ArcGIS. 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 longitudinal spacing and bankfull width (Richards, 1976a; O'Neill and Abrahams, 1984; Gregory et al., 1994; Knighton, 2014). 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 longitudinal spacing 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).
1.2 Study objectives
The studies that used wavelet analysis in the geomorphological field consist in extracting components of a given spatial series (e.g., w(x), v(x)), but they are not specifically designed to identify pseudoperiodic components in a univariate case, let alone in a multivariate case. For this reason, we introduce an automatic procedure called wavelet ridge extraction defined by Lilly and Olhede (2011) and used in this study to extract the longitudinal spacing of the alternating MUs.
The objective is to extract some quantitative properties of these alternating morphological units, such as the mean and the median of their longitudinal spacing, with a continuous vision of the topography instead of a discrete classification. For that, we focus on two numerical criteria computed at reach scale: the distribution of spacings between morphological units (mean, median, etc.) and the evaluation of correlations between all geometrical and flow variables. We use in this work four classical variables (velocity, hydraulic radius, bottom shear stress, and the local channel direction angle) because they respond directly to morphodynamic processes (flow convergence routing or meander migration), and they are independent hydraulic degrees of freedom.
In this study, we focus mainly on alternating alluvial channels, especially pool–riffle sequences, even though the method presented here could be used to analyze any morphology characterized by alternating topographic forms. We first present the dataset of six river reaches in France used for this analysis (Sect. 2). In Sect. 3, we present the wavelet ridge extraction method to identify pool–riffle sequences in the univariate and multivariate cases with the four variables. Section 4 presents results and compares them with the BDT developed by O'Neill and Abrahams (1984) to determine whether they yield the same results in terms of spacing. We choose this method instead of threshold methods because the latter require ad hoc thresholding/parameter range definition from independent calibration data, which was not possible in our case.
The rationale behind this approach is to provide a continuous description of geometric and flow patterns along a reach using the wavelength, a description that could be subsequently used to create a synthetic river as in the RiverBuilder (Brown et al., 2014).
Six reaches of small French rivers are used in this study (Navratil, 2005; Navratil et al., 2006): the Graulade at St Sylvain Montaigut (1), the Semme at Droux (2), the Olivet at Beaumont (3), the Ozanne at Tirzay lès Bonneval (4), the Avenelles at BoissyleChâtel Les Avenelles (5), and the Orgeval at BoissyleChâtel Le Theil (6) (Fig. 2). These reaches contain mainly pool–riffle sequences. They have slopes between 0.002 and 0.013 m m^{−1} (estimated from the talweg elevation, which is the lowest point in the section), mobile gravel beds, stable banks, and welldefined 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), and their length ranges from 155 to 495 m (Table 2). All reaches are located at or near the stream gauging stations of the French national hydrometric network. Longterm (about 20 years) hydrological records are available for most reaches. The bankfull widths vary from 4 to 12 m, with an average value of about 9 m.
Cross sections were surveyed along the river reaches at the level of hydraulic controls and morphological breaks 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. Cross sections and water surface profile measurements were surveyed in 2002–2004, covering the main channel and floodplain and using an electronic, digital, totalstation theodolite. Water surface profiles were measured at different flow discharges (Navratil et al., 2006). Using this dataset, we solely rely on measurements at the lowest surveyed discharge in the development of the method because it is the discharge through which we can visualize the variability of the bathymetry (alternating morphological units). We select four spatial series:

velocity v(x);

hydraulic radius ${R}_{\mathrm{h}}\left(x\right)=\frac{A\left(x\right)}{P\left(x\right)}\approx \frac{A\left(x\right)}{w\left(x\right)}$ with A(x) the crosssectional area and P(x) the wetted perimeter;

bed shear stress ${\mathit{\tau}}_{\mathrm{b}}\left(x\right)=\left(\mathit{\rho}g\right){n}^{\mathrm{2}}v\left(x{)}^{\mathrm{2}}{R}_{\mathrm{h}}\right(x{)}^{\mathrm{1}/\mathrm{3}}$ with ρ water density (1000 kg m^{−3}) and n Manning's roughness coefficient;

local channel direction angle (planform) θ(x).
All descriptors are derived from in situ observations taken from Navratil et al. (2006), except the calibrated estimates of Manning’s roughness coefficient n. These values were estimated by Navratil et al. (2006) using a 1D open channel steady and step backwater model, FLUVIA (Baume and Poirson, 1984). However, we use these values to compute the bed shear stress τ_{b}(x), along the reach: even if it partly relies on calibration, it is a more robust way of computing τ_{b} here than through the finite differentiation of the total head function $\frac{v(x{)}^{\mathrm{2}}}{\mathrm{2}g}+z\left(x\right)$ between adjacent cross sections to get the energy slope J, given the typical number and spacing of surveyed cross sections for each reach in the dataset.
The fourth variable chosen is related to the channel planform: we define θ(x) as the local angular deviation of the channel direction from a lowerfrequency curve. There are many possible definitions of this lowfrequency behavior, such as parametric splines or Bezier curves; in order to avoid overparameterization, we define this lowfrequency planform as a constant curvature curve, i.e., the bestfitting arccircle (Fig. 3), a choice suitable for all six reaches studied. Since θ is signed, it is expected to have a pseudoperiodicity, which is approximately twice as slow as other 1D variables. Indeed, a large positive value of θ indicates a counterclockwise deviation from the lowfrequency direction, while a large negative value of the same amplitude indicates a clockwise deviation. From a hydraulic perspective, both deviations have the same effect since they are symmetrical with respect to the lowfrequency direction. For this reason, we chose to analyze the variable cos (θ(x)).
Classical mathematical methods, such as Fourier analysis, extract the wavelengths in the frequency domain for stationary signals, but can also be used for nonstationary signals using an “evolutive” methodology based on spectral estimators (Thomson, 1982; Pasternack and Hinnov, 2003). Wavelet transform standardly does the same for nonstationary signals: analyzing a signal basically consists in looking for the local similarity between the signal and a given waveform (the wavelet). In this paper, we use the continuous wavelet transform with the Morlet wavelet (Gabor, 1946) (Fig. 4) applied to spatial series instead of time series, so periods and frequencies in time series are replaced by wavelengths (in m) and wavenumbers (in rad m^{−1}). The choice of the Morlet wavelet is justified by the analytical properties in its derivation and its flexibility due to the exponential form (see Appendix B).
The wavelet transform uses a whole family of “daughter” wavelets generated by scaling and translating the mother wavelet ψ; the value of the transform at location x and scale s is the scalar product of the signal and this daughter wavelet ψ_{s,x}.
Wavelet analysis is prevalent 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 FoufoulaGeorgiou, 1993; Kumar, 1996), geophysics (e.g., Ng and Chan, 2012; Grinsted et al., 2004), hydrology (e.g., Rossi et al., 2011; Schaefli et al., 2007; Nourani et al., 2014), and geomorphology (Gangodagamage et al., 2007; Lashermes and FoufoulaGeorgiou, 2007; McKean et al., 2009). In the literature of the alternating bedforms' identification, McKean et al. (2009) used a derivative of a Gaussian wavelet (DOG) of order 6 to investigate the spatial patterns (pools and riffles) of channel morphology and salmon spawning using a 1D elevation profile of the channel bed morphology.
In this study, we use another application of the wavelet analysis called the wavelet ridge extraction method (Mallat, 1999; Lilly and Olhede, 2010). This analysis is based on the existence of special 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; Ozkurt and Savaci, 2005). Along such a curve, the signal can be approximated by a single component modulated in both amplitude and frequency. So, the rationale behind the method is that the existence of alternating morphological units along a reach (such as pool–riffle sequences) could be translated into a pseudoperiodicity in geometric and flow variables. Hence, identifying these bedforms requires identification of a local wavenumber K(x) and phase Φ(x) for each variable, a task that can be performed by wavelet analysis and especially wavelet ridge extraction (Mallat, 1999; Lilly and Olhede, 2010).
3.1 Wavelet analysis and ridge extraction
Few methods in the literature have been trying to identify river characteristics with wavelets. For example, Gangodagamage et al. (2007) used wavelet transform modulus maxima (WTMM, Muzy et al., 1993) in a fractal analysis to extract multiscale statistical properties of a corridor width. Procedures such as the WTMM consist in extracting components of the signal, but they are not specifically designed to identify pseudoperiodic components in a univariate case, let alone in a multivariate case.
In the present study, we test a new wavelet ridge analysis on spatial series with the Morlet mother basis function represented in Fig. 4. Its expression is
with ψ the mother wavelet function that depends on the dimensionless “position” parameter η and β 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 ψ.
with x the translation offset, 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 i𝒲[f](x,s) depending on 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 ϕ is the phase (argument) at position x with the scale s.
where Im is the imaginary part of ln (W[f(x)](x,s)). To maintain the nomenclature in the spatial definition and facilitate the extraction of wavelengths, we choose the angular wavenumber (in rad m^{−1}) $k=\frac{\mathrm{2}\mathit{\pi}}{\mathit{\lambda}}$ 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, and
Thus, the wavelet transform of the function f(x) is defined in the space wavenumber as
Except for the channel angle, all input variables are always positive and may substantially vary in magnitude, so we perform the wavelet transform on the Neperian logarithm of these variables. The whole analysis is performed in a simple Scilab script, using the functions that compute the wavelet transform 𝒲[f] and are provided by Torrence and Compo (1998) (https://atoc.colorado.edu/research/wavelets/, 11 June 2020). We followed the procedure in Appendix B to compute $\frac{\partial \mathit{\varphi}}{\partial x}$ and extract the curves that satisfy Eqs. (12) and (13).
The complex wavelet transform can be classically visualized using a scalogram, i.e., a colored map of the modulus R(x,k) in the (x,k) plane (Fig. 5, bottom). The wavelet analysis neglects parts of the signal at both extremities of the series: this is the cone of influence (Torrence and Compo, 1998) that is the region of the wavelet spectrum in which edge effects become important. However, as explained previously, the complex transform also yields a phase ϕ(x,k) in rad (Eq. 8); it can also be plotted in the same plane (Fig. 5, top). In our study, we search for space/wavenumber curves mainly using phase information, i.e., search for phase ridges as opposed to amplitude ridges (Lilly and Olhede, 2010).
In Sect. 3.2, we give a rigorous definition of wavelet ridge points and curves in a univariate case (i.e., a single spatial series). Then, in Sect. 3.3, we generalize the definition to the multivariate, i.e., when the series consists of several correlated variables.
3.2 Univariate case
In the univariate case, we choose a single variable f (velocity, hydraulic radius, bed shear stress, or local channel direction angle). For the wavelet ψ(η), the ridge point of 𝒲[f](x,k) is a space–wavenumber pair (x,k) satisfying the phase ridge point conditions (Lilly and Olhede, 2010):
or, according to the definition of the phase (Eq. 8):
This condition states that the rate of change of the transform phase at scale k exactly matches k at location x; from this condition, the instantaneous frequency of the signal can be derived (Lilly and Olhede, 2008, 2010). The sets of points satisfying the condition form a parametric curve (ridge curve) noted (x,K(x)) implicitly defined by
This property is illustrated in Fig. 5, where a ridge curve is superposed both on the scalogram and on the phase map.
There may be several curves that verify Eq. (14); in practice, we choose curves that cross the domain of the wavelet transform (from one cone of influence to another) continuously and belong to the region where the maximum power of the wavelet is. This curve K(x) also represents the local wavenumber, which is defined on a support ℓ<L named assessed length, 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. B1a.
In the end, we can extract the wavelength function of pool–riffle sequences, which corresponds to a pseudoperiod function of the signal f, and which is
Also, the shape's amplitude A_{m}, with which pools and riffles vary, is corrected by a coefficient $\sqrt{\frac{\mathrm{1}}{\mathit{\alpha}K\left(x\right)}}$. This correction comes from the inversion of the direct transformation equation (Eq. 11) which holds the coefficient $\sqrt{\mathit{\alpha}K\left(x\right)}$.
The signal is locally similar to a sinusoid f_{mod} of wavenumber K in rad m^{−1}, which models the variability of f. We can define the pseudoperiodic variable as presented in Fig. 6 with
In the example below (Fig. 6), the modeled velocity function follows the variability of the observed velocity; it is a pseudoperiodic, continuous function that approximates the firstorder variability of this hydraulic parameter across pool–riffle sequences. The statistics of the K(x) function can be translated into statistics of longitudinal spacings of alternating bedforms, e.g., mean spacing ${\mathit{\lambda}}_{\mathrm{mean}}^{*}$, median spacing ${\mathit{\lambda}}_{\mathrm{median}}^{*}$ or spacing standard deviation σ(λ^{*}). In Fig. 6 we would find ${\mathit{\lambda}}_{\mathrm{mean}}^{*}\approx \mathrm{8.7}$, ${\mathit{\lambda}}_{\mathrm{median}}^{*}\approx \mathrm{9.12}$ and $\mathit{\sigma}\left({\mathit{\lambda}}^{*}\right)\approx \mathrm{0.79}$ if we were to analyze velocity only. The pseudoperiodicity of v_{mod} results in the identification of six pools (white) and seven riffles (gray).
In the next section, we extend the definition of phase ridge points and ridges to the case where several variables are sampled along the reach, all of them potentially correlated and embedding information about the pseudoperiodicity of the channel's hydraulic behavior.
3.3 Multivariate case
The multivariate case is the extension of the univariate case to a set of N realvalued signals. We use the coevolution of more than one variable to extract the wavelength of the reach and therefore identify the pool–riffle sequences. We start by computing the wavelet transform for each variable i=1 … N and extract their phase functions ϕ_{i}(x,k). According to the previous section, univariate ridge curves K_{i}(x) would be defined by
But then the local wavenumber would be specific to a given variable. On the other hand, the multivariate case requires one to determine a common wavenumber between all the variables such that
The identification of a “master” ridge point/curve is now a minimization problem. We define it as a local minimum of the squared norm of the vector $\left({\left(\right)}_{\frac{\partial {\mathit{\varphi}}_{\mathrm{1}}}{\partial x}}(x,k),,k,,,{\left(\right)}_{\frac{\partial {\mathit{\varphi}}_{\mathrm{2}}}{\partial x}},(x,k)\right)k,\phantom{\rule{0.125em}{0ex}}\mathrm{\dots},\phantom{\rule{0.125em}{0ex}}{\left(\right)}_{\frac{\partial {\mathit{\varphi}}_{N}}{\partial x}}(x,k)k$:
This minimum is calculated by searching for the wavenumbers and positions where the derivatives (Eq. 22) of this quantity satisfy
and
The procedure is applied to a set of variables [v, R_{h}, τ_{b}, θ] and the goal is to seek for the common wavenumber between all these variables. In Fig. 7 we illustrate the result of this procedure applied to the Olivet River for all four variables. A unique wavenumber is extracted, which represents a coevolution of all these variables.
As a result, the phase shift of every variable is calculated by
This ridge curve K(x) is common between all variables, yet Φ_{i} varies according to each variable. Therefore, each one can be represented as a pseudoperiodic function f_{i,mod} with the pair (K(x), Φ_{i}(x)).
In our case, after calculating the phase and amplitude, we modeled each variable as in Eq. (24) and represented them in Fig. 8.
The amplitude shape of the modeled variable is calculated by the same way in the univariate case:
The results in Fig. 8 show that a common pseudoperiod has been successfully identified and allows consistent pseudoperiodic representation of all four variables.
Figure 9 shows the correlations between these variables: it can be seen that the pseudoperiodic reconstruction preserves the correlation structure between the three flow variables; an anticorrelated hydraulic radius with bed shear stress and velocity and a strong correlation between bed shear stress and velocity. However, concerning the angle, the results show a small phase shift, which is corrected in the following x positions. But generally, a deviation (clockwise or counterclockwise) from the average direction of the channel (i.e., cos (θ) much smaller than 1) is associated with a low hydraulic radius and large values of τ_{b} and v, a consistent characterization of a riffle. These pseudoperiodic functions give us an identification of reach features: pools (in white) and riffles (in gray).
As already mentioned in Sect. 3.2 (univariate analysis), the statistics of the K(x) function can be translated into statistics of local wavelength $\mathit{\lambda}\left(x\right)=\frac{\mathrm{2}\mathit{\pi}}{K\left(x\right)}$, which can, in turn, be interpreted as statistics of longitudinal spacings of alternating bedforms, e.g., mean spacing ${\mathit{\lambda}}_{\mathrm{mean}}^{*}$, median spacing ${\mathit{\lambda}}_{\mathrm{median}}^{*}$ or spacing standard deviation σ(λ^{*}). In the example of the Olivet River (Fig. 9) ${\mathit{\lambda}}_{\mathrm{mean}}^{*}\approx \mathrm{8.16}$, ${\mathit{\lambda}}_{\mathrm{median}}^{*}\approx \mathrm{8.62}$ and $\mathit{\sigma}\left({\mathit{\lambda}}^{*}\right)\approx \mathrm{0.70}$. The pseudoperiodicity of the set [v_{mod}, R_{h,mod}, τ_{b,mod}, cos (θ)_{mod}] results in the identification of five pools and five riffles.
In this section, we present the results of the analysis on the six reaches presented in Sect. 2. We compare the univariate to the multivariate approach and also the multivariate to the benchmark method. First, the methods are compared in terms of the statistics (mean, median, etc.) they yield. Second, we present the benchmark method called the BDT and compare its results for the six reaches with the multivariate case.
4.1 Univariate vs. multivariate
First, both approaches are employed on all reaches to extract statistics such as the mean, median, and standard deviation wavelengths of morphological units (pool–riffle sequences). The wavelet method extracts the wavelength for an assessed length ℓ (which is the K(x) support in Figs. 6 and 9) that is generally small compared to the total length of the reach. Consequently, we have results that are valuable only for the lengths shown in Table 3. In this table, we give the values of the assessed lengths for each approach, including the variables used in it. These values generally depend on the number of alternating bedforms and also on the total length of the reach. The greater the number of alternating bedforms and the reach length are, the greater the assessed length is.
Moreover, the multivariate approach takes into account all the variables and therefore looks for a single pseudoperiodicity between the four variables, and then we are going to have a pseudoperiodicity that represents the reach and not the chosen variable.
Figure S1 in the Supplement shows the comparison between the univariate and multivariate results for the six reaches (from 1 to 6) using the four variables (velocity, hydraulic radius, bed shear stress, and cosine of local channel direction angle. Table 4 gives some statistics on both approaches. Longitudinal spacing is calculated using the wavelengths extracted automatically by the wavelet ridge method from K(x).
We compare the methods in terms of longitudinal spacing (λ^{*}). In each reach, there seems to be one variable which drives the wavelength identified in the multivariate approach:

in the Graulade River, the longitudinal spacing identified using the multivariate approach matches closely the one associated with the hydraulic radius (in the mean and the median with a deviation of 0.05w_{bf}) and also with the local channel direction angle (in the median with a deviation of 0.06w_{bf});

in the Semme River, it matches those of the local channel direction angle (in the mean and the median with a deviation of 0.14w_{bf} and 0.12w_{bf} consecutively);

in the Olivet River, it matches the bed shear stress (in the mean with a deviation of 0.25w_{bf}) and the velocity (in the median with a deviation of 0.5w_{bf});

in the Ozanne River, it matches those of the hydraulic radius and the velocity (in the mean and the median with a deviation less than 0.6w_{bf});

in the Avenelles, it matches those of the velocity, hydraulic radius, and bed shear stress (in the mean with a deviation less than 0.15w_{bf});

in the Orgeval River, it matches those of the hydraulic radius (in the mean with a deviation of 0.28w_{bf} and the median with 0.06w_{bf}) and also with the local channel direction angle (in the mean with a deviation of 0.23w_{bf} and in the median with 0.11w_{bf}).
Consequently, the multivariate estimates of λ^{*} compare with univariate estimates in a similar way:

the distribution of λ^{*} in the multivariate case is included in the envelope of univariate distributions;

the dispersion of this multivariate distribution, measured by σ(λ^{*}), is always close to the minimum value that can be achieved by any of the univariate distributions.
Hence, the multivariate method improves the identification of the wavelength: it is less sensitive to a local highfrequency variation of a given variable if this variation is not associated with a variation of the other variables. However, there is no direct way of validating the estimates from these raw results: a way of doing so would be to build a synthetic, equivalent periodic geometry parameterized by the identified wavelength to verify that it yields, for example, a similar reachaverage rating curve. This will be the subject of further work.
In the following section, we compare the wavelet method with a benchmark method using talweg elevation.
4.2 Comparison with benchmark method
In this section, we compare our method's results with a selected benchmark method from the literature (i.e., BDT). This method shows good results in the identification of these bedforms according to some studies (e.g., Frothingham and Brown, 2002; Krueger and Frothingham, 2007).
The technique of O'Neill and Abrahams (1984) (BDT) 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 SD 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), Avenelles (5), and Orgeval (6) reaches, we find the same results. We choose to check one tolerance value for each reach with T=SD. This method identifies pool and riffle positions by assigning a crest as a riffle and a bottom as a pool, and therefore the computation of the wavelengths becomes a little difficult. So, we choose to calculate a series of pool–pool and riffle–riffle spacings, their medians, and their standard deviations and then calculate their averages.
This procedure is applied to all rivers, and the results are depicted in Fig. 10. Table 5 presents statistics of the BDT and displays a comparison between these two types of morphological unit identification and mostly the identification of an average wavelength of the reach.
Figure 10 shows the BDT results on all reaches; this method relies only on topography to determine the positions of pools and riffles. Moreover, it also uses a threshold T (tolerance), but the technique does not need a calibration reach or field investigation to know how to set this threshold. In this figure, round points are pools or riffles, and we can calculate the wavelengths and longitudinal spacing of each reach using positions of these points.
The work of the wavelet analysis is done on the assessed length ℓ. However, the BDT method works on the total length of the reach. This comparison is made to determine how effectively the wavelength extracted by the wavelet analysis can represent the entire reach even if an entire part is left unassessed.
For the wavelet method (Fig. 9), the wavelength extraction is among its objectives, while the BDT does not directly calculate the wavelength. It is computed by averaging the pooltopool and riffletoriffle distances. To compare these two methods, we use only the longitudinal spacing (λ^{*}) as a criterion.
In Table 5, we present the results of the BDT on the total length L and the assessed length ℓ of all reaches. By using the total length L, the longitudinal spacings found with the BDT are close to the ones found with the wavelet analysis (deviation less than 1 times the bankfull width for the median) in all the reaches except the Olivet (deviation of 4w_{bf}). Over the assessed length ℓ, we find very similar results with deviations less than 1 times the bankfull width. However, the shortening of the length (ℓ<L) reduces the number of pools and riffles identified – Graulade (1) and Avenelles (5) – and therefore introduces bias. This indicates that a length greater than two cycles (pool–riffle) is always required to produce a pseudoperiodicity of the reach by both methods, a condition which is not fulfilled for all reaches of our dataset. But for the other rivers except for the Olivet (3) and Orgeval (6) reaches, there is not much improvement if we replace the total length with the assessed one. In this comparison, we found that the wavelengths extracted by the multivariate wavelet analysis are generally included in the variance intervals of the wavelengths found by the BDT. This conclusion is verified in all reaches except the Olivet River (3), where there is a big difference between the longitudinal spacings found by the BDT and by wavelets. This difference is due to the choice of the tolerance value, which is low in our case to the point of not filtering out the highfrequency variability of bed elevation and therefore gives a lower periodicity compared to the wavelets.
In this study, we consider the BDT method to be a benchmark method. We do not consider a specific method to be the “true” one, and we only apply these methods to have a general idea on the uncertainties in the identification of morphological units. This method was chosen not because it is the “best” method for pool–riffle identification, but because it does not use thresholds (except for the tolerance T, which does not depend on field data). It means that it does not require a preliminary calibration of thresholds on velocity, hydraulic radius, etc. on an independent reach (e.g., Wyrick and Pasternack, 2014; Hauer et al., 2009). These thresholds vary from one reach to another and according to the characteristics of each river. For this reason, we did not compare our method with threshold methods on this dataset. In contrast, the results of the longitudinal spacing intervals will be compared with the literature.
For a long time, researchers have found a common interval of longitudinal spacings that vary between 5 and 7 times the channel width (Leopold et al., 1964; Keller, 1972; Richards, 1976a; Gregory et al., 1994). Keller (1972) found that the median is less and varies between 3 and 5 times the channel width. O'Neill and Abrahams (1984), using the BDT method, found the same results but with a median close to 3 times the channel width, and this value can vary according to the tolerance T. Carling and Orr (2000) found lower values than before at about 3w. Recent studies (e.g., Wyrick and Pasternack, 2014) have calculated the longitudinal spacing of six morphological units using 2D identification methods. The averages of these pool and riffle spacings are, respectively, 3.3 and 4.3 times the channel width, which is less than the commonly accepted values of 5–7w.
In this study, the longitudinal spacing varies in the mean and the median from ∼1.8 to 8.6 times the bankfull width, supporting the conclusion of Carling and Orr (2000) that pools are spaced approximately 3 to 7 times the channel width. However, the quoted longitudinal spacing relationships should be considered in the context that the bankfull width and spacing distance are inherently variable even for short length reaches. To illustrate this inherent variability, we found the example of Keller and Melhorn (1978), where the pool–pool spacing values ranged from 1.5 to 23.3 channel widths, with an overall mean of 5.9 (Gregory et al., 1994; Knighton, 2014). This variability in longitudinal spacing is probably related to a short assessed length, a small number of cross sections surveyed, or other factors such as geology, bank characteristics (cohesion), the grain size of the river bed, or artificial channel modifications.
We worked with a dataset that contains cross sections spaced 0.46 to 2.9 times the bankfull width. Other studies have used much shorter spacings (e.g., Pasternack and Arroyo, 2018; Legleiter, 2014) to identify morphological units. Of course, the larger the number of cross sections, the more robust the identified correlations are. Also, we worked with irregularly spaced cross sections, which normally lead to biases in the results. Despite this, the “biased” placement does not impair the overall methodology. This methodology has provided good results in terms of longitudinal spacings, and therefore it can be applied for shorter crosssection spacings to clearly identify these alternate morphologies. The short lengths we found raise questions about the naturality of the rivers. In our case, the rivers are subject to artificial modifications (e.g., bridges, weirs) and rehabilitations, which will have a significant impact on the hydromorphological parameters (width, depth, meandering, etc.). This can have a very important impact on the identification of pseudoperiods.
The wavelet ridge analysis is powerful in identifying pseudoperiods, amplitude, and phase while preserving the correlations between parameters. We can thus identify alternating morphological units more objectively in terms of frequency/wavenumber.
This wavelength can be used to represent the variability of the bathymetry in hydraulic models in cases where we do not have full access to the geometry of the channel (e.g., remote sensing data such as the overcoming Surface Water and Ocean Topography Mission) and the morphology can be modeled by pseudoperiodic functions. Furthermore, it can be implemented in synthetic geometry generators (e.g., River Builder, Pasternack and Zhang, 2020) where the bathymetry and sinuosity wavelengths extracted by the wavelets can be used to model meandering rivers with alternating morphologies. Ultimately, hydraulic modeling will be the true test of the potential of a pseudoperiodic equivalent geometry (e.g., for simulating a reachaverage rating curve).
On the other hand, it presents drawbacks compared to other methods. First, the cone of influence ignores a large part of the river and sometimes biases the results in the case of small total lengths (the Graulade (1) and Semme (2) reaches). Similarly for reach length and number of morphological units as for the number of cross sections, the larger it is, the more robust the results are, and the smaller the relative portion of “unassessed length” is. Still, the method remains a powerful tool for nonstationary analysis. Another problem is the amplitude, which is sometimes overestimated in some regions of the topography. We visualized this in several cases in our study since we used the Neperian logarithm to avoid negative values, and therefore the inverse function (exponential) gives slightly larger values. However, this does not bias the identified wavelength of the reach.
In this study, we present an automatic procedure based on wavelet ridge extraction to identify some characteristics of alternating morphological units (MUs), such as their longitudinal spacing and amplitude. The method does not rely on any a priori thresholds to identify MU sequences. It was applied to six rivers with a maximum length of 500 m. We chose to work with classical hydromorphological variables (velocity, hydraulic radius, bed shear stress) in addition to the planform channel direction angle that evaluates the impact of river sinuosity in the determination of the wavelength.
As a result, identified wavelengths are consistent with the values of the literature (mean in 3–7w_{bf}). The use of a multivariate approach yields more robust results than the univariate approaches by ensuring a consistent covariance of flow variables in the pseudoperiodic behavior.
Given the short length of several reaches, the relatively small number of cross sections for each reach, and the possible impacts of artificial modifications, this paper is mainly a proofofconcept of the wavelet approach. We foresee many perspectives for this work, such as the possibility of extending the work to other rivers with other types of MUs or other longer reaches with a large number of cross sections.
A(x)  Crosssectional area (m^{2}) 
A_{i,mod}  Signal amplitude of the shape of the modeled variable number i 
A_{m}  Signal amplitude of the shape 
cos (θ)  Cosine of the local channel direction angle 
cos (θ)_{mod}  Modeled cosine of the local channel direction angle 
f  Space series function (m) 
f_{i}  Measured space series function number i 
f_{i,mod}  Modeled space series function number i with multivariate wavelet analysis 
f_{mod}  Modeled variable with the univariate wavelet analysis 
g  Acceleration due to gravity and its value is 9.81 m s^{−2} 
k(x)  Wavenumber (rad m^{−1}) 
K(x)  Local wavenumber that corresponds to the maximum variance of the signal (rad m^{−1}) 
K_{i}(x)  Local wavenumber that corresponds to the maximum variance of the variable i (rad m^{−1}) 
ℓ  K(x) support (m) 
L  Total reach length (m) 
N  Number of total chosen variables 
n  Manning's roughness coefficient 
P(x)  Wetted perimeter (m) 
Q_{min}  Minimum discharge modeled (m^{3} s^{−1}) 
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) 
R_{h,mod}  Modeled hydraulic radius (m) 
s  Dilation or scale factor 
S  Reach slope (m m^{−1}) 
SD  The standard deviation of the bed elevation difference series (m) 
T  Tolerance value used in BDT method (m) 
v(x)  Velocity (m s^{−1}) 
v_{mod}  Modeled velocity (m s^{−1}) 
v_{obs}  Measured velocity (m s^{−1}) 
w  Reach width (m) 
w_{bf}  Reach bankfull width (m) 
w_{m}  Mean bankfull width (m) 
x  Translation factor in the wavelet transform or the abscissa position (m) 
y=y_{max}  Water depth measured from the talweg elevation $y={z}_{\mathrm{ws}}z$ (m) 
y_{m}  Mean depth (m) 
z  Measured bed elevation 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 is taken to be 6 recommended by Torrence and Compo (1998) 
λ  Reach wavelength (m) 
λ^{*}  Typical pool (riffle) spacing or dimensionless reach wavelength 
ρ  Water density (997 kg m^{−3}) 
θ  local channel direction angle in rad 
σ(w)  The standard deviation of the width along the reach (m) 
σ( )  Standard deviation 
τ_{b}(x)  Bed shear stress in the x abscissa (Pa) 
τ_{b,mod}(x)  Modeled bed shear stress in the x abscissa (Pa) 
Φ  Corresponding phase at the position x and the wavenumber K with $\mathrm{\Phi}\left(x\right)=\mathit{\varphi}(x,K(x\left)\right)$ (rad) 
Φ_{i}  Phase at the position x and the wavenumber K for the variable number i 
ϕ(x,s) or ϕ(x,k) or ϕ  Phase or argument at a position x and with a scale s or wavenumber k (rad) 
ϕ_{i}  The phase of the variable number i 
ψ  Mother wavelet function 
ψ_{s,x}  Daughter wavelet function 
η  Dimensionless position parameter 
ℂ  Complex numbers 
ℝ  Real numbers 
${\mathbb{R}}_{+}^{*}$  Positive real numbers 
𝒲[f](x,s)  Continuous wavelet transform of f(x) with the wavelet ψ 
B1 The univariate case
The conjugate form of the mother wavelet is
Its derivative depending on the mute variable η is
In Sect. 3.1, η is a mute integration variable, and x appears only in the argument αk(ξ−x) of the function ψ^{*}. By applying the derivation formula of a composite function, the derivative of the wavelet transform is expressed by
On the other hand, we have
Finally,
The previous expression numerically avoids the derivative of the function ϕ(x,k), which varies quickly for large wavenumbers.
B2 The multivariate case
In the multivariate case, we should resolve Eq. (20), which contains three derivatives to compute. The first one is already done in the univariate case:
The second one is the computation of $\frac{{\partial}^{\mathrm{2}}{\mathit{\varphi}}_{i}(x,k)}{\partial k\partial x}$:
For that, we should develop $\frac{\partial}{\partial k}\left(\frac{\mathcal{W}\left[x{f}_{i}\left(x\right)\right](x,k)}{\mathcal{W}\left[{f}_{i}\left(x\right)\right](x,k)}\right)$:
We calculate each derivative:
We find a general formulation with p=0 … N−2:
The third one is the computation of $\frac{{\partial}^{\mathrm{3}}{\mathit{\varphi}}_{i}(x,k)}{\partial {k}^{\mathrm{2}}\partial x}$:
with
and
The script is available upon request from the authors.
The supplement related to this article is available online at: https://doi.org/10.5194/hess2435132020supplement.
NLM, MM, and RM designed the study. ON provided the dataset. MM and NLM built the methodology and analyzed the results. MM wrote the paper with contributions from NLM, PR, and RM.
The authors declare that they have no conflict of interest.
The authors acknowledge financial support from the French National Space Agency (Centre National d'Etudes Spatiales, CNES) and Sorbonne University through the PhD grant of Mounir Mahdade. We would also like to thank Gregory B. Pasternack, the second anonymous reviewer, and the editor for the valuable suggestions and comments that greatly helped improve this paper.
This research has been supported by the Centre National d'Etudes Spatiales (CNES; grant no. CNES5100018525) and Sorbonne University (grant no. C18/1749).
This paper was edited by Theresa Blume and reviewed by Gregory Pasternack and one anonymous referee.
Alfieri, L., Feyen, L., Salamon, P., Thielen, J., Bianchi, A., Dottori, F., and Burek, P.: Modelling the socioeconomic impact of river floods in Europe, Nat. Hazards Earth Syst. Sci., 16, 1401–1411, https://doi.org/10.5194/nhess1614012016, 2016. a
Baume, J.P. and Poirson, M.: Modélisation numérique d'un écoulement permanent dans un réseau hydraulique maillé à surface libre, en régime fluvial, La Houille Blanche, 1–2, 95–100, 1984. a
Box, G. E. and Jenkins, G. M.: Time series analysis: Forecasting and control HoldenDay, San Francisco, California, 1976. a
Brown, R. A. and Pasternack, G. B.: Bed and width oscillations form coherent patterns in a partially confined, regulated gravel–cobblebedded river adjusting to anthropogenic disturbances, Earth Surf. Dynam., 5, 1–20, https://doi.org/10.5194/esurf512017, 2017. a, b
Brown, R. A., Pasternack, G. B., and Wallender, W. W.: Synthetic river valleys: Creating prescribed topography for form–process inquiry and river rehabilitation design, Geomorphology, 214, 40–55, 2014. a
Carling, P. A. and Orr, H. G.: Morphology of rifflepool sequences in the River Severn, England, Earth Surf. Proc. Land., 25, 369–384, 2000. a, b, c, d, e
Carmona, R. A., Hwang, W. L., and Torrésani, B.: Multiridge detection and timefrequency reconstruction, IEEE T. Sig. Process., 47, 480–492, 1999. a
Clifford, N. J.: Formation of riffle–pool sequences: field evidence for an autogenetic process, Sediment. Geol., 85, 39–51, 1993. a
Dury, G.: Osagetype underfitness on the River Severn near Shrewsbury, Shropshire, England, Background to palaeohydrology, Wiley, Chichester, 399–412, 1983. a
Frothingham, K. M. and Brown, N.: Objective identification of pools and riffles in a humanmodified stream system, Mid. Stat. Geogr., 35, 52–60, 2002. a
Gabor, D.: Theory of communication. Part 1: The analysis of information, J. Inst. Elect. Eng. Pt. III, 93, 429–441, 1946. a
Gangodagamage, C., Barnes, E., and FoufoulaGeorgiou, E.: Scaling in river corridor widths depicts organization in valley morphology, Geomorphology, 91, 198–215, 2007. a, b
Gregory, K., Gurnell, A., Hill, C., and Tooth, S.: Stability of the poolriffle sequence in changing river channels, River Res. Appl., 9, 35–43, 1994. a, b, c, d
Grinsted, A., Moore, J. C., and Jevrejeva, S.: Application of the cross wavelet transform and wavelet coherence to geophysical time series, Nonlin. Processes Geophys., 11, 561–566, https://doi.org/10.5194/npg115612004, 2004. a
Harvey, A.: Some aspects of the relations between channel characteristics and riffle spacing in meandering streams, Am. J. Sci., 275, 470–478, 1975. a, b
Hauer, C., Mandlburger, G., and Habersack, H.: Hydraulically related hydromorphological units: description based on a new conceptual mesohabitat evaluation model (MEM) using LiDAR data as geometric input, River Res. Appl., 25, 29–47, 2009. a, b, c
Hauer, C., Unfer, G., Tritthart, M., Formann, E., and Habersack, H.: Variability of mesohabitat characteristics in rifflepool reaches: Testing an integrative evaluation concept (FGC) for MEMapplication, River Res. Appl., 27, 403–430, 2011. a
Higuchi, H., Lewalle, J., and Crane, P.: On the structure of a twodimensional wake behind a pair of flat plates, Phys. Fluids, 6, 297–305, 1994. a
Hogan, D. L.: Channel morphology of unlogged, logged, and debris torrented streams in the Queen Charlotte Islands, Ministry of Forests and Lands, British Columbia, Canada, 1986. a
Jowett, I. G.: A method for objectively identifying pool, run, and riffle habitats from physical measurements, N. Zeal. J. Mar. Freshw. Res., 27, 241–248, 1993. a, b
Katul, G. G. and Parlange, M. B.: Analysis of land surface heat fluxes using the orthonormal wavelet approach, Water Resour. Res., 31, 2743–2749, 1995a. a
Katul, G. G. and Parlange, M. B.: The spatial structure of turbulence at production wavenumbers using orthonormal wavelets, Bound.Lay. Meteorol., 75, 81–108, 1995b. a
Katul, G. G., Parlange, M. B., and Chu, C. R.: Intermittency, local isotropy, and nonGaussian statistics in atmospheric surface layer turbulence, Phys. Fluids, 6, 2480–2492, 1994. a
Keller, E.: Development of alluvial stream channels: a fivestage model, Geol. Soc. Am. Bull., 83, 1531–1536, 1972. a, b
Keller, E. and Melhorn, W.: Bedforms and fluvial processes in alluvial stream channels: selected observations, edited by: Morisawa, M., Fluvial Geomorphology, St. Univ. New York, Binghamton, USA, 253–283, 1973. a
Keller, E. and Melhorn, W.: Rhythmic spacing and origin of pools and riffles, Geol. Soc. Am. Bull., 89, 723–730, 1978. a, b
Knighton, A.: Asymmetry of river channel cross‐sections: Part I. Quantitative indices, Earth Surf. Proc. Land., 6, 581–588, 1981. a, b, c
Knighton, D.: Fluvial forms and processes: a new perspective, Routledge, New York, USA, 2014. a, b
Kondolf, G. M.: Geomorphological stream channel classification in aquatic habitat restoration: uses and limitations, Aquat. Conserv.: Mar. Freshw. Ecosyst. 5, 127–141, 1995. a
Krueger, A. and Frothingham, K.: Application and comparison of geomorphological and hydrological pool and riffle quantification methods, Geogr. Bull., 48, 85–95, 2007. a, b, c
Kumar, P.: Role of coherent structures in the stochasticdynamic variability of precipitation, J. Geophys. Res.Atmos., 101, 26393–26404, 1996. a
Kumar, P. and FoufoulaGeorgiou, E.: A multicomponent decomposition of spatial rainfall fields: 1. Segregation of largeand smallscale features using wavelet transforms, Water Resour. Res., 29, 2515–2532, 1993. a
Lashermes, B. and FoufoulaGeorgiou, E.: Area and width functions of river networks: New results on multifractal properties, Water Resour. Res., 43, W09405, https://doi.org/10.1029/2006WR005329, 2007. a
Legleiter, C. J.: A geostatistical framework for quantifying the reachscale spatial structure of river morphology: 1. Variogram models, related metrics, and relation to channel form, Geomorphology, 205, 65–84, 2014. a, b
Leopold, L., Wolman, M., and Miller, J.: Fluvial processes in geomorphology, Freeman, San Francisco, 1964. a, b, c
Lilly, J. M. and Olhede, S. C.: Higherorder properties of analytic wavelets, IEEE T. Sig. Process., 57, 146–160, 2008. a
Lilly, J. M. and Olhede, S. C.: On the analytic wavelet transform, IEEE T. Inform. Theory, 56, 4135–4156, 2010. a, b, c, d, e, f
Lilly, J. M. and Olhede, S. C.: Analysis of modulated multivariate oscillations, IEEE T. Sig. Process., 60, 600–612, 2011. a
Mallat, S.: A wavelet tour of signal processing, Elsevier, San Diego, CA, USA, 1999. a, b
McKean, J., Nagel, D., Tonina, D., Bailey, P., Wright, C. W., Bohn, C., and Nayegandhi, A.: Remote sensing of channels and riparian zones with a narrowbeam aquaticterrestrial LIDAR, Remote Sens., 1, 1065–1096, 2009. a, b
Milne, J.: Bedmaterial size and the rifflepool sequence, Sedimentology, 29, 267–278, 1982. a
Milne, J. and Sear, D.: Modelling river channel topography using GIS, Int. J. Geogr. Inform. Sci., 11, 499–519, 1997. a, b
Montgomery, D. R., Buffington, J. M., Smith, R. D., Schmidt, K. M., and Pess, G.: Pool spacing in forest channels, Water Resour. Res., 31, 1097–1105, 1995. a, b
Muzy, J.F., Bacry, E., and Arneodo, A.: Multifractal formalism for fractal signals: The structurefunction approach versus the wavelettransform modulusmaxima method, Phys. Rev. E, 47, 875–884, 1993. a
Navratil, O.: Débit de pleins bords et géométrie hydraulique: une description synthétique de la morphologie des cours d'eau pour relier le bassin versant et les habitats aquatiques, Institut National Polytechnique de Grenoble, Grenoble, France, 2005. a
Navratil, O., Albert, M., Herouin, E., and Gresillon, J.: Determination of bankfull discharge magnitude and frequency: comparison of methods on 16 gravelbed river reaches, Earth Surf. Proc. Land., 31, 1345–1363, 2006. a, b, c, d, e, f
Neal, J. C., Odoni, N. A., Trigg, M. A., Freer, J. E., GarciaPintado, J., Mason, D. C., Wood, M., and Bates, P. D.: Efficient incorporation of channel crosssection geometry uncertainty into regional and global scale flood inundation models, J. Hydrol., 529, 169–183, 2015. a
Ng, E. K. and Chan, J. C.: Geophysical applications of partial wavelet coherence and multiple wavelet coherence, J. Atmos. Ocean. Tech., 29, 1845–1853, 2012. a
Nordin, C. F.: Statistical properties of dune profiles, in: vol. 562, US Government Printing Office, Washington, USA, 1971. a
Nourani, V., Baghanam, A. H., Adamowski, J., and Kisi, O.: Applications of hybrid wavelet – Artificial Intelligence models in hydrology: A review, J. Hydrol., 514, 358–377, 2014. a
O'Neill, M. P. and Abrahams, A. D.: Objective identification of pools and riffles, Water Resour. Res., 20, 921–926, 1984. a, b, c, d, e, f, g
Ozkurt, N. and Savaci, F. A.: Determination of wavelet ridges of nonstationary signals by singular value decomposition, IEEE T. Circuit. Syst. II, 52, 480–485, 2005. a
Pasternack, G. and Arroyo, R.: RiverBuilder: River Generation for Given Data Sets, R package version 0.1, availabe at: https://cran.rproject.org/web/packages/RiverBuilder/index.html, (last access: 11 June 2020), 2018. a
Pasternack, G. and Zhang, M.: River Builder User's Manual For Version 1.0.0, University of California, Davis, CA, available at: https://github.com/RiverBuilder/RiverBuilder, last access: 11 June 2020. a, b
Pasternack, G. B. and Brown, R. A.: 20 Ecohydraulic Design of RifflePool Relief and Morphological Unit Geometry in Support of Regulated GravelBed River Rehabilitation, in: Ecohydraulics, Wiley Online Library, p. 337, 2013. a
Pasternack, G. B. and Hinnov, L. A.: Hydrometeorological controls on water level in a vegetated Chesapeake Bay tidal freshwater delta, Estuar. Coast. Shelf Sci., 58, 367–387, 2003. a
Richards, K.: The morphology of rifflepool sequences, Earth Surf. Proc. Land., 1, 71–88, 1976a. a, b, c, d, e
Richards, K.: Channel width and the rifflepool sequence, Geol. Soc. Am. Bull., 87, 883–890, 1976b. a, b
Rodríguez, J. F., García, C. M., and García, M. H.: Threedimensional flow in centered poolriffle sequences, Water Resour. Res., 49, 202–215, 2013. a
Rosgen, D. L.: The crossvane, wweir and jhook vane structures … their description, design and application for stream stabilization and river restoration, in: Designing successful stream and wetland restoration projects, editec by: Hays, D. F., Proceedings of the 2001 Wetlands Engineering and River Restoration, 27–31 August 2001, Reno, NV, USA, Am. Soc. Civil Eng., Reston, VA, USA, 1–22, https://doi.org/10.1061/40581(2001)72, 2001. a
Rossi, A., Massei, N., and Laignel, B.: A synthesis of the timescale variability of commonly used climate indices using continuous wavelet transform, Global Planet. Change, 78, 1–13, 2011. a
Saleh, F., Ducharne, A., Flipo, N., Oudin, L., and Ledoux, E.: Impact of river bed morphology on discharge and water levels simulated by a 1D Saint–Venant hydraulic model at regional scale, J. Hydrol., 476, 169–177, 2013. a
Schaefli, B., Maraun, D., and Holschneider, M.: What drives high flow events in the Swiss Alps? Recent developments in wavelet spectral analysis and their application to hydrology, Adv. Water Resour., 30, 2511–2525, 2007. a
Schneider, K. and Vasilyev, O. V.: Wavelet methods in computational fluid dynamics, Annu. Rev. Fluid Mech., 42, 473–503, 2010. a
Schweizer, S., Borsuk, M. E., Jowett, I., and Reichert, P.: Predicting joint frequency distributions of depth and velocity for instream habitat assessment, River Res. Appl., 23, 287–302, 2007. a, b
Thompson, D. M.: Random controls on semirhythmic spacing of pools and riffles in constrictiondominated rivers, Earth Surf. Proc. Land., 26, 1195–1212, 2001. a, b
Thomson, D. J.: Spectrum estimation and harmonic analysis, Proc. IEEE, 70, 1055–1096, 1982. a
Torrence, C. and Compo, G. P.: A practical guide to wavelet analysis, B. Am. Meteorol. Soc., 79, 61–78, 1998. a, b, c, d
Trigg, M. A., Wilson, M. D., Bates, P. D., Horritt, M. S., Alsdorf, D. E., Forsberg, B. R., and Vega, M. C.: Amazon flood wave hydraulics, J. Hydrol., 374, 92–105, 2009. a
Wadeson, R.: A geomorphological approach to the identification and classification of instream flow environments, S. Afr. J. Aquat. Sci., 20, 38–61, 1994. a
Wheaton, J. M., Pasternack, G. B., and Merz, J. E.: Spawning habitat rehabilitation – I. Conceptual approach and methods, Int. J. River Basin Manage., 2, 3–20, 2004. a
Wyrick, J. R. and Pasternack, G.: Geospatial organization of fluvial landforms in a gravel–cobble river: beyond the riffle–pool couplet, Geomorphology, 213, 48–65, 2014. a, b, c, d, e
Wyrick, J. R., Senter, A. E., and Pasternack, G. B.: Revealing the natural complexity of fluvial morphology through 2D hydrodynamic delineation of river landforms, Geomorphology, 210, 14–22, 2014. a, b, c, d
Yang, C. T.: Formation of Riffles and Pools, Water Resour. Res., 7, 1567–1574, https://doi.org/10.1029/WR007i006p01567, 1971. a, b
 Abstract
 Introduction
 Dataset and study reaches
 Wavelet method
 Results
 Discussion
 Conclusions
 Appendix A: List of symbols
 Appendix B: Mathematical calculus for the wavelet transform
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Dataset and study reaches
 Wavelet method
 Results
 Discussion
 Conclusions
 Appendix A: List of symbols
 Appendix B: Mathematical calculus for the wavelet transform
 Code availability
 Author contributions
 Competing interests
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement