Automatic identification of alternating morphological units in river channels using wavelet analysis and ridge extraction

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 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 in the downstream direction), allowing parameterizations of main cross-sectional parameters with large-scale 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 high-frequency 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 pseudo-periodicity 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 pseudo-periodicity 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.


Introduction
Hydraulic modeling is based on the description of river morphology (cross-sectional 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 cross-sectional 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 cross-sectional direction yields a variation in flow parameters and is known to occur Published by Copernicus Publications on behalf of the European Geosciences Union. for many morphological channel types, each type being characterized by typical morphological units (MUs), e.g., pools, riffles, steps, point bars, and meanders.

State-of-the-art methods for a quantitative assessment of morphological variability within a reach
Morphological units are topographic forms that shape the river corridor (Wadeson, 1994;. 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;. 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 cross-section 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 cross-section 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 zero-crossing 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 cross-sectional area.
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 depth-averaged velocity, water depth, and bottom shear stress to describe and quantify six different hydro-morphological 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 field-surveyed 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 . So to overcome this and take into account the lateral variation of rivers,  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 quasi-periodic 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, Figure 1. Different views of pool-riffle sequences. (a) Plan view pattern that includes bankfull width w bf , floodplain extent, talweg line, velocity v, pools and riffles, and channel direction (planform); (b) cross-sectional view of a pool 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 = y max = z ws − z, with z the bed elevation, z ws the water surface elevation, and y m the mean water depth; (c) cross-sectional view of a riffle 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. 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 River-Builder (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 pseudo-period. A few methods are developed to extract this pseudo-period 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 continu-ous 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 andMelhorn, 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

Variables MUs References
Control-point method Energy gradient Pools and riffles Yang (1971) Zero-crossing method Bed topography Pools and riffles Richards (1976a), Milne (1982) Areal difference asymmetry Cross-section area Pools and riffles Knighton (1981)  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 pool-to-pool 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 mean-annual flood and a 5-year recurrence interval (Thompson, 2001). Recently,  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).

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 pseudo-periodic 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 River-Builder (Brown et al., 2014).

Dataset and study reaches
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 Boissy-le-Châtel Les Avenelles (5), and the Orgeval at Boissy-le-Châ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 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), 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. Long-term (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, total-station 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: the cross-sectional area and P (x) the wetted perimeter; ρ water density (1000 kg m −3 ) and n Manning's roughness coefficient; 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 v(x) 2 2g + z(x) 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 lower-frequency curve. There are many possible definitions of this low-frequency behavior, such as parametric splines or Bezier curves; in order to avoid over-parameterization, we define this low-frequency planform as a constant curvature curve, i.e., the best-fitting arccircle (Fig. 3), a choice suitable for all six reaches studied. Since θ is signed, it is expected to have a pseudo-periodicity, which is approximately twice as slow as other 1D variables. Indeed, a large positive value of θ indicates a counterclockwise deviation from the low-frequency 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 low-frequency direction. For this reason, we chose to analyze the variable cos(θ (x)).  Table 2. Characteristics of the six reaches and their catchment. The bankfull width w bf is taken from the study of Navratil et al. (2006), and the average width w m and the standard deviation σ (w) are calculated for minimum discharge Q min .

Wavelet method
Classical mathematical methods, such as Fourier analysis, extract the wavelengths in the frequency domain for stationary signals, but can also be used for non-stationary signals using an "evolutive" methodology based on spectral estimators (Thomson, 1982;Pasternack and Hinnov, 2003). Wavelet transform standardly does the same for non-stationary 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) ap-plied 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 . Figure 3. Definition of θ , the local angular deviation of the channel direction from a lower-frequency behavior. Here this low-frequency planform is defined as an arc-circle (illustration on the Olivet River reach). It is worth noting that θ is signed: at the location indicated in the figure, θ is negative. 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 Foufoula-Georgiou, 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 Foufoula-Georgiou, 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 pseudo-periodicity 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).

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 pseudo-periodic 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 iW[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 = 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, 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 W[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 ∂φ ∂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.

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 W[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 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 poolriffle sequences, which corresponds to a pseudo-period 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 1 αK(x) . This correction comes from the inversion of the direct transformation equation (Eq. 11) which holds the coefficient 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 pseudo-periodic 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 pseudo-periodic, continuous function that approximates the first-order 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 al-ternating bedforms, e.g., mean spacing λ * mean , median spacing λ * median or spacing standard deviation σ (λ * ). In Fig. 6 we would find λ * mean ≈ 8.7, λ * median ≈ 9.12 and σ (λ * ) ≈ 0.79 if we were to analyze velocity only. The pseudo-periodicity 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 pseudo-periodicity of the channel's hydraulic behavior.

Multivariate case
The multivariate case is the extension of the univariate case to a set of N real-valued 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 This minimum is calculated by searching for the wavenumbers and positions where the derivatives (Eq. 22) of this quantity satisfy 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 pseudo-periodic 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 pseudo-period 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 pseudo-periodic reconstruction preserves the correlation structure between the three flow variables; an anti-correlated 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 pseudo-periodic 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 λ(x) = 2π K(x) , which can, in turn, be interpreted as statistics of longitudinal spacings of alternating bedforms, e.g., mean spacing λ * mean , median spacing λ * median or spacing standard deviation σ (λ * ). In the example of the Olivet River (Fig. 9) λ * mean ≈ 8.16, λ * median ≈ 8.62 and σ (λ * ) ≈ 0.70. The pseudo-periodicity of the set [v mod , R h,mod , τ b,mod , cos(θ ) mod ] results in the identification of five pools and five riffles.

Results
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.

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 Figure 8. Variation of the modeled function f i,mod , which represents the pseudo-periodic variable (in red) for the velocity, the hydraulic radius, the bed shear stress, and the local channel direction angle of the Olivet River compared to the observed ones (in blue). 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 pseudo-periodicity 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 high-frequency 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. Table 4. Summary of results for all reaches in the univariate case using the velocity, hydraulic radius, bed shear stress, or local channel direction angle and in the multivariate case using all these four variables. For each variable, we compute the mean, median, and standard deviation σ of the wavelength and the longitudinal spacing. This one (λ * ) is calculated by λ w bf , and w bf is taken from Table 2. Bold values represent the median or the mean wavelengths that are similar in both the univariate case and the multivariate case in a reach.

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 re-sults. 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 Figure 10. Results of the BDT method using a tolerance equal to the standard deviation on the total length (blue) and the assessed one (red) for all reaches (1 to 6). Round points are pools or riffles: pools are high, and riffles are low points. 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 poolto-pool and riffle-to-riffle 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 Table 5. Results of BDT and the multivariate wavelet methods for all reaches. λ meth,mean is the mean wavelength using one of the methods (for BDT, it used on the total length and the assessed one ). λ meth,median is the median, and λ * meth,mean is the mean longitudinal spacing, λ * meth,median is the median, and σ (λ meth,median ) is the standard deviation of the longitudinal spacing. (-) means that we find only one longitudinal spacing, which is the mean and the median, so there is no standard deviation.

Reaches
(1) Graulade (2)   (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 pseudo-periodicity 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 high-frequency variability of bed elevation and therefore gives a lower periodicity compared to the wavelets.

Discussion
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., 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.,  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 cross-section 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 hydro-morphological parameters (width, depth, meandering, etc.). This can have a very important impact on the identification of pseudo-periods.
The wavelet ridge analysis is powerful in identifying pseudo-periods, 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 pseudo-periodic 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 morpholo-gies. Ultimately, hydraulic modeling will be the true test of the potential of a pseudo-periodic equivalent geometry (e.g., for simulating a reach-average 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 non-stationary 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.

Conclusions
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 hydro-morphological 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 pseudo-periodic 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 proof-of-concept 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)
Cross-sectional 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 ) Local wavenumber that corresponds to the maximum variance of the signal (rad m −1 ) Local wavenumber that corresponds to the maximum variance of the variable i (rad m −1 ) K(x) support ( Translation factor in the wavelet transform or the abscissa position (m) y = y max Water depth measured from the talweg elevation y = z 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 (x) = φ(x, K(x)) (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 The conjugate form of the mother wavelet is Its derivative depending on the mute variable η is ψ * (η) = −π − 1 4 (iβ + η)e −iβη− η 2 2 (B2) = −(iβ + η)ψ * (η).
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 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 ∂ 2 φ i (x,k) ∂k∂x : For that, we should develop ∂ ∂k W[xf i (x)](x,k) We calculate each derivative: We find a general formulation with p = 0 . . . N − 2: The third one is the computation of ∂ 3 φ i (x,k) ∂k 2 ∂x :