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, higher-frequency 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.

1 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 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 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; 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 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).

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_{\mathrm{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.

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.

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)

Table 1 Review of some methods of morphological units' identification (variable used and MU types).

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 (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 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 pseudo-periodicity 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 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 continuous wavelet transform (CWT) to calculate the wavelength λ and the dimensionless wavelength spacing λ^* (longitudinal spacing) which is

$$\begin{array}{}\text{(1)}& {\mathit{\lambda }}^{*}={\displaystyle \frac{\mathit{\lambda }}{w_{\mathrm{bf}}}},\end{array}$$

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 pool-to-pool spacing/bankfull width and the channel slope expressed as a percentage based on a negative power function of slope S:

$$\begin{array}{}\text{(2)}& {\mathit{\lambda }}^{*}=\mathrm{8.2513}{S}^{-\mathrm{0.9799}}.\end{array}$$

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, 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 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 RiverBuilder (Brown et al., 2014).

2 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⁻¹ (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.

Figure 2 The locations of the study reaches in France.

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.

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:

1. velocity v(x);

2. 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 cross-sectional area and P(x) the wetted perimeter;

3. 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⁻³) and n Manning's roughness coefficient;

4. 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 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 arc-circle (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)).

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.

3 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) 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⁻¹). 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).

Figure 4 The wavelet Morlet mother function; the plot shows the real and imaginary parts of the wavelets in the space domain (distance).

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

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

$$\begin{array}{}\text{(3)}& \mathit{\psi }\left(\mathit{\eta }\right)={\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{4}}}{e}^{i\mathit{\beta }\mathit{\eta }}{e}^{-\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{2}}},\end{array}$$

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

$$\begin{array}{}\text{(4)}& {\mathit{\psi }}_{\mathrm{s},x}\left(\mathit{\eta }\right)={\displaystyle \frac{\mathrm{1}}{\sqrt{s}}}\mathit{\psi }\left({\displaystyle \frac{\mathit{\eta }-x}{s}}\right),x\in \mathbb{R},s>\mathrm{0},\end{array}$$

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

$$\begin{array}{}\text{(5)}& \begin{array}{rl}& \mathcal{W}\left[f\right]:\mathbb{R}\times {\mathbb{R}}_{+}^{*}\to \mathbb{C},\\ & (x,s)\mapsto {\displaystyle \frac{\mathrm{1}}{\sqrt{s}}}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right){\mathit{\psi }}^{*}\left({\displaystyle \frac{\mathit{\eta }-x}{s}}\right)\mathrm{d}\mathit{\eta }.\end{array}\end{array}$$

(^*) indicates the complex conjugate. This complex function can also be written as

$$\begin{array}{}\text{(6)}& \mathcal{W}\left[f\right](x,s)=R(x,s)e^{i\mathit{\varphi }(x,s)},\end{array}$$

where R is the absolute value (modulus) and ϕ is the phase (argument) at position x with the scale s.

$$\begin{array}{}\text{(7)}& {\displaystyle }& {\displaystyle }R(x,s)=\left|\mathcal{W}\right[f\left]\right(x,s\left)\right|\text{(8)}& {\displaystyle }& {\displaystyle }\mathit{\varphi }(x,s)=\mathrm{Im}\left(\ln \mathcal{W}\right[f\left(x\right)\left]\right(x,s\left)\right)\end{array}$$

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⁻¹) $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

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }\mathit{\alpha }={\displaystyle \frac{\mathrm{2}}{\mathit{\beta }+\sqrt{\mathrm{2}+{\mathit{\beta }}^{\mathrm{2}}}}},\text{(10)}& {\displaystyle }& {\displaystyle }s={\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }k}}.\end{array}$$

Thus, the wavelet transform of the function f(x) is defined in the space wavenumber as

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& \mathcal{W}\left[f\right]:\mathbb{R}\times {\mathbb{R}}_{+}^{*}\to \mathbb{C},\\ & (x,k)\mapsto \sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right){\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }.\end{array}\end{array}$$

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

Figure 5 (a) The phase function from which we get the function K(x); (b) the power of the wavelet with the region where there is maximum variability depicted by the black curve K(x) (ridge curve). These two figures are represented in a wavenumber/distance space for the Olivet River, and the wavelet transform is performed on the logarithm of the velocity. The part of the figure with low opacity shows the cone of influence, which is ignored in this study (edge effects are more important for short wavelengths than for long wavelengths).

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):

$$\begin{array}{}\text{(12)}& {\displaystyle \frac{\partial }{\partial x}}\mathrm{Im}\left(\ln \mathcal{W}\left[f\right(x\left)\right](x,k)\right)-k=\mathrm{0},\end{array}$$

or, according to the definition of the phase (Eq. 8):

$$\begin{array}{}\text{(13)}& {\left.{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right|}_{(x,k)}-k=\mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(14)}& {\left.{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right|}_{(x,K(x\left)\right)}-K\left(x\right)=\mathrm{0}.\end{array}$$

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.

$$\begin{array}{}\text{(15)}& \mathrm{\Phi }\left(x\right)=\mathit{\varphi }(x,K(x\left)\right)\end{array}$$

In the end, we can extract the wavelength function of pool–riffle sequences, which corresponds to a pseudo-period function of the signal f, and which is

$$\begin{array}{}\text{(16)}& \mathit{\lambda }\left(x\right)={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{K\left(x\right)}}.\end{array}$$

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)}$.

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{A}_{\mathrm{m}}\left(x\right)& =\left|\mathcal{W}\right[f\left]\right(x,K\left(x\right)\left)\right|\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }K\left(x\right)}}}\\ & =R(x,K(x\left)\right)\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }K\left(x\right)}}}\end{array}\end{array}$$

The signal is locally similar to a sinusoid f_mod of wavenumber K in rad m⁻¹, which models the variability of f. We can define the pseudo-periodic variable as presented in Fig. 6 with

$$\begin{array}{}\text{(18)}& f_{\mathrm{mod}}\left(x\right)=A_{\mathrm{m}}\left(x\right)\cos \left(\mathrm{\Phi }\right(x\left)\right)=A_{\mathrm{m}}\left(x\right)\cos \left(\mathit{\varphi }\right(x,K\left(x\right)\left)\right).\end{array}$$

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 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 pseudo-periodicity of v_mod results in the identification of six pools (white) and seven riffles (gray).

Figure 6 (a) Variation of the modeled function f_mod which represents the pseudo-periodic variable (e.g., the velocity of the Olivet River) compared to the observed one. This pseudo-periodicity results in the identification of pools (white) and riffles (gray) in panel (b). The unstudied part is due to the cone of influence of the wavelet method.

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.

3.3 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

$$\begin{array}{}\text{(19)}& {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right|}_{(x,K_i(x\left)\right)}-K_i\left(x\right)=\mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(20)}& {\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right|}_{(x,K(x\left)\right)}-K\left(x\right)\approx \mathrm{0}\forall i.\end{array}$$

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.\frac{\partial {\mathit{\varphi }}_{\mathrm{1}}}{\partial x}\right|}_{(x,k)}-k,{\left.\frac{\partial {\mathit{\varphi }}_{\mathrm{2}}}{\partial x}\right|}_{(x,k)}-k,\mathrm{\dots },{\left.\frac{\partial {\mathit{\varphi }}_N}{\partial x}\right|}_{(x,k)}-k\right)$:

$$\begin{array}{}\text{(21)}& E(x,k)=\sum _{i=\mathrm{1}}^N{\left({\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right|}_{(x,k)}-k\right)}^{\mathrm{2}}.\end{array}$$

This minimum is calculated by searching for the wavenumbers and positions where the derivatives (Eq. 22) of this quantity satisfy

$${\displaystyle \frac{\partial E(x,k)}{\partial k}}=\sum _{i=\mathrm{1}}^N\left({\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_i}{\partial k\partial x}}\right|}_{(x,k)}-\mathrm{1}\right)\left({\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right|}_{(x,k)}-k\right)=\mathrm{0}$$

and

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{\displaystyle \frac{{\partial }^{\mathrm{2}}E(x,k)}{\partial k^{\mathrm{2}}}}& =\sum _{i=\mathrm{1}}^N\left[{\left.{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\varphi }}_i}{\partial k^{\mathrm{2}}\partial x}}\right|}_{(x,k)}\left({\left.{\displaystyle \frac{\partial {\mathit{\varphi }}_i}{\partial x}}\right|}_{(x,k)}-k\right)\right.\\ & \left.+{\left({\left.{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_i}{\partial k\partial x}}\right|}_{(x,k)}-\mathrm{1}\right)}^{\mathrm{2}}\right]>\mathrm{0}.\end{array}\end{array}$$

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 co-evolution of all these variables.

Figure 7 Power of the wavelet of the four variables: velocity, hydraulic radius, bed shear stress, and local channel direction angle. The black curve K(x) is the extracted ridge curve of the Olivet River in the multivariate case.

As a result, the phase shift of every variable is calculated by

$$\begin{array}{}\text{(23)}& {\mathrm{\Phi }}_i\left(x\right)={\mathit{\varphi }}_i(x,K(x\left)\right).\end{array}$$

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.

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

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{f}_{i,\mathrm{mod}}\left(x\right)& =A_{i,\mathrm{m}}\left(x\right)\cos \left({\mathrm{\Phi }}_i\left(x\right)\right)\\ & =A_{i,\mathrm{m}}\left(x\right)\cos \left({\mathit{\varphi }}_i(x,K(x\left)\right)\right)\end{array}\end{array}$$

The amplitude shape of the modeled variable is calculated by the same way in the univariate case:

$$\begin{array}{}\text{(25)}& A_{i,\mathrm{m}}\left(x\right)=\left|\mathcal{W}\left[f_i\right](x,K(x\left)\right)\right|\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathit{\alpha }K\left(x\right)}}}.\end{array}$$

The results in Fig. 8 show that a common pseudo-period has been successfully identified and allows consistent pseudo-periodic 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).

Figure 9 (a) Correlation between the modeled functions f_i,mod, which represents the pseudo-periodic variables (velocity in red, hydraulic radius in blue, bed shear stress in green, and local channel direction angle in black) of the Olivet River. These pseudo-periodic functions result in the identification of pools (white) and riffles (gray) in panel (b).

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 pseudo-periodicity of the set [v_mod, R_h,mod, τ_b,mod, cos (θ)_mod] results in the identification of five pools and five riffles.

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

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.

Table 3 The assessed length provided by the wavelet analysis 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.

Moreover, the multivariate approach takes into account all the variables and therefore looks for a single pseudo-periodicity 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).

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 $\frac{\mathit{\lambda }}{w_{\mathrm{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.

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 reach-average 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 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.

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 ${\mathit{\lambda }}_{\mathrm{meth},\mathrm{mean}}^{*}$ is the mean longitudinal spacing, ${\mathit{\lambda }}_{\mathrm{meth},\mathrm{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.

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 pool-to-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 (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.

5 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., 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 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 morphologies. 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.

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

Appendix A: List of symbols

A(x)	Cross-sectional area (m2)
Ai,mod	Signal amplitude of the shape of the modeled variable number i
Am	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)
fi	Measured space series function number i
fi,mod	Modeled space series function number i with multivariate wavelet analysis
fmod	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)
Ki(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)
Qmin	Minimum discharge modeled (m3 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
Rh	Hydraulic radius (m)
Rh,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)
vmod	Modeled velocity (m s−1)
vobs	Measured velocity (m s−1)
w	Reach width (m)
wbf	Reach bankfull width (m)
wm	Mean bankfull width (m)
x	Translation factor in the wavelet transform or the abscissa position (m)
y=ymax	Water depth measured from the talweg elevation $y=z_{\mathrm{ws}}-z$ (m)
ym	Mean depth (m)
z	Measured bed elevation or talweg elevation (m)
zws	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 ψ

Appendix B: Mathematical calculus for the wavelet transform

B1 The univariate case

The conjugate form of the mother wavelet is

$$\begin{array}{}\text{(B1)}& \mathit{\psi }\left(\mathit{\eta }\right)={\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{4}}}{e}^{-i\mathit{\beta }\mathit{\eta }-\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{2}}}.\end{array}$$

Its derivative depending on the mute variable η is

$$\begin{array}{}\text{(B2)}& {\displaystyle }{{\mathit{\psi }}^{\prime }}^{*}\left(\mathit{\eta }\right)& {\displaystyle }=-{\mathit{\pi }}^{-\frac{\mathrm{1}}{\mathrm{4}}}(i\mathit{\beta }+\mathit{\eta })e^{-i\mathit{\beta }\mathit{\eta }-\frac{{\mathit{\eta }}^{\mathrm{2}}}{\mathrm{2}}}\text{(B3)}& {\displaystyle }& {\displaystyle }=-(i\mathit{\beta }+\mathit{\eta }){\mathit{\psi }}^{*}\left(\mathit{\eta }\right).\end{array}$$

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

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }\mathcal{W}\left[f\right(x\left)\right](x,k)=\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right){\displaystyle \frac{\partial }{\partial x}}\left[{\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\right]\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right)\cdot (-\mathit{\alpha }k)\cdot {{\mathit{\psi }}^{\prime }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left(\mathit{\alpha }k\right)\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right)\cdot (i\mathit{\beta }+\mathit{\alpha }k(\mathit{\eta }-x\left)\right)\cdot {\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left(\mathit{\alpha }k\right)\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\left[\right(i\mathit{\beta }-\mathit{\alpha }kx\left)f\right(\mathit{\eta })+\mathit{\alpha }k\mathit{\eta }f(\mathit{\eta }\left)\right]\cdot {\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left(\mathit{\alpha }k\right)(i\mathit{\beta }-\mathit{\alpha }kx)\mathcal{W}\left[f\right(x\left)\right](x,k)+\left(\mathit{\alpha }k{)}^{\mathrm{2}}\mathcal{W}\right[xf\left(x\right)\left]\right(x,k).\end{array}$$

On the other hand, we have

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}\mathcal{W}\left[f\right(x\left)\right](x,k)& {\displaystyle }={\displaystyle \frac{\partial }{\partial x}}\left(R(x,k)e^{i\mathit{\varphi }(x,k)}\right)\\ {\displaystyle }& {\displaystyle }=\left[{\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial R}{\partial x}}+i{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\right]R{e}^{i\mathit{\varphi }},\end{array}$$

$$\begin{array}{ll}{\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }Re\left(\ln \mathcal{W}\right[f\left(x\right)\left]\right(x,k\left)\right)={\displaystyle \frac{\mathrm{1}}{R}}{\displaystyle \frac{\partial R}{\partial x}}\\ {\displaystyle }& {\displaystyle }=Re\left({\displaystyle \frac{\mathrm{1}}{\mathcal{W}\left[f\right(x\left)\right](x,k)}}{\displaystyle \frac{\partial }{\partial x}}\mathcal{W}g\left[f\right(x\left)\right](x,k)\right),\\ {\displaystyle \frac{\partial }{\partial x}}& {\displaystyle }\mathrm{Im}\left(\ln \mathcal{W}\right[f\left(x\right)\left]\right(x,k\left)\right)={\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}\\ {\displaystyle }& {\displaystyle }=\mathrm{Im}\left({\displaystyle \frac{\mathrm{1}}{\mathcal{W}\left[f\right(x\left)\right](x,k)}}{\displaystyle \frac{\partial }{\partial x}}\mathcal{W}\left[f\right(x\left)\right](x,k)\right).\end{array}$$

Finally,

$$\begin{array}{}\text{(B4)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}& =\mathrm{Im}\left\{\left(\mathit{\alpha }k\right)(i\mathit{\beta }-\mathit{\alpha }kx)+(\mathit{\alpha }k{)}^{\mathrm{2}}{\displaystyle \frac{\mathcal{W}\left[xf\right(x\left)\right](x,k)}{\mathcal{W}\left[f\right(x\left)\right](x,k)}}\right\},\\ {\displaystyle \frac{\partial \mathit{\varphi }}{\partial x}}& =\left(\mathit{\alpha }k\right)\mathit{\beta }+(\mathit{\alpha }k{)}^{\mathrm{2}}\mathrm{Im}\left\{{\displaystyle \frac{\mathcal{W}\left[xf\right(x\left)\right](x,k)}{\mathcal{W}\left[f\right(x\left)\right](x,k)}}\right\}.\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(B5)}& {\displaystyle \frac{\partial {\mathit{\varphi }}_i(x,k)}{\partial x}}=\left(\mathit{\alpha }k\right)\mathit{\beta }+(\mathit{\alpha }k{)}^{\mathrm{2}}\mathrm{Im}\left\{{\displaystyle \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\}.\end{array}$$

The second one is the computation of $\frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_i(x,k)}{\partial k\partial x}$:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\varphi }}_i(x,k)}{\partial k\partial x}}& =\left(\mathit{\alpha }\mathit{\beta }\right)+\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}kIm\left\{{\displaystyle \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\}\\ & +(\mathit{\alpha }k{)}^{\mathrm{2}}\mathrm{Im}\left\{{\displaystyle \frac{\partial }{\partial k}}\left({\displaystyle \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)\right\}.\end{array}\end{array}$$

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)$:

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{\displaystyle \frac{\partial }{\partial k}}& \left({\displaystyle \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)={\displaystyle \frac{\mathrm{1}}{\mathcal{W}\left[f_i\left(x\right)\right](x,k)}}\\ & {\displaystyle \frac{\partial \mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)}{\partial k}}-{\displaystyle \frac{\mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)}{\mathcal{W}\left[f_i\left(x\right)\right](x,k{)}^{\mathrm{2}}}}\\ & {\displaystyle \frac{\partial \mathcal{W}\left[f_i\left(x\right)\right](x,k)}{\partial k}}.\end{array}\end{array}$$

We calculate each derivative:

$$\begin{array}{ll}{\displaystyle \frac{\partial \mathcal{W}\left[f_i\left(x\right)\right](x,k)}{\partial k}}& {\displaystyle }=\left({\displaystyle \frac{\mathrm{1}}{\sqrt{k}}}{\displaystyle \frac{\partial \sqrt{k}}{\partial k}}\right)\mathcal{W}\left[f_i\left(x\right)\right](x,k)\\ {\displaystyle }& {\displaystyle }+\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}f\left(\mathit{\eta }\right){\displaystyle \frac{\partial }{\partial k}}\left[{\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\right]\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}\right)\mathcal{W}\left[f_i\left(x\right)\right](x,k)+\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\\ {\displaystyle }& {\displaystyle }f\left(\mathit{\eta }\right)\mathit{\alpha }(\mathit{\eta }-x){{\mathit{\psi }}^{\prime }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}\right)\mathcal{W}\left[f_i\left(x\right)\right](x,k)+\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\\ {\displaystyle }& {\displaystyle }f\left(\mathit{\eta }\right)\mathit{\alpha }(\mathit{\eta }-x)(i\mathit{\beta }+\mathit{\alpha }k(\mathit{\eta }-x\left)\right){\mathit{\psi }}^{*}\\ {\displaystyle }& {\displaystyle }\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}\right)\mathcal{W}\left[f_i\left(x\right)\right](x,k)\\ {\displaystyle }& {\displaystyle }+\sqrt{\mathit{\alpha }k}\underset{-\mathrm{\infty }}{\overset{+\mathrm{\infty }}{\int }}\left[\left(i\mathit{\beta }-{\mathit{\alpha }}^{\mathrm{2}}k{x}^{\mathrm{2}}\right)\right.\\ {\displaystyle }& {\displaystyle }+\left(-i\mathit{\beta }\mathit{\alpha }+\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}kx\right)\mathit{\eta }\\ {\displaystyle }& {\displaystyle }\left.-\left({\mathit{\alpha }}^{\mathrm{2}}k\right){\mathit{\eta }}^{\mathrm{2}}\right]f\left(\mathit{\eta }\right){\mathit{\psi }}^{*}\left(\mathit{\alpha }k\right(\mathit{\eta }-x\left)\right)\mathrm{d}\mathit{\eta }\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}+i\mathit{\beta }\mathit{\alpha }x-{\mathit{\alpha }}^{\mathrm{2}}k{x}^{\mathrm{2}}\right)\mathcal{W}\left[f_i\left(x\right)\right](x,k)\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}kx-i\mathit{\beta }\mathit{\alpha }\right)\mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^{\mathrm{2}}k\right)\mathcal{W}\left[x^{\mathrm{2}}{f}_i\left(x\right)\right](x,k).\end{array}$$

We find a general formulation with p=0 … N−2:

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathcal{W}\left[f_i\left(x\right)\right](x,k)}{\partial k}}& =\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}+i\mathit{\beta }\mathit{\alpha }x-\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}k{x}^{\mathrm{2}}\right)\\ & \mathcal{W}\left[x^p{f}_i\left(x\right)\right](x,k)\\ & +\left(\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}kx-i\mathit{\beta }\mathit{\alpha }\right)\mathcal{W}\left[x^{p+\mathrm{1}}{f}_i\left(x\right)\right](x,k)\\ & -\left({\mathit{\alpha }}^{\mathrm{2}}k\right)\mathcal{W}\left[x^{p+\mathrm{2}}{f}_i\left(x\right)\right](x,k).\end{array}\end{array}$$

The third one is the computation of $\frac{{\partial }^{\mathrm{3}}{\mathit{\varphi }}_i(x,k)}{\partial k^{\mathrm{2}}\partial x}$:

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\varphi }}_i(x,k)}{\partial k^{\mathrm{2}}\partial x}}& =\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}\mathrm{Im}\left\{{\displaystyle \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\}\\ & +\mathrm{4}{\mathit{\alpha }}^{\mathrm{2}}kIm\left\{{\displaystyle \frac{\partial }{\partial k}}\left({\displaystyle \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)\right\}\\ & +(\mathit{\alpha }k{)}^{\mathrm{2}}\mathrm{Im}\left\{{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial k^{\mathrm{2}}}}\left({\displaystyle \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)\right\},\end{array}\end{array}$$

with

$$\begin{array}{ll}{\displaystyle \frac{{\partial }^{\mathrm{2}}}{\partial k^{\mathrm{2}}}}\left({\displaystyle \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)& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathcal{W}\left[f_i\left(x\right)\right](x,k)}}\\ {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)}{\partial k^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)}{\mathcal{W}\left[f_i\left(x\right)\right](x,k{)}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{W}\left[f_i\left(x\right)\right](x,k)}{\partial k^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{2}}{\mathcal{W}\left[f_i\left(x\right)\right](x,k{)}^{\mathrm{2}}}}\\ {\displaystyle }& {\displaystyle \frac{\partial \mathcal{W}\left[f_i\left(x\right)\right](x,k)}{\partial k}}\\ {\displaystyle }& {\displaystyle \frac{\partial \mathcal{W}\left[x{f}_i\left(x\right)\right](x,k)}{\partial k}}\end{array}$$

and

$$\begin{array}{ll}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{W}\left[x^p{f}_i\left(x\right)\right](x,k)}{\partial k^{\mathrm{2}}}}& {\displaystyle }=\left(-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}{k}^{\mathrm{2}}}}-{\mathit{\alpha }}^{\mathrm{2}}{x}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\mathcal{W}\left[x^p{f}_i\left(x\right)\right](x,k)+\left(\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}x\right)\\ {\displaystyle }& {\displaystyle }\mathcal{W}\left[x^{p+\mathrm{1}}{f}_i\left(x\right)\right](x,k)-\left({\mathit{\alpha }}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\mathcal{W}\left[x^{p+\mathrm{2}}{f}_i\left(x\right)\right](x,k)\\ {\displaystyle }& {\displaystyle }+\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}k}}+i\mathit{\beta }\mathit{\alpha }x-{\mathit{\alpha }}^{\mathrm{2}}k{x}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle \frac{\partial \mathcal{W}\left[x^p{f}_i\left(x\right)\right](x,k)}{\partial k}}\\ {\displaystyle }& {\displaystyle }+\left(\mathrm{2}{\mathit{\alpha }}^{\mathrm{2}}kx-i\mathit{\beta }\mathit{\alpha }\right)\\ {\displaystyle }& {\displaystyle \frac{\partial \mathcal{W}\left[x^{p+\mathrm{1}}{f}_i\left(x\right)\right](x,k)}{\partial k}}\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^{\mathrm{2}}k\right){\displaystyle \frac{\partial \mathcal{W}\left[x^{p+\mathrm{2}}{f}_i\left(x\right)\right](x,k)}{\partial k}}.\end{array}$$

Figure B1 Steps of determining the local wavenumber K(x) using the wavelet univariate ridge analysis of the velocity of the Olivet (3) reach, represented in the four panels. (a) The phase function ϕ(x,k); (b) the power’s cone of influence of the wavelet to characterize the region where there is a maximum variability of the velocity in the Neperian logarithm; (c) the function $\frac{\partial \mathit{\varphi }}{\partial x}$; (d) the function k.